14
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Unfortunately I don’t know the answer for BWA-MEM (since it differs from BWA!) but pretty much all other tools are described here: sequencing.qcfail.com/articles/… $\endgroup$
    – Konrad Rudolph
    May 31, 2017 at 13:06
  • $\begingroup$ Maybe this page can help 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. $\endgroup$
    – dariober
    Jun 1, 2017 at 13:44

1 Answer 1

Reset to default
3
$\begingroup$

Yes, there bwa-mem was published as a preprint

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 (lines 22 - 44), but maybe the only solution is really to go though the source code.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Karel Břinda
    May 31, 2017 at 0:47
  • 1
    $\begingroup$ Have you checked the source code? line 22 - 44. $\endgroup$
    – Kamil S Jaron
    May 31, 2017 at 5:33
  • 5
    $\begingroup$ @KamilSJaron Wow, that’s terribly hard to understand. The actual code unfortunately isn’t better. :-( At any rate, could you update your answer to include this more prominently? $\endgroup$
    – Konrad Rudolph
    May 31, 2017 at 13:01
  • $\begingroup$ @ KamilSJaron I did and it is still not completely clear for me even in the easier case of single-end reads. $\endgroup$
    – Karel Břinda
    May 31, 2017 at 13:15
  • $\begingroup$ 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. $\endgroup$
    – Kamil S Jaron
    May 31, 2017 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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