qualityscores.cpp

   1 /*
   2  *  qualityscores.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 7/12/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "qualityscores.h"
  11
  12 /**************************************************************************************************/
  13
  14 QualityScores::QualityScores(){
  15         try {
  16                 m = MothurOut::getInstance();
  17                 seqName = "";
  18                 seqLength = -1;
  19
  20         }
  21         catch(exception& e) {
  22                 m->errorOut(e, "QualityScores", "QualityScores");
  23                 exit(1);
  24         }
  25 }
  26
  27 /**************************************************************************************************/
  28
  29 QualityScores::QualityScores(ifstream& qFile){
  30         try {
  31
  32                 m = MothurOut::getInstance();
  33
  34                 seqName = "";
  35                 int score;
  36
  37                 qFile >> seqName;
  38                 m->getline(qFile);
  39
  40                 if (seqName == "")      {
  41                         m->mothurOut("Error reading quality file, name blank at position, " + toString(qFile.tellg()));
  42                         m->mothurOutEndLine();
  43                 }
  44                 else{
  45                         seqName = seqName.substr(1);
  46                 }
  47
  48                 string qScoreString = m->getline(qFile);
  49
  50                 istringstream qScoreStringStream(qScoreString);
  51                 while(!qScoreStringStream.eof()){
  52                         qScoreStringStream >> score;
  53                         qScores.push_back(score);
  54                 }
  55                 qScores.pop_back();
  56                 seqLength = qScores.size();
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "QualityScores", "QualityScores");
  60                 exit(1);
  61         }
  62
  63 }
  64
  65 /**************************************************************************************************/
  66
  67 string QualityScores::getName(){
  68
  69         try {
  70                 return seqName;
  71         }
  72         catch(exception& e) {
  73                 m->errorOut(e, "QualityScores", "getName");
  74                 exit(1);
  75         }
  76 }
  77
  78 /**************************************************************************************************/
  79
  80 void QualityScores::printQScores(ofstream& qFile){
  81         try {
  82
  83                 double aveQScore = calculateAverage();
  84
  85                 qFile << '>' << seqName << '\t' << aveQScore << endl;
  86
  87                 for(int i=0;i<seqLength;i++){
  88                         qFile << qScores[i] << ' ';
  89                 }
  90                 qFile << endl;
  91
  92         }
  93         catch(exception& e) {
  94                 m->errorOut(e, "QualityScores", "printQScores");
  95                 exit(1);
  96         }
  97 }
  98
  99 /**************************************************************************************************/
 100
 101 void QualityScores::trimQScores(int start, int end){
 102         try {
 103                 vector<int> hold;
 104
 105                 if(end == -1){
 106                         hold = vector<int>(qScores.begin()+start, qScores.end());
 107                         qScores = hold;
 108                 }
 109                 if(start == -1){
 110                         if(qScores.size() > end){
 111                                 hold = vector<int>(qScores.begin(), qScores.begin()+end);
 112                                 qScores = hold;
 113                         }
 114                 }
 115
 116                 seqLength = qScores.size();
 117         }
 118         catch(exception& e) {
 119                 m->errorOut(e, "QualityScores", "trimQScores");
 120                 exit(1);
 121         }
 122 }
 123
 124 /**************************************************************************************************/
 125
 126 void QualityScores::flipQScores(){
 127         try {
 128
 129                 vector<int> temp = qScores;
 130                 for(int i=0;i<seqLength;i++){
 131                         qScores[seqLength - i - 1] = temp[i];
 132                 }
 133
 134         }
 135         catch(exception& e) {
 136                 m->errorOut(e, "QualityScores", "flipQScores");
 137                 exit(1);
 138         }
 139 }
 140
 141 /**************************************************************************************************/
 142
 143 bool QualityScores::stripQualThreshold(Sequence& sequence, double qThreshold){
 144         try {
 145                 string rawSequence = sequence.getUnaligned();
 146                 int seqLength = sequence.getNumBases();
 147
 148                 if(seqName != sequence.getName()){
 149                         m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
 150                         m->mothurOutEndLine();
 151                 }
 152
 153                 int end;
 154                 for(int i=0;i<seqLength;i++){
 155                         end = i;
 156                         if(qScores[i] < qThreshold){
 157                                 break;
 158                         }
 159                 }
 160
 161                 //every score passed
 162                 if (end == (seqLength-1)) { end = seqLength; }
 163
 164                 sequence.setUnaligned(rawSequence.substr(0,end));
 165                 trimQScores(-1, end);
 166
 167                 return 1;
 168         }
 169         catch(exception& e) {
 170                 m->errorOut(e, "QualityScores", "flipQScores");
 171                 exit(1);
 172         }
 173
 174 }
 175
 176 /**************************************************************************************************/
 177
 178 bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
 179         try {
 180                 string rawSequence = sequence.getUnaligned();
 181                 int seqLength = sequence.getNumBases();
 182
 183                 if(seqName != sequence.getName()){
 184                         m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
 185                         m->mothurOutEndLine();
 186                 }
 187
 188                 int end = -1;
 189                 double rollingSum = 0.0000;
 190
 191                 for(int i=0;i<seqLength;i++){
 192
 193                         rollingSum += (double)qScores[i];
 194
 195                         if(rollingSum / (double)(i+1) < qThreshold){
 196                                 end = i;
 197                                 break;
 198                         }
 199                 }
 200
 201                 if(end == -1){  end = seqLength;        }
 202
 203
 204                 sequence.setUnaligned(rawSequence.substr(0,end));
 205                 trimQScores(-1, end);
 206
 207
 208                 return 1;
 209         }
 210         catch(exception& e) {
 211                 m->errorOut(e, "QualityScores", "flipQScores");
 212                 exit(1);
 213         }
 214
 215 }
 216
 217 /**************************************************************************************************/
 218
 219 bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
 220         try {
 221                 string rawSequence = sequence.getUnaligned();
 222                 int seqLength = sequence.getNumBases();
 223
 224                 if(seqName != sequence.getName()){
 225                         m->mothurOut("sequence name mismatch between fasta: " + sequence.getName() + " and qual file: " + seqName);
 226                         m->mothurOutEndLine();
 227                 }
 228
 229                 int end = windowSize;
 230                 int start = 0;
 231
 232                 if(seqLength < windowSize) {    return 0;       }
 233
 234                 while(start < seqLength){
 235                         double windowSum = 0.0000;
 236
 237                         for(int i=start;i<end;i++){
 238                                 windowSum += qScores[i];
 239                         }
 240                         double windowAverage = windowSum / (double)(end-start);
 241
 242                         if(windowAverage < qThreshold){
 243                                 end = end - stepSize;
 244                                 break;
 245                         }
 246                         start += stepSize;
 247                         end = start + windowSize;
 248                         if(end >= seqLength){   end = seqLength - 1;    }
 249                 }
 250
 251                 if(end == -1){  end = seqLength;        }
 252
 253                 sequence.setUnaligned(rawSequence.substr(0,end));
 254                 trimQScores(-1, end);
 255
 256                 return 1;
 257         }
 258         catch(exception& e) {
 259                 m->errorOut(e, "QualityScores", "flipQScores");
 260                 exit(1);
 261         }
 262
 263 }
 264
 265 /**************************************************************************************************/
 266
 267 double QualityScores::calculateAverage(){
 268
 269         double aveQScore = 0.0000;
 270
 271         for(int i=0;i<seqLength;i++){
 272                 aveQScore += (double) qScores[i];
 273         }
 274         aveQScore /= (double) seqLength;
 275
 276         return aveQScore;
 277 }
 278
 279 /**************************************************************************************************/
 280
 281 bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
 282         try {
 283                 string rawSequence = sequence.getUnaligned();
 284                 bool success = 0;       //guilty until proven innocent
 285
 286                 if(seqName != sequence.getName())       {
 287                         m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
 288                         m->mothurOutEndLine();
 289                 }
 290
 291                 double aveQScore = calculateAverage();
 292
 293                 if(aveQScore >= qAverage)       {       success = 1;    }
 294                 else                                            {       success = 0;    }
 295
 296                 return success;
 297         }
 298         catch(exception& e) {
 299                 m->errorOut(e, "QualityScores", "cullQualAverage");
 300                 exit(1);
 301         }
 302 }
 303
 304 /**************************************************************************************************/
 305
 306 void QualityScores::updateQScoreErrorMap(map<char, vector<int> >& qualErrorMap, string errorSeq, int start, int stop, int weight){
 307         try {
 308
 309                 int seqLength = errorSeq.size();
 310
 311                 int qIndex = start - 1;
 312
 313                 for(int i=0;i<seqLength;i++){
 314
 315                         if(errorSeq[i] == 'm')          {       qualErrorMap['m'][qScores[qIndex]] += weight;   }
 316                         else if(errorSeq[i] == 's')     {       qualErrorMap['s'][qScores[qIndex]] += weight;   }
 317                         else if(errorSeq[i] == 'i')     {       qualErrorMap['i'][qScores[qIndex]] += weight;   }
 318                         else if(errorSeq[i] == 'a')     {       qualErrorMap['a'][qScores[qIndex]] += weight;   }
 319                         else if(errorSeq[i] == 'd')     {       /*      there are no qScores for deletions      */              }
 320
 321                         if(errorSeq[i] != 'd')          {       qIndex++;       }
 322                         if(qIndex > stop){      break;  }
 323                 }
 324         }
 325         catch(exception& e) {
 326                 m->errorOut(e, "QualityScores", "updateQScoreErrorMap");
 327                 exit(1);
 328         }
 329 }
 330
 331 /**************************************************************************************************/
 332
 333 void QualityScores::updateForwardMap(vector<vector<int> >& forwardMap, int start, int stop, int weight){
 334         try {
 335
 336                 int index = 0;
 337                 for(int i=start-1;i<stop;i++){
 338                         forwardMap[index++][qScores[i]] += weight;
 339                 }
 340
 341         }
 342         catch(exception& e) {
 343                 m->errorOut(e, "QualityScores", "updateForwardMap");
 344                 exit(1);
 345         }
 346 }
 347
 348 /**************************************************************************************************/
 349
 350 void QualityScores::updateReverseMap(vector<vector<int> >& reverseMap, int start, int stop, int weight){
 351         try {
 352
 353                 int index = 0;
 354                 for(int i=stop-1;i>=start;i--){
 355                         reverseMap[index++][qScores[i]] += weight;
 356                 }
 357
 358         }
 359         catch(exception& e) {
 360                 m->errorOut(e, "QualityScores", "updateForwardMap");
 361                 exit(1);
 362         }
 363 }
 364
 365 /**************************************************************************************************/