Given a BAM file, I'd like to calculate the sum, over all reads, of the mapping quality and the number of mapped bases (i.e. number of M's in the CIGAR string).

For example, given two reads like this:

read1: mapq=40, mapped_bases=10
read2: mapq=20, mapped_bases=50

The number I'm looking for is:

(40*10) + (20*50) = 1400

One option here is to stream the file and get the two quantities line-by-line. But this seems likely to be slow. Thoughts welcome.


BAM files are stored line by line, so your only option is to process it that way. In python with pysam:

import pysam

bam = pysam.AlignmentFile("something.bam")
tot = 0
for read in bam.fetch():
    if read.is_unmapped:
    tl = 0
    for op, l in read.cigartuples:
        if op in [0, 7, 8]:
            tl += l
    tot += tl * read.mapping_quality
print("the total is {}".format(tot))

One possibility uses bamutils (http://ngsutils.org/modules/bamutils/removeclipping/) to remove the clipped bases, then awk to get the sum of the products, something like:

bamutils removeclipping val.bam val_rc.bam
samtools view val_rc.bam | awk '{print length($10)*$5}' | awk '{sum+=$0} END { print "mapq*len sum =",sum}'

the bamutils step is rather slow, but it will do the job.

  • $\begingroup$ Why use three separate awk calls instead of samtools view val_rc.bam | awk '{sum+=length($10)*$5} END { print "mapq*len sum =",sum}' ? $\endgroup$
    – terdon
    Nov 3 '17 at 8:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.