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