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                         if (isalpha(q[i]) && isalpha(s[i])) { front = i; break;  }
  63                 }
  64
  65                 int back = 0;
  66                 for (int i = q.length(); i >= 0; i--) {
  67                         if (isalpha(q[i]) && isalpha(s[i])) { back = i; break;  }
  68                 }
  69
  70                 trim[front] = back;
  71
  72         }
  73         catch(exception& e) {
  74                 errorOut(e, "DeCalculator", "trimSeqs");
  75                 exit(1);
  76         }
  77 }
  78 //***************************************************************************************************************
  79 //find the window breaks for each sequence - this is so you can move ahead by bases.
  80 vector<int>  DeCalculator::findWindows(Sequence* query, int front, int back, int& size, int increment) {
  81         try {
  82
  83                 vector<int> win;
  84
  85                 int cutoff = back - front;  //back - front
  86
  87                 //if window is set to default
  88                 if (size == 0) {  if (cutoff > 1200) {  size = 300; }
  89                                                         else{  size = (cutoff / 4); }  }
  90                 else if (size > (cutoff / 4)) {
  91                                 mothurOut("You have selected to large a window size for sequence " + query->getName() + ".  I will choose an appropriate window size."); mothurOutEndLine();
  92                                 size = (cutoff / 4);
  93                 }
  94
  95                 string seq = query->getAligned().substr(front, cutoff);
  96
  97                 //count bases
  98                 int numBases = 0;
  99                 for (int l = 0; l < seq.length(); l++) {  if (isalpha(seq[l])) { numBases++; }  }
 100
 101                 //save start of seq
 102                 win.push_back(front);
 103
 104                 //move ahead increment bases at a time until all bases are in a window
 105                 int countBases = 0;
 106                 int totalBases = 0;  //used to eliminate window of blanks at end of sequence
 107
 108                 seq = query->getAligned();
 109                 for (int m = front; m < (back - size) ; m++) {
 110
 111                         //count number of bases you see
 112                         if (isalpha(seq[m])) { countBases++; totalBases++;  }
 113
 114                         //if you have seen enough bases to make a new window
 115                         if (countBases >= increment) {
 116                                 win.push_back(m);  //save spot in alignment
 117                                 countBases = 0;                         //reset bases you've seen in this window
 118                         }
 119
 120                         //no need to continue if all your bases are in a window
 121                         if (totalBases == numBases) {   break;  }
 122                 }
 123
 124                 return win;
 125
 126         }
 127         catch(exception& e) {
 128                 errorOut(e, "DeCalculator", "findWindows");
 129                 exit(1);
 130         }
 131 }
 132
 133 //***************************************************************************************************************
 134 vector<float> DeCalculator::calcObserved(Sequence* query, Sequence subject, vector<int> window, int size) {
 135         try {
 136
 137                 vector<float> temp;
 138 //cout << "query length = " << query->getAligned().length() << '\t' << " subject length = " << subject.getAligned().length() << endl;
 139                 for (int m = 0; m < window.size(); m++) {
 140
 141                         string seqFrag = query->getAligned().substr(window[m], size);
 142                         string seqFragsub = subject.getAligned().substr(window[m], size);
 143         //cout << "start point = " << window[m] << " end point = " << window[m]+size << endl;
 144                         int diff = 0;
 145                         for (int b = 0; b < seqFrag.length(); b++) {
 146
 147                                 if (seqFrag[b] != seqFragsub[b]) { diff++; }
 148                         }
 149
 150                         //percentage of mismatched bases
 151                         float dist;
 152                         dist = diff / (float) seqFrag.length() * 100;
 153
 154                         temp.push_back(dist);
 155                 }
 156
 157                 return temp;
 158         }
 159         catch(exception& e) {
 160                 errorOut(e, "DeCalculator", "calcObserved");
 161                 exit(1);
 162         }
 163 }
 164 //***************************************************************************************************************
 165 float DeCalculator::calcDist(Sequence* query, Sequence subject, int front, int back) {
 166         try {
 167
 168                 //so you only look at the trimmed part of the sequence
 169                 int cutoff = back - front;
 170
 171                 //from first startpoint with length back-front
 172                 string seqFrag = query->getAligned().substr(front, cutoff);
 173                 string seqFragsub = subject.getAligned().substr(front, cutoff);
 174
 175                 int diff = 0;
 176                 for (int b = 0; b < seqFrag.length(); b++) {
 177                         if (seqFrag[b] != seqFragsub[b]) { diff++; }
 178                 }
 179
 180                 //percentage of mismatched bases
 181                 float dist = diff / (float) seqFrag.length() * 100;
 182
 183                 return dist;
 184         }
 185         catch(exception& e) {
 186                 errorOut(e, "DeCalculator", "calcDist");
 187                 exit(1);
 188         }
 189 }
 190
 191 //***************************************************************************************************************
 192 vector<float> DeCalculator::calcExpected(vector<float> qav, float coef) {
 193         try {
 194
 195                 //for each window
 196                 vector<float> queryExpected;
 197
 198                 for (int m = 0; m < qav.size(); m++) {
 199
 200                         float expected = qav[m] * coef;
 201
 202                         queryExpected.push_back(expected);
 203                 }
 204
 205                 return queryExpected;
 206
 207         }
 208         catch(exception& e) {
 209                 errorOut(e, "DeCalculator", "calcExpected");
 210                 exit(1);
 211         }
 212 }
 213 //***************************************************************************************************************
 214 float DeCalculator::calcDE(vector<float> obs, vector<float> exp) {
 215         try {
 216
 217                 //for each window
 218                 float sum = 0.0;  //sum = sum from 1 to m of (oi-ei)^2
 219                 for (int m = 0; m < obs.size(); m++) {          sum += ((obs[m] - exp[m]) * (obs[m] - exp[m]));         }
 220
 221                 float de = sqrt((sum / (obs.size() - 1)));
 222
 223                 return de;
 224         }
 225         catch(exception& e) {
 226                 errorOut(e, "DeCalculator", "calcDE");
 227                 exit(1);
 228         }
 229 }
 230
 231 //***************************************************************************************************************
 232
 233 vector<float> DeCalculator::calcFreq(vector<Sequence*> seqs, string filename) {
 234         try {
 235
 236                 vector<float> prob;
 237                 string freqfile = getRootName(filename) + "prob";
 238                 ofstream outFreq;
 239
 240                 openOutputFile(freqfile, outFreq);
 241
 242                 //at each position in the sequence
 243                 for (int i = 0; i < seqs[0]->getAligned().length(); i++) {
 244
 245                         vector<int> freq;   freq.resize(4,0);
 246                         int gaps = 0;
 247
 248                         //find the frequency of each nucleotide
 249                         for (int j = 0; j < seqs.size(); j++) {
 250
 251                                 char value = seqs[j]->getAligned()[i];
 252
 253                                 if(toupper(value) == 'A')                                                                       {       freq[0]++;      }
 254                                 else if(toupper(value) == 'T' || toupper(value) == 'U')         {       freq[1]++;      }
 255                                 else if(toupper(value) == 'G')                                                          {       freq[2]++;      }
 256                                 else if(toupper(value) == 'C')                                                          {       freq[3]++;      }
 257                                 else { gaps++; }
 258                         }
 259
 260                         //find base with highest frequency
 261                         int highest = 0;
 262                         for (int m = 0; m < freq.size(); m++) {   if (freq[m] > highest) {  highest = freq[m];  }               }
 263
 264                         float highFreq;
 265                         //subtract gaps to "ignore them"
 266                         if ( (seqs.size() - gaps) == 0 ) {  highFreq = 1.0;  }
 267                         else { highFreq = highest / (float) (seqs.size() - gaps);        }
 268
 269                         float Pi;
 270                         Pi =  (highFreq - 0.25) / 0.75;
 271
 272                         //cannot have probability less than 0.
 273                         if (Pi < 0) { Pi = 0.0; }
 274
 275                         //saves this for later
 276                         outFreq << i+1 << '\t' << Pi << endl;
 277
 278                         prob.push_back(Pi);
 279                 }
 280
 281                 outFreq.close();
 282
 283                 return prob;
 284
 285         }
 286         catch(exception& e) {
 287                 errorOut(e, "DeCalculator", "calcFreq");
 288                 exit(1);
 289         }
 290 }
 291 //***************************************************************************************************************
 292 vector<float>  DeCalculator::findQav(vector<int> window, int size, vector<float> probabilityProfile) {
 293         try {
 294                 vector<float>  averages;
 295
 296                 //for each window find average
 297                 for (int m = 0; m < window.size(); m++) {
 298
 299                         float average = 0.0;
 300
 301                         //while you are in the window for this sequence
 302                         int count = 0;
 303                         for (int j = window[m]; j < (window[m]+size); j++) {
 304
 305                                 //is this a spot that is included in the mask
 306                                 if (h.count(j) > 0) {
 307                                         average += probabilityProfile[j];
 308                                         count++;
 309                                 }
 310                         }
 311
 312                         average = average / count;
 313
 314                         //save this windows average
 315                         averages.push_back(average);
 316                 }
 317
 318                 return averages;
 319         }
 320         catch(exception& e) {
 321                 errorOut(e, "DeCalculator", "findQav");
 322                 exit(1);
 323         }
 324 }
 325
 326 //***************************************************************************************************************
 327 vector< vector<float> > DeCalculator::getQuantiles(vector<Sequence*> seqs, vector<int> windowSizesTemplate, int window, vector<float> probProfile, int increment, int start, int end) {
 328         try {
 329                 vector< vector<float> > quan;
 330
 331                 //percentage of mismatched pairs 1 to 100
 332                 quan.resize(100);
 333
 334
 335                 //for each sequence
 336                 for(int i = start; i < end; i++){
 337
 338                         mothurOut("Processing template sequence " + toString(i)); mothurOutEndLine();
 339                         Sequence* query = seqs[i];
 340
 341                         //compare to every other sequence in template
 342                         for(int j = 0; j < i; j++){
 343
 344                                 Sequence subject = *(seqs[j]);
 345
 346                                 map<int, int> trim;
 347                                 map<int, int>::iterator it;
 348
 349                                 trimSeqs(query, subject, trim);
 350
 351                                 it = trim.begin();
 352                                 int front = it->first; int back = it->second;
 353
 354                                 //reset window for each new comparison
 355                                 windowSizesTemplate[i] = window;
 356
 357                                 vector<int> win = findWindows(query, front, back, windowSizesTemplate[i], increment);
 358
 359                                 vector<float> obsi = calcObserved(query, subject, win, windowSizesTemplate[i]);
 360
 361                                 vector<float> q = findQav(win, windowSizesTemplate[i], probProfile);
 362
 363                                 float alpha = getCoef(obsi, q);
 364
 365                                 vector<float> exp = calcExpected(q, alpha);
 366
 367                                 float de = calcDE(obsi, exp);
 368
 369                                 float dist = calcDist(query, subject, front, back);
 370
 371                                 dist = ceil(dist);
 372
 373                                 //dist-1 because vector indexes start at 0.
 374                                 quan[dist-1].push_back(de);
 375
 376                         }
 377                 }
 378
 379                 return quan;
 380
 381         }
 382         catch(exception& e) {
 383                 errorOut(e, "DeCalculator", "findQav");
 384                 exit(1);
 385         }
 386 }
 387
 388 //***************************************************************************************************************
 389 float DeCalculator::getCoef(vector<float> obs, vector<float> qav) {
 390         try {
 391
 392                 //find average prob for this seqs windows
 393                 float probAverage = 0.0;
 394                 for (int j = 0; j < qav.size(); j++) {   probAverage += qav[j]; }
 395                 probAverage = probAverage / (float) qav.size();
 396
 397                 //find observed average
 398                 float obsAverage = 0.0;
 399                 for (int j = 0; j < obs.size(); j++) {   obsAverage += obs[j];  }
 400                 obsAverage = obsAverage / (float) obs.size();
 401 //cout << "sum ai / m = " << probAverage << endl;
 402 //cout << "sum oi / m = " << obsAverage << endl;
 403                 float coef = obsAverage / probAverage;
 404
 405                 return coef;
 406         }
 407         catch(exception& e) {
 408                 errorOut(e, "DeCalculator", "getCoef");
 409                 exit(1);
 410         }
 411 }
 412 //***************************************************************************************************************
 413
 414