6 extern unsigned char qualRounds[];
7 extern unsigned char solToPhred[];
9 /// Translate a Phred-encoded ASCII character into a Phred quality
10 static inline uint8_t phredCharToPhredQual(char c) {
11 return ((uint8_t)c >= 33 ? ((uint8_t)c - 33) : 0);
15 * Convert a Solexa-scaled quality value into a Phred-scale quality
18 * p = probability that base is miscalled
19 * Qphred = -10 * log10 (p)
20 * Qsolexa = -10 * log10 (p / (1 - p))
21 * See: http://en.wikipedia.org/wiki/FASTQ_format
24 static inline uint8_t solexaToPhred(int sol) {
26 if(sol < -10) return 0;
27 return solToPhred[sol+10];
30 class SimplePhredPenalty {
32 static uint8_t mmPenalty (uint8_t qual) {
35 static uint8_t delPenalty(uint8_t qual) {
38 static uint8_t insPenalty(uint8_t qual_left, uint8_t qual_right) {
39 return std::max(qual_left, qual_right);
43 class MaqPhredPenalty {
45 static uint8_t mmPenalty (uint8_t qual) {
46 return qualRounds[qual];
48 static uint8_t delPenalty(uint8_t qual) {
49 return qualRounds[qual];
51 static uint8_t insPenalty(uint8_t qual_left, uint8_t qual_right) {
52 return qualRounds[std::max(qual_left, qual_right)];
56 static inline uint8_t mmPenalty(bool maq, uint8_t qual) {
58 return MaqPhredPenalty::mmPenalty(qual);
60 return SimplePhredPenalty::mmPenalty(qual);
64 static inline uint8_t delPenalty(bool maq, uint8_t qual) {
66 return MaqPhredPenalty::delPenalty(qual);
68 return SimplePhredPenalty::delPenalty(qual);
72 static inline uint8_t insPenalty(bool maq, uint8_t qual_left, uint8_t qual_right) {
74 return MaqPhredPenalty::insPenalty(qual_left, qual_right);
76 return SimplePhredPenalty::insPenalty(qual_left, qual_right);
81 * Take an ASCII-encoded quality value and convert it to a Phred33
84 inline static char charToPhred33(char c, bool solQuals, bool phred64Quals) {
86 cerr << "Saw a space but expected an ASCII-encoded quality value." << endl
87 << "Are quality values formatted as integers? If so, try --integer-quals." << endl;
91 // Convert solexa-scaled chars to phred
92 // http://maq.sourceforge.net/fastq.shtml
93 char cc = solexaToPhred((int)c - 64) + 33;
95 cerr << "Saw ASCII character "
97 << " but expected 64-based Solexa qual (converts to " << (int)cc << ")." << endl
98 << "Try not specifying --solexa-quals." << endl;
103 else if(phred64Quals) {
105 cerr << "Saw ASCII character "
107 << " but expected 64-based Phred qual." << endl
108 << "Try not specifying --solexa1.3-quals/--phred64-quals." << endl;
111 // Convert to 33-based phred
115 // Keep the phred quality
117 cerr << "Saw ASCII character "
119 << " but expected 33-based Phred qual." << endl;
127 * Take an integer quality value and convert it to a Phred33 ASCII
130 inline static char intToPhred33(int iQ, bool solQuals) {
133 // Convert from solexa quality to phred
134 // quality and translate to ASCII
135 // http://maq.sourceforge.net/qual.shtml
136 pQ = solexaToPhred((int)iQ) + 33;
138 // Keep the phred quality and translate
140 pQ = (iQ <= 93 ? iQ : 93) + 33;
143 cerr << "Saw negative Phred quality " << ((int)pQ-33) << "." << endl;
151 * Fill the q[] array with the penalties that are determined by
152 * subtracting the quality values of the alternate basecalls from
153 * the quality of the primary basecall.
155 inline static uint8_t penaltiesAt(size_t off, uint8_t *q,
157 const String<char>& qual,
158 const String<Dna5>* altQry,
159 const String<char>* altQual)
161 uint8_t primQ = qual[off]; // qual of primary call
162 uint8_t bestPenalty = primQ - 33;
163 q[0] = q[1] = q[2] = q[3] = bestPenalty;
164 for(int i = 0; i < alts; i++) {
165 uint8_t altQ = altQual[i][off]; // qual of alt call
166 if(altQ == 33) break; // no alt call
167 assert_leq(altQ, primQ);
168 if(primQ - altQ < bestPenalty) {
169 bestPenalty = primQ - altQ;
172 int altC = (int)altQry[i][off];
174 q[altC] = primQ - altQ;
180 * Fill the q[] array with the penalties that are determined by
181 * subtracting the quality values of the alternate basecalls from
182 * the quality of the primary basecall.
184 inline static uint8_t loPenaltyAt(size_t off, int alts,
185 const String<char>& qual,
186 const String<char>* altQual)
188 uint8_t primQ = qual[off]; // qual of primary call
189 uint8_t bestPenalty = primQ - 33;
190 for(int i = 0; i < alts; i++) {
191 uint8_t altQ = altQual[i][off]; // qual of alt call
192 if(altQ == 33) break; // no more alt calls at this position
193 assert_leq(altQ, primQ);
194 if(primQ - altQ < bestPenalty) {
195 bestPenalty = primQ - altQ;