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 /***********************************************************************/ //int num, int match, int misMatch, , string mode, Database* dataLeft, Database* dataRight
  13 Maligner::Maligner(vector<Sequence*> temp, int match, int misMatch, float div, int ms, int minCov) : db(temp), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov) {
  14                         //numWanted(num),  , searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight)
  15
  16                         m = MothurOut::getInstance();
  17
  18 }
  19 /***********************************************************************/
  20 string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
  21         try {
  22
  23                 outputResults.clear();
  24
  25                 //make copy so trimming doesn't destroy query from calling class - remember to deallocate
  26                 query = new Sequence(q->getName(), q->getAligned());
  27
  28                 string chimera;
  29
  30                 //copy refSeqs so that filter does not effect original
  31                 for(int i = 0; i < db.size(); i++) {
  32                         Sequence* newSeq = new Sequence(db[i]->getName(), db[i]->getAligned());
  33                         refSeqs.push_back(newSeq);
  34                 }
  35
  36                 refSeqs = minCoverageFilter(refSeqs);
  37
  38                 if (refSeqs.size() < 2)  {
  39                         for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
  40                         percentIdenticalQueryChimera = 0.0;
  41                         return "unknown";
  42                 }
  43
  44                 int chimeraPenalty = computeChimeraPenalty();
  45
  46                 //fills outputResults
  47                 chimera = chimeraMaligner(chimeraPenalty, decalc);
  48
  49                 if (m->control_pressed) { return chimera;  }
  50
  51                 //free memory
  52                 delete query;
  53
  54                 for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
  55
  56                 return chimera;
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "Maligner", "getResults");
  60                 exit(1);
  61         }
  62 }
  63 /***********************************************************************/
  64 string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
  65         try {
  66
  67                 string chimera;
  68
  69                 //trims seqs to first non gap char in all seqs and last non gap char in all seqs
  70                 spotMap = decalc->trimSeqs(query, refSeqs);
  71
  72                 //you trimmed the whole sequence, skip
  73                 if (query->getAligned() == "") { return "no"; }
  74
  75                 vector<Sequence*> temp = refSeqs;
  76                 temp.push_back(query);
  77
  78                 verticalFilter(temp);
  79
  80                 //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl << refSeqs[i]->getAligned() << endl; }
  81
  82                 vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
  83
  84                 if (m->control_pressed) { return chimera;  }
  85
  86                 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
  87
  88                 vector<trace_struct> trace = extractHighestPath(matrix);
  89
  90                 //cout << "traces\n";
  91                 //for(int i=0;i<trace.size();i++){
  92                 //      cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
  93                 //}
  94
  95                 if (trace.size() > 1) {         chimera = "yes";        }
  96                 else { chimera = "no";  return chimera; }
  97
  98                 int traceStart = trace[0].col;
  99                 int traceEnd = trace[trace.size()-1].oldCol;
 100                 string queryInRange = query->getAligned();
 101                 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
 102
 103                 if (m->control_pressed) { return chimera;  }
 104
 105                 //save output results
 106                 for (int i = 0; i < trace.size(); i++) {
 107                         int regionStart = trace[i].col;
 108                         int regionEnd = trace[i].oldCol;
 109                         int seqIndex = trace[i].row;
 110
 111                         results temp;
 112
 113                         temp.parent = refSeqs[seqIndex]->getName();
 114                         temp.parentAligned = db[seqIndex]->getAligned();
 115                         temp.nastRegionStart = spotMap[regionStart];
 116                         temp.nastRegionEnd = spotMap[regionEnd];
 117                         temp.regionStart = regionStart;
 118                         temp.regionEnd = regionEnd;
 119
 120                         string parentInRange = refSeqs[seqIndex]->getAligned();
 121                         parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
 122
 123                         temp.queryToParent = computePercentID(queryInRange, parentInRange);
 124                         temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
 125
 126                         string queryInRegion = query->getAligned();
 127                         queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
 128
 129                         string parentInRegion = refSeqs[seqIndex]->getAligned();
 130                         parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
 131
 132                         temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
 133
 134                         //cout << temp.parent << '\t' << "NAST:" << temp.nastRegionStart << '-' << temp.nastRegionEnd << " G:" << temp.queryToParent << " L:" << temp.queryToParentLocal << endl;
 135                         outputResults.push_back(temp);
 136                 }
 137
 138                 return chimera;
 139         }
 140         catch(exception& e) {
 141                 m->errorOut(e, "Maligner", "chimeraMaligner");
 142                 exit(1);
 143         }
 144 }
 145 /***********************************************************************/
 146 //removes top matches that do not have minimum coverage with query.
 147 vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
 148         try {
 149                 vector<Sequence*> newRefs;
 150
 151                 string queryAligned = query->getAligned();
 152
 153                 for (int i = 0; i < ref.size(); i++) {
 154
 155                         string refAligned = ref[i]->getAligned();
 156
 157                         int numBases = 0;
 158                         int numCovered = 0;
 159
 160                         //calculate coverage
 161                         for (int j = 0; j < queryAligned.length(); j++) {
 162
 163                                 if (isalpha(queryAligned[j])) {
 164                                         numBases++;
 165
 166                                         if (isalpha(refAligned[j])) {
 167                                                 numCovered++;
 168                                         }
 169                                 }
 170                         }
 171
 172                         int coverage = ((numCovered/(float)numBases)*100);
 173
 174                         //if coverage above minimum
 175                         if (coverage > minCoverage) {
 176                                 newRefs.push_back(ref[i]);
 177                         }else {
 178                                 delete ref[i];
 179                         }
 180                 }
 181
 182                 return newRefs;
 183         }
 184         catch(exception& e) {
 185                 m->errorOut(e, "Maligner", "minCoverageFilter");
 186                 exit(1);
 187         }
 188 }
 189 /***********************************************************************/
 190 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
 191 int Maligner::computeChimeraPenalty() {
 192         try {
 193
 194                 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
 195
 196 //              if(numAllowable < 1){   numAllowable = 1;       }
 197
 198                 int penalty = int(numAllowable + 1) * misMatchPenalty;
 199
 200                 return penalty;
 201
 202         }
 203         catch(exception& e) {
 204                 m->errorOut(e, "Maligner", "computeChimeraPenalty");
 205                 exit(1);
 206         }
 207 }
 208 /***********************************************************************/
 209 //this is a vertical filter
 210 void Maligner::verticalFilter(vector<Sequence*> seqs) {
 211         try {
 212                 vector<int> gaps;       gaps.resize(query->getAligned().length(), 0);
 213
 214                 string filterString = (string(query->getAligned().length(), '1'));
 215
 216                 //for each sequence
 217                 for (int i = 0; i < seqs.size(); i++) {
 218
 219                         string seqAligned = seqs[i]->getAligned();
 220
 221                         for (int j = 0; j < seqAligned.length(); j++) {
 222                                 //if this spot is a gap
 223                                 if ((seqAligned[j] == '-') || (seqAligned[j] == '.'))   {       gaps[j]++;      }
 224                         }
 225                 }
 226
 227                 //zero out spot where all sequences have blanks
 228                 int numColRemoved = 0;
 229                 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
 230                         if(gaps[i] == seqs.size())      {       filterString[i] = '0';  numColRemoved++;  }
 231                 }
 232
 233                 map<int, int> newMap;
 234                 //for each sequence
 235                 for (int i = 0; i < seqs.size(); i++) {
 236
 237                         string seqAligned = seqs[i]->getAligned();
 238                         string newAligned = "";
 239                         int count = 0;
 240
 241                         for (int j = 0; j < seqAligned.length(); j++) {
 242                                 //if this spot is not a gap
 243                                 if (filterString[j] == '1') {
 244                                         newAligned += seqAligned[j];
 245                                         newMap[count] = spotMap[j];
 246                                         count++;
 247                                 }
 248                         }
 249
 250                         seqs[i]->setAligned(newAligned);
 251                 }
 252
 253                 spotMap = newMap;
 254         }
 255         catch(exception& e) {
 256                 m->errorOut(e, "Maligner", "verticalFilter");
 257                 exit(1);
 258         }
 259 }
 260 //***************************************************************************************************************
 261 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
 262         try{
 263
 264                 vector< vector<score_struct> > m(rows);
 265
 266                 for (int i = 0; i < rows; i++) {
 267                         for (int j = 0; j < cols; j++) {
 268
 269                                 //initialize each cell
 270                                 score_struct temp;
 271                                 temp.prev = -1;
 272                                 temp.score = -9999999;
 273                                 temp.col = j;
 274                                 temp.row = i;
 275
 276                                 m[i].push_back(temp);
 277                         }
 278                 }
 279
 280                 return m;
 281         }
 282         catch(exception& e) {
 283                 m->errorOut(e, "Maligner", "buildScoreMatrix");
 284                 exit(1);
 285         }
 286 }
 287
 288 //***************************************************************************************************************
 289
 290 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
 291         try{
 292
 293                 //get matrix dimensions
 294                 int numCols = query->getAligned().length();
 295                 int numRows = seqs.size();
 296
 297                 //initialize first col
 298                 string queryAligned = query->getAligned();
 299                 for (int i = 0; i < numRows; i++) {
 300                         string subjectAligned = seqs[i]->getAligned();
 301
 302                         //are you both gaps?
 303                         if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
 304                                 ms[i][0].score = 0;
 305 //                              ms[i][0].mismatches = 0;
 306                         }else if (queryAligned[0] == subjectAligned[0])  { //|| subjectAligned[0] == 'N')
 307                                 ms[i][0].score = matchScore;
 308 //                              ms[i][0].mismatches = 0;
 309                         }else{
 310                                 ms[i][0].score = 0;
 311 //                              ms[i][0].mismatches = 1;
 312                         }
 313                 }
 314
 315                 //fill rest of matrix
 316                 for (int j = 1; j < numCols; j++) {  //iterate through matrix columns
 317
 318                         for (int i = 0; i < numRows; i++) {  //iterate through matrix rows
 319
 320                                 string subjectAligned = seqs[i]->getAligned();
 321
 322                                 int matchMisMatchScore = 0;
 323                                 //are you both gaps?
 324                                 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
 325                                         //leave the same
 326                                 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
 327                                         //matchMisMatchScore = matchScore;
 328                                         //leave the same
 329                                 }else if (queryAligned[j] == subjectAligned[j]) {
 330                                         matchMisMatchScore = matchScore;
 331 //                                      ms[i][j].mismatches = ms[i][j-1].mismatches;
 332                                 }else if (queryAligned[j] != subjectAligned[j]) {
 333                                         matchMisMatchScore = misMatchPenalty;
 334 //                                      ms[i][j].mismatches = ms[i][j-1].mismatches + 1;
 335                                 }
 336
 337                                 //compute score based on previous columns scores
 338                                 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
 339
 340                                         int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
 341
 342                                         //you are not at yourself
 343                                         if (prevIndex != i) {   sumScore += penalty;    }
 344                                         if (sumScore < 0)       {       sumScore = 0;                   }
 345
 346                                         if (sumScore > ms[i][j].score) {
 347                                                 ms[i][j].score = sumScore;
 348                                                 ms[i][j].prev = prevIndex;
 349                                         }
 350                                 }
 351                         }
 352                 }
 353
 354         /*      for(int i=0;i<numRows;i++){
 355                         cout << seqs[i]->getName();
 356                         for(int j=0;j<numCols;j++){
 357                                 cout << '\t' << ms[i][j].mismatches;
 358                         }
 359                         cout << endl;
 360                 }
 361                 cout << endl;*/
 362                 /*cout << numRows << '\t' << numCols << endl;
 363                 for(int i=0;i<numRows;i++){
 364                         cout << seqs[i]->getName() << endl << seqs[i]->getAligned() << endl << endl;
 365                         if ((seqs[i]->getName() == "S000003470") || (seqs[i]->getName() == "S000383265") || (seqs[i]->getName() == "7000004128191054")) {
 366                         for(int j=0;j<numCols;j++){
 367                                 cout << '\t' << ms[i][j].score;
 368                         }
 369                         cout << endl;
 370                         }
 371                 }
 372                 cout << endl;*/
 373                 /*for(int i=0;i<numRows;i++){
 374                         cout << seqs[i]->getName();
 375                         for(int j=0;j<numCols;j++){
 376                                 cout << '\t' << ms[i][j].prev;
 377                         }
 378                         cout << endl;
 379                 }*/
 380
 381
 382         }
 383         catch(exception& e) {
 384                 m->errorOut(e, "Maligner", "fillScoreMatrix");
 385                 exit(1);
 386         }
 387 }
 388
 389 //***************************************************************************************************************
 390
 391 vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
 392         try {
 393
 394
 395                 //get matrix dimensions
 396                 int numCols = query->getAligned().length();
 397                 int numRows = ms.size();
 398
 399
 400                 //find highest score scoring matrix
 401                 vector<score_struct> highestStruct;
 402                 int highestScore = 0;
 403
 404                 for (int i = 0; i < numRows; i++) {
 405                         for (int j = 0; j < numCols; j++) {
 406                                 if (ms[i][j].score > highestScore) {
 407                                         highestScore = ms[i][j].score;
 408                                         highestStruct.resize(0);
 409                                         highestStruct.push_back(ms[i][j]);
 410                                 }
 411                                 else if(ms[i][j].score == highestScore){
 412                                         highestStruct.push_back(ms[i][j]);
 413                                 }
 414                         }
 415                 }
 416
 417                 //cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl;
 418
 419                 vector<trace_struct> maxTrace;
 420                 double maxPercentIdenticalQueryAntiChimera = 0;
 421
 422                 for(int i=0;i<highestStruct.size();i++){
 423
 424                         vector<score_struct> path;
 425
 426                         int rowIndex = highestStruct[i].row;
 427                         int pos = highestStruct[i].col;
 428                         int score = highestStruct[i].score;
 429
 430                         while (pos >= 0 && score > 0) {
 431                                 score_struct temp = ms[rowIndex][pos];
 432                                 score = temp.score;
 433
 434                                 if (score > 0) {        path.push_back(temp);   }
 435
 436                                 rowIndex = temp.prev;
 437                                 pos--;
 438                         }
 439
 440                         reverse(path.begin(), path.end());
 441
 442                         vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
 443
 444                         //cout << "traces\n";
 445                         //for(int j=0;j<trace.size();j++){
 446                         //      cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
 447                         //}
 448
 449                         int traceStart = path[0].col;
 450                         int traceEnd = path[path.size()-1].col;
 451 //                      cout << "traceStart/End\t" << traceStart << '\t' << traceEnd << endl;
 452
 453                         string queryInRange = query->getAligned();
 454                         queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
 455 //                      cout << "here" << endl;
 456                         string chimeraSeq = constructChimericSeq(trace, refSeqs);
 457                         string antiChimeraSeq = constructAntiChimericSeq(trace, refSeqs);
 458
 459                         percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
 460                         double percentIdenticalQueryAntiChimera = computePercentID(queryInRange, antiChimeraSeq);
 461 //                      cout << i << '\t' << percentIdenticalQueryChimera << '\t' << percentIdenticalQueryAntiChimera << endl;
 462
 463                         if(percentIdenticalQueryAntiChimera > maxPercentIdenticalQueryAntiChimera){
 464                                 maxPercentIdenticalQueryAntiChimera = percentIdenticalQueryAntiChimera;
 465                                 maxTrace = trace;
 466                         }
 467                 }
 468 //              cout << maxPercentIdenticalQueryAntiChimera << endl;
 469                 return maxTrace;
 470
 471
 472         }
 473         catch(exception& e) {
 474                 m->errorOut(e, "Maligner", "extractHighestPath");
 475                 exit(1);
 476         }
 477 }
 478
 479 //***************************************************************************************************************
 480
 481 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
 482         try {
 483                 vector<trace_struct> trace;
 484
 485                 int region_index = path[0].row;
 486                 int region_start = path[0].col;
 487
 488                 for (int i = 1; i < path.size(); i++) {
 489
 490                         int next_region_index = path[i].row;
 491
 492                         if (next_region_index != region_index) {
 493
 494                                 // add trace region
 495                                 int col_index = path[i].col;
 496                                 trace_struct temp;
 497                                 temp.col = region_start;
 498                                 temp.oldCol = col_index-1;
 499                                 temp.row = region_index;
 500
 501                                 trace.push_back(temp);
 502
 503                                 region_index = path[i].row;
 504                                 region_start = col_index;
 505                         }
 506                 }
 507
 508                 // get last one
 509                 trace_struct temp;
 510                 temp.col = region_start;
 511                 temp.oldCol = path[path.size()-1].col;
 512                 temp.row = region_index;
 513                 trace.push_back(temp);
 514
 515 //              cout << endl;
 516 //              cout << trace.size() << endl;
 517 //              for(int i=0;i<trace.size();i++){
 518 //                      cout << seqs[trace[i].row]->getName() << endl;
 519 //              }
 520 //              cout << endl;
 521
 522                 return trace;
 523
 524         }
 525         catch(exception& e) {
 526                 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
 527                 exit(1);
 528         }
 529 }
 530
 531 //***************************************************************************************************************
 532
 533 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
 534         try {
 535                 string chimera = "";
 536
 537                 for (int i = 0; i < trace.size(); i++) {
 538 //                      cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
 539
 540                         string seqAlign = seqs[trace[i].row]->getAligned();
 541                         seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
 542                         chimera += seqAlign;
 543                 }
 544
 545                 return chimera;
 546         }
 547         catch(exception& e) {
 548                 m->errorOut(e, "Maligner", "constructChimericSeq");
 549                 exit(1);
 550         }
 551 }
 552
 553 //***************************************************************************************************************
 554
 555 string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
 556         try {
 557                 string antiChimera = "";
 558
 559                 for (int i = 0; i < trace.size(); i++) {
 560 //                      cout << i << '\t' << (trace.size() - i - 1) << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
 561
 562                         int oppositeIndex = trace.size() - i - 1;
 563
 564                         string seqAlign = seqs[trace[oppositeIndex].row]->getAligned();
 565                         seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
 566                         antiChimera += seqAlign;
 567                 }
 568
 569                 return antiChimera;
 570         }
 571         catch(exception& e) {
 572                 m->errorOut(e, "Maligner", "constructChimericSeq");
 573                 exit(1);
 574         }
 575 }
 576
 577 //***************************************************************************************************************
 578 float Maligner::computePercentID(string queryAlign, string chimera) {
 579         try {
 580
 581                 if (queryAlign.length() != chimera.length()) {
 582                         m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
 583                         m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine();
 584                         m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
 585                         return -1.0;
 586                 }
 587
 588
 589                 int numBases = 0;
 590                 int numIdentical = 0;
 591
 592                 for (int i = 0; i < queryAlign.length(); i++) {
 593                         if ((isalpha(queryAlign[i])) || (isalpha(chimera[i])))  {
 594                                 numBases++;
 595                                 if (queryAlign[i] == chimera[i]) {
 596                                         numIdentical++;
 597                                 }
 598                         }
 599                 }
 600
 601                 if (numBases == 0) { return 0; }
 602
 603                 float percentIdentical = (numIdentical/(float)numBases) * 100;
 604
 605                 return percentIdentical;
 606
 607         }
 608         catch(exception& e) {
 609                 m->errorOut(e, "Maligner", "computePercentID");
 610                 exit(1);
 611         }
 612 }
 613 //***************************************************************************************************************