maligner.cpp

   1 /*
   2  *  maligner.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/23/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "maligner.h"
  11
  12 /***********************************************************************/
  13 Maligner::Maligner(vector<Sequence*> temp, int num, int match, int misMatch, float div, int minCov) :
  14                 db(temp), numWanted(num), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minCoverage(minCov) {}
  15 /***********************************************************************/
  16 string Maligner::getResults(Sequence* q) {
  17         try {
  18
  19                 //make copy so trimming doesn't destroy query from calling class - remember to deallocate
  20                 query = new Sequence(q->getName(), q->getAligned());
  21
  22                 string chimera;
  23
  24                 decalc = new DeCalculator();
  25
  26                 //find closest seqs to query in template - returns copies of seqs so trim does not destroy - remember to deallocate
  27                 refSeqs = decalc->findClosest(query, db, numWanted);
  28
  29                 refSeqs = minCoverageFilter(refSeqs);
  30
  31                 if (refSeqs.size() < 2)  {
  32                         for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
  33                         percentIdenticalQueryChimera = 0.0;
  34                         return "unknown";
  35                 }
  36
  37                 int chimeraPenalty = computeChimeraPenalty();
  38
  39                 //trims seqs to first non gap char in all seqs and last non gap char in all seqs
  40                 decalc->trimSeqs(query, refSeqs);
  41
  42                 vector<Sequence*> temp = refSeqs;
  43                 temp.push_back(query);
  44
  45                 verticalFilter(temp);
  46
  47                 vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
  48
  49                 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
  50
  51                 vector<score_struct> path = extractHighestPath(matrix);
  52
  53                 vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
  54
  55                 if (trace.size() > 1) {         chimera = "yes";        }
  56                 else { chimera = "no";  }
  57
  58                 int traceStart = path[0].col;
  59                 int traceEnd = path[path.size()-1].col;
  60
  61                 string queryInRange = query->getAligned();
  62                 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
  63
  64                 string chimeraSeq = constructChimericSeq(trace, refSeqs);
  65
  66                 percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
  67
  68                 delete decalc;
  69
  70                 //save output results
  71                 for (int i = 0; i < trace.size(); i++) {
  72                         int regionStart = trace[i].col;
  73                         int regionEnd = trace[i].oldCol;
  74                         int seqIndex = trace[i].row;
  75
  76                         results temp;
  77
  78                         temp.parent = refSeqs[seqIndex]->getName();
  79                         temp.regionStart = regionStart;
  80                         temp.regionEnd = regionEnd;
  81
  82                         string parentInRange = refSeqs[seqIndex]->getAligned();
  83                         parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
  84
  85                         temp.queryToParent = computePercentID(queryInRange, parentInRange);
  86                         temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
  87
  88                         string queryInRegion = query->getAligned();
  89                         queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
  90
  91                         string parentInRegion = refSeqs[seqIndex]->getAligned();
  92                         parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
  93
  94                         temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
  95
  96                         outputResults.push_back(temp);
  97                 }
  98
  99                 //free memory
 100                 delete query;
 101                 for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
 102
 103                 return chimera;
 104         }
 105         catch(exception& e) {
 106                 errorOut(e, "Maligner", "getResults");
 107                 exit(1);
 108         }
 109 }
 110 /***********************************************************************/
 111 //removes top matches that do not have minimum coverage with query.
 112 vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
 113         try {
 114                 vector<Sequence*> newRefs;
 115
 116                 string queryAligned = query->getAligned();
 117
 118                 for (int i = 0; i < ref.size(); i++) {
 119
 120                         string refAligned = ref[i]->getAligned();
 121
 122                         int numBases = 0;
 123                         int numCovered = 0;
 124
 125                         //calculate coverage
 126                         for (int j = 0; j < queryAligned.length(); j++) {
 127
 128                                 if (isalpha(queryAligned[j])) {
 129                                         numBases++;
 130
 131                                         if (isalpha(refAligned[j])) {
 132                                                 numCovered++;
 133                                         }
 134                                 }
 135                         }
 136
 137                         int coverage = ((numCovered/(float)numBases)*100);
 138
 139                         //if coverage above minimum
 140                         if (coverage > minCoverage) {
 141                                 newRefs.push_back(ref[i]);
 142                         }
 143                 }
 144
 145                 return newRefs;
 146         }
 147         catch(exception& e) {
 148                 errorOut(e, "Maligner", "minCoverageFilter");
 149                 exit(1);
 150         }
 151 }
 152 /***********************************************************************/
 153 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
 154 int Maligner::computeChimeraPenalty() {
 155         try {
 156
 157                 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
 158
 159                 int penalty = int(numAllowable + 1) * misMatchPenalty;
 160
 161                 return penalty;
 162
 163         }
 164         catch(exception& e) {
 165                 errorOut(e, "Maligner", "computeChimeraPenalty");
 166                 exit(1);
 167         }
 168 }
 169 /***********************************************************************/
 170 //this is a vertical filter
 171 void Maligner::verticalFilter(vector<Sequence*> seqs) {
 172         try {
 173                 vector<int> gaps;       gaps.resize(query->getAligned().length(), 0);
 174
 175                 string filterString = (string(query->getAligned().length(), '1'));
 176
 177                 //for each sequence
 178                 for (int i = 0; i < seqs.size(); i++) {
 179
 180                         string seqAligned = seqs[i]->getAligned();
 181
 182                         for (int j = 0; j < seqAligned.length(); j++) {
 183                                 //if this spot is a gap
 184                                 if ((seqAligned[j] == '-') || (seqAligned[j] == '.'))   {       gaps[j]++;      }
 185                         }
 186                 }
 187
 188                 //zero out spot where all sequences have blanks
 189                 int numColRemoved = 0;
 190                 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
 191                         if(gaps[i] == seqs.size())      {       filterString[i] = '0';  numColRemoved++;  }
 192                 }
 193
 194                 //for each sequence
 195                 for (int i = 0; i < seqs.size(); i++) {
 196
 197                         string seqAligned = seqs[i]->getAligned();
 198                         string newAligned = "";
 199
 200                         for (int j = 0; j < seqAligned.length(); j++) {
 201                                 //if this spot is not a gap
 202                                 if (filterString[j] == '1') { newAligned += seqAligned[j]; }
 203                         }
 204
 205                         seqs[i]->setAligned(newAligned);
 206                 }
 207
 208
 209         }
 210         catch(exception& e) {
 211                 errorOut(e, "Maligner", "verticalFilter");
 212                 exit(1);
 213         }
 214 }
 215 //***************************************************************************************************************
 216 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
 217         try{
 218
 219                 vector< vector<score_struct> > m; m.resize(rows);
 220
 221                 for (int i = 0; i < m.size(); i++) {
 222                         for (int j = 0; j < cols; j++) {
 223
 224                                 //initialize each cell
 225                                 score_struct temp;
 226                                 temp.prev = -1;
 227                                 temp.score = -9999999;
 228                                 temp.col = j;
 229                                 temp.row = i;
 230
 231                                 m[i].push_back(temp);
 232                         }
 233                 }
 234
 235                 return m;
 236         }
 237         catch(exception& e) {
 238                 errorOut(e, "Maligner", "buildScoreMatrix");
 239                 exit(1);
 240         }
 241 }
 242 //***************************************************************************************************************
 243 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& m, vector<Sequence*> seqs, int penalty) {
 244         try{
 245
 246                 //get matrix dimensions
 247                 int numCols = query->getAligned().length();
 248                 int numRows = seqs.size();
 249
 250                 //initialize first col
 251                 string queryAligned = query->getAligned();
 252                 for (int i = 0; i < numRows; i++) {
 253                         string subjectAligned = seqs[i]->getAligned();
 254
 255                         //are you both gaps?
 256                         if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
 257                                 m[i][0].score = 0;
 258                         }else if (queryAligned[0] == subjectAligned[0]) {
 259                                 m[i][0].score = matchScore;
 260                         }else{
 261                                 m[i][0].score = 0;
 262                         }
 263                 }
 264
 265                 //fill rest of matrix
 266                 for (int j = 1; j < numCols; j++) {  //iterate through matrix columns
 267
 268                         for (int i = 0; i < numRows; i++) {  //iterate through matrix rows
 269
 270                                 string subjectAligned = seqs[i]->getAligned();
 271
 272                                 int matchMisMatchScore = 0;
 273                                 //are you both gaps?
 274                                 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
 275                                         //leave the same
 276                                 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
 277                                         //leave the same
 278                                 }else if (queryAligned[j] == subjectAligned[j]) {
 279                                         matchMisMatchScore = matchScore;
 280                                 }else if (queryAligned[j] != subjectAligned[j]) {
 281                                         matchMisMatchScore = misMatchPenalty;
 282                                 }
 283
 284                                 //compute score based on previous columns scores
 285                                 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
 286
 287                                         int sumScore = matchMisMatchScore + m[prevIndex][j-1].score;
 288
 289                                         //you are not at yourself
 290                                         if (prevIndex != i) {   sumScore += penalty;    }
 291                                         if (sumScore < 0)       {       sumScore = 0;                   }
 292
 293                                         if (sumScore > m[i][j].score) {
 294                                                 m[i][j].score = sumScore;
 295                                                 m[i][j].prev = prevIndex;
 296                                         }
 297                                 }
 298                         }
 299                 }
 300
 301         }
 302         catch(exception& e) {
 303                 errorOut(e, "Maligner", "fillScoreMatrix");
 304                 exit(1);
 305         }
 306 }
 307 //***************************************************************************************************************
 308 vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > m) {
 309         try {
 310
 311                 //get matrix dimensions
 312                 int numCols = query->getAligned().length();
 313                 int numRows = m.size();
 314
 315
 316                 //find highest score scoring matrix
 317                 score_struct highestStruct;
 318                 int highestScore = 0;
 319
 320                 for (int i = 0; i < numRows; i++) {
 321                         for (int j = 0; j < numCols; j++) {
 322                                 if (m[i][j].score > highestScore) {
 323                                         highestScore = m[i][j].score;
 324                                         highestStruct = m[i][j];
 325                                 }
 326                         }
 327                 }
 328
 329                 vector<score_struct> path;
 330
 331                 int rowIndex = highestStruct.row;
 332                 int pos = highestStruct.col;
 333                 int score = highestStruct.score;
 334
 335                 while (pos >= 0 && score > 0) {
 336                         score_struct temp = m[rowIndex][pos];
 337                         score = temp.score;
 338
 339                         if (score > 0) {        path.push_back(temp);   }
 340
 341                         rowIndex = temp.prev;
 342                         pos--;
 343                 }
 344
 345                 reverse(path.begin(), path.end());
 346
 347                 return path;
 348
 349         }
 350         catch(exception& e) {
 351                 errorOut(e, "Maligner", "extractHighestPath");
 352                 exit(1);
 353         }
 354 }
 355 //***************************************************************************************************************
 356 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
 357         try {
 358                 vector<trace_struct> trace;
 359
 360                 int region_index = path[0].row;
 361                 int region_start = path[0].col;
 362
 363                 for (int i = 1; i < path.size(); i++) {
 364
 365                         int next_region_index = path[i].row;
 366
 367                         if (next_region_index != region_index) {
 368
 369                                 // add trace region
 370                                 int col_index = path[i].col;
 371                                 trace_struct temp;
 372                                 temp.col = region_start;
 373                                 temp.oldCol = col_index-1;
 374                                 temp.row = region_index;
 375
 376                                 trace.push_back(temp);
 377
 378                                 region_index = path[i].row;
 379                                 region_start = col_index;
 380                         }
 381                 }
 382
 383                 // get last one
 384                 trace_struct temp;
 385                 temp.col = region_start;
 386                 temp.oldCol = path[path.size()-1].col;
 387                 temp.row = region_index;
 388                 trace.push_back(temp);
 389
 390                 return trace;
 391
 392         }
 393         catch(exception& e) {
 394                 errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
 395                 exit(1);
 396         }
 397 }
 398 //***************************************************************************************************************
 399 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
 400         try {
 401                 string chimera = "";
 402
 403                 for (int i = 0; i < trace.size(); i++) {
 404                         string seqAlign = seqs[trace[i].row]->getAligned();
 405                         seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
 406                         chimera += seqAlign;
 407                 }
 408
 409                 return chimera;
 410         }
 411         catch(exception& e) {
 412                 errorOut(e, "Maligner", "constructChimericSeq");
 413                 exit(1);
 414         }
 415 }
 416 //***************************************************************************************************************
 417 float Maligner::computePercentID(string queryAlign, string chimera) {
 418         try {
 419
 420                 if (queryAlign.length() != chimera.length()) {
 421                         mothurOut("Error, alignment strings are of different lengths: "); mothurOutEndLine();
 422                         mothurOut(toString(queryAlign.length())); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine(); mothurOutEndLine();
 423                         mothurOut(toString(chimera.length())); mothurOutEndLine();
 424                         return -1.0;
 425                 }
 426
 427
 428                 int numBases = 0;
 429                 int numIdentical = 0;
 430
 431                 for (int i = 0; i < queryAlign.length(); i++) {
 432                         if ((isalpha(queryAlign[i])) || (isalpha(chimera[i])))  {
 433                                 numBases++;
 434                                 if (queryAlign[i] == chimera[i]) {
 435                                         numIdentical++;
 436                                 }
 437                         }
 438                 }
 439
 440                 if (numBases == 0) { return 0; }
 441
 442                 float percentIdentical = (numIdentical/(float)numBases) * 100;
 443
 444                 return percentIdentical;
 445
 446         }
 447         catch(exception& e) {
 448                 errorOut(e, "Maligner", "computePercentID");
 449                 exit(1);
 450         }
 451 }
 452 //***************************************************************************************************************
 453