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
  19                         m = MothurOut::getInstance();
  20
  21 //                      cout << matchScore << '\t' << misMatchPenalty << endl;
  22 //
  23 //                      matchScore = 1;
  24 //                      misMatchPenalty = -1;
  25
  26                 }
  27
  28 /***********************************************************************/
  29 string Maligner::getResults(Sequence* q, DeCalculator* decalc) {
  30         try {
  31
  32                 outputResults.clear();
  33
  34                 //make copy so trimming doesn't destroy query from calling class - remember to deallocate
  35                 query = new Sequence(q->getName(), q->getAligned());
  36
  37                 string chimera;
  38
  39                 if (searchMethod == "distance") {
  40                         //find closest seqs to query in template - returns copies of seqs so trim does not destroy - remember to deallocate
  41                         refSeqs = decalc->findClosest(query, db, numWanted, indexes);
  42                 }else if (searchMethod == "blast")  {
  43                         refSeqs = getBlastSeqs(query, numWanted); //fills indexes
  44                 }else if (searchMethod == "kmer") {
  45                         refSeqs = getKmerSeqs(query, numWanted); //fills indexes
  46                 }else { m->mothurOut("not valid search."); exit(1);  } //should never get here
  47
  48                 if (m->control_pressed) { return chimera;  }
  49
  50                 refSeqs = minCoverageFilter(refSeqs);
  51
  52
  53
  54                 if (refSeqs.size() < 2)  {
  55                         for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
  56                         percentIdenticalQueryChimera = 0.0;
  57                         return "unknown";
  58                 }
  59
  60                 int chimeraPenalty = computeChimeraPenalty();
  61                 //cout << "chimeraPenalty = " << chimeraPenalty << endl;
  62                 //fills outputResults
  63                 chimera = chimeraMaligner(chimeraPenalty, decalc);
  64
  65                 if (m->control_pressed) { return chimera;  }
  66
  67                 //free memory
  68                 delete query;
  69
  70                 for (int i = 0; i < refSeqs.size(); i++) {  delete refSeqs[i];  }
  71
  72                 return chimera;
  73         }
  74         catch(exception& e) {
  75                 m->errorOut(e, "Maligner", "getResults");
  76                 exit(1);
  77         }
  78 }
  79 /***********************************************************************/
  80 string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) {
  81         try {
  82
  83                 string chimera;
  84
  85                 //trims seqs to first non gap char in all seqs and last non gap char in all seqs
  86                 spotMap = decalc->trimSeqs(query, refSeqs);
  87
  88                 //you trimmed the whole sequence, skip
  89                 if (query->getAligned() == "") { return "no"; }
  90
  91                 vector<Sequence*> temp = refSeqs;
  92                 temp.push_back(query);
  93
  94                 verticalFilter(temp);
  95
  96                 //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl << refSeqs[i]->getAligned() << endl; }
  97
  98                 vector< vector<score_struct> > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes
  99
 100                 if (m->control_pressed) { return chimera;  }
 101
 102                 fillScoreMatrix(matrix, refSeqs, chimeraPenalty);
 103
 104                 vector<trace_struct> trace = extractHighestPath(matrix);
 105
 106                 //cout << "traces\n";
 107                 //for(int i=0;i<trace.size();i++){
 108                 //      cout << trace[i].col << '\t' << trace[i].oldCol << '\t' << refSeqs[trace[i].row]->getName() << endl;
 109                 //}
 110
 111                 if (trace.size() > 1) {         chimera = "yes";        }
 112                 else { chimera = "no";  return chimera; }
 113
 114                 int traceStart = trace[0].col;
 115                 int traceEnd = trace[trace.size()-1].oldCol;
 116                 string queryInRange = query->getAligned();
 117                 queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
 118
 119                 if (m->control_pressed) { return chimera;  }
 120
 121                 //save output results
 122                 for (int i = 0; i < trace.size(); i++) {
 123                         int regionStart = trace[i].col;
 124                         int regionEnd = trace[i].oldCol;
 125                         int seqIndex = trace[i].row;
 126
 127                         results temp;
 128
 129                         temp.parent = refSeqs[seqIndex]->getName();
 130                         temp.parentAligned = db[indexes[seqIndex]]->getAligned();
 131                         temp.nastRegionStart = spotMap[regionStart];
 132                         temp.nastRegionEnd = spotMap[regionEnd];
 133                         temp.regionStart = regionStart;
 134                         temp.regionEnd = regionEnd;
 135
 136                         string parentInRange = refSeqs[seqIndex]->getAligned();
 137                         parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1));
 138
 139                         temp.queryToParent = computePercentID(queryInRange, parentInRange);
 140                         temp.divR = (percentIdenticalQueryChimera / temp.queryToParent);
 141
 142                         string queryInRegion = query->getAligned();
 143                         queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1));
 144
 145                         string parentInRegion = refSeqs[seqIndex]->getAligned();
 146                         parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1));
 147
 148                         temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion);
 149
 150                         outputResults.push_back(temp);
 151                 }
 152
 153                 return chimera;
 154         }
 155         catch(exception& e) {
 156                 m->errorOut(e, "Maligner", "chimeraMaligner");
 157                 exit(1);
 158         }
 159 }
 160 /***********************************************************************/
 161 //removes top matches that do not have minimum coverage with query.
 162 vector<Sequence*> Maligner::minCoverageFilter(vector<Sequence*> ref){
 163         try {
 164                 vector<Sequence*> newRefs;
 165
 166                 string queryAligned = query->getAligned();
 167
 168                 for (int i = 0; i < ref.size(); i++) {
 169
 170                         string refAligned = ref[i]->getAligned();
 171
 172                         int numBases = 0;
 173                         int numCovered = 0;
 174
 175                         //calculate coverage
 176                         for (int j = 0; j < queryAligned.length(); j++) {
 177
 178                                 if (isalpha(queryAligned[j])) {
 179                                         numBases++;
 180
 181                                         if (isalpha(refAligned[j])) {
 182                                                 numCovered++;
 183                                         }
 184                                 }
 185                         }
 186
 187                         int coverage = ((numCovered/(float)numBases)*100);
 188
 189                         //if coverage above minimum
 190                         if (coverage > minCoverage) {
 191                                 newRefs.push_back(ref[i]);
 192                         }else {
 193                                 delete ref[i];
 194                         }
 195                 }
 196
 197                 return newRefs;
 198         }
 199         catch(exception& e) {
 200                 m->errorOut(e, "Maligner", "minCoverageFilter");
 201                 exit(1);
 202         }
 203 }
 204 /***********************************************************************/
 205 // a breakpoint should yield fewer mismatches than this number with respect to the best parent sequence.
 206 int Maligner::computeChimeraPenalty() {
 207         try {
 208
 209                 int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases());
 210
 211 //              if(numAllowable < 1){   numAllowable = 1;       }
 212
 213                 int penalty = int(numAllowable + 1) * misMatchPenalty;
 214
 215                 return penalty;
 216
 217         }
 218         catch(exception& e) {
 219                 m->errorOut(e, "Maligner", "computeChimeraPenalty");
 220                 exit(1);
 221         }
 222 }
 223 /***********************************************************************/
 224 //this is a vertical filter
 225 void Maligner::verticalFilter(vector<Sequence*> seqs) {
 226         try {
 227                 vector<int> gaps;       gaps.resize(query->getAligned().length(), 0);
 228
 229                 string filterString = (string(query->getAligned().length(), '1'));
 230
 231                 //for each sequence
 232                 for (int i = 0; i < seqs.size(); i++) {
 233
 234                         string seqAligned = seqs[i]->getAligned();
 235
 236                         for (int j = 0; j < seqAligned.length(); j++) {
 237                                 //if this spot is a gap
 238                                 if ((seqAligned[j] == '-') || (seqAligned[j] == '.'))   {       gaps[j]++;      }
 239                         }
 240                 }
 241
 242                 //zero out spot where all sequences have blanks
 243                 int numColRemoved = 0;
 244                 for(int i = 0; i < seqs[0]->getAligned().length(); i++){
 245                         if(gaps[i] == seqs.size())      {       filterString[i] = '0';  numColRemoved++;  }
 246                 }
 247
 248                 map<int, int> newMap;
 249                 //for each sequence
 250                 for (int i = 0; i < seqs.size(); i++) {
 251
 252                         string seqAligned = seqs[i]->getAligned();
 253                         string newAligned = "";
 254                         int count = 0;
 255
 256                         for (int j = 0; j < seqAligned.length(); j++) {
 257                                 //if this spot is not a gap
 258                                 if (filterString[j] == '1') {
 259                                         newAligned += seqAligned[j];
 260                                         newMap[count] = spotMap[j];
 261                                         count++;
 262                                 }
 263                         }
 264
 265                         seqs[i]->setAligned(newAligned);
 266                 }
 267
 268                 spotMap = newMap;
 269         }
 270         catch(exception& e) {
 271                 m->errorOut(e, "Maligner", "verticalFilter");
 272                 exit(1);
 273         }
 274 }
 275 //***************************************************************************************************************
 276 vector< vector<score_struct> > Maligner::buildScoreMatrix(int cols, int rows) {
 277         try{
 278
 279                 vector< vector<score_struct> > m(rows);
 280
 281                 for (int i = 0; i < rows; i++) {
 282                         for (int j = 0; j < cols; j++) {
 283
 284                                 //initialize each cell
 285                                 score_struct temp;
 286                                 temp.prev = -1;
 287                                 temp.score = -9999999;
 288                                 temp.col = j;
 289                                 temp.row = i;
 290
 291                                 m[i].push_back(temp);
 292                         }
 293                 }
 294
 295                 return m;
 296         }
 297         catch(exception& e) {
 298                 m->errorOut(e, "Maligner", "buildScoreMatrix");
 299                 exit(1);
 300         }
 301 }
 302
 303 //***************************************************************************************************************
 304
 305 void Maligner::fillScoreMatrix(vector<vector<score_struct> >& ms, vector<Sequence*> seqs, int penalty) {
 306         try{
 307
 308                 //get matrix dimensions
 309                 int numCols = query->getAligned().length();
 310                 int numRows = seqs.size();
 311
 312                 //initialize first col
 313                 string queryAligned = query->getAligned();
 314                 for (int i = 0; i < numRows; i++) {
 315                         string subjectAligned = seqs[i]->getAligned();
 316
 317                         //are you both gaps?
 318                         if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) {
 319                                 ms[i][0].score = 0;
 320 //                              ms[i][0].mismatches = 0;
 321                         }else if (queryAligned[0] == subjectAligned[0])  { //|| subjectAligned[0] == 'N')
 322                                 ms[i][0].score = matchScore;
 323 //                              ms[i][0].mismatches = 0;
 324                         }else{
 325                                 ms[i][0].score = 0;
 326 //                              ms[i][0].mismatches = 1;
 327                         }
 328                 }
 329
 330                 //fill rest of matrix
 331                 for (int j = 1; j < numCols; j++) {  //iterate through matrix columns
 332
 333                         for (int i = 0; i < numRows; i++) {  //iterate through matrix rows
 334
 335                                 string subjectAligned = seqs[i]->getAligned();
 336
 337                                 int matchMisMatchScore = 0;
 338                                 //are you both gaps?
 339                                 if ((!isalpha(queryAligned[j])) && (!isalpha(subjectAligned[j]))) {
 340                                         //leave the same
 341                                 }else if ((toupper(queryAligned[j]) == 'N') || (toupper(subjectAligned[j]) == 'N')) {
 342                                         //matchMisMatchScore = matchScore;
 343                                         //leave the same
 344                                 }else if (queryAligned[j] == subjectAligned[j]) {
 345                                         matchMisMatchScore = matchScore;
 346 //                                      ms[i][j].mismatches = ms[i][j-1].mismatches;
 347                                 }else if (queryAligned[j] != subjectAligned[j]) {
 348                                         matchMisMatchScore = misMatchPenalty;
 349 //                                      ms[i][j].mismatches = ms[i][j-1].mismatches + 1;
 350                                 }
 351
 352                                 //compute score based on previous columns scores
 353                                 for (int prevIndex = 0; prevIndex < numRows; prevIndex++) { //iterate through rows
 354
 355                                         int sumScore = matchMisMatchScore + ms[prevIndex][j-1].score;
 356
 357                                         //you are not at yourself
 358                                         if (prevIndex != i) {   sumScore += penalty;    }
 359                                         if (sumScore < 0)       {       sumScore = 0;                   }
 360
 361                                         if (sumScore > ms[i][j].score) {
 362                                                 ms[i][j].score = sumScore;
 363                                                 ms[i][j].prev = prevIndex;
 364                                         }
 365                                 }
 366                         }
 367                 }
 368
 369         /*      for(int i=0;i<numRows;i++){
 370                         cout << seqs[i]->getName();
 371                         for(int j=0;j<numCols;j++){
 372                                 cout << '\t' << ms[i][j].mismatches;
 373                         }
 374                         cout << endl;
 375                 }
 376                 cout << endl;*/
 377                 /*cout << numRows << '\t' << numCols << endl;
 378                 for(int i=0;i<numRows;i++){
 379                         cout << seqs[i]->getName() << endl << seqs[i]->getAligned() << endl << endl;
 380                         if ((seqs[i]->getName() == "S000003470") || (seqs[i]->getName() == "S000383265") || (seqs[i]->getName() == "7000004128191054")) {
 381                         for(int j=0;j<numCols;j++){
 382                                 cout << '\t' << ms[i][j].score;
 383                         }
 384                         cout << endl;
 385                         }
 386                 }
 387                 cout << endl;*/
 388                 /*for(int i=0;i<numRows;i++){
 389                         cout << seqs[i]->getName();
 390                         for(int j=0;j<numCols;j++){
 391                                 cout << '\t' << ms[i][j].prev;
 392                         }
 393                         cout << endl;
 394                 }*/
 395
 396
 397         }
 398         catch(exception& e) {
 399                 m->errorOut(e, "Maligner", "fillScoreMatrix");
 400                 exit(1);
 401         }
 402 }
 403
 404 //***************************************************************************************************************
 405
 406 vector<trace_struct> Maligner::extractHighestPath(vector<vector<score_struct> > ms) {
 407         try {
 408
 409
 410                 //get matrix dimensions
 411                 int numCols = query->getAligned().length();
 412                 int numRows = ms.size();
 413
 414
 415                 //find highest score scoring matrix
 416                 vector<score_struct> highestStruct;
 417                 int highestScore = 0;
 418
 419                 for (int i = 0; i < numRows; i++) {
 420                         for (int j = 0; j < numCols; j++) {
 421                                 if (ms[i][j].score > highestScore) {
 422                                         highestScore = ms[i][j].score;
 423                                         highestStruct.resize(0);
 424                                         highestStruct.push_back(ms[i][j]);
 425                                 }
 426                                 else if(ms[i][j].score == highestScore){
 427                                         highestStruct.push_back(ms[i][j]);
 428                                 }
 429                         }
 430                 }
 431
 432                 //cout << endl << highestScore << '\t' << highestStruct.size() << '\t' << highestStruct[0].row << endl;
 433
 434                 vector<trace_struct> maxTrace;
 435                 double maxPercentIdenticalQueryAntiChimera = 0;
 436
 437                 for(int i=0;i<highestStruct.size();i++){
 438
 439                         vector<score_struct> path;
 440
 441                         int rowIndex = highestStruct[i].row;
 442                         int pos = highestStruct[i].col;
 443                         int score = highestStruct[i].score;
 444
 445                         while (pos >= 0 && score > 0) {
 446                                 score_struct temp = ms[rowIndex][pos];
 447                                 score = temp.score;
 448
 449                                 if (score > 0) {        path.push_back(temp);   }
 450
 451                                 rowIndex = temp.prev;
 452                                 pos--;
 453                         }
 454
 455                         reverse(path.begin(), path.end());
 456
 457                         vector<trace_struct> trace = mapTraceRegionsToAlignment(path, refSeqs);
 458
 459                         //cout << "traces\n";
 460                         //for(int j=0;j<trace.size();j++){
 461                         //      cout << trace[j].col << '\t' << trace[j].oldCol << '\t' << refSeqs[trace[j].row]->getName() << endl;
 462                         //}
 463
 464                         int traceStart = path[0].col;
 465                         int traceEnd = path[path.size()-1].col;
 466 //                      cout << "traceStart/End\t" << traceStart << '\t' << traceEnd << endl;
 467
 468                         string queryInRange = query->getAligned();
 469                         queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1));
 470 //                      cout << "here" << endl;
 471                         string chimeraSeq = constructChimericSeq(trace, refSeqs);
 472                         string antiChimeraSeq = constructAntiChimericSeq(trace, refSeqs);
 473
 474                         percentIdenticalQueryChimera = computePercentID(queryInRange, chimeraSeq);
 475                         double percentIdenticalQueryAntiChimera = computePercentID(queryInRange, antiChimeraSeq);
 476 //                      cout << i << '\t' << percentIdenticalQueryChimera << '\t' << percentIdenticalQueryAntiChimera << endl;
 477
 478                         if(percentIdenticalQueryAntiChimera > maxPercentIdenticalQueryAntiChimera){
 479                                 maxPercentIdenticalQueryAntiChimera = percentIdenticalQueryAntiChimera;
 480                                 maxTrace = trace;
 481                         }
 482                 }
 483 //              cout << maxPercentIdenticalQueryAntiChimera << endl;
 484                 return maxTrace;
 485
 486
 487         }
 488         catch(exception& e) {
 489                 m->errorOut(e, "Maligner", "extractHighestPath");
 490                 exit(1);
 491         }
 492 }
 493
 494 //***************************************************************************************************************
 495
 496 vector<trace_struct> Maligner::mapTraceRegionsToAlignment(vector<score_struct> path, vector<Sequence*> seqs) {
 497         try {
 498                 vector<trace_struct> trace;
 499
 500                 int region_index = path[0].row;
 501                 int region_start = path[0].col;
 502
 503                 for (int i = 1; i < path.size(); i++) {
 504
 505                         int next_region_index = path[i].row;
 506
 507                         if (next_region_index != region_index) {
 508
 509                                 // add trace region
 510                                 int col_index = path[i].col;
 511                                 trace_struct temp;
 512                                 temp.col = region_start;
 513                                 temp.oldCol = col_index-1;
 514                                 temp.row = region_index;
 515
 516                                 trace.push_back(temp);
 517
 518                                 region_index = path[i].row;
 519                                 region_start = col_index;
 520                         }
 521                 }
 522
 523                 // get last one
 524                 trace_struct temp;
 525                 temp.col = region_start;
 526                 temp.oldCol = path[path.size()-1].col;
 527                 temp.row = region_index;
 528                 trace.push_back(temp);
 529
 530 //              cout << endl;
 531 //              cout << trace.size() << endl;
 532 //              for(int i=0;i<trace.size();i++){
 533 //                      cout << seqs[trace[i].row]->getName() << endl;
 534 //              }
 535 //              cout << endl;
 536
 537                 return trace;
 538
 539         }
 540         catch(exception& e) {
 541                 m->errorOut(e, "Maligner", "mapTraceRegionsToAlignment");
 542                 exit(1);
 543         }
 544 }
 545
 546 //***************************************************************************************************************
 547
 548 string Maligner::constructChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
 549         try {
 550                 string chimera = "";
 551
 552                 for (int i = 0; i < trace.size(); i++) {
 553 //                      cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
 554
 555                         string seqAlign = seqs[trace[i].row]->getAligned();
 556                         seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
 557                         chimera += seqAlign;
 558                 }
 559
 560                 return chimera;
 561         }
 562         catch(exception& e) {
 563                 m->errorOut(e, "Maligner", "constructChimericSeq");
 564                 exit(1);
 565         }
 566 }
 567
 568 //***************************************************************************************************************
 569
 570 string Maligner::constructAntiChimericSeq(vector<trace_struct> trace, vector<Sequence*> seqs) {
 571         try {
 572                 string antiChimera = "";
 573
 574                 for (int i = 0; i < trace.size(); i++) {
 575 //                      cout << i << '\t' << (trace.size() - i - 1) << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl;
 576
 577                         int oppositeIndex = trace.size() - i - 1;
 578
 579                         string seqAlign = seqs[trace[oppositeIndex].row]->getAligned();
 580                         seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1));
 581                         antiChimera += seqAlign;
 582                 }
 583
 584                 return antiChimera;
 585         }
 586         catch(exception& e) {
 587                 m->errorOut(e, "Maligner", "constructChimericSeq");
 588                 exit(1);
 589         }
 590 }
 591
 592 //***************************************************************************************************************
 593 float Maligner::computePercentID(string queryAlign, string chimera) {
 594         try {
 595
 596                 if (queryAlign.length() != chimera.length()) {
 597                         m->mothurOut("Error, alignment strings are of different lengths: "); m->mothurOutEndLine();
 598                         m->mothurOut(toString(queryAlign.length())); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOutEndLine();
 599                         m->mothurOut(toString(chimera.length())); m->mothurOutEndLine();
 600                         return -1.0;
 601                 }
 602
 603
 604                 int numBases = 0;
 605                 int numIdentical = 0;
 606
 607                 for (int i = 0; i < queryAlign.length(); i++) {
 608                         if ((isalpha(queryAlign[i])) || (isalpha(chimera[i])))  {
 609                                 numBases++;
 610                                 if (queryAlign[i] == chimera[i]) {
 611                                         numIdentical++;
 612                                 }
 613                         }
 614                 }
 615
 616                 if (numBases == 0) { return 0; }
 617
 618                 float percentIdentical = (numIdentical/(float)numBases) * 100;
 619
 620                 return percentIdentical;
 621
 622         }
 623         catch(exception& e) {
 624                 m->errorOut(e, "Maligner", "computePercentID");
 625                 exit(1);
 626         }
 627 }
 628 //***************************************************************************************************************
 629 vector<Sequence*> Maligner::getBlastSeqs(Sequence* q, int num) {
 630         try {
 631                 indexes.clear();
 632                 vector<Sequence*> refResults;
 633
 634                 //get parts of query
 635                 string queryUnAligned = q->getUnaligned();
 636                 string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
 637                 string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
 638
 639                 Sequence* queryLeft = new Sequence(q->getName()+"left", leftQuery);
 640                 Sequence* queryRight = new Sequence(q->getName()+"right", rightQuery);
 641
 642                 vector<int> tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1);
 643                 vector<float> leftScores = databaseLeft->getSearchScores();
 644                 vector<int> tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1);
 645                 vector<float> rightScores = databaseLeft->getSearchScores();
 646
 647                 //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; }
 648
 649                 vector<int> smaller;
 650                 vector<float> smallerScores;
 651                 vector<int> larger;
 652                 vector<float> largerScores;
 653
 654                 if (tempIndexesRight.size() < tempIndexesLeft.size()) { smaller = tempIndexesRight; smallerScores = rightScores; larger = tempIndexesLeft; largerScores = leftScores; }
 655                 else { smaller = tempIndexesLeft; smallerScores = leftScores; larger = tempIndexesRight; largerScores = rightScores; }
 656
 657                 //for (int i = 0; i < smaller.size(); i++) { cout << "smaller = " << smaller[i] << '\t' << smallerScores[i] << endl; }
 658                 //cout << endl;
 659                 //for (int i = 0; i < larger.size(); i++) { cout << "larger = " << larger[i] << '\t' << largerScores[i] << endl; }
 660
 661                 //merge results
 662                 map<int, int> seen;
 663                 map<int, int>::iterator it;
 664                 float lastSmaller = smallerScores[0];
 665                 float lastLarger = largerScores[0];
 666                 int lasti = 0;
 667                 vector<int> mergedResults;
 668                 for (int i = 0; i < smaller.size(); i++) {
 669                         //add left if you havent already
 670                         it = seen.find(smaller[i]);
 671                         if (it == seen.end()) {
 672                                 mergedResults.push_back(smaller[i]);
 673                                 seen[smaller[i]] = smaller[i];
 674                                 lastSmaller = smallerScores[i];
 675                         }
 676
 677                         //add right if you havent already
 678                         it = seen.find(larger[i]);
 679                         if (it == seen.end()) {
 680                                 mergedResults.push_back(larger[i]);
 681                                 seen[larger[i]] = larger[i];
 682                                 lastLarger = largerScores[i];
 683                         }
 684
 685                         lasti = i;
 686                         //if (mergedResults.size() > num) { break; }
 687                 }
 688
 689                 //save lasti for smaller ties below
 690                 /*lasti++;
 691                 int iSmaller = lasti;
 692
 693                 if (!(mergedResults.size() > num)) { //if we still need more results.
 694                         for (int i = smaller.size(); i < larger.size(); i++) {
 695                                 it = seen.find(larger[i]);
 696                                 if (it == seen.end()) {
 697                                         mergedResults.push_back(larger[i]);
 698                                         seen[larger[i]] = larger[i];
 699                                         lastLarger = largerScores[i];
 700                                 }
 701
 702                                 lasti = i;
 703                                 if (mergedResults.size() > num) {  break; }
 704                         }
 705                 }
 706
 707
 708                 //add in any ties from smaller
 709                 while (iSmaller < smaller.size()) {
 710                         if (smallerScores[iSmaller] == lastSmaller) {
 711                                 it = seen.find(smaller[iSmaller]);
 712
 713                                 if (it == seen.end()) {
 714                                         mergedResults.push_back(smaller[iSmaller]);
 715                                         seen[smaller[iSmaller]] = smaller[iSmaller];
 716                                 }
 717                         }
 718                         else { break; }
 719                         iSmaller++;
 720                 }
 721
 722                 lasti++;
 723                 //add in any ties from larger
 724                 while (lasti < larger.size()) {
 725                         if (largerScores[lasti] == lastLarger) { //is it a tie
 726                                 it = seen.find(larger[lasti]);
 727
 728                                 if (it == seen.end()) {  //we haven't already seen it
 729                                         mergedResults.push_back(larger[lasti]);
 730                                         seen[larger[lasti]] = larger[lasti];
 731                                 }
 732                         }
 733                         else { break; }
 734                         lasti++;
 735                 }
 736                 */
 737
 738                 for (int i = smaller.size(); i < larger.size(); i++) {
 739                         //add right if you havent already
 740                         it = seen.find(larger[i]);
 741                         if (it == seen.end()) {
 742                                 mergedResults.push_back(larger[i]);
 743                                 seen[larger[i]] = larger[i];
 744                                 lastLarger = largerScores[i];
 745                         }
 746                 }
 747                 numWanted = mergedResults.size();
 748
 749                 if (mergedResults.size() < numWanted) { numWanted = mergedResults.size(); }
 750 //cout << q->getName() << " merged results size = " << mergedResults.size() << '\t' << "numwanted = " << numWanted <<  endl;
 751                 for (int i = 0; i < numWanted; i++) {
 752 //cout << db[mergedResults[i]]->getName()  << '\t' << mergedResults[i] << endl;
 753
 754                         if (db[mergedResults[i]]->getName() != q->getName()) {
 755                                 Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
 756                                 refResults.push_back(temp);
 757                                 indexes.push_back(mergedResults[i]);
 758                                 //cout << db[mergedResults[i]]->getName() << endl;
 759                         }
 760
 761 //cout << mergedResults[i] << endl;
 762                 }
 763 //cout << "done " << q->getName()  << endl;
 764                 delete queryRight;
 765                 delete queryLeft;
 766
 767                 return refResults;
 768         }
 769         catch(exception& e) {
 770                 m->errorOut(e, "Maligner", "getBlastSeqs");
 771                 exit(1);
 772         }
 773 }
 774 //***************************************************************************************************************
 775 vector<Sequence*> Maligner::getKmerSeqs(Sequence* q, int num) {
 776         try {
 777                 indexes.clear();
 778
 779                 //get parts of query
 780                 string queryUnAligned = q->getUnaligned();
 781                 string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence
 782                 string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence
 783
 784                 Sequence* queryLeft = new Sequence(q->getName(), leftQuery);
 785                 Sequence* queryRight = new Sequence(q->getName(), rightQuery);
 786
 787                 vector<int> tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, numWanted);
 788                 vector<int> tempIndexesRight = databaseRight->findClosestSequences(queryRight, numWanted);
 789                 vector<float> scoresLeft = databaseLeft->getSearchScores();
 790                 vector<float> scoresRight = databaseRight->getSearchScores();
 791
 792                 //merge results
 793                 map<int, int> seen;
 794                 map<int, int>::iterator it;
 795                 float lastRight = scoresRight[0];
 796                 float lastLeft = scoresLeft[0];
 797                 //int lasti = 0;
 798                 vector<int> mergedResults;
 799                 for (int i = 0; i < tempIndexesLeft.size(); i++) {
 800                         //add left if you havent already
 801                         it = seen.find(tempIndexesLeft[i]);
 802                         if (it == seen.end()) {
 803                                 mergedResults.push_back(tempIndexesLeft[i]);
 804                                 seen[tempIndexesLeft[i]] = tempIndexesLeft[i];
 805                                 lastLeft = scoresLeft[i];
 806                         }
 807
 808                         //add right if you havent already
 809                         it = seen.find(tempIndexesRight[i]);
 810                         if (it == seen.end()) {
 811                                 mergedResults.push_back(tempIndexesRight[i]);
 812                                 seen[tempIndexesRight[i]] = tempIndexesRight[i];
 813                                 lastRight = scoresRight[i];
 814                         }
 815
 816                         //if (mergedResults.size() > numWanted) { lasti = i; break; } //you have enough results
 817                 }
 818
 819                 //add in sequences with same distance as last sequence added
 820                 /*lasti++;
 821                 int i = lasti;
 822                 while (i < tempIndexesLeft.size()) {
 823                         if (scoresLeft[i] == lastLeft) {
 824                                 it = seen.find(tempIndexesLeft[i]);
 825
 826                                 if (it == seen.end()) {
 827                                         mergedResults.push_back(tempIndexesLeft[i]);
 828                                         seen[tempIndexesLeft[i]] = tempIndexesLeft[i];
 829                                 }
 830                         }
 831                         else { break; }
 832                         i++;
 833                 }
 834
 835                 //              cout << "lastRight\t" << lastRight << endl;
 836                 //add in sequences with same distance as last sequence added
 837                 i = lasti;
 838                 while (i < tempIndexesRight.size()) {
 839                         if (scoresRight[i] == lastRight) {
 840                                 it = seen.find(tempIndexesRight[i]);
 841
 842                                 if (it == seen.end()) {
 843                                         mergedResults.push_back(tempIndexesRight[i]);
 844                                         seen[tempIndexesRight[i]] = tempIndexesRight[i];
 845                                 }
 846                         }
 847                         else { break; }
 848                         i++;
 849                 }*/
 850
 851                 numWanted = mergedResults.size();
 852
 853                 if (numWanted > mergedResults.size()) {
 854                         //m->mothurOut("numwanted is larger than the number of template sequences, adjusting numwanted."); m->mothurOutEndLine();
 855                         numWanted = mergedResults.size();
 856                 }
 857
 858
 859 //cout << q->getName() << endl;
 860                 vector<Sequence*> refResults;
 861                 for (int i = 0; i < numWanted; i++) {
 862 //cout << db[mergedResults[i]]->getName() << endl;
 863                         if (db[mergedResults[i]]->getName() != q->getName()) {
 864                                 Sequence* temp = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned());
 865                                 refResults.push_back(temp);
 866                                 indexes.push_back(mergedResults[i]);
 867                         }
 868                 }
 869 //cout << endl;
 870                 delete queryRight;
 871                 delete queryLeft;
 872
 873                 return refResults;
 874         }
 875         catch(exception& e) {
 876                 m->errorOut(e, "Maligner", "getKmerSeqs");
 877                 exit(1);
 878         }
 879 }
 880 //***************************************************************************************************************
 881