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                         }
 182                 }
 183
 184                 return newRefs;
 185         }
 186         catch(exception& e) {
 187                 m->errorOut(e, "Maligner", "minCoverageFilter");
 188                 exit(1);
 189         }
 190 }
 191 /***********************************************************************/
 192 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
 193 int Maligner::computeChimeraPenalty() {
 194         try {
 195
 196                 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
 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; m.resize(rows);
 265
 266                 for (int i = 0; i < m.size(); 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 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
 289         try{
 290
 291                 //get matrix dimensions
 292                 int numCols = query->getAligned().length();
 293                 int numRows = seqs.size();
 294
 295                 //initialize first col
 296                 string queryAligned = query->getAligned();
 297                 for (int i = 0; i < numRows; i++) {
 298                         string subjectAligned = seqs[i]->getAligned();
 299
 300                         //are you both gaps?
 301                         if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
 302                                 ms[i][0].score = 0;
 303                         }else if (queryAligned[0] == subjectAligned[0]) {
 304                                 ms[i][0].score = matchScore;
 305                         }else{
 306                                 ms[i][0].score = 0;
 307                         }
 308                 }
 309
 310                 //fill rest of matrix
 311                 for (int j = 1; j < numCols; j++) {  //iterate through matrix columns
 312
 313                         for (int i = 0; i < numRows; i++) {  //iterate through matrix rows
 314
 315                                 string subjectAligned = seqs[i]->getAligned();
 316
 317                                 int matchMisMatchScore = 0;
 318                                 //are you both gaps?
 319                                 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
 320                                         //leave the same
 321                                 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
 322                                         //leave the same
 323                                 }else if (queryAligned[j] == subjectAligned[j]) {
 324                                         matchMisMatchScore = matchScore;
 325                                 }else if (queryAligned[j] != subjectAligned[j]) {
 326                                         matchMisMatchScore = misMatchPenalty;
 327                                 }
 328
 329                                 //compute score based on previous columns scores
 330                                 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
 331
 332                                         int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
 333
 334                                         //you are not at yourself
 335                                         if (prevIndex != i) {   sumScore += penalty;    }
 336                                         if (sumScore < 0)       {       sumScore = 0;                   }
 337
 338                                         if (sumScore > ms[i][j].score) {
 339                                                 ms[i][j].score = sumScore;
 340                                                 ms[i][j].prev = prevIndex;
 341                                         }
 342                                 }
 343                         }
 344                 }
 345
 346         }
 347         catch(exception& e) {
 348                 m->errorOut(e, "Maligner", "fillScoreMatrix");
 349                 exit(1);
 350         }
 351 }
 352 //***************************************************************************************************************
 353 vector<score_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
 354         try {
 355
 356                 //get matrix dimensions
 357                 int numCols = query->getAligned().length();
 358                 int numRows = ms.size();
 359
 360
 361                 //find highest score scoring matrix
 362                 score_struct highestStruct;
 363                 int highestScore = 0;
 364
 365                 for (int i = 0; i < numRows; i++) {
 366                         for (int j = 0; j < numCols; j++) {
 367                                 if (ms[i][j].score > highestScore) {
 368                                         highestScore = ms[i][j].score;
 369                                         highestStruct = ms[i][j];
 370                                 }
 371                         }
 372                 }
 373
 374                 vector<score_struct> path;
 375
 376                 int rowIndex = highestStruct.row;
 377                 int pos = highestStruct.col;
 378                 int score = highestStruct.score;
 379
 380                 while (pos >= 0 && score > 0) {
 381                         score_struct temp = ms[rowIndex][pos];
 382                         score = temp.score;
 383
 384                         if (score > 0) {        path.push_back(temp);   }
 385
 386                         rowIndex = temp.prev;
 387                         pos--;
 388                 }
 389
 390                 reverse(path.begin(), path.end());
 391
 392                 return path;
 393
 394         }
 395         catch(exception& e) {
 396                 m->errorOut(e, "Maligner", "extractHighestPath");
 397                 exit(1);
 398         }
 399 }
 400 //***************************************************************************************************************
 401 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
 402         try {
 403                 vector<trace_struct> trace;
 404
 405                 int region_index = path[0].row;
 406                 int region_start = path[0].col;
 407
 408                 for (int i = 1; i < path.size(); i++) {
 409
 410                         int next_region_index = path[i].row;
 411
 412                         if (next_region_index != region_index) {
 413
 414                                 // add trace region
 415                                 int col_index = path[i].col;
 416                                 trace_struct temp;
 417                                 temp.col = region_start;
 418                                 temp.oldCol = col_index-1;
 419                                 temp.row = region_index;
 420
 421                                 trace.push_back(temp);
 422
 423                                 region_index = path[i].row;
 424                                 region_start = col_index;
 425                         }
 426                 }
 427
 428                 // get last one
 429                 trace_struct temp;
 430                 temp.col = region_start;
 431                 temp.oldCol = path[path.size()-1].col;
 432                 temp.row = region_index;
 433                 trace.push_back(temp);
 434
 435                 return trace;
 436
 437         }
 438         catch(exception& e) {
 439                 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
 440                 exit(1);
 441         }
 442 }
 443 //***************************************************************************************************************
 444 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
 445         try {
 446                 string chimera = "";
 447
 448                 for (int i = 0; i < trace.size(); i++) {
 449                         string seqAlign = seqs[trace[i].row]->getAligned();
 450                         seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
 451                         chimera += seqAlign;
 452                 }
 453
 454                 return chimera;
 455         }
 456         catch(exception& e) {
 457                 m->errorOut(e, "Maligner", "constructChimericSeq");
 458                 exit(1);
 459         }
 460 }
 461 //***************************************************************************************************************
 462 float Maligner::computePercentID(string queryAlign, string chimera) {
 463         try {
 464
 465                 if (queryAlign.length() != chimera.length()) {
 466                         m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
 467                         m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine();
 468                         m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
 469                         return -1.0;
 470                 }
 471
 472
 473                 int numBases = 0;
 474                 int numIdentical = 0;
 475
 476                 for (int i = 0; i < queryAlign.length(); i++) {
 477                         if ((isalpha(queryAlign[i])) || (isalpha(chimera[i])))  {
 478                                 numBases++;
 479                                 if (queryAlign[i] == chimera[i]) {
 480                                         numIdentical++;
 481                                 }
 482                         }
 483                 }
 484
 485                 if (numBases == 0) { return 0; }
 486
 487                 float percentIdentical = (numIdentical/(float)numBases) * 100;
 488
 489                 return percentIdentical;
 490
 491         }
 492         catch(exception& e) {
 493                 m->errorOut(e, "Maligner", "computePercentID");
 494                 exit(1);
 495         }
 496 }
 497 //***************************************************************************************************************
 498 vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
 499         try {
 500                 indexes.clear();
 501                 vector<Sequence*> refResults;
 502
 503                 //get parts of query
 504                 string queryUnAligned = q->getUnaligned();
 505                 string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
 506                 string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
 507
 508                 Sequence* queryLeft = new Sequence(q->getName()+"left", leftQuery);
 509                 Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery);
 510
 511                 vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
 512                 vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
 513
 514                 //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; }
 515
 516                 vector<int> smaller;
 517                 vector<int> larger;
 518
 519                 if (tempIndexesRight.size() < tempIndexesLeft.size()) { smaller = tempIndexesRight;  larger = tempIndexesLeft; }
 520                 else { smaller = tempIndexesLeft; larger = tempIndexesRight; }
 521
 522                 //merge results
 523                 map<int, int> seen;
 524                 map<int, int>::iterator it;
 525
 526                 vector<int> mergedResults;
 527                 for (int i = 0; i < smaller.size(); i++) {
 528                         //add left if you havent already
 529                         it = seen.find(smaller[i]);
 530                         if (it == seen.end()) {
 531                                 mergedResults.push_back(smaller[i]);
 532                                 seen[smaller[i]] = smaller[i];
 533                         }
 534
 535                         //add right if you havent already
 536                         it = seen.find(larger[i]);
 537                         if (it == seen.end()) {
 538                                 mergedResults.push_back(larger[i]);
 539                                 seen[larger[i]] = larger[i];
 540                         }
 541                 }
 542
 543                 for (int i = smaller.size(); i < larger.size(); i++) {
 544                         it = seen.find(larger[i]);
 545                         if (it == seen.end()) {
 546                                 mergedResults.push_back(larger[i]);
 547                                 seen[larger[i]] = larger[i];
 548                         }
 549                 }
 550
 551                 if (mergedResults.size() < numWanted) { numWanted = mergedResults.size(); }
 552 //cout << q->getName() <<  endl;
 553                 for (int i = 0; i < numWanted; i++) {
 554 //cout << db[mergedResults[i]]->getName() << endl;
 555                         if (db[mergedResults[i]]->getName() != q->getName()) {
 556                                 Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
 557                                 refResults.push_back(temp);
 558                                 indexes.push_back(mergedResults[i]);
 559                         }
 560 //cout << mergedResults[i] << endl;
 561                 }
 562 //cout << endl;
 563                 delete queryRight;
 564                 delete queryLeft;
 565
 566                 return refResults;
 567         }
 568         catch(exception& e) {
 569                 m->errorOut(e, "Maligner", "getBlastSeqs");
 570                 exit(1);
 571         }
 572 }
 573 //***************************************************************************************************************
 574 vector<Sequence*> Maligner::getKmerSeqs(Sequence* q, int num) {
 575         try {
 576                 indexes.clear();
 577
 578                 //get parts of query
 579                 string queryUnAligned = q->getUnaligned();
 580                 string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
 581                 string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
 582
 583                 Sequence* queryLeft = new Sequence(q->getName(), leftQuery);
 584                 Sequence* queryRight = new Sequence(q->getName(), rightQuery);
 585
 586                 vector<int> tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted);
 587                 vector<int> tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted);
 588
 589                 //merge results
 590                 map<int, int> seen;
 591                 map<int, int>::iterator it;
 592
 593                 vector<int> mergedResults;
 594                 for (int i = 0; i < tempIndexesLeft.size(); i++) {
 595                         //add left if you havent already
 596                         it = seen.find(tempIndexesLeft[i]);
 597                         if (it == seen.end()) {
 598                                 mergedResults.push_back(tempIndexesLeft[i]);
 599                                 seen[tempIndexesLeft[i]] = tempIndexesLeft[i];
 600                         }
 601
 602                         //add right if you havent already
 603                         it = seen.find(tempIndexesRight[i]);
 604                         if (it == seen.end()) {
 605                                 mergedResults.push_back(tempIndexesRight[i]);
 606                                 seen[tempIndexesRight[i]] = tempIndexesRight[i];
 607                         }
 608                 }
 609
 610 //cout << q->getName() << endl;
 611                 vector<Sequence*> refResults;
 612                 for (int i = 0; i < numWanted; i++) {
 613 //cout << db[mergedResults[i]]->getName() << endl;
 614                         Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
 615                         refResults.push_back(temp);
 616                         indexes.push_back(mergedResults[i]);
 617                 }
 618 //cout << endl;
 619                 delete queryRight;
 620                 delete queryLeft;
 621
 622                 return refResults;
 623         }
 624         catch(exception& e) {
 625                 m->errorOut(e, "Maligner", "getBlastSeqs");
 626                 exit(1);
 627         }
 628 }
 629 //***************************************************************************************************************
 630