How to count number of mapped paired-end and single-end rna-seq reads

[repinementer]

Does any one know know, how to count number of mapped paired-end and single-end rna-seq reads using BAM files.
It seems samtools idx stats does not give exactly mapped reads information ? Any suggestion will be appreciated!

[rdeborja]

Try using samtools flagstat.

[repinementer]

that gives the no.of mapped loci but not mapped reads.

[rdeborja]

It generates a summary of reads based on the SAM FLAG in column 2 of the BAM file:

4255310402 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
4252238423 + 0 mapped (99.93%:nan%)
4255310402 + 0 paired in sequencing
2102851470 + 0 read1
2152458932 + 0 read2
362406042 + 0 properly paired (8.52%:nan%)
4217472878 + 0 with itself and mate mapped
34765545 + 0 singletons (0.82%:nan%)
3616654841 + 0 with mate mapped to a different chr
3273787 + 0 with mate mapped to a different chr (mapQ>=5)

[repinementer]

99.93% mapping ? I think it is not referring 99.93% of your reads are mapped. 100% mapping is not possible or at least too good be true.

[rdeborja]

Yes, 99.93% read mapping, although it doesn't include the quality of the mapping. You'll have to look that up in the BAM file independently.

[repinementer]

If you look at any published studies (2010-12), you will typically see 80-90% but not ~100%. What thats tells ? Tophat always reports 100%. Something wrong isn't it ?

[dariober]

If you look at any published studies (2010-12), you will typically see 80-90% but not ~100%. What thats tells ? Tophat always reports 100%. Something wrong isn't it ?

There's nothing wrong there. If I remember correctly Tophat produces bam files containing only the mapped reads (accepted_hits.bam). The unmapped reads are written to a separate file I think. That's the reason why the bam files have 100% mapped reads (in fact it shoud be 100% not ~99%).

Dario

[kopi-o]

I find it helpful to use bam_stat.py from RSeQC or Picard's CollectAlignmentSummaryMetrics to get the number of reads that mapped one or more times (which you don't get from flagstat)