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                 int score;
  35                 seqName = getSequenceName(qFile); m->gobble(qFile);
  36
  37         if (m->debug) { m->mothurOut("[DEBUG]: name = '" + seqName + "'\n.");  }
  38
  39                 if (!m->control_pressed) {
  40             string qScoreString = m->getline(qFile); m->gobble(qFile);
  41
  42             if (m->debug) { m->mothurOut("[DEBUG]: scores = '" + qScoreString + "'\n.");  }
  43
  44             while(qFile.peek() != '>' && qFile.peek() != EOF){
  45                 if (m->control_pressed) { break; }
  46                 string temp = m->getline(qFile); m->gobble(qFile);
  47                 //if (m->debug) { m->mothurOut("[DEBUG]: scores = '" + temp + "'\n.");  }
  48                 qScoreString +=  ' ' + temp;
  49             }
  50             //cout << "done reading " << endl;
  51             istringstream qScoreStringStream(qScoreString);
  52             int count = 0;
  53             while(!qScoreStringStream.eof()){
  54                 if (m->control_pressed) { break; }
  55                 string temp;
  56                 qScoreStringStream >> temp;  m->gobble(qScoreStringStream);
  57
  58                 //if (m->debug) { m->mothurOut("[DEBUG]: score " + toString(qScores.size()) + " = '" + temp + "'\n.");  }
  59
  60                 //check temp to make sure its a number
  61                 if (!m->isContainingOnlyDigits(temp)) { m->mothurOut("[ERROR]: In sequence " + seqName + "'s quality scores, expected a number and got " + temp + ", setting score to 0."); m->mothurOutEndLine(); temp = "0"; }
  62                 convert(temp, score);
  63
  64                 //cout << count << '\t' << score << endl;
  65                 qScores.push_back(score);
  66                 count++;
  67             }
  68         }
  69
  70                 seqLength = qScores.size();
  71                 //cout << "seqlength = " << seqLength  << endl;
  72
  73         }
  74         catch(exception& e) {
  75                 m->errorOut(e, "QualityScores", "QualityScores");
  76                 exit(1);
  77         }
  78
  79 }
  80 //********************************************************************************************************************
  81 string QualityScores::getSequenceName(ifstream& qFile) {
  82         try {
  83                 string name = "";
  84
  85         qFile >> name;
  86         m->getline(qFile);
  87
  88                 if (name.length() != 0) {
  89
  90                         name = name.substr(1);
  91
  92             m->checkName(name);
  93
  94         }else{ m->mothurOut("Error in reading your qfile, at position " + toString(qFile.tellg()) + ". Blank name."); m->mothurOutEndLine(); m->control_pressed = true;  }
  95
  96                 return name;
  97         }
  98         catch(exception& e) {
  99                 m->errorOut(e, "QualityScores", "getSequenceName");
 100                 exit(1);
 101         }
 102 }
 103 //********************************************************************************************************************
 104 void QualityScores::setName(string name) {
 105         try {
 106
 107         m->checkName(name);
 108         seqName = name;
 109         }
 110         catch(exception& e) {
 111                 m->errorOut(e, "QualityScores", "setName");
 112                 exit(1);
 113         }
 114 }
 115 /**************************************************************************************************/
 116
 117 string QualityScores::getName(){
 118
 119         try {
 120                 return seqName;
 121         }
 122         catch(exception& e) {
 123                 m->errorOut(e, "QualityScores", "getName");
 124                 exit(1);
 125         }
 126 }
 127
 128 /**************************************************************************************************/
 129
 130 void QualityScores::printQScores(ofstream& qFile){
 131         try {
 132
 133                 double aveQScore = calculateAverage();
 134
 135                 qFile << '>' << seqName << '\t' << aveQScore << endl;
 136
 137                 for(int i=0;i<seqLength;i++){
 138                         qFile << qScores[i] << ' ';
 139                 }
 140                 qFile << endl;
 141
 142         }
 143         catch(exception& e) {
 144                 m->errorOut(e, "QualityScores", "printQScores");
 145                 exit(1);
 146         }
 147 }
 148
 149 /**************************************************************************************************/
 150
 151 void QualityScores::trimQScores(int start, int end){
 152         try {
 153                 vector<int> hold;
 154
 155
 156                 //cout << seqName << '\t' << start << '\t' << end << '\t' << qScores.size() << endl;
 157                 //for (int i = 0; i < qScores.size(); i++) { cout << qScores[i] << end; }
 158                 if(end == -1){
 159                         hold = vector<int>(qScores.begin()+start, qScores.end());
 160                         qScores = hold;
 161                 }
 162                 if(start == -1){
 163                         if(qScores.size() > end){
 164                                 hold = vector<int>(qScores.begin(), qScores.begin()+end);
 165                                 qScores = hold;
 166                         }
 167                 }
 168
 169                 seqLength = qScores.size();
 170         }
 171         catch(exception& e) {
 172                 m->errorOut(e, "QualityScores", "trimQScores");
 173                 exit(1);
 174         }
 175 }
 176
 177 /**************************************************************************************************/
 178
 179 void QualityScores::flipQScores(){
 180         try {
 181
 182                 vector<int> temp = qScores;
 183                 for(int i=0;i<seqLength;i++){
 184                         qScores[seqLength - i - 1] = temp[i];
 185                 }
 186
 187         }
 188         catch(exception& e) {
 189                 m->errorOut(e, "QualityScores", "flipQScores");
 190                 exit(1);
 191         }
 192 }
 193
 194 /**************************************************************************************************/
 195
 196 bool QualityScores::stripQualThreshold(Sequence& sequence, double qThreshold){
 197         try {
 198                 string rawSequence = sequence.getUnaligned();
 199                 int seqLength = sequence.getNumBases();
 200
 201                 if(seqName != sequence.getName()){
 202                         m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
 203                         m->mothurOutEndLine();  m->control_pressed = true;
 204                 }
 205
 206                 int end;
 207                 for(int i=0;i<seqLength;i++){
 208                         end = i;
 209                         if(qScores[i] < qThreshold){
 210                                 break;
 211                         }
 212                 }
 213
 214                 //every score passed
 215                 if (end == (seqLength-1)) { end = seqLength; }
 216
 217                 sequence.setUnaligned(rawSequence.substr(0,end));
 218                 trimQScores(-1, end);
 219
 220                 return 1;
 221         }
 222         catch(exception& e) {
 223                 m->errorOut(e, "QualityScores", "flipQScores");
 224                 exit(1);
 225         }
 226
 227 }
 228
 229 /**************************************************************************************************/
 230
 231 bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
 232         try {
 233                 string rawSequence = sequence.getUnaligned();
 234                 int seqLength = sequence.getNumBases();
 235
 236                 if(seqName != sequence.getName()){
 237                         m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
 238                         m->mothurOutEndLine();
 239                 }
 240
 241                 int end = -1;
 242                 double rollingSum = 0.0000;
 243
 244                 for(int i=0;i<seqLength;i++){
 245
 246                         rollingSum += (double)qScores[i];
 247
 248                         if(rollingSum / (double)(i+1) < qThreshold){
 249                                 end = i;
 250                                 break;
 251                         }
 252                 }
 253
 254                 if(end == -1){  end = seqLength;        }
 255
 256
 257                 sequence.setUnaligned(rawSequence.substr(0,end));
 258                 trimQScores(-1, end);
 259
 260
 261                 return 1;
 262         }
 263         catch(exception& e) {
 264                 m->errorOut(e, "QualityScores", "flipQScores");
 265                 exit(1);
 266         }
 267
 268 }
 269
 270 /**************************************************************************************************/
 271
 272 bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
 273         try {
 274                 string rawSequence = sequence.getUnaligned();
 275                 int seqLength = sequence.getNumBases();
 276
 277                 if(seqName != sequence.getName()){
 278                         m->mothurOut("sequence name mismatch between fasta: " + sequence.getName() + " and qual file: " + seqName);
 279                         m->mothurOutEndLine();
 280                 }
 281
 282                 int end = windowSize;
 283                 int start = 0;
 284
 285                 if(seqLength < windowSize) {    return 0;       }
 286
 287                 while((start+windowSize) < seqLength){
 288                         double windowSum = 0.0000;
 289
 290                         for(int i=start;i<end;i++){
 291                                 windowSum += qScores[i];
 292                         }
 293                         double windowAverage = windowSum / (double)(end-start);
 294
 295                         if(windowAverage < qThreshold){
 296                                 end = end - stepSize;
 297                                 break;
 298                         }
 299
 300                         start += stepSize;
 301                         end = start + windowSize;
 302
 303                         if(end >= seqLength){   end = seqLength;        }
 304
 305                 }
 306
 307                 if(end == -1){  end = seqLength;        }
 308
 309                 //failed first window
 310                 if (end < windowSize) { return 0; }
 311
 312                 sequence.setUnaligned(rawSequence.substr(0,end));
 313                 trimQScores(-1, end);
 314
 315                 return 1;
 316         }
 317         catch(exception& e) {
 318                 m->errorOut(e, "QualityScores", "stripQualWindowAverage");
 319                 exit(1);
 320         }
 321
 322 }
 323
 324 /**************************************************************************************************/
 325
 326 double QualityScores::calculateAverage(){
 327
 328         double aveQScore = 0.0000;
 329
 330         for(int i=0;i<seqLength;i++){
 331                 aveQScore += (double) qScores[i];
 332         }
 333         aveQScore /= (double) seqLength;
 334
 335         return aveQScore;
 336 }
 337
 338 /**************************************************************************************************/
 339
 340 bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
 341         try {
 342                 string rawSequence = sequence.getUnaligned();
 343                 bool success = 0;       //guilty until proven innocent
 344
 345                 if(seqName != sequence.getName())       {
 346                         m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
 347                         m->mothurOutEndLine();
 348                 }
 349
 350                 double aveQScore = calculateAverage();
 351
 352                 if(aveQScore >= qAverage)       {       success = 1;    }
 353                 else                                            {       success = 0;    }
 354
 355                 return success;
 356         }
 357         catch(exception& e) {
 358                 m->errorOut(e, "QualityScores", "cullQualAverage");
 359                 exit(1);
 360         }
 361 }
 362
 363 /**************************************************************************************************/
 364
 365 void QualityScores::updateQScoreErrorMap(map<char, vector<int> >& qualErrorMap, string errorSeq, int start, int stop, int weight){
 366         try {
 367
 368                 int seqLength = errorSeq.size();
 369
 370                 int qIndex = start - 1;
 371
 372                 for(int i=0;i<seqLength;i++){
 373
 374                         if(errorSeq[i] == 'm')          {       qualErrorMap['m'][qScores[qIndex]] += weight;   }
 375                         else if(errorSeq[i] == 's')     {       qualErrorMap['s'][qScores[qIndex]] += weight;   }
 376                         else if(errorSeq[i] == 'i')     {       qualErrorMap['i'][qScores[qIndex]] += weight;   }
 377                         else if(errorSeq[i] == 'a')     {       qualErrorMap['a'][qScores[qIndex]] += weight;   /*if(qScores[qIndex] != 0){     cout << qIndex << '\t';         }*/     }
 378                         else if(errorSeq[i] == 'd')     {       /*      there are no qScores for deletions      */              }
 379
 380                         if(errorSeq[i] != 'd')          {       qIndex++;       }
 381
 382                         if(qIndex > stop){      break;  }
 383                 }
 384         }
 385         catch(exception& e) {
 386                 m->errorOut(e, "QualityScores", "updateQScoreErrorMap");
 387                 exit(1);
 388         }
 389 }
 390
 391 /**************************************************************************************************/
 392
 393 void QualityScores::updateForwardMap(vector<vector<int> >& forwardMap, int start, int stop, int weight){
 394         try {
 395
 396                 int index = 0;
 397                 for(int i=start-1;i<stop;i++){
 398                         forwardMap[index++][qScores[i]] += weight;
 399                 }
 400
 401         }
 402         catch(exception& e) {
 403                 m->errorOut(e, "QualityScores", "updateForwardMap");
 404                 exit(1);
 405         }
 406 }
 407
 408 /**************************************************************************************************/
 409
 410 void QualityScores::updateReverseMap(vector<vector<int> >& reverseMap, int start, int stop, int weight){
 411         try {
 412
 413                 int index = 0;
 414                 for(int i=stop-1;i>=start-1;i--){
 415                         reverseMap[index++][qScores[i]] += weight;
 416                 }
 417
 418         }
 419         catch(exception& e) {
 420                 m->errorOut(e, "QualityScores", "updateReverseMap");
 421                 exit(1);
 422         }
 423 }
 424
 425 /**************************************************************************************************/