Answers.Body: <p>Yes, there bwa-mem was published as a <a href="https://arxiv.org/pdf/1303.3997.pdf" rel="nofollow noreferrer">preprint</a></p>
  
  <blockquote>
    <p>BWA-MEM’s seed extension differs from the standard seed extension in two aspects. Firstly, suppose at a certain extension step we come to reference
    position x with the best extension score achieved at query position y.</p>
    
    <p>...</p>
    
    <p>Secondly, while extending a seed, BWA-MEM tries to keep track of the
    best extension score reaching the end of the query sequence</p>
  </blockquote>
  
  <p>And there is a description of the scoring algorithm directly in the source code of <a href="https://github.com/lh3/bwa/blob/master/bwamem.c" rel="nofollow noreferrer">bwa-mem</a> (lines 22 - 44), but maybe the only solution is really to go though the source code.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-31T19:44:04.327Z
Answers.Text: Well I had no intention to actually explain the score (since I really do not know and the question explicitly asked about resources). I just knew about the preprint and I also got the idea to look at the source code, where I found those 22 lines of maths that seemed to explain the scoring.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-31T00:47:54.873Z
Answers.Text: Thank you for your answer. However, the question is more about assigning mapping qualities. Even reads with a very high alignment score can have a mapping quality equal to zero.
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-31T05:33:00.883Z
Answers.Text: Have you checked the source code? line 22 - 44.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-31T13:01:26.537Z
Answers.Text: @KamilSJaron Wow, that’s terribly hard to understand. The [actual code](https://github.com/lh3/bwa/blob/master/bwamem.c#L945-L969) unfortunately isn’t better. :-( At any rate, could you update your answer to include this more prominently?
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-31T13:15:28.960Z
Answers.Text: @ KamilSJaron I did and it is still not completely clear for me even in the easier case of single-end reads.
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T22:15:53.533Z
Answers.LastActivityDate: 2017-06-02T18:39:54.083Z
Answers.LastEditDate: 2017-06-02T18:39:54.083Z
Answers.UserDisplayName: Kamil S Jaron
Answers.Body: <p>The MAPQ score appears to be calculated in the mem_approx_mapq_se function. I tried porting it to Python and inserting comments using ChatGPT 4, though I haven't tested it to ensure the output is identical to the original.</p>
  <pre><code>import math
  
  # Constants
  MEM_MAPQ_COEFFICIENT = 30.0
  MEM_MAPQ_MAXIMUM = 60
  
  def calculate_mapping_quality(alignment_options, alignment_region):
      # Determine the suboptimal alignment score, which is the second-best SW alignment score if available,
      # otherwise calculated from minimum seed length and match score.
      suboptimal_alignment_score = (alignment_region.suboptimal_score if alignment_region.suboptimal_score 
                                    else alignment_options.minimum_seed_length * alignment_options.match_score)
  
      # Update the suboptimal alignment score to the higher of the two: current suboptimal score or tandem hit score.
      suboptimal_alignment_score = (alignment_region.tandem_hit_score if alignment_region.tandem_hit_score &gt; suboptimal_alignment_score 
                                    else suboptimal_alignment_score)
  
      # Early exit with zero mapping quality if the suboptimal score is greater than or equal to the alignment score.
      if suboptimal_alignment_score &gt;= alignment_region.best_alignment_score:
          return 0
  
      # Calculate the maximum alignment length between the query and reference sequences.
      alignment_length = max(alignment_region.query_end - alignment_region.query_start, 
                             alignment_region.reference_end - alignment_region.reference_start)
  
      # Calculate the identity of the alignment.
      identity = (1.0 - (alignment_length * alignment_options.match_score - alignment_region.best_alignment_score) 
                  / (alignment_options.match_score + alignment_options.mismatch_penalty) / alignment_length)
  
      # Compute the mapping quality based on the alignment score, identity, and other factors.
      if alignment_region.best_alignment_score == 0:
          mapping_quality = 0
      elif alignment_options.mapQ_coef_length &gt; 0:
          # Adjust mapping quality based on alignment length and identity.
          length_factor = (1.0 if alignment_length &lt; alignment_options.mapQ_coef_length 
                           else alignment_options.mapQ_coef_fac / math.log(alignment_length))
          length_factor *= identity * identity
          mapping_quality = int(6.02 * (alignment_region.best_alignment_score - suboptimal_alignment_score) 
                                / alignment_options.match_score * length_factor * length_factor + 0.499)
      else:
          # Use the constant MEM_MAPQ_COEFFICIENT for mapping quality calculation.
          mapping_quality = int(MEM_MAPQ_COEFFICIENT * (1.0 - suboptimal_alignment_score / alignment_region.best_alignment_score) 
                                * math.log(alignment_region.seed_coverage) + 0.499)
          mapping_quality = int(mapping_quality * identity * identity + 0.499) if identity &lt; 0.95 else mapping_quality
  
      # Penalize the mapping quality for the presence of suboptimal hits.
      if alignment_region.suboptimal_hit_count &gt; 0:
          mapping_quality -= int(4.343 * math.log(alignment_region.suboptimal_hit_count + 1) + 0.499)
  
      # Ensure the mapping quality does not exceed the maximum allowed value.
      mapping_quality = min(max(mapping_quality, 0), MEM_MAPQ_MAXIMUM)
  
      # Adjust the mapping quality based on the fraction of repetitive regions in the alignment.
      mapping_quality = int(mapping_quality * (1.0 - alignment_region.repetitive_region_fraction) + 0.499)
  
      # Return the final mapping quality value.
      return mapping_quality
  
  <span class="math-container">```</span>
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-11-09T08:52:32.393Z
Answers.LastActivityDate: 2023-11-09T08:52:32.393Z
Answers.UserDisplayName: Autodidactyle
Body: <p>Is there any resource (paper, blogpost, Github gist, etc.) describing the BWA-MEM algorithm for assigning mapping qualities? I vaguely remember that I have somewhere seen a formula for SE reads, which looked like</p>
  
  <p>$C * (s_1 - s_2) / s_1,$</p>
  
  <p>where $s_1$ and $s_2$ denoted the alignment scores of two best alignments and <code>C</code> was some constant.</p>
  
  <p>I believe that a reimplementation of this algorithm in some scripting language could be very useful for the bioinfo community. For instance, I sometimes test various mapping methods and some of them tend to find good alignments, but fail in assigning appropriate qualities. Therefore, I would like to re-assign all the mapping qualities in a SAM file with the BWA-MEM algorithm.</p>
  
  <p>Btw. This algorithm must already have been implemented outside BWA, see the BWA-MEM paper: </p>
  
  <blockquote>
    <p>GEM does not compute mapping quality. Its
    mapping quality is estimated with a BWA-like algorithm with suboptimal
    alignments available.</p>
  </blockquote>
  
  <p>Unfortunately, the <a href="https://github.com/lh3/mem-paper/" rel="noreferrer">BWA-MEM paper repo</a> contains only the resulting <code>.eval</code> files.</p>
  
  <p><strong>Update:</strong> The question is <em>not</em> about the algorithm for computing alignment scores. Mapping qualities and alignment scores are two different things:</p>
  
  <ul>
  <li>Alignment score quantifies the similarity between two sequences (e.g., a read and a reference sequence)</li>
  <li>Mapping quality (MAQ) quantifies the probability that a read is aligned to a wrong position.</li>
  </ul>
  
  <p>Even alignments with high scores can have a very low mapping quality.</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-01T13:44:37.613Z
Comments.Text: Maybe this page can help http://genome.sph.umich.edu/wiki/Mapping_Quality_Scores . Once you know the best and alternative positions where a read can align (or even the best and second best only?) it's not too difficult to implement I guess.
Comments.UserDisplayName: dariober
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-31T13:06:04.893Z
Comments.Text: Unfortunately I don’t know the answer for BWA-MEM (since it differs from BWA!) but pretty much all other tools are described here: https://sequencing.qcfail.com/articles/mapq-values-are-really-useful-but-their-implementation-is-a-mess/
Comments.UserDisplayName: Konrad Rudolph
Tags: <phylogenetics><bwa><read-mapping>
Title: How does the BWA-MEM algorithm assign its mapping qualities?
ThreadText: # How does the BWA-MEM algorithm assign its mapping qualities?
  
  (Asked by: username_0 on 2017-05-30)
  
  Is there any resource (paper, blogpost, Github gist, etc.) describing the BWA-MEM algorithm for assigning mapping qualities? I vaguely remember that I have somewhere seen a formula for SE reads, which looked like
  
  $C \* (s\_1 - s\_2) / s\_1,$
  
  where $s\_1$ and $s\_2$ denoted the alignment scores of two best alignments and `C` was some constant.
  
  I believe that a reimplementation of this algorithm in some scripting language could be very useful for the bioinfo community. For instance, I sometimes test various mapping methods and some of them tend to find good alignments, but fail in assigning appropriate qualities. Therefore, I would like to re-assign all the mapping qualities in a SAM file with the BWA-MEM algorithm.
  
  Btw. This algorithm must already have been implemented outside BWA, see the BWA-MEM paper:
  
  > GEM does not compute mapping quality. Its
  > mapping quality is estimated with a BWA-like algorithm with suboptimal
  > alignments available.
  
  Unfortunately, the [BWA-MEM paper repo](https://github.com/lh3/mem-paper/) contains only the resulting `.eval` files.
  
  **Update:** The question is *not* about the algorithm for computing alignment scores. Mapping qualities and alignment scores are two different things:
  
  * Alignment score quantifies the similarity between two sequences (e.g., a read and a reference sequence)
  * Mapping quality (MAQ) quantifies the probability that a read is aligned to a wrong position.
  
  Even alignments with high scores can have a very low mapping quality.
  
  **Tags:** phylogenetics, bwa, read-mapping
  
  > **Comments:**
  > **username_2**: Unfortunately I don’t know the answer for BWA-MEM (since it differs from BWA!) but pretty much all other tools are described here: https://sequencing.qcfail.com/articles/mapq-values-are-really-useful-but-their-implementation-is-a-mess/  
  > **username_1**: Maybe this page can help http://genome.sph.umich.edu/wiki/Mapping_Quality_Scores . Once you know the best and alternative positions where a read can align (or even the best and second best only?) it's not too difficult to implement I guess.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_3 on 2017-05-30. Score: 3
  
  Yes, there bwa-mem was published as a [preprint](https://arxiv.org/pdf/1303.3997.pdf)
  
  > BWA-MEM’s seed extension differs from the standard seed extension in two aspects. Firstly, suppose at a certain extension step we come to reference
  > position x with the best extension score achieved at query position y.
  >
  > ...
  >
  > Secondly, while extending a seed, BWA-MEM tries to keep track of the
  > best extension score reaching the end of the query sequence
  
  And there is a description of the scoring algorithm directly in the source code of [bwa-mem](https://github.com/lh3/bwa/blob/master/bwamem.c) (lines 22 - 44), but maybe the only solution is really to go though the source code.
  
  > **Comments:**
  > **username_0**: Thank you for your answer. However, the question is more about assigning mapping qualities. Even reads with a very high alignment score can have a mapping quality equal to zero.  
  > **username_3**: Have you checked the source code? line 22 - 44.  
  > **username_2**: @username_3 Wow, that’s terribly hard to understand. The [actual code](https://github.com/lh3/bwa/blob/master/bwamem.c#L945-L969) unfortunately isn’t better. :-( At any rate, could you update your answer to include this more prominently?  
  > **username_0**: @ username_3 I did and it is still not completely clear for me even in the easier case of single-end reads.  
  > **username_3**: Well I had no intention to actually explain the score (since I really do not know and the question explicitly asked about resources). I just knew about the preprint and I also got the idea to look at the source code, where I found those 22 lines of maths that seemed to explain the scoring.  
  
  ---
  
  ### Answer by username_4 on 2023-11-09. Score: 0
  
  The MAPQ score appears to be calculated in the mem\_approx\_mapq\_se function. I tried porting it to Python and inserting comments using ChatGPT 4, though I haven't tested it to ensure the output is identical to the original.
  
  ```
  import math
  
  # Constants
  MEM_MAPQ_COEFFICIENT = 30.0
  MEM_MAPQ_MAXIMUM = 60
  
  def calculate_mapping_quality(alignment_options, alignment_region):
      # Determine the suboptimal alignment score, which is the second-best SW alignment score if available,
      # otherwise calculated from minimum seed length and match score.
      suboptimal_alignment_score = (alignment_region.suboptimal_score if alignment_region.suboptimal_score 
                                    else alignment_options.minimum_seed_length * alignment_options.match_score)
  
      # Update the suboptimal alignment score to the higher of the two: current suboptimal score or tandem hit score.
      suboptimal_alignment_score = (alignment_region.tandem_hit_score if alignment_region.tandem_hit_score > suboptimal_alignment_score 
                                    else suboptimal_alignment_score)
  
      # Early exit with zero mapping quality if the suboptimal score is greater than or equal to the alignment score.
      if suboptimal_alignment_score >= alignment_region.best_alignment_score:
          return 0
  
      # Calculate the maximum alignment length between the query and reference sequences.
      alignment_length = max(alignment_region.query_end - alignment_region.query_start, 
                             alignment_region.reference_end - alignment_region.reference_start)
  
      # Calculate the identity of the alignment.
      identity = (1.0 - (alignment_length * alignment_options.match_score - alignment_region.best_alignment_score) 
                  / (alignment_options.match_score + alignment_options.mismatch_penalty) / alignment_length)
  
      # Compute the mapping quality based on the alignment score, identity, and other factors.
      if alignment_region.best_alignment_score == 0:
          mapping_quality = 0
      elif alignment_options.mapQ_coef_length > 0:
          # Adjust mapping quality based on alignment length and identity.
          length_factor = (1.0 if alignment_length < alignment_options.mapQ_coef_length 
                           else alignment_options.mapQ_coef_fac / math.log(alignment_length))
          length_factor *= identity * identity
          mapping_quality = int(6.02 * (alignment_region.best_alignment_score - suboptimal_alignment_score) 
                                / alignment_options.match_score * length_factor * length_factor + 0.499)
      else:
          # Use the constant MEM_MAPQ_COEFFICIENT for mapping quality calculation.
          mapping_quality = int(MEM_MAPQ_COEFFICIENT * (1.0 - suboptimal_alignment_score / alignment_region.best_alignment_score) 
                                * math.log(alignment_region.seed_coverage) + 0.499)
          mapping_quality = int(mapping_quality * identity * identity + 0.499) if identity < 0.95 else mapping_quality
  
      # Penalize the mapping quality for the presence of suboptimal hits.
      if alignment_region.suboptimal_hit_count > 0:
          mapping_quality -= int(4.343 * math.log(alignment_region.suboptimal_hit_count + 1) + 0.499)
  
      # Ensure the mapping quality does not exceed the maximum allowed value.
      mapping_quality = min(max(mapping_quality, 0), MEM_MAPQ_MAXIMUM)
  
      # Adjust the mapping quality based on the fraction of repetitive regions in the alignment.
      mapping_quality = int(mapping_quality * (1.0 - alignment_region.repetitive_region_fraction) + 0.499)
  
      # Return the final mapping quality value.
      return mapping_quality
  
  ```
  
  ```
  
  
  ---
  
  