SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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

Reply
 
Thread Tools
Old 01-04-2013, 11:25 PM   #1
repinementer
Member
 
Location: asia

Join Date: Dec 2009
Posts: 80
Default 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!
repinementer is offline   Reply With Quote
Old 01-05-2013, 03:44 AM   #2
rdeborja
Member
 
Location: Toronto

Join Date: Aug 2008
Posts: 42
Default

Try using samtools flagstat.
rdeborja is offline   Reply With Quote
Old 01-05-2013, 08:38 AM   #3
repinementer
Member
 
Location: asia

Join Date: Dec 2009
Posts: 80
Default

that gives the no.of mapped loci but not mapped reads.
repinementer is offline   Reply With Quote
Old 01-05-2013, 09:46 AM   #4
rdeborja
Member
 
Location: Toronto

Join Date: Aug 2008
Posts: 42
Default

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)
rdeborja is offline   Reply With Quote
Old 01-05-2013, 03:09 PM   #5
repinementer
Member
 
Location: asia

Join Date: Dec 2009
Posts: 80
Default

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.
repinementer is offline   Reply With Quote
Old 01-05-2013, 03:34 PM   #6
rdeborja
Member
 
Location: Toronto

Join Date: Aug 2008
Posts: 42
Default

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.
rdeborja is offline   Reply With Quote
Old 01-06-2013, 01:34 AM   #7
repinementer
Member
 
Location: asia

Join Date: Dec 2009
Posts: 80
Default

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 ?
repinementer is offline   Reply With Quote
Old 01-06-2013, 01:53 AM   #8
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 311
Default

Quote:
Originally Posted by repinementer View Post
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
dariober is offline   Reply With Quote
Old 01-06-2013, 06:06 AM   #9
kopi-o
Senior Member
 
Location: Stockholm, Sweden

Join Date: Feb 2008
Posts: 319
Default

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)
kopi-o is offline   Reply With Quote
Reply

Tags
rna-seq

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:14 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO