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