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

repinementer · 01-04-2013, 11:25 PM

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 · 01-05-2013, 03:44 AM

Try using samtools flagstat.

repinementer · 01-05-2013, 08:38 AM

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

rdeborja · 01-05-2013, 09:46 AM

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 · 01-05-2013, 03:09 PM

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 · 01-05-2013, 03:34 PM

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 · 01-06-2013, 01:34 AM

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 · 01-06-2013, 01:53 AM

Quote:

Originally Posted by 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 ?

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 · 01-06-2013, 06:06 AM

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)

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Extract mapped reads at certain position from paired end RNA-seq data	RedMary	Bioinformatics	2	09-04-2012 07:16 AM
paired-end reads mapped to genome.. gene with only one direction of paired-end reads?	danwiththeplan	Bioinformatics	2	09-22-2011 03:06 AM
How to count Paired-End reads for RNA-Seq read counts	epistatic	RNA Sequencing	0	08-31-2011 03:44 PM
RNA-seq: Replicates, single-end, paired-end story	pasta	Bioinformatics	2	07-05-2011 12:51 AM

01-04-2013, 11:25 PM	#1
repinementer Member Location: asia Join Date: Dec 2009 Posts: 80	How to count number of mapped paired-end and single-end rna-seq reads 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!

01-05-2013, 03:34 PM	#6
rdeborja Member Location: Toronto Join Date: Aug 2008 Posts: 42	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. Last edited by rdeborja; 01-05-2013 at 03:42 PM.

01-05-2013, 03:44 AM	#2
rdeborja Member Location: Toronto Join Date: Aug 2008 Posts: 42	Try using samtools flagstat.

01-05-2013, 08:38 AM	#3
repinementer Member Location: asia Join Date: Dec 2009 Posts: 80	that gives the no.of mapped loci but not mapped reads.

01-05-2013, 09:46 AM	#4
rdeborja Member Location: Toronto Join Date: Aug 2008 Posts: 42	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)

01-05-2013, 03:09 PM	#5
repinementer Member Location: asia Join Date: Dec 2009 Posts: 80	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.

01-06-2013, 01:34 AM	#7
repinementer Member Location: asia Join Date: Dec 2009 Posts: 80	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 ?

01-06-2013, 06:06 AM	#9
kopi-o Senior Member Location: Stockholm, Sweden Join Date: Feb 2008 Posts: 319	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)