Answers.Body: <p>One possibility uses bamutils (<a href="http://ngsutils.org/modules/bamutils/removeclipping/" rel="nofollow noreferrer">http://ngsutils.org/modules/bamutils/removeclipping/</a>) to remove the clipped bases, then awk to get the sum of the products, something like:</p>
  
  <pre><code>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}'
  </code></pre>
  
  <p>the bamutils step is rather slow, but it will do the job.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-03T08:44:19.850Z
Answers.Text: Why use three separate awk calls instead of `samtools view val_rc.bam | awk '{sum+=length($10)*$5} END { print "mapq*len sum =",sum}'` ?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-03T04:13:52.200Z
Answers.LastActivityDate: 2017-11-03T06:00:07.460Z
Answers.LastEditDate: 2017-11-03T06:00:07.460Z
Answers.UserDisplayName: roblanf
Answers.Body: <p>BAM files are stored line by line, so your only option is to process it that way. In python with <code>pysam</code>:</p>
  
  <pre class="lang-py prettyprint-override"><code>import pysam
  
  bam = pysam.AlignmentFile("something.bam")
  tot = 0
  for read in bam.fetch():
      if read.is_unmapped:
          continue
      tl = 0
      for op, l in read.cigartuples:
          if op in [0, 7, 8]:
              tl += l
      tot += tl * read.mapping_quality
  bam.close()
  print("the total is {}".format(tot))
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-03T08:21:24.693Z
Answers.LastActivityDate: 2017-11-03T08:21:24.693Z
Answers.UserDisplayName: Devon Ryan
Body: <p>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).</p>
  
  <p>For example, given two reads like this:</p>
  
  <pre><code>read1: mapq=40, mapped_bases=10
  read2: mapq=20, mapped_bases=50
  </code></pre>
  
  <p>The number I'm looking for is: </p>
  
  <pre><code>(40*10) + (20*50) = 1400
  </code></pre>
  
  <p>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.</p>
  
Tags: <sam><samtools><mapq>
Title: sum of products of mapq and mapped bases for each read in a from a BAM file
ThreadText: # sum of products of mapq and mapped bases for each read in a from a BAM file
  
  (Asked by: username_0 on 2017-11-03)
  
  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.
  
  **Tags:** sam, samtools, mapq
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2017-11-03. Score: 3 | Accepted answer
  
  BAM files are stored line by line, so your only option is to process it that way. In python with `pysam`:
  
  ```py
  import pysam
  
  bam = pysam.AlignmentFile("something.bam")
  tot = 0
  for read in bam.fetch():
      if read.is_unmapped:
          continue
      tl = 0
      for op, l in read.cigartuples:
          if op in [0, 7, 8]:
              tl += l
      tot += tl * read.mapping_quality
  bam.close()
  print("the total is {}".format(tot))
  
  ```
  
  
  ---
  
  ### Answer by username_0 on 2017-11-03. Score: 0
  
  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.
  
  > **Comments:**
  > **username_1**: Why use three separate awk calls instead of `samtools view val_rc.bam | awk '{sum+=length($10)*$5} END { print "mapq*len sum =",sum}'` ?  
  
  ---
  
  