Revision: 29 Committed: Tue Aug 2 21:24:54 2011 UTC (9 years, 2 months ago) by gpertea File size: 3832 byte(s) Log Message: ```adding tophat source work ```
Line File contents
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_*/