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