1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
|
/*++
Module Name:
mapq.h
Abstract:
Support functions for mapping quality
Authors:
Bill Bolosky, December, 2012
Environment:
User mode service.
Revision History:
--*/
#pragma once
#include "directions.h"
void initializeMapqTables();
double mapqToProbability(int mapq); // The probability of a match for the given MAPQ
inline int computeMAPQ(
double probabilityOfAllCandidates,
double probabilityOfBestCandidate,
int score,
int popularSeedsSkipped)
{
probabilityOfAllCandidates = __max(probabilityOfAllCandidates, probabilityOfBestCandidate); // You'd think this wouldn't be necessary, but floating point limited precision causes it to be.
_ASSERT(probabilityOfBestCandidate >= 0.0);
// Special case for MAPQ 70, which we generate only if there is no evidence of a mismatch at all.
// cheese is off, so no special casing MAPQ 70. If you want to turn is back on, return these three lines and then change the two instance of 70 to 69 in the
// next set below (baseMAPQ =, twice).
//
// if (probabilityOfAllCandidates == probabilityOfBestCandidate && popularSeedsSkipped == 0 && score < 5) {
// return 70;
// }
double correctnessProbability = probabilityOfBestCandidate / probabilityOfAllCandidates;
int baseMAPQ;
if (correctnessProbability >= 1) {
baseMAPQ = 70;
} else {
baseMAPQ = __min(70, (int)(-10 * log10(1 - correctnessProbability)));
}
//
// Apply a penalty based on the number of overly popular seeds in the read
//
baseMAPQ = __max(0, baseMAPQ - __max(0, popularSeedsSkipped-10) / 2);
#ifdef TRACE_ALIGNER
printf("computeMAPQ called at %u: score %d, pThis %g, pAll %g, result %d\n",
location, score, probabilityOfBestCandidate, probabilityOfAllCandidates, baseMAPQ);
#endif
return baseMAPQ;
}
|
