1 |
#ifndef QUAL_H_ |
2 |
#define QUAL_H_ |
3 |
|
4 |
/* NOTE: This file was written by Ben Langmead, and is borrowed from Bowtie */ |
5 |
|
6 |
#include <stdexcept> |
7 |
#include <iostream> |
8 |
#include "assert_helpers.h" |
9 |
|
10 |
using namespace std; |
11 |
|
12 |
extern unsigned char qualRounds[]; |
13 |
extern unsigned char solToPhred[]; |
14 |
|
15 |
/// Translate a Phred-encoded ASCII character into a Phred quality |
16 |
static inline uint8_t phredCharToPhredQual(char c) { |
17 |
return ((uint8_t)c >= 33 ? ((uint8_t)c - 33) : 0); |
18 |
} |
19 |
|
20 |
/** |
21 |
* Convert a Solexa-scaled quality value into a Phred-scale quality |
22 |
* value. |
23 |
* |
24 |
* p = probability that base is miscalled |
25 |
* Qphred = -10 * log10 (p) |
26 |
* Qsolexa = -10 * log10 (p / (1 - p)) |
27 |
* See: http://en.wikipedia.org/wiki/FASTQ_format |
28 |
* |
29 |
*/ |
30 |
static inline uint8_t solexaToPhred(int sol) { |
31 |
assert_lt(sol, 256); |
32 |
if(sol < -10) return 0; |
33 |
return solToPhred[sol+10]; |
34 |
} |
35 |
|
36 |
class SimplePhredPenalty { |
37 |
public: |
38 |
static uint8_t mmPenalty (uint8_t qual) { |
39 |
return qual; |
40 |
} |
41 |
static uint8_t delPenalty(uint8_t qual) { |
42 |
return qual; |
43 |
} |
44 |
static uint8_t insPenalty(uint8_t qual_left, uint8_t qual_right) { |
45 |
return std::max(qual_left, qual_right); |
46 |
} |
47 |
}; |
48 |
|
49 |
class MaqPhredPenalty { |
50 |
public: |
51 |
static uint8_t mmPenalty (uint8_t qual) { |
52 |
return qualRounds[qual]; |
53 |
} |
54 |
static uint8_t delPenalty(uint8_t qual) { |
55 |
return qualRounds[qual]; |
56 |
} |
57 |
static uint8_t insPenalty(uint8_t qual_left, uint8_t qual_right) { |
58 |
return qualRounds[std::max(qual_left, qual_right)]; |
59 |
} |
60 |
}; |
61 |
|
62 |
static inline uint8_t mmPenalty(bool maq, uint8_t qual) { |
63 |
if(maq) { |
64 |
return MaqPhredPenalty::mmPenalty(qual); |
65 |
} else { |
66 |
return SimplePhredPenalty::mmPenalty(qual); |
67 |
} |
68 |
} |
69 |
|
70 |
static inline uint8_t delPenalty(bool maq, uint8_t qual) { |
71 |
if(maq) { |
72 |
return MaqPhredPenalty::delPenalty(qual); |
73 |
} else { |
74 |
return SimplePhredPenalty::delPenalty(qual); |
75 |
} |
76 |
} |
77 |
|
78 |
static inline uint8_t insPenalty(bool maq, uint8_t qual_left, uint8_t qual_right) { |
79 |
if(maq) { |
80 |
return MaqPhredPenalty::insPenalty(qual_left, qual_right); |
81 |
} else { |
82 |
return SimplePhredPenalty::insPenalty(qual_left, qual_right); |
83 |
} |
84 |
} |
85 |
|
86 |
/** |
87 |
* Take an ASCII-encoded quality value and convert it to a Phred33 |
88 |
* ASCII char. |
89 |
*/ |
90 |
inline static char charToPhred33(char c, bool solQuals, bool phred64Quals) { |
91 |
if(c == ' ') { |
92 |
cerr << "Saw a space but expected an ASCII-encoded quality value." << endl |
93 |
<< "Are quality values formatted as integers? If so, try --integer-quals." << endl; |
94 |
throw 1; |
95 |
} |
96 |
if (solQuals) { |
97 |
// Convert solexa-scaled chars to phred |
98 |
// http://maq.sourceforge.net/fastq.shtml |
99 |
char cc = solexaToPhred((int)c - 64) + 33; |
100 |
if (cc < 33) { |
101 |
cerr << "Saw ASCII character " |
102 |
<< ((int)c) |
103 |
<< " but expected 64-based Solexa qual (converts to " << (int)cc << ")." << endl |
104 |
<< "Try not specifying --solexa-quals." << endl; |
105 |
throw 1; |
106 |
} |
107 |
c = cc; |
108 |
} |
109 |
else if(phred64Quals) { |
110 |
if (c < 64) { |
111 |
cerr << "Saw ASCII character " |
112 |
<< ((int)c) |
113 |
<< " but expected 64-based Phred qual." << endl |
114 |
<< "Try not specifying --solexa1.3-quals/--phred64-quals." << endl; |
115 |
throw 1; |
116 |
} |
117 |
// Convert to 33-based phred |
118 |
c -= (64-33); |
119 |
} |
120 |
else { |
121 |
// Keep the phred quality |
122 |
if (c < 33) { |
123 |
cerr << "Saw ASCII character " |
124 |
<< ((int)c) |
125 |
<< " but expected 33-based Phred qual." << endl; |
126 |
throw 1; |
127 |
} |
128 |
} |
129 |
return c; |
130 |
} |
131 |
|
132 |
/** |
133 |
* Take an integer quality value and convert it to a Phred33 ASCII |
134 |
* char. |
135 |
*/ |
136 |
inline static char intToPhred33(int iQ, bool solQuals) { |
137 |
int pQ; |
138 |
if (solQuals) { |
139 |
// Convert from solexa quality to phred |
140 |
// quality and translate to ASCII |
141 |
// http://maq.sourceforge.net/qual.shtml |
142 |
pQ = solexaToPhred((int)iQ) + 33; |
143 |
} else { |
144 |
// Keep the phred quality and translate |
145 |
// to ASCII |
146 |
pQ = (iQ <= 93 ? iQ : 93) + 33; |
147 |
} |
148 |
if (pQ < 33) { |
149 |
cerr << "Saw negative Phred quality " << ((int)pQ-33) << "." << endl; |
150 |
throw 1; |
151 |
} |
152 |
assert_geq(pQ, 0); |
153 |
return (int)pQ; |
154 |
} |
155 |
|
156 |
#endif /*QUAL_H_*/ |