<A href="<a href=" http:="" samtools.sourceforge.net="" "="" rel="nofollow">http://samtools.sourceforge.net/">samtools's flagstat command can tell you basic information like number of reads aligned from looking at the flags in the sam/bam output. Use:
samtools flagstat bwaOutput.bam
$ bwa aln $ref $input | bwa samse $ref - $input | samtools view -S -b -o $output - $ samtools index $output $ samtools flagstat $output
'samtools import' doesn't seem to be in their latest documentation. I think I used to use it too, but 'view' is now the prefered command for converting .sam to.bam
You need to make a .bam from bwa's .sam, and you run flagstat on that. You don't have to have sorted it.
Thanks @swbarnes2. Commands are from script developed circa summer 2010 and clearly need to be updated.
updated as suggested.
It is slightly dumb that flagstat won't work on .sam files, but really, there's pretty much no reason not to pipe from bwa sampe/samse straight into samtools view to convert on the fly to .bam format. .sam format just takes up so much space, and there's very little that you can do on .sam files that you can't do on.bam files.
This is the long way 'round, I know, but I usually extract the unmapped reads using a perl script which compares the alignment file to the input read file. I then count the number of lines (wc -l) on the input read file and unmapped read file (of course, dividing by 4 for fastq files).
The script I use is from this seqanswers post: http://seqanswers.com/forums/showthread.php?t=6847
Hope you find this useful!
Log in to answer this question.
Are you talking about percentage of total reads aligned? Alignment identity? Alignment score is a raw score based on number of matches/mismathces/gaps. I don't think it's represented as a percentage.
Was talking about percentage of total reads aligned. Sorry for the confusion then.