mapped reads statistics (from .bam)

[thedamian]

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

[dpryan]

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).

[Michael.Ante]

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.

[dpryan]

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.

[jwfoley]

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).