Meaning of BWA-MEM MAPQ scores

Question

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.

Answer 1

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.

Answer 2

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.

Another way of looking at it is to consider the decade to mean the number of 9s in the [mapping] accuracy, e.g.:

Q10: 90%
Q20: 99%
Q30: 99.9%
Q40: 99.99%

With the formula adjusted to make all the numbers in between a smooth transition.

With this transformation, 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.