BWA - number of sequences aligned

[Cbon]

Hello everyone,

I am quite a newbie in bioinformatics, and I have some difficulties using bwa software.

I successfully aligned my dataset to my reference genome. Now, I want to see the impact of the settings on my results, but I don't know how to obtain the number of reads aligned in the .sam file.

Is there a useful trick to obtain the number of reads aligned ?

I have found no clue on the bwa manual, neither in the samtool package and I begin to feel quite helpless ...

Thank you for your answers !

[dawe]

Hello everyone,

Is there a useful trick to obtain the number of reads aligned ?

Try

samtools flagstat [I]file.bam[/I]

Although I don't know if it works on sam files (in case you should convert in BAM, which I believe should be the default...).

d

[racetrout]

If you are using UNIX and GNU awk,

gawk '!/^[@#]/ && (!and($2,4){print}' bwa.sam| wc -l

should give you the number of reads aligned - providing BWA is setting the flag for unaligned reads. If you have many reads and BWA writes out only the aligned reads,

wc -l bwa.sam

is a very good simple approximate of aligned reads (it counts the header of course, but the ~5 lines of the header compared to the millions of reads are negligible).

[Cbon]

Thank you both for your answers !

My files are already converted in .bam, so I hope I can use the flagstat command. If it doesn't work, I will try the code.

thank you again !

[saml]

If you are using UNIX and GNU awk,

gawk '!/^[@#]/ && (!and($2,4){print}' bwa.sam| wc -l

should give you the number of reads aligned - providing BWA is setting the flag for unaligned reads. If you have many reads and BWA writes out only the aligned reads,

wc -l bwa.sam

is a very good simple approximate of aligned reads (it counts the header of course, but the ~5 lines of the header compared to the millions of reads are negligible).

Yes, but removing the lines is easy too. Skipping 2 header lines (think I had that), and starting on the 3rd line, can be done with the tail command and adding a '+' before the number of lines, to denote which line to start from:

tail -n +3 bwa.sam | gawk '!/^[@#]/ && (!and($2,4){print}' | wc -l

... or for the case with only aligned reads coming from bwa:

tail -n +3 bwa.sam | wc -l

[swbarnes2]

bwa by default includes all reads, even unaligned ones. You can tell that they are unmapped, from the binary tag.

So the flagstat command will distinguish between the mapped an unmapped reads.

You can also use samtools view -bF0x4 all.bam > mapped .bam to get a bam which has all the unmapped reads filtered out.

[nilshomer]

Use "samtools view -c -S <in.sam>" and add the "-f/-F" options as desired.

[solonas13]

I personally find it---at the best case---unacceptable the fact that a short-read aligner does not have an extra line of code to report the total number of aligned reads... Especially when it outputs (on the standard error) things like:

`[infer_isize] (25, 50, 75) percentile: (137, 179, 227)'

Anyhow, isn't there a closing parenthesis missing from this one?

tail -n +3 bwa.sam | gawk '!/^[@#]/ && (!and($2,4){print}' | wc -l

Thanks