11
$\begingroup$

Does anyone know what the MAPQ values produced by BWA-MEM mean?

I'm looking for something similar to what Keith Bradnam discovered for Tophat v 1.4.1, where he realized that:

0 = maps to 5 or more locations

1 = maps to 3-4 locations

3 = maps to 2 locations

255 = unique mapping

I'm familiar with the notion that MAPQ should be theoretically be related to the probability of an "incorrect" alignment (10^(-MAPQ/10)), although this is vaguely-specified enough that aligners in practice tend to actually just use something like the above.

Note that this is different than the BWA MAPQ scoring interpretation, because BWA-MEM gives MAPQ scores in the range $[0,60]$, rather than $[0,37]$ as has been established for BWA.

$\endgroup$
13
$\begingroup$

First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway.

Strictly speaking, you have two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the phred scaled probability of the alignment/base being wrong. It often amuses me that we question mapQ but take baseQ for granted. BaseQ is also scaled and discretized differently; even fewer people know how Illumina/pacbio/nanopore/historical sequencers estimates baseQ.

On the second question, Section 2 of the MAQ supplementary explains the theoretical aspects of mapQ, which is still correct today. Briefly, mapping quality consists of three components: 1) the probability of contamination; 2) the effect of mapping heuristics and 3) the error due to the repetitiveness of the reference. Only 3) can be modeled theoretically.

In case of bwa-mem, if we assume the matching score is 1, type-3 error is estimated with:

$$10/\log10\cdot[\log4\cdot(S_1-S_2)-\log n_{\rm sub}]$$

where $S_1$ is the best alignment score, $S_2$ is the second best and $n_{\rm sub}$ is the number of suboptimal alignments. Factor $\log 4$ comes from the scoring matrix. Factor $10/\log10$ is the Phred scale. This equation assumes gap-free alignment and is very close to Section 2.5.2 in the MAQ supplementary. It is ok-ish for short reads, but often overestimates for long reads. I am not aware of a practical approach in general cases. In addition to this method, you can estimate mapQ by read simulation: just try to find a function such that it fits empirical mapQ. Some have tried machine learning, too.

$\endgroup$
  • 1
    $\begingroup$ This has an ambiguous definition, at least because it doesn't define what happens when there are two equally similar and optimal matches in different parts of the reference index. If "second best" is allowed to have the same alignment score as "best", then this score reduces to $$10/\log10\cdot[-\log n_{\rm sub}]$$ in such a case. $\endgroup$ – gringer Sep 6 '17 at 14:20
  • $\begingroup$ Beautiful. Do you have a reference for the type-3 error estimation expression for bwa-mem? $\endgroup$ – ijoseph Sep 6 '17 at 17:14
  • $\begingroup$ @gringer Read the MAQ supplementary. It is an approximate on the condition that $S_1\gg S_2$. Furthermore, you do need to consider type-2 errors and pairing. The real estimate is full of ad hoc. ijoseph, that would the MAQ paper. $\endgroup$ – user172818 Sep 7 '17 at 12:02
4
$\begingroup$

Those numbers are not arbitrarily picked (well... maybe 255/60/40 is arbitrarily picked).

To convert from log10 Q values like these (also used for error rates in FASTQ files) to probabilities, divide the number by 10, negate it, then raise 10 to the power of the result.

A MAPQ score of 3 corresponds to a probability of an incorrect mapping of almost exactly 50%, which is as confident as you can be about a specific mapping if you know that it maps equally well to two locations:

> 10^-0.3
[1] 0.5011872

A MAPQ score of 2 is an incorrect mapping probability of 63% (i.e. close to 3 different locations), and a MAPQ of 1 is an incorrect mapping probability of 79% (close to 5 different locations). These fit with the numbers that you have described.

However, as you point out, BWA doesn't do this; it seems to ignore the number of mappable locations when considering the MAPQ score. I explore the reason why this is a problem in my question about multi-mapped reads.

I notice that there's a bit of an attempt to answer this in a similar previous question:

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

See the answer by user172818 for more insight into the details of BWA-MEM MAPQ scores.

$\endgroup$
  • $\begingroup$ "Similar" question is generous, that's basically a duplicate of my question that I didn't do due diligence to find! $\endgroup$ – ijoseph Aug 31 '17 at 21:09
  • 1
    $\begingroup$ So it basically comes down to understanding this uncommented C code. Maybe I'll try contacting Heng Li directly. github.com/lh3/bwa/blob/master/bwamem.c#L952-L976 $\endgroup$ – ijoseph Aug 31 '17 at 21:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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