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