decalc.cpp

   1 /*
   2  *  decalc.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 7/22/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "decalc.h"
  11
  12 //***************************************************************************************************************
  13 void DeCalculator::setMask(string m) {
  14         try {
  15                 seqMask = m;
  16
  17                 if (seqMask.length() != 0) {
  18                         //whereever there is a base in the mask, save that value is query and subject
  19                         for (int i = 0; i < seqMask.length(); i++) {
  20                                 if (isalpha(seqMask[i])) {
  21                                         h.insert(i);
  22                                 }
  23                         }
  24                 }else {
  25                         for (int i = 0; i < alignLength; i++) {   h.insert(i);  }
  26                 }
  27         }
  28         catch(exception& e) {
  29                 errorOut(e, "DeCalculator", "setMask");
  30                 exit(1);
  31         }
  32 }
  33 //***************************************************************************************************************
  34 void DeCalculator::runMask(Sequence* seq) {
  35         try{
  36
  37                 string q = seq->getAligned();
  38                 string tempQuery = "";
  39
  40                 //whereever there is a base in the mask, save that value is query and subject
  41                 set<int>::iterator setit;
  42                 for ( setit=h.begin() ; setit != h.end(); setit++ )  {
  43                         tempQuery += q[*setit];
  44                 }
  45
  46                 //save masked values
  47                 seq->setAligned(tempQuery);
  48                 seq->setUnaligned(tempQuery);
  49         }
  50         catch(exception& e) {
  51                 errorOut(e, "DeCalculator", "runMask");
  52                 exit(1);
  53         }
  54 }
  55 //***************************************************************************************************************
  56 //num is query's spot in querySeqs
  57 void DeCalculator::trimSeqs(Sequence* query, Sequence* subject, map<int, int>& trim) {
  58         try {
  59
  60                 string q = query->getAligned();
  61                 string s = subject->getAligned();
  62
  63                 int front = 0;
  64                 for (int i = 0; i < q.length(); i++) {
  65 //cout << "query = " << q[i] << " subject = " << s[i] << endl;
  66                         if (isalpha(q[i]) && isalpha(s[i])) { front = i; break;  }
  67                 }
  68 //cout << endl << endl;
  69                 int back = 0;
  70                 for (int i = q.length(); i >= 0; i--) {
  71 //cout << "query = " << q[i] << " subject = " << s[i] << endl;
  72                         if (isalpha(q[i]) && isalpha(s[i])) { back = i; break;  }
  73                 }
  74
  75                 trim[front] = back;
  76
  77         }
  78         catch(exception& e) {
  79                 errorOut(e, "DeCalculator", "trimSeqs");
  80                 exit(1);
  81         }
  82 }
  83 //***************************************************************************************************************
  84 //find the window breaks for each sequence - this is so you can move ahead by bases.
  85 vector<int>  DeCalculator::findWindows(Sequence* query, int front, int back, int& size, int increment) {
  86         try {
  87
  88                 vector<int> win;
  89
  90                 int cutoff = back - front;  //back - front
  91
  92                 //if window is set to default
  93                 if (size == 0) {  if (cutoff > 1200) {  size = 300; }
  94                                                         else{  size = (cutoff / 4); }  }
  95                 else if (size > (cutoff / 4)) {
  96                                 mothurOut("You have selected to large a window size for sequence " + query->getName() + ".  I will choose an appropriate window size."); mothurOutEndLine();
  97                                 size = (cutoff / 4);
  98                 }
  99
 100         /*      string seq = query->getAligned().substr(front, cutoff);
 101
 102                 //count bases
 103                 int numBases = 0;
 104                 for (int l = 0; l < seq.length(); l++) {  if (isalpha(seq[l])) { numBases++; }  }
 105 //cout << "num Bases = " << numBases << endl;
 106                 //save start of seq
 107                 win.push_back(front);
 108 //cout << front << '\t';
 109                 //move ahead increment bases at a time until all bases are in a window
 110                 int countBases = 0;
 111                 int totalBases = 0;  //used to eliminate window of blanks at end of sequence
 112
 113                 seq = query->getAligned();
 114                 for (int m = front; m < (back - size) ; m++) {
 115
 116                         //count number of bases you see
 117                         if (isalpha(seq[m])) { countBases++;  }
 118
 119                         //if you have seen enough bases to make a new window
 120                         if (countBases >= increment) {
 121                                 //total bases is the number of bases in a window already.
 122                                 totalBases += countBases;
 123 //cout << "total bases = " << totalBases << endl;
 124                                 win.push_back(m);  //save spot in alignment
 125 //cout << m << '\t';
 126                                 countBases = 0;                         //reset bases you've seen in this window
 127                         }
 128
 129                         //no need to continue if all your bases are in a window
 130                         if (totalBases == numBases) {   break;  }
 131                 }
 132
 133
 134                 //get last window if needed
 135                 if (totalBases < numBases) {   win.push_back(back-size);  }
 136 //cout << endl << endl;
 137 */
 138
 139                 //this follows wigeon, but we may want to consider that it chops off the end values if the sequence cannot be evenly divided into steps
 140                 for (int m = front; m < (back - size) ; m+=increment) {  win.push_back(m);  }
 141
 142
 143
 144                 return win;
 145
 146         }
 147         catch(exception& e) {
 148                 errorOut(e, "DeCalculator", "findWindows");
 149                 exit(1);
 150         }
 151 }
 152
 153 //***************************************************************************************************************
 154 vector<float> DeCalculator::calcObserved(Sequence* query, Sequence* subject, vector<int> window, int size) {
 155         try {
 156
 157                 vector<float> temp;
 158 //cout << "query length = " << query->getAligned().length() << '\t' << " subject length = " << subject.getAligned().length() << endl;
 159                 for (int m = 0; m < window.size(); m++) {
 160
 161                         string seqFrag = query->getAligned().substr(window[m], size);
 162                         string seqFragsub = subject->getAligned().substr(window[m], size);
 163         //cout << "start point = " << window[m] << " end point = " << window[m]+size << endl;
 164                         int diff = 0;
 165                         for (int b = 0; b < seqFrag.length(); b++) {
 166
 167                                 if (seqFrag[b] != seqFragsub[b]) { diff++; }
 168                         }
 169
 170                         //percentage of mismatched bases
 171                         float dist;
 172                         dist = diff / (float) seqFrag.length() * 100;
 173
 174                         temp.push_back(dist);
 175                 }
 176
 177                 return temp;
 178         }
 179         catch(exception& e) {
 180                 errorOut(e, "DeCalculator", "calcObserved");
 181                 exit(1);
 182         }
 183 }
 184 //***************************************************************************************************************
 185 float DeCalculator::calcDist(Sequence* query, Sequence* subject, int front, int back) {
 186         try {
 187
 188                 //so you only look at the trimmed part of the sequence
 189                 int cutoff = back - front;
 190                 int gaps = 0;
 191
 192                 //from first startpoint with length back-front
 193                 string seqFrag = query->getAligned().substr(front, cutoff);
 194                 string seqFragsub = subject->getAligned().substr(front, cutoff);
 195
 196                 int diff = 0;
 197                 for (int b = 0; b < seqFrag.length(); b++) {
 198                         //ignore gaps
 199                         if((!isalpha(seqFrag[b])) && (!isalpha(seqFragsub[b]))) { gaps++; }
 200                         if (seqFrag[b] != seqFragsub[b]) { diff++; }
 201                 }
 202
 203                 //percentage of mismatched bases
 204                 float dist = diff / (float) (seqFrag.length()-gaps) * 100;
 205
 206                 return dist;
 207         }
 208         catch(exception& e) {
 209                 errorOut(e, "DeCalculator", "calcDist");
 210                 exit(1);
 211         }
 212 }
 213
 214 //***************************************************************************************************************
 215 vector<float> DeCalculator::calcExpected(vector<float> qav, float coef) {
 216         try {
 217
 218                 //for each window
 219                 vector<float> queryExpected;
 220
 221                 for (int m = 0; m < qav.size(); m++) {
 222
 223                         float expected = qav[m] * coef;
 224
 225                         queryExpected.push_back(expected);
 226                 }
 227
 228                 return queryExpected;
 229
 230         }
 231         catch(exception& e) {
 232                 errorOut(e, "DeCalculator", "calcExpected");
 233                 exit(1);
 234         }
 235 }
 236 //***************************************************************************************************************
 237 float DeCalculator::calcDE(vector<float> obs, vector<float> exp) {
 238         try {
 239
 240                 //for each window
 241                 float sum = 0.0;  //sum = sum from 1 to m of (oi-ei)^2
 242                 for (int m = 0; m < obs.size(); m++) {          sum += ((obs[m] - exp[m]) * (obs[m] - exp[m]));         }
 243
 244                 float de = sqrt((sum / (obs.size() - 1)));
 245
 246                 return de;
 247         }
 248         catch(exception& e) {
 249                 errorOut(e, "DeCalculator", "calcDE");
 250                 exit(1);
 251         }
 252 }
 253
 254 //***************************************************************************************************************
 255
 256 vector<float> DeCalculator::calcFreq(vector<Sequence*> seqs, string filename) {
 257         try {
 258
 259                 vector<float> prob;
 260                 string freqfile = getRootName(filename) + "freq";
 261                 ofstream outFreq;
 262
 263                 openOutputFile(freqfile, outFreq);
 264
 265                 string length = toString(seqs.size());  //if there are 5000 seqs in the template then set precision to 3
 266                 int precision = length.length() - 1;
 267
 268                 //format output
 269                 outFreq.setf(ios::fixed, ios::floatfield); outFreq.setf(ios::showpoint);
 270
 271                 //at each position in the sequence
 272                 for (int i = 0; i < seqs[0]->getAligned().length(); i++) {
 273
 274                         vector<int> freq;   freq.resize(4,0);
 275                         int gaps = 0;
 276
 277                         //find the frequency of each nucleotide
 278                         for (int j = 0; j < seqs.size(); j++) {
 279
 280                                 char value = seqs[j]->getAligned()[i];
 281
 282                                 if(toupper(value) == 'A')                                                                       {       freq[0]++;      }
 283                                 else if(toupper(value) == 'T' || toupper(value) == 'U')         {       freq[1]++;      }
 284                                 else if(toupper(value) == 'G')                                                          {       freq[2]++;      }
 285                                 else if(toupper(value) == 'C')                                                          {       freq[3]++;      }
 286                                 else { gaps++; }
 287                         }
 288
 289                         //find base with highest frequency
 290                         int highest = 0;
 291                         for (int m = 0; m < freq.size(); m++) {   if (freq[m] > highest) {  highest = freq[m];  }               }
 292
 293                         float highFreq = highest / (float) (seqs.size());
 294
 295                         float Pi;
 296                         Pi =  (highFreq - 0.25) / 0.75;
 297
 298                         //cannot have probability less than 0.
 299                         if (Pi < 0) { Pi = 0.0; }
 300
 301                         //saves this for later
 302                         outFreq << setprecision(precision) << i << '\t' << highFreq << endl;
 303
 304                         if (h.count(i) > 0) {
 305                                 prob.push_back(Pi);
 306                         }
 307                 }
 308
 309                 outFreq.close();
 310
 311                 return prob;
 312
 313         }
 314         catch(exception& e) {
 315                 errorOut(e, "DeCalculator", "calcFreq");
 316                 exit(1);
 317         }
 318 }
 319 //***************************************************************************************************************
 320 vector<float>  DeCalculator::findQav(vector<int> window, int size, vector<float> probabilityProfile) {
 321         try {
 322                 vector<float>  averages;
 323
 324                 //for each window find average
 325                 for (int m = 0; m < window.size(); m++) {
 326
 327                         float average = 0.0;
 328
 329                         //while you are in the window for this sequence
 330                         int count = 0;
 331                         for (int j = window[m]; j < (window[m]+size); j++) {
 332                                 average += probabilityProfile[j];
 333                                 count++;
 334                         }
 335
 336                         average = average / count;
 337
 338                         //save this windows average
 339                         averages.push_back(average);
 340                 }
 341
 342                 return averages;
 343         }
 344         catch(exception& e) {
 345                 errorOut(e, "DeCalculator", "findQav");
 346                 exit(1);
 347         }
 348 }
 349
 350 //***************************************************************************************************************
 351 vector< vector<float> > DeCalculator::getQuantiles(vector<Sequence*> seqs, vector<int> windowSizesTemplate, int window, vector<float> probProfile, int increment, int start, int end) {
 352         try {
 353                 vector< vector<float> > quan;
 354
 355                 //percentage of mismatched pairs 1 to 100
 356                 quan.resize(100);
 357
 358
 359                 //for each sequence
 360                 for(int i = start; i < end; i++){
 361
 362                         mothurOut("Processing template sequence " + toString(i)); mothurOutEndLine();
 363                         Sequence* query = seqs[i];
 364
 365                         //compare to every other sequence in template
 366                         for(int j = 0; j < i; j++){
 367
 368                                 Sequence* subject = seqs[j];
 369
 370                                 map<int, int> trim;
 371                                 map<int, int>::iterator it;
 372
 373                                 trimSeqs(query, subject, trim);
 374
 375                                 it = trim.begin();
 376                                 int front = it->first; int back = it->second;
 377
 378                                 //reset window for each new comparison
 379                                 windowSizesTemplate[i] = window;
 380
 381                                 vector<int> win = findWindows(query, front, back, windowSizesTemplate[i], increment);
 382
 383                                 vector<float> obsi = calcObserved(query, subject, win, windowSizesTemplate[i]);
 384
 385                                 vector<float> q = findQav(win, windowSizesTemplate[i], probProfile);
 386
 387                                 float alpha = getCoef(obsi, q);
 388
 389                                 vector<float> exp = calcExpected(q, alpha);
 390
 391                                 float de = calcDE(obsi, exp);
 392
 393                                 float dist = calcDist(query, subject, front, back);
 394
 395                                 dist = ceil(dist);
 396
 397                                 //dist-1 because vector indexes start at 0.
 398                                 quan[dist-1].push_back(de);
 399
 400                         }
 401                 }
 402
 403                 return quan;
 404
 405         }
 406         catch(exception& e) {
 407                 errorOut(e, "DeCalculator", "findQav");
 408                 exit(1);
 409         }
 410 }
 411 //***************************************************************************************************************
 412 void DeCalculator::removeObviousOutliers(vector< vector<float> >& quantiles) {
 413         try {
 414
 415
 416                 for (int i = 0; i < quantiles.size(); i++) {
 417
 418                         //find mean of this quantile score
 419                         sort(quantiles[i].begin(), quantiles[i].end());
 420
 421                         float average = quantiles[i][int(quantiles[i].size() * 0.5)];
 422 cout << i << "\taverage = " << average << "\tquantiles[i].size = " << quantiles[i].size() << endl;
 423                         vector<float> newQuanI;
 424                         //look at each value in quantiles to see if it is an outlier
 425                         for (int j = 0; j < quantiles[i].size(); j++) {
 426
 427                                 float highCutOff, lowCutOff;
 428
 429                                 //99%
 430                                 highCutOff = sqrt(((quantiles[i][j] - average + 3) * (quantiles[i][j] - average + 3)) / (float)(quantiles[i].size() - 1));
 431
 432                                 //1%
 433                                 lowCutOff = sqrt(((quantiles[i][j] - average - 3) * (quantiles[i][j] - average + 3)) / (float)(quantiles[i].size() - 1));
 434 //cout << "high = " << highCutOff << " low = " << lowCutOff << " de = " << quantiles[i][j] << endl;
 435                                 //if this is below the highcutff and above the lowcutoff
 436                                 if ((quantiles[i][j] < highCutOff) && (quantiles[i][j] > lowCutOff)) {
 437
 438                                         newQuanI.push_back(quantiles[i][j]);
 439
 440                                 }else { cout << "removed outlier:  high = " << highCutOff << " low = " << lowCutOff << " de = " << quantiles[i][j] << endl;      }
 441                         }
 442
 443                         quantiles[i] = newQuanI;
 444
 445                 }
 446
 447         }
 448         catch(exception& e) {
 449                 errorOut(e, "DeCalculator", "removeObviousOutliers");
 450                 exit(1);
 451         }
 452 }
 453
 454 //***************************************************************************************************************
 455 float DeCalculator::findAverage(vector<float> myVector) {
 456         try{
 457
 458                 float total = 0.0;
 459                 for (int i = 0; i < myVector.size(); i++)  {  total += myVector[i];  }
 460
 461                 float average = total / (float) myVector.size();
 462
 463                 return average;
 464
 465         }
 466         catch(exception& e) {
 467                 errorOut(e, "DeCalculator", "findAverage");
 468                 exit(1);
 469         }
 470 }
 471
 472 //***************************************************************************************************************
 473 float DeCalculator::getCoef(vector<float> obs, vector<float> qav) {
 474         try {
 475
 476                 //find average prob for this seqs windows
 477                 float probAverage = findAverage(qav);
 478
 479                 //find observed average
 480                 float obsAverage = findAverage(obs);
 481
 482                 float coef = obsAverage / probAverage;
 483
 484                 return coef;
 485         }
 486         catch(exception& e) {
 487                 errorOut(e, "DeCalculator", "getCoef");
 488                 exit(1);
 489         }
 490 }
 491 //***************************************************************************************************************
 492
 493