qual.h

   1 #ifndef QUAL_H_
   2 #define QUAL_H_
   3
   4 #include <stdexcept>
   5
   6 extern unsigned char qualRounds[];
   7 extern unsigned char solToPhred[];
   8
   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);
  12 }
  13
  14 /**
  15  * Convert a Solexa-scaled quality value into a Phred-scale quality
  16  * value.
  17  *
  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
  22  *
  23  */
  24 static inline uint8_t solexaToPhred(int sol) {
  25         assert_lt(sol, 256);
  26         if(sol < -10) return 0;
  27         return solToPhred[sol+10];
  28 }
  29
  30 class SimplePhredPenalty {
  31 public:
  32         static uint8_t mmPenalty (uint8_t qual) {
  33                 return qual;
  34         }
  35         static uint8_t delPenalty(uint8_t qual) {
  36                 return qual;
  37         }
  38         static uint8_t insPenalty(uint8_t qual_left, uint8_t qual_right) {
  39                 return std::max(qual_left, qual_right);
  40         }
  41 };
  42
  43 class MaqPhredPenalty {
  44 public:
  45         static uint8_t mmPenalty (uint8_t qual) {
  46                 return qualRounds[qual];
  47         }
  48         static uint8_t delPenalty(uint8_t qual) {
  49                 return qualRounds[qual];
  50         }
  51         static uint8_t insPenalty(uint8_t qual_left, uint8_t qual_right) {
  52                 return qualRounds[std::max(qual_left, qual_right)];
  53         }
  54 };
  55
  56 static inline uint8_t mmPenalty(bool maq, uint8_t qual) {
  57         if(maq) {
  58                 return MaqPhredPenalty::mmPenalty(qual);
  59         } else {
  60                 return SimplePhredPenalty::mmPenalty(qual);
  61         }
  62 }
  63
  64 static inline uint8_t delPenalty(bool maq, uint8_t qual) {
  65         if(maq) {
  66                 return MaqPhredPenalty::delPenalty(qual);
  67         } else {
  68                 return SimplePhredPenalty::delPenalty(qual);
  69         }
  70 }
  71
  72 static inline uint8_t insPenalty(bool maq, uint8_t qual_left, uint8_t qual_right) {
  73         if(maq) {
  74                 return MaqPhredPenalty::insPenalty(qual_left, qual_right);
  75         } else {
  76                 return SimplePhredPenalty::insPenalty(qual_left, qual_right);
  77         }
  78 }
  79
  80 /**
  81  * Take an ASCII-encoded quality value and convert it to a Phred33
  82  * ASCII char.
  83  */
  84 inline static char charToPhred33(char c, bool solQuals, bool phred64Quals) {
  85         if(c == ' ') {
  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;
  88                 throw 1;
  89         }
  90         if (solQuals) {
  91                 // Convert solexa-scaled chars to phred
  92                 // http://maq.sourceforge.net/fastq.shtml
  93                 char cc = solexaToPhred((int)c - 64) + 33;
  94                 if (cc < 33) {
  95                         cerr << "Saw ASCII character "
  96                              << ((int)c)
  97                              << " but expected 64-based Solexa qual (converts to " << (int)cc << ")." << endl
  98                              << "Try not specifying --solexa-quals." << endl;
  99                         throw 1;
 100                 }
 101                 c = cc;
 102         }
 103         else if(phred64Quals) {
 104                 if (c < 64) {
 105                         cerr << "Saw ASCII character "
 106                              << ((int)c)
 107                              << " but expected 64-based Phred qual." << endl
 108                              << "Try not specifying --solexa1.3-quals/--phred64-quals." << endl;
 109                         throw 1;
 110                 }
 111                 // Convert to 33-based phred
 112                 c -= (64-33);
 113         }
 114         else {
 115                 // Keep the phred quality
 116                 if (c < 33) {
 117                         cerr << "Saw ASCII character "
 118                              << ((int)c)
 119                              << " but expected 33-based Phred qual." << endl;
 120                         throw 1;
 121                 }
 122         }
 123         return c;
 124 }
 125
 126 /**
 127  * Take an integer quality value and convert it to a Phred33 ASCII
 128  * char.
 129  */
 130 inline static char intToPhred33(int iQ, bool solQuals) {
 131         int pQ;
 132         if (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;
 137         } else {
 138                 // Keep the phred quality and translate
 139                 // to ASCII
 140                 pQ = (iQ <= 93 ? iQ : 93) + 33;
 141         }
 142         if (pQ < 33) {
 143                 cerr << "Saw negative Phred quality " << ((int)pQ-33) << "." << endl;
 144                 throw 1;
 145         }
 146         assert_geq(pQ, 0);
 147         return (int)pQ;
 148 }
 149
 150 /**
 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.
 154  */
 155 inline static uint8_t penaltiesAt(size_t off, uint8_t *q,
 156                                   int alts,
 157                                   const String<char>& qual,
 158                                   const String<Dna5>* altQry,
 159                                   const String<char>* altQual)
 160 {
 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;
 170                 }
 171                 // Get the base
 172                 int altC = (int)altQry[i][off];
 173                 assert_lt(altC, 4);
 174                 q[altC] = primQ - altQ;
 175         }
 176         return bestPenalty;
 177 }
 178
 179 /**
 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.
 183  */
 184 inline static uint8_t loPenaltyAt(size_t off, int alts,
 185                                   const String<char>& qual,
 186                                   const String<char>* altQual)
 187 {
 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;
 196                 }
 197         }
 198         return bestPenalty;
 199 }
 200
 201 #endif /*QUAL_H_*/