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

CC BY-SA 3.0
2017-05-31T19:44:04.327Z
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.
Kamil S Jaron
CC BY-SA 3.0
2017-05-31T00:47:54.873Z
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.
Karel Břinda
CC BY-SA 3.0
2017-05-31T05:33:00.883Z
Have you checked the source code? line 22 - 44.
Kamil S Jaron
CC BY-SA 3.0
2017-05-31T13:01:26.537Z
@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?
Konrad Rudolph
CC BY-SA 3.0
2017-05-31T13:15:28.960Z
@ KamilSJaron I did and it is still not completely clear for me even in the easier case of single-end reads.
Karel Břinda
CC BY-SA 3.0
2017-05-30T22:15:53.533Z
2017-06-02T18:39:54.083Z
2017-06-02T18:39:54.083Z
Kamil S Jaron
<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>

CC BY-SA 4.0
2023-11-09T08:52:32.393Z
2023-11-09T08:52:32.393Z
Autodidactyle
<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>

CC BY-SA 3.0
2017-06-01T13:44:37.613Z
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.
dariober
CC BY-SA 3.0
2017-05-31T13:06:04.893Z
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/
Konrad Rudolph
<phylogenetics><bwa><read-mapping>
How does the BWA-MEM algorithm assign its mapping qualities?
# 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**: @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?  
> **username_0**: @ KamilSJaron 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

```


---

