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