maligner.cpp

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