SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Extract mean mapped base phred quality by mapped location from BAM file ugm6hr Bioinformatics 2 10-05-2012 03:52 AM
Filter paired & mapped reads from SAM/BAM file pirates.genome Bioinformatics 2 06-19-2012 01:21 AM
How to get mapped reads ID from bam files using samtools? fabrice Bioinformatics 15 08-20-2011 06:03 PM
Extract perfectly mapped reads from SAM/BAM file Graham Etherington Bioinformatics 2 07-21-2011 07:27 AM
SAM/BAM how to retrieve all mapped reads ? Sylphide Bioinformatics 2 05-31-2011 06:57 AM

Reply
 
Thread Tools
Old 04-07-2014, 01:26 AM   #1
thedamian
Member
 
Location: Barcelona

Join Date: Feb 2012
Posts: 49
Default mapped reads statistics (from .bam)

Hi,

Do you know if it is possible to obtain such statistics from the .bam file:

Total Unique - total reads aligning to only a single location in the reference
U0 - # of unique hits with 0 mismatches
U1 - # of unique hits with 1 mismatches
U2 - # of unique hits with 2 mismatches
R012 - # of hits aligning to multiple location in the reference with 0,1 or 2 mismatches
NM - # of reads without a match to the reference

(I know samtools flagstat but it is not what I need)

Thanks in advance
thedamian is offline   Reply With Quote
Old 04-07-2014, 03:46 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Yes, it's possible, but something like samtools flagstat won't produce that. You could code something easily enough, though. Also, think long and hard about what you mean when you write "unique". What you likely really want are reads whose probability of incorrect alignment is sufficiently low (i.e., their MAPQ is sufficiently high).
dpryan is offline   Reply With Quote
Old 04-07-2014, 07:06 AM   #3
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 123
Default

You can use awk or perl to count your frequencies with an associative array. If you have the computational resources, you may count the reads occurring only once in your alignment. This can lead to problems, dpryan already mentioned.
Moreover, you should have a look at the fields in the sam-file's alignment section, which were set by your aligner.
The number of mismatches can be extracted from the CIGAR-string or directly found in the NM-field. Some aligners are clipping the reads; the clipping is not represented in the NM-field.

Just have a look at http://samtools.sourceforge.net/SAM1.pdf for further information.
Michael.Ante is offline   Reply With Quote
Old 04-07-2014, 07:11 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

A warning regarding the CIGAR field is that mismatches are typically not indicated unless they're indels. A few aligners, such as BBMap, can be told to include this information in the cigar string, but you're otherwise stuck parsing the MD auxiliary tag (remember that the NM tag indicates the edit distance, which can include indels...whether those should be included in the counts depends on ones goals).

BTW, I'd recommend using the pysam module in python, since you then don't have to reinvent the wheel in terms of parsing the alignments.
dpryan is offline   Reply With Quote
Old 04-07-2014, 06:24 PM   #5
jwfoley
Senior Member
 
Location: Stanford

Join Date: Jun 2009
Posts: 181
Default

Quote:
Originally Posted by dpryan View Post
Also, think long and hard about what you mean when you write "unique". What you likely really want are reads whose probability of incorrect alignment is sufficiently low (i.e., their MAPQ is sufficiently high).
Yes, counting mismatches is not generally a useful thing to do, and "unique" best hits aren't nearly stringent enough; if you want to count confident alignments, the command is "samtools view -cq 10 somefile.bam" (or a more stringent number, but 10 works fine for me).
jwfoley is offline   Reply With Quote
Reply

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 02:19 AM.


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