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 vector<int>  DeCalculator::findWindows(Sequence* query, int front, int back, int& size, int increment) {
  85         try {
  86
  87                 vector<int> win;
  88
  89                 int cutoff = back - front;  //back - front
  90
  91                 //if window is set to default
  92                 if (size == 0) {  if (cutoff > 1200) {  size = 300; }
  93                                                         else{  size = (cutoff / 4); }  }
  94                 else if (size > (cutoff / 4)) {
  95                                 mothurOut("You have selected to large a window size for sequence " + query->getName() + ".  I will choose an appropriate window size."); mothurOutEndLine();
  96                                 size = (cutoff / 4);
  97                 }
  98
  99         /*      string seq = query->getAligned().substr(front, cutoff);
 100
 101                 //count bases
 102                 int numBases = 0;
 103                 for (int l = 0; l < seq.length(); l++) {  if (isalpha(seq[l])) { numBases++; }  }
 104 //cout << "num Bases = " << numBases << endl;
 105                 //save start of seq
 106                 win.push_back(front);
 107 //cout << front << '\t';
 108                 //move ahead increment bases at a time until all bases are in a window
 109                 int countBases = 0;
 110                 int totalBases = 0;  //used to eliminate window of blanks at end of sequence
 111
 112                 seq = query->getAligned();
 113                 for (int m = front; m < (back - size) ; m++) {
 114
 115                         //count number of bases you see
 116                         if (isalpha(seq[m])) { countBases++;  }
 117
 118                         //if you have seen enough bases to make a new window
 119                         if (countBases >= increment) {
 120                                 //total bases is the number of bases in a window already.
 121                                 totalBases += countBases;
 122 //cout << "total bases = " << totalBases << endl;
 123                                 win.push_back(m);  //save spot in alignment
 124 //cout << m << '\t';
 125                                 countBases = 0;                         //reset bases you've seen in this window
 126                         }
 127
 128                         //no need to continue if all your bases are in a window
 129                         if (totalBases == numBases) {   break;  }
 130                 }
 131
 132
 133                 //get last window if needed
 134                 if (totalBases < numBases) {   win.push_back(back-size);  }
 135 //cout << endl << endl;
 136 */
 137
 138                 //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
 139                 for (int m = front; m < (back - size) ; m+=increment) {  win.push_back(m);  }
 140
 141
 142
 143                 return win;
 144
 145         }
 146         catch(exception& e) {
 147                 errorOut(e, "DeCalculator", "findWindows");
 148                 exit(1);
 149         }
 150 }
 151
 152 //***************************************************************************************************************
 153 vector<float> DeCalculator::calcObserved(Sequence* query, Sequence* subject, vector<int> window, int size) {
 154         try {
 155
 156                 vector<float> temp;
 157
 158                 for (int m = 0; m < window.size(); m++) {
 159
 160                         string seqFrag = query->getAligned().substr(window[m], size);
 161                         string seqFragsub = subject->getAligned().substr(window[m], size);
 162
 163                         int diff = 0;
 164                         for (int b = 0; b < seqFrag.length(); b++) {
 165                                 //if at least one is a base and they are not equal
 166                                 if( (isalpha(seqFrag[b]) || isalpha(seqFragsub[b])) && (seqFrag[b] != seqFragsub[b]) ) { diff++; }
 167                         }
 168
 169                         //percentage of mismatched bases
 170                         float dist;
 171
 172                         //if the whole fragment is 0 distance = 0
 173
 174                         dist = diff / (float) (seqFrag.length()) * 100;
 175
 176                         temp.push_back(dist);
 177                 }
 178
 179                 return temp;
 180         }
 181         catch(exception& e) {
 182                 errorOut(e, "DeCalculator", "calcObserved");
 183                 exit(1);
 184         }
 185 }
 186 //***************************************************************************************************************
 187 float DeCalculator::calcDist(Sequence* query, Sequence* subject, int front, int back) {
 188         try {
 189
 190                 //so you only look at the trimmed part of the sequence
 191                 int cutoff = back - front;
 192                 int gaps = 0;
 193
 194                 //from first startpoint with length back-front
 195                 string seqFrag = query->getAligned().substr(front, cutoff);
 196                 string seqFragsub = subject->getAligned().substr(front, cutoff);
 197
 198                 int diff = 0;
 199                 for (int b = 0; b < seqFrag.length(); b++) {
 200                         //ignore gaps
 201                         if((!isalpha(seqFrag[b])) && (!isalpha(seqFragsub[b]))) { gaps++; }
 202                         if (seqFrag[b] != seqFragsub[b]) { diff++; }
 203                 }
 204
 205                 //if the whole fragment is 0 distance = 0
 206                 if ((seqFrag.length()-gaps) == 0)  { return 0.0; }
 207
 208                 //percentage of mismatched bases
 209                 float dist = diff / (float) (seqFrag.length()-gaps) * 100;
 210
 211                 return dist;
 212         }
 213         catch(exception& e) {
 214                 errorOut(e, "DeCalculator", "calcDist");
 215                 exit(1);
 216         }
 217 }
 218
 219 //***************************************************************************************************************
 220 vector<float> DeCalculator::calcExpected(vector<float> qav, float coef) {
 221         try {
 222
 223                 //for each window
 224                 vector<float> queryExpected;
 225
 226                 for (int m = 0; m < qav.size(); m++) {
 227
 228                         float expected = qav[m] * coef;
 229
 230                         queryExpected.push_back(expected);
 231                 }
 232
 233                 return queryExpected;
 234
 235         }
 236         catch(exception& e) {
 237                 errorOut(e, "DeCalculator", "calcExpected");
 238                 exit(1);
 239         }
 240 }
 241 //***************************************************************************************************************
 242 float DeCalculator::calcDE(vector<float> obs, vector<float> exp) {
 243         try {
 244
 245                 //for each window
 246                 float sum = 0.0;  //sum = sum from 1 to m of (oi-ei)^2
 247                 for (int m = 0; m < obs.size(); m++) {          sum += ((obs[m] - exp[m]) * (obs[m] - exp[m]));         }
 248
 249                 float de = sqrt((sum / (obs.size() - 1)));
 250
 251                 return de;
 252         }
 253         catch(exception& e) {
 254                 errorOut(e, "DeCalculator", "calcDE");
 255                 exit(1);
 256         }
 257 }
 258
 259 //***************************************************************************************************************
 260
 261 vector<float> DeCalculator::calcFreq(vector<Sequence*> seqs, string filename) {
 262         try {
 263
 264                 vector<float> prob;
 265                 string freqfile = getRootName(filename) + "freq";
 266                 ofstream outFreq;
 267
 268                 openOutputFile(freqfile, outFreq);
 269
 270                 string length = toString(seqs.size());  //if there are 5000 seqs in the template then set precision to 3
 271                 int precision = length.length() - 1;
 272
 273                 //format output
 274                 outFreq.setf(ios::fixed, ios::floatfield); outFreq.setf(ios::showpoint);
 275
 276                 //at each position in the sequence
 277                 for (int i = 0; i < seqs[0]->getAligned().length(); i++) {
 278
 279                         vector<int> freq;   freq.resize(4,0);
 280                         int gaps = 0;
 281
 282                         //find the frequency of each nucleotide
 283                         for (int j = 0; j < seqs.size(); j++) {
 284
 285                                 char value = seqs[j]->getAligned()[i];
 286
 287                                 if(toupper(value) == 'A')                                                                       {       freq[0]++;      }
 288                                 else if(toupper(value) == 'T' || toupper(value) == 'U')         {       freq[1]++;      }
 289                                 else if(toupper(value) == 'G')                                                          {       freq[2]++;      }
 290                                 else if(toupper(value) == 'C')                                                          {       freq[3]++;      }
 291                                 else { gaps++; }
 292                         }
 293
 294                         //find base with highest frequency
 295                         int highest = 0;
 296                         for (int m = 0; m < freq.size(); m++) {   if (freq[m] > highest) {  highest = freq[m];  }               }
 297
 298                         float highFreq = highest / (float) (seqs.size());
 299
 300                         float Pi;
 301                         Pi =  (highFreq - 0.25) / 0.75;
 302
 303                         //cannot have probability less than 0.
 304                         if (Pi < 0) { Pi = 0.0; }
 305
 306                         //saves this for later
 307                         outFreq << setprecision(precision) << i << '\t' << highFreq << endl;
 308
 309                         if (h.count(i) > 0) {
 310                                 prob.push_back(Pi);
 311                         }
 312                 }
 313
 314                 outFreq.close();
 315
 316                 return prob;
 317
 318         }
 319         catch(exception& e) {
 320                 errorOut(e, "DeCalculator", "calcFreq");
 321                 exit(1);
 322         }
 323 }
 324 //***************************************************************************************************************
 325 vector<float>  DeCalculator::findQav(vector<int> window, int size, vector<float> probabilityProfile) {
 326         try {
 327                 vector<float>  averages;
 328
 329                 //for each window find average
 330                 for (int m = 0; m < window.size(); m++) {
 331
 332                         float average = 0.0;
 333
 334                         //while you are in the window for this sequence
 335                         int count = 0;
 336                         for (int j = window[m]; j < (window[m]+size); j++) {
 337                                 average += probabilityProfile[j];
 338                                 count++;
 339                         }
 340
 341                         average = average / count;
 342
 343                         //save this windows average
 344                         averages.push_back(average);
 345                 }
 346
 347                 return averages;
 348         }
 349         catch(exception& e) {
 350                 errorOut(e, "DeCalculator", "findQav");
 351                 exit(1);
 352         }
 353 }
 354 //********************************************************************************************************************
 355 //sorts lowest to highest
 356 inline bool compareQuanMembers(quanMember left, quanMember right){
 357         return (left.score < right.score);
 358 }
 359 //***************************************************************************************************************
 360 //seqs have already been masked
 361 vector< vector<quanMember> > DeCalculator::getQuantiles(vector<Sequence*> seqs, vector<int> windowSizesTemplate, int window, vector<float> probProfile, int increment, int start, int end, vector<float>& highestDE) {
 362         try {
 363                 vector< vector<quanMember> > quan;
 364
 365                 //percentage of mismatched pairs 1 to 100
 366                 quan.resize(100);
 367
 368                 //for each sequence
 369                 for(int i = start; i < end; i++){
 370
 371                         mothurOut("Processing template sequence " + toString(i)); mothurOutEndLine();
 372                         Sequence* query = new Sequence(seqs[i]->getName(), seqs[i]->getAligned());
 373
 374                         //compare to every other sequence in template
 375                         for(int j = 0; j < i; j++){
 376
 377                                 Sequence* subject = new Sequence(seqs[j]->getName(), seqs[j]->getAligned());
 378
 379                                 map<int, int> trim;
 380                                 map<int, int>::iterator it;
 381
 382                                 trimSeqs(query, subject, trim);
 383
 384                                 it = trim.begin();
 385                                 int front = it->first; int back = it->second;
 386
 387                                 //reset window for each new comparison
 388                                 windowSizesTemplate[i] = window;
 389
 390                                 vector<int> win = findWindows(query, front, back, windowSizesTemplate[i], increment);
 391
 392                                 vector<float> obsi = calcObserved(query, subject, win, windowSizesTemplate[i]);
 393
 394                                 vector<float> q = findQav(win, windowSizesTemplate[i], probProfile);
 395
 396                                 float alpha = getCoef(obsi, q);
 397
 398                                 vector<float> exp = calcExpected(q, alpha);
 399
 400                                 float de = calcDE(obsi, exp);
 401
 402                                 float dist = calcDist(query, subject, front, back);
 403
 404                                 dist = ceil(dist);
 405
 406                                 quanMember newScore(de, i, j);
 407
 408                                 //dist-1 because vector indexes start at 0.
 409                                 quan[dist-1].push_back(newScore);
 410
 411                                 //save highestDE
 412                                 if(de > highestDE[i]) { highestDE[i] = de;  }
 413                                 if(de > highestDE[j]) { highestDE[j] = de;  }
 414
 415                                 delete subject;
 416
 417                         }
 418
 419                         delete query;
 420                 }
 421
 422
 423                 return quan;
 424
 425         }
 426         catch(exception& e) {
 427                 errorOut(e, "DeCalculator", "findQav");
 428                 exit(1);
 429         }
 430 }
 431
 432 //***************************************************************************************************************
 433 vector< vector<float> > DeCalculator::removeObviousOutliers(vector< vector<quanMember> >& quantiles, int num) {
 434         try {
 435                 vector< vector<float> > quan;
 436                 quan.resize(100);
 437
 438                 /*vector<quanMember> contributions;
 439                 vector<int> seen;  //seen[0] is the number of outliers that template seqs[0] was part of.
 440                 seen.resize(num,0);
 441
 442                 //find contributions
 443                 for (int i = 0; i < quantiles.size(); i++) {
 444
 445                         //find mean of this quantile score
 446                         sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
 447
 448                         float high = quantiles[i][int(quantiles[i].size() * 0.99)].score;
 449                         float low =  quantiles[i][int(quantiles[i].size() * 0.01)].score;
 450
 451                         //look at each value in quantiles to see if it is an outlier
 452                         for (int j = 0; j < quantiles[i].size(); j++) {
 453
 454                                 //is this score between 1 and 99%
 455                                 if ((quantiles[i][j].score > low) && (quantiles[i][j].score < high)) {
 456
 457                                 }else {
 458                                         //add to contributions
 459                                         contributions.push_back(quantiles[i][j]);
 460                                         seen[quantiles[i][j].member1]++;
 461                                         seen[quantiles[i][j].member2]++;
 462                                 }
 463                         }
 464                 }
 465
 466                 //find contributer with most offending score related to it
 467                 int largestContrib = findLargestContrib(seen);
 468
 469                 //while you still have guys to eliminate
 470                 while (contributions.size() > 0) {
 471
 472                         mothurOut("Removing scores contributed by sequence " + toString(largestContrib) + " in your template file."); mothurOutEndLine();
 473
 474                         //remove from quantiles all scores that were made using this largestContrib
 475                         for (int i = 0; i < quantiles.size(); i++) {
 476 //cout << "before remove " << quantiles[i].size() << '\t';
 477                                 removeContrib(largestContrib, quantiles[i]);
 478 //cout << "after remove " << quantiles[i].size() << endl;
 479                         }
 480 //cout << count << " largest contrib = " << largestContrib << endl;  count++;
 481                         //remove from contributions all scores that were made using this largestContrib
 482                         removeContrib(largestContrib, contributions);
 483
 484                         //"erase" largestContrib
 485                         seen[largestContrib] = -1;
 486
 487                         //get next largestContrib
 488                         largestContrib = findLargestContrib(seen);
 489                 }
 490                 */
 491
 492                 vector<int> marked = returnObviousOutliers(quantiles, num);
 493                 vector<int> copyMarked = marked;
 494
 495                 //find the 99% of marked
 496                 sort(copyMarked.begin(), copyMarked.end());
 497                 int high = copyMarked[int(copyMarked.size() * 0.99)];
 498 cout << "high = " << high << endl;
 499
 500                 for(int i = 0; i < marked.size(); i++) {
 501                         if (marked[i] > high) {
 502                                 mothurOut("Removing scores contributed by sequence " + toString(marked[i]) + " in your template file."); mothurOutEndLine();
 503                                 for (int i = 0; i < quantiles.size(); i++) {
 504                                         removeContrib(marked[i], quantiles[i]);
 505                                 }
 506                         }
 507
 508                 }
 509
 510
 511                 //adjust quantiles
 512                 for (int i = 0; i < quantiles.size(); i++) {
 513                         vector<float> temp;
 514
 515                         if (quantiles[i].size() == 0) {
 516                                 //in case this is not a distance found in your template files
 517                                 for (int g = 0; g < 6; g++) {
 518                                         temp.push_back(0.0);
 519                                 }
 520                         }else{
 521
 522                                 sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
 523
 524                                 //save 10%
 525                                 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.10)].score);
 526                                 //save 25%
 527                                 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.25)].score);
 528                                 //save 50%
 529                                 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.5)].score);
 530                                 //save 75%
 531                                 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.75)].score);
 532                                 //save 95%
 533                                 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.95)].score);
 534                                 //save 99%
 535                                 temp.push_back(quantiles[i][int(quantiles[i].size() * 0.99)].score);
 536
 537                         }
 538
 539                         quan[i] = temp;
 540
 541                 }
 542
 543                 return quan;
 544         }
 545         catch(exception& e) {
 546                 errorOut(e, "DeCalculator", "removeObviousOutliers");
 547                 exit(1);
 548         }
 549 }
 550 //***************************************************************************************************************
 551 //follows Mallard algorythn in paper referenced from mallard class
 552 vector<int> DeCalculator::returnObviousOutliers(vector< vector<quanMember> > quantiles, int num) {
 553         try {
 554                 vector< vector<float> > quan;
 555                 quan.resize(100);
 556
 557                 map<quanMember*, float> contributions;  //map of quanMember to distance from high or low - how bad is it.
 558                 vector<int> marked;  //marked[0] is the penalty of template seqs[0]. the higher the penalty the more likely the sequence is chimeric
 559                 marked.resize(num,0);
 560
 561                 //find contributions
 562                 for (int i = 0; i < quantiles.size(); i++) {
 563
 564                         //find mean of this quantile score
 565                         sort(quantiles[i].begin(), quantiles[i].end(), compareQuanMembers);
 566
 567                         float high = quantiles[i][int(quantiles[i].size() * 0.99)].score;
 568                         float low =  quantiles[i][int(quantiles[i].size() * 0.01)].score;
 569
 570                         //look at each value in quantiles to see if it is an outlier
 571                         for (int j = 0; j < quantiles[i].size(); j++) {
 572
 573                                 //is this score between 1 and 99%
 574                                 if ((quantiles[i][j].score > low) && (quantiles[i][j].score < high)) {
 575
 576                                 }else {
 577                                         float dist;
 578                                         //add to contributions
 579                                         if (quantiles[i][j].score < low) {
 580                                                 dist = low - quantiles[i][j].score;
 581                                                 contributions[&(quantiles[i][j])] = dist;
 582                                         }else { //you are higher than high
 583                                                 dist = quantiles[i][j].score - high;
 584                                                 contributions[&(quantiles[i][j])] = dist;
 585                                         }
 586                                 }
 587                         }
 588                 }
 589
 590                 //find contributer with most offending score related to it
 591                 vector<quanMember> outliers = sortContrib(contributions);
 592
 593                 //go through the outliers marking the potential chimeras
 594                 for (int i = 0; i < outliers.size(); i++) {
 595
 596                         //who is responsible for this outlying score?
 597                         //if member1 has greater score mark him
 598                         //if member2 has greater score mark her
 599                         //if they are the same mark both
 600                         if (marked[outliers[i].member1] > marked[outliers[i].member2])                  {       marked[outliers[i].member1]++;  }
 601                         else if (marked[outliers[i].member2] > marked[outliers[i].member1])             {       marked[outliers[i].member2]++;  }
 602                         else if (marked[outliers[i].member2] == marked[outliers[i].member1])    {       marked[outliers[i].member2]++;  marked[outliers[i].member1]++;  }
 603                 }
 604
 605                 return marked;
 606         }
 607         catch(exception& e) {
 608                 errorOut(e, "DeCalculator", "removeObviousOutliers");
 609                 exit(1);
 610         }
 611 }
 612 //***************************************************************************************************************
 613 vector<quanMember> DeCalculator::sortContrib(map<quanMember*, float> quan) {
 614         try{
 615
 616                 vector<quanMember> newQuan;
 617                 map<quanMember*, float>::iterator it;
 618
 619                 while (quan.size() > 0) {
 620
 621                          map<quanMember*, float>::iterator largest = quan.begin();
 622
 623                         //find biggest member
 624                         for (it = quan.begin(); it != quan.end(); it++) {
 625                                 if (it->second > largest->second) {  largest = it;  }
 626                         }
 627
 628                         //save it
 629                         newQuan.push_back(*(largest->first));
 630
 631                         //erase from quan
 632                         quan.erase(largest);
 633                 }
 634
 635                 return newQuan;
 636
 637         }
 638         catch(exception& e) {
 639                 errorOut(e, "DeCalculator", "sortContrib");
 640                 exit(1);
 641         }
 642 }
 643
 644 //***************************************************************************************************************
 645 int DeCalculator::findLargestContrib(vector<int> seen) {
 646         try{
 647
 648                 int largest = 0;
 649
 650                 int largestContribs;
 651
 652                 for (int i = 0; i < seen.size(); i++)  {
 653
 654                         if (seen[i] > largest) {
 655                                 largestContribs = i;
 656                                 largest = seen[i];
 657                         }
 658                 }
 659
 660                 return largestContribs;
 661
 662         }
 663         catch(exception& e) {
 664                 errorOut(e, "DeCalculator", "findLargestContribs");
 665                 exit(1);
 666         }
 667 }
 668 //***************************************************************************************************************
 669 void DeCalculator::removeContrib(int bad, vector<quanMember>& quan) {
 670         try{
 671
 672                 vector<quanMember> newQuan;
 673                 for (int i = 0; i < quan.size(); i++)  {
 674                         if ((quan[i].member1 != bad) && (quan[i].member2 != bad) ) {
 675                                 newQuan.push_back(quan[i]);
 676                         }
 677                 }
 678
 679                 quan = newQuan;
 680
 681         }
 682         catch(exception& e) {
 683                 errorOut(e, "DeCalculator", "removeContrib");
 684                 exit(1);
 685         }
 686 }
 687
 688 //***************************************************************************************************************
 689 float DeCalculator::findAverage(vector<float> myVector) {
 690         try{
 691
 692                 float total = 0.0;
 693                 for (int i = 0; i < myVector.size(); i++)  {  total += myVector[i];  }
 694
 695                 float average = total / (float) myVector.size();
 696
 697                 return average;
 698
 699         }
 700         catch(exception& e) {
 701                 errorOut(e, "DeCalculator", "findAverage");
 702                 exit(1);
 703         }
 704 }
 705
 706 //***************************************************************************************************************
 707 float DeCalculator::getCoef(vector<float> obs, vector<float> qav) {
 708         try {
 709
 710                 //find average prob for this seqs windows
 711                 float probAverage = findAverage(qav);
 712
 713                 //find observed average
 714                 float obsAverage = findAverage(obs);
 715
 716                 float coef = obsAverage / probAverage;
 717
 718                 return coef;
 719         }
 720         catch(exception& e) {
 721                 errorOut(e, "DeCalculator", "getCoef");
 722                 exit(1);
 723         }
 724 }
 725 //***************************************************************************************************************
 726
 727