X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=chimeraslayer.cpp;h=5d714fd1ee88ff18a00619e96ae253847883a151;hb=532c39ca689a9ad7fdff7a81671b2128c71f9671;hp=1336876966f110dc25e0e18eee8bd13ec4e8d6fd;hpb=f1e2f581f91d88e9d6c7a69e495b36031b29e4ab;p=mothur.git diff --git a/chimeraslayer.cpp b/chimeraslayer.cpp index 1336876..5d714fd 100644 --- a/chimeraslayer.cpp +++ b/chimeraslayer.cpp @@ -13,7 +13,7 @@ #include "blastdb.hpp" //*************************************************************************************************************** -ChimeraSlayer::ChimeraSlayer(string file, string temp, string mode, int k, int ms, int mms, int win, float div, +ChimeraSlayer::ChimeraSlayer(string file, string temp, bool trim, string mode, int k, int ms, int mms, int win, float div, int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int numw, bool r) : Chimera() { try { fastafile = file; @@ -33,9 +33,9 @@ int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int num increment = inc; numWanted = numw; realign = r; + trimChimera = trim; + numNoParents = 0; - decalc = new DeCalculator(); - doPrep(); } catch(exception& e) { @@ -44,24 +44,75 @@ int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int num } } //*************************************************************************************************************** -int ChimeraSlayer::doPrep() { +ChimeraSlayer::ChimeraSlayer(string file, string temp, bool trim, map& prior, string mode, int k, int ms, int mms, int win, float div, + int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int numw, bool r) : Chimera() { try { + fastafile = file; templateSeqs = readSeqs(fastafile); + templateFileName = temp; + searchMethod = mode; + kmerSize = k; + match = ms; + misMatch = mms; + window = win; + divR = div; + minSim = minsim; + minCov = mincov; + minBS = minbs; + minSNP = minsnp; + parents = par; + iters = it; + increment = inc; + numWanted = numw; + realign = r; + trimChimera = trim; + priority = prior; + numNoParents = 0; - //read in all query seqs - vector tempQuerySeqs = readSeqs(fastafile); + createFilter(templateSeqs, 0.0); //just removed columns where all seqs have a gap - vector temp = templateSeqs; - for (int i = 0; i < tempQuerySeqs.size(); i++) { temp.push_back(tempQuerySeqs[i]); } + if (searchMethod == "distance") { + createFilter(templateSeqs, 0.0); //just removed columns where all seqs have a gap + + //run filter on template copying templateSeqs into filteredTemplateSeqs + for (int i = 0; i < templateSeqs.size(); i++) { + if (m->control_pressed) { break; } + + Sequence* newSeq = new Sequence(templateSeqs[i]->getName(), templateSeqs[i]->getAligned()); + runFilter(newSeq); + filteredTemplateSeqs.push_back(newSeq); + } + } + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "ChimeraSlayer"); + exit(1); + } +} +//*************************************************************************************************************** +int ChimeraSlayer::doPrep() { + try { + if (searchMethod == "distance") { + //read in all query seqs + vector tempQuerySeqs = readSeqs(fastafile); - createFilter(temp, 0.0); //just removed columns where all seqs have a gap + vector temp = templateSeqs; + for (int i = 0; i < tempQuerySeqs.size(); i++) { temp.push_back(tempQuerySeqs[i]); } - for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i]; } + createFilter(temp, 0.0); //just removed columns where all seqs have a gap - if (m->control_pressed) { return 0; } + for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i]; } - //run filter on template - for (int i = 0; i < templateSeqs.size(); i++) { if (m->control_pressed) { return 0; } runFilter(templateSeqs[i]); } + if (m->control_pressed) { return 0; } + //run filter on template copying templateSeqs into filteredTemplateSeqs + for (int i = 0; i < templateSeqs.size(); i++) { + if (m->control_pressed) { return 0; } + + Sequence* newSeq = new Sequence(templateSeqs[i]->getName(), templateSeqs[i]->getAligned()); + runFilter(newSeq); + filteredTemplateSeqs.push_back(newSeq); + } + } string kmerDBNameLeft; string kmerDBNameRight; @@ -164,7 +215,10 @@ int ChimeraSlayer::doPrep() { }else if (searchMethod == "blast") { //generate blastdb - databaseLeft = new BlastDB(-2.0, -1.0, match, misMatch); + databaseLeft = new BlastDB(m->getRootName(m->getSimpleName(fastafile)), -1.0, -1.0, 1, -3); + + if (m->control_pressed) { return 0; } + for (int i = 0; i < templateSeqs.size(); i++) { databaseLeft->addSequence(*templateSeqs[i]); } databaseLeft->generateDB(); databaseLeft->setNumSeqs(templateSeqs.size()); @@ -178,11 +232,122 @@ int ChimeraSlayer::doPrep() { exit(1); } } +//*************************************************************************************************************** +vector ChimeraSlayer::getTemplate(Sequence q, vector& userTemplateFiltered) { + try { + + //when template=self, the query file is sorted from most abundance to least abundant + //userTemplate grows as the query file is processed by adding sequences that are not chimeric and more abundant + vector userTemplate; + + int myAbund = priority[q.getName()]; + + for (int i = 0; i < templateSeqs.size(); i++) { + + if (m->control_pressed) { return userTemplate; } + + //have I reached a sequence with the same abundance as myself? + if (!(priority[templateSeqs[i]->getName()] > myAbund)) { break; } + + //if its am not chimeric add it + if (chimericSeqs.count(templateSeqs[i]->getName()) == 0) { + userTemplate.push_back(templateSeqs[i]); + if (searchMethod == "distance") { userTemplateFiltered.push_back(filteredTemplateSeqs[i]); } + } + } + + string kmerDBNameLeft; + string kmerDBNameRight; + + //generate the kmerdb to pass to maligner + if (searchMethod == "kmer") { + string templatePath = m->hasPath(templateFileName); + string rightTemplateFileName = templatePath + "right." + m->getRootName(m->getSimpleName(templateFileName)); + databaseRight = new KmerDB(rightTemplateFileName, kmerSize); + + string leftTemplateFileName = templatePath + "left." + m->getRootName(m->getSimpleName(templateFileName)); + databaseLeft = new KmerDB(leftTemplateFileName, kmerSize); +#ifdef USE_MPI + for (int i = 0; i < userTemplate.size(); i++) { + + if (m->control_pressed) { return userTemplate; } + + string leftFrag = userTemplate[i]->getUnaligned(); + leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33)); + + Sequence leftTemp(userTemplate[i]->getName(), leftFrag); + databaseLeft->addSequence(leftTemp); + } + databaseLeft->generateDB(); + databaseLeft->setNumSeqs(userTemplate.size()); + + for (int i = 0; i < userTemplate.size(); i++) { + if (m->control_pressed) { return userTemplate; } + + string rightFrag = userTemplate[i]->getUnaligned(); + rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66)); + + Sequence rightTemp(userTemplate[i]->getName(), rightFrag); + databaseRight->addSequence(rightTemp); + } + databaseRight->generateDB(); + databaseRight->setNumSeqs(userTemplate.size()); + +#else + + + for (int i = 0; i < userTemplate.size(); i++) { + + if (m->control_pressed) { return userTemplate; } + + string leftFrag = userTemplate[i]->getUnaligned(); + leftFrag = leftFrag.substr(0, int(leftFrag.length() * 0.33)); + + Sequence leftTemp(userTemplate[i]->getName(), leftFrag); + databaseLeft->addSequence(leftTemp); + } + databaseLeft->generateDB(); + databaseLeft->setNumSeqs(userTemplate.size()); + + for (int i = 0; i < userTemplate.size(); i++) { + if (m->control_pressed) { return userTemplate; } + + string rightFrag = userTemplate[i]->getUnaligned(); + rightFrag = rightFrag.substr(int(rightFrag.length() * 0.66)); + + Sequence rightTemp(userTemplate[i]->getName(), rightFrag); + databaseRight->addSequence(rightTemp); + } + databaseRight->generateDB(); + databaseRight->setNumSeqs(userTemplate.size()); +#endif + }else if (searchMethod == "blast") { + + //generate blastdb + databaseLeft = new BlastDB(m->getRootName(m->getSimpleName(templateFileName)), -1.0, -1.0, 1, -3); + + if (m->control_pressed) { return userTemplate; } + + for (int i = 0; i < userTemplate.size(); i++) { if (m->control_pressed) { return userTemplate; } databaseLeft->addSequence(*userTemplate[i]); } + databaseLeft->generateDB(); + databaseLeft->setNumSeqs(userTemplate.size()); + } + + return userTemplate; + + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "getTemplate"); + exit(1); + } +} + //*************************************************************************************************************** ChimeraSlayer::~ChimeraSlayer() { - delete decalc; - if (searchMethod == "kmer") { delete databaseRight; delete databaseLeft; } - else if (searchMethod == "blast") { delete databaseLeft; } + if (templateFileName != "self") { + if (searchMethod == "kmer") { delete databaseRight; delete databaseLeft; } + else if (searchMethod == "blast") { delete databaseLeft; } + } } //*************************************************************************************************************** void ChimeraSlayer::printHeader(ostream& out) { @@ -193,8 +358,11 @@ void ChimeraSlayer::printHeader(ostream& out) { out << "Name\tLeftParent\tRightParent\tDivQLAQRB\tPerIDQLAQRB\tBootStrapA\tDivQLBQRA\tPerIDQLBQRA\tBootStrapB\tFlag\tLeftWindow\tRightWindow\n"; } //*************************************************************************************************************** -int ChimeraSlayer::print(ostream& out, ostream& outAcc) { +Sequence ChimeraSlayer::print(ostream& out, ostream& outAcc) { try { + Sequence trim; + if (trimChimera) { trim.setName(trimQuery.getName()); trim.setAligned(trimQuery.getAligned()); } + if (chimeraFlags == "yes") { string chimeraFlag = "no"; if( (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR) @@ -204,16 +372,137 @@ int ChimeraSlayer::print(ostream& out, ostream& outAcc) { if (chimeraFlag == "yes") { if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) { - m->mothurOut(querySeq->getName() + "\tyes"); m->mothurOutEndLine(); - outAcc << querySeq->getName() << endl; + m->mothurOut(querySeq.getName() + "\tyes"); m->mothurOutEndLine(); + outAcc << querySeq.getName() << endl; + + if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); } + + if (trimChimera) { + int lengthLeft = chimeraResults[0].winLEnd - chimeraResults[0].winLStart; + int lengthRight = chimeraResults[0].winREnd - chimeraResults[0].winRStart; + + string newAligned = trim.getAligned(); + + if (lengthLeft > lengthRight) { //trim right + for (int i = (chimeraResults[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else { //trim left + for (int i = 0; i < chimeraResults[0].winLEnd; i++) { newAligned[i] = '.'; } + } + trim.setAligned(newAligned); + } } } printBlock(chimeraResults[0], chimeraFlag, out); out << endl; - }else { out << querySeq->getName() << "\tno" << endl; } + }else { + out << querySeq.getName() << "\tno" << endl; + } - return 0; + return trim; + + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "print"); + exit(1); + } +} +//*************************************************************************************************************** +Sequence ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftPiece, data_results rightPiece) { + try { + Sequence trim; + + if (trimChimera) { + string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned(); + trim.setName(leftPiece.trimQuery.getName()); trim.setAligned(aligned); + } + + if ((leftPiece.flag == "yes") || (rightPiece.flag == "yes")) { + + string chimeraFlag = "no"; + if (leftPiece.flag == "yes") { + + if( (leftPiece.results[0].bsa >= minBS && leftPiece.results[0].divr_qla_qrb >= divR) + || + (leftPiece.results[0].bsb >= minBS && leftPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; } + } + + if (rightPiece.flag == "yes") { + if ( (rightPiece.results[0].bsa >= minBS && rightPiece.results[0].divr_qla_qrb >= divR) + || + (rightPiece.results[0].bsb >= minBS && rightPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; } + } + + bool rightChimeric = false; + bool leftChimeric = false; + + if (chimeraFlag == "yes") { + //which peice is chimeric or are both + if (rightPiece.flag == "yes") { if ((rightPiece.results[0].bsa >= minBS) || (rightPiece.results[0].bsb >= minBS)) { rightChimeric = true; } } + if (leftPiece.flag == "yes") { if ((leftPiece.results[0].bsa >= minBS) || (leftPiece.results[0].bsb >= minBS)) { leftChimeric = true; } } + + if (rightChimeric || leftChimeric) { + m->mothurOut(querySeq.getName() + "\tyes"); m->mothurOutEndLine(); + outAcc << querySeq.getName() << endl; + + if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); } + + if (trimChimera) { + string newAligned = trim.getAligned(); + + //right side is fine so keep that + if ((leftChimeric) && (!rightChimeric)) { + for (int i = 0; i < leftPiece.results[0].winREnd; i++) { newAligned[i] = '.'; } + }else if ((!leftChimeric) && (rightChimeric)) { //leftside is fine so keep that + for (int i = (rightPiece.results[0].winLStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else { //both sides are chimeric, keep longest piece + + int lengthLeftLeft = leftPiece.results[0].winLEnd - leftPiece.results[0].winLStart; + int lengthLeftRight = leftPiece.results[0].winREnd - leftPiece.results[0].winRStart; + + int longest = 1; // leftleft = 1, leftright = 2, rightleft = 3 rightright = 4 + int length = lengthLeftLeft; + if (lengthLeftLeft < lengthLeftRight) { longest = 2; length = lengthLeftRight; } + + int lengthRightLeft = rightPiece.results[0].winLEnd - rightPiece.results[0].winLStart; + int lengthRightRight = rightPiece.results[0].winREnd - rightPiece.results[0].winRStart; + + if (lengthRightLeft > length) { longest = 3; length = lengthRightLeft; } + if (lengthRightRight > length) { longest = 4; } + + if (longest == 1) { //leftleft + for (int i = (leftPiece.results[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else if (longest == 2) { //leftright + //get rid of leftleft + for (int i = (leftPiece.results[0].winLStart-1); i < (leftPiece.results[0].winLEnd-1); i++) { newAligned[i] = '.'; } + //get rid of right + for (int i = (rightPiece.results[0].winLStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else if (longest == 3) { //rightleft + //get rid of left + for (int i = 0; i < leftPiece.results[0].winREnd; i++) { newAligned[i] = '.'; } + //get rid of rightright + for (int i = (rightPiece.results[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else { //rightright + //get rid of left + for (int i = 0; i < leftPiece.results[0].winREnd; i++) { newAligned[i] = '.'; } + //get rid of rightleft + for (int i = (rightPiece.results[0].winLStart-1); i < (rightPiece.results[0].winLEnd-1); i++) { newAligned[i] = '.'; } + } + } + + trim.setAligned(newAligned); + } + + } + } + + printBlock(leftPiece, rightPiece, leftChimeric, rightChimeric, chimeraFlag, out); + out << endl; + }else { + out << querySeq.getName() << "\tno" << endl; + } + + return trim; } catch(exception& e) { @@ -221,15 +510,156 @@ int ChimeraSlayer::print(ostream& out, ostream& outAcc) { exit(1); } } + #ifdef USE_MPI //*************************************************************************************************************** -int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { +Sequence ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results leftPiece, data_results rightPiece) { try { MPI_Status status; bool results = false; string outAccString = ""; string outputString = ""; + Sequence trim; + + if (trimChimera) { + string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned(); + trim.setName(leftPiece.trimQuery.getName()); trim.setAligned(aligned); + } + + + if ((leftPiece.flag == "yes") || (rightPiece.flag == "yes")) { + + string chimeraFlag = "no"; + if (leftPiece.flag == "yes") { + + if( (leftPiece.results[0].bsa >= minBS && leftPiece.results[0].divr_qla_qrb >= divR) + || + (leftPiece.results[0].bsb >= minBS && leftPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; } + } + + if (rightPiece.flag == "yes") { + if ( (rightPiece.results[0].bsa >= minBS && rightPiece.results[0].divr_qla_qrb >= divR) + || + (rightPiece.results[0].bsb >= minBS && rightPiece.results[0].divr_qlb_qra >= divR) ) { chimeraFlag = "yes"; } + } + + bool rightChimeric = false; + bool leftChimeric = false; + + cout << endl; + + if (chimeraFlag == "yes") { + //which peice is chimeric or are both + if (rightPiece.flag == "yes") { if ((rightPiece.results[0].bsa >= minBS) || (rightPiece.results[0].bsb >= minBS)) { rightChimeric = true; } } + if (leftPiece.flag == "yes") { if ((leftPiece.results[0].bsa >= minBS) || (leftPiece.results[0].bsb >= minBS)) { leftChimeric = true; } } + + if (rightChimeric || leftChimeric) { + cout << querySeq.getName() << "\tyes" << endl; + outAccString += querySeq.getName() + "\n"; + results = true; + + if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); } + + //write to accnos file + int length = outAccString.length(); + char* buf2 = new char[length]; + memcpy(buf2, outAccString.c_str(), length); + + MPI_File_write_shared(outAcc, buf2, length, MPI_CHAR, &status); + delete buf2; + + if (trimChimera) { + string newAligned = trim.getAligned(); + + //right side is fine so keep that + if ((leftChimeric) && (!rightChimeric)) { + for (int i = 0; i < leftPiece.results[0].winREnd; i++) { newAligned[i] = '.'; } + }else if ((!leftChimeric) && (rightChimeric)) { //leftside is fine so keep that + for (int i = (rightPiece.results[0].winLStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else { //both sides are chimeric, keep longest piece + + int lengthLeftLeft = leftPiece.results[0].winLEnd - leftPiece.results[0].winLStart; + int lengthLeftRight = leftPiece.results[0].winREnd - leftPiece.results[0].winRStart; + + int longest = 1; // leftleft = 1, leftright = 2, rightleft = 3 rightright = 4 + int length = lengthLeftLeft; + if (lengthLeftLeft < lengthLeftRight) { longest = 2; length = lengthLeftRight; } + + int lengthRightLeft = rightPiece.results[0].winLEnd - rightPiece.results[0].winLStart; + int lengthRightRight = rightPiece.results[0].winREnd - rightPiece.results[0].winRStart; + + if (lengthRightLeft > length) { longest = 3; length = lengthRightLeft; } + if (lengthRightRight > length) { longest = 4; } + + if (longest == 1) { //leftleft + for (int i = (leftPiece.results[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else if (longest == 2) { //leftright + //get rid of leftleft + for (int i = (leftPiece.results[0].winLStart-1); i < (leftPiece.results[0].winLEnd-1); i++) { newAligned[i] = '.'; } + //get rid of right + for (int i = (rightPiece.results[0].winLStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else if (longest == 3) { //rightleft + //get rid of left + for (int i = 0; i < leftPiece.results[0].winREnd; i++) { newAligned[i] = '.'; } + //get rid of rightright + for (int i = (rightPiece.results[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else { //rightright + //get rid of left + for (int i = 0; i < leftPiece.results[0].winREnd; i++) { newAligned[i] = '.'; } + //get rid of rightleft + for (int i = (rightPiece.results[0].winLStart-1); i < (rightPiece.results[0].winLEnd-1); i++) { newAligned[i] = '.'; } + } + } + + trim.setAligned(newAligned); + } + + } + } + + outputString = getBlock(leftPiece, rightPiece, leftChimeric, rightChimeric, chimeraFlag); + outputString += "\n"; + + //write to output file + int length = outputString.length(); + char* buf = new char[length]; + memcpy(buf, outputString.c_str(), length); + + MPI_File_write_shared(out, buf, length, MPI_CHAR, &status); + delete buf; + + }else { + outputString += querySeq.getName() + "\tno\n"; + + //write to output file + int length = outputString.length(); + char* buf = new char[length]; + memcpy(buf, outputString.c_str(), length); + + MPI_File_write_shared(out, buf, length, MPI_CHAR, &status); + delete buf; + } + + + return trim; + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "print"); + exit(1); + } +} +//*************************************************************************************************************** +Sequence ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { + try { + MPI_Status status; + bool results = false; + string outAccString = ""; + string outputString = ""; + + Sequence trim; + if (trimChimera) { trim.setName(trimQuery.getName()); trim.setAligned(trimQuery.getAligned()); } + if (chimeraFlags == "yes") { string chimeraFlag = "no"; if( (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR) @@ -239,45 +669,59 @@ int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { if (chimeraFlag == "yes") { if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) { - cout << querySeq->getName() << "\tyes" << endl; - outAccString += querySeq->getName() + "\n"; + cout << querySeq.getName() << "\tyes" << endl; + outAccString += querySeq.getName() + "\n"; results = true; + if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); } + //write to accnos file int length = outAccString.length(); char* buf2 = new char[length]; memcpy(buf2, outAccString.c_str(), length); - + MPI_File_write_shared(outAcc, buf2, length, MPI_CHAR, &status); delete buf2; + + if (trimChimera) { + int lengthLeft = chimeraResults[0].winLEnd - chimeraResults[0].winLStart; + int lengthRight = chimeraResults[0].winREnd - chimeraResults[0].winRStart; + + string newAligned = trim.getAligned(); + if (lengthLeft > lengthRight) { //trim right + for (int i = (chimeraResults[0].winRStart-1); i < newAligned.length(); i++) { newAligned[i] = '.'; } + }else { //trim left + for (int i = 0; i < (chimeraResults[0].winLEnd-1); i++) { newAligned[i] = '.'; } + } + trim.setAligned(newAligned); + } } } outputString = getBlock(chimeraResults[0], chimeraFlag); outputString += "\n"; - //cout << outputString << endl; + //write to output file int length = outputString.length(); char* buf = new char[length]; memcpy(buf, outputString.c_str(), length); - + MPI_File_write_shared(out, buf, length, MPI_CHAR, &status); delete buf; - + }else { - outputString += querySeq->getName() + "\tno\n"; - //cout << outputString << endl; + outputString += querySeq.getName() + "\tno\n"; + //write to output file int length = outputString.length(); char* buf = new char[length]; memcpy(buf, outputString.c_str(), length); - + MPI_File_write_shared(out, buf, length, MPI_CHAR, &status); delete buf; } - - return results; + return trim; } catch(exception& e) { m->errorOut(e, "ChimeraSlayer", "print"); @@ -289,109 +733,135 @@ int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { //*************************************************************************************************************** int ChimeraSlayer::getChimeras(Sequence* query) { try { + + trimQuery.setName(query->getName()); trimQuery.setAligned(query->getAligned()); + printResults.trimQuery = trimQuery; + chimeraFlags = "no"; - - //filter query - spotMap = runFilter(query); + printResults.flag = "no"; - querySeq = query; + querySeq = *query; + + //you must create a template + vector thisTemplate; + vector thisFilteredTemplate; + if (templateFileName != "self") { thisTemplate = templateSeqs; thisFilteredTemplate = filteredTemplateSeqs; } + else { thisTemplate = getTemplate(*query, thisFilteredTemplate); } //fills this template and creates the databases - //referenceSeqs, numWanted, matchScore, misMatchPenalty, divR, minSimilarity - Maligner maligner(templateSeqs, numWanted, match, misMatch, divR, minSim, minCov, searchMethod, databaseLeft, databaseRight); - Slayer slayer(window, increment, minSim, divR, iters, minSNP); - if (m->control_pressed) { return 0; } - string chimeraFlag = maligner.getResults(query, decalc); + if (thisTemplate.size() == 0) { return 0; } //not chimeric + + //moved this out of maligner - 4/29/11 + vector refSeqs = getRefSeqs(*query, thisTemplate, thisFilteredTemplate); + + Maligner maligner(refSeqs, match, misMatch, divR, minSim, minCov); + Slayer slayer(window, increment, minSim, divR, iters, minSNP, minBS); + + if (templateFileName == "self") { + if (searchMethod == "kmer") { delete databaseRight; delete databaseLeft; } + else if (searchMethod == "blast") { delete databaseLeft; } + } + + if (m->control_pressed) { return 0; } + + string chimeraFlag = maligner.getResults(*query, decalc); + if (m->control_pressed) { return 0; } + vector Results = maligner.getOutput(); + + //for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } + + if (chimeraFlag == "yes") { - //found in testing realigning only made things worse - if (realign) { - ChimeraReAligner realigner(templateSeqs, match, misMatch); - realigner.reAlign(query, Results); - } + if (realign) { + vector parents; + for (int i = 0; i < Results.size(); i++) { + parents.push_back(Results[i].parentAligned); + } + + ChimeraReAligner realigner; + realigner.reAlign(query, parents); - if (chimeraFlag == "yes") { + } +// cout << query->getAligned() << endl; //get sequence that were given from maligner results - vector seqs; + vector seqs; map removeDups; map::iterator itDup; map parentNameSeq; map::iterator itSeq; for (int j = 0; j < Results.size(); j++) { + float dist = (Results[j].regionEnd - Results[j].regionStart + 1) * Results[j].queryToParentLocal; //only add if you are not a duplicate - itDup = removeDups.find(Results[j].parent); - if (itDup == removeDups.end()) { //this is not duplicate - removeDups[Results[j].parent] = dist; - parentNameSeq[Results[j].parent] = Results[j].parentAligned; - }else if (dist > itDup->second) { //is this a stronger number for this parent - removeDups[Results[j].parent] = dist; - parentNameSeq[Results[j].parent] = Results[j].parentAligned; +// cout << Results[j].parent << '\t' << Results[j].regionEnd << '\t' << Results[j].regionStart << '\t' << Results[j].regionEnd - Results[j].regionStart +1 << '\t' << Results[j].queryToParentLocal << '\t' << dist << endl; + + + if(Results[j].queryToParentLocal >= 90){ //local match has to be over 90% similarity + + itDup = removeDups.find(Results[j].parent); + if (itDup == removeDups.end()) { //this is not duplicate + removeDups[Results[j].parent] = dist; + parentNameSeq[Results[j].parent] = Results[j].parentAligned; + }else if (dist > itDup->second) { //is this a stronger number for this parent + removeDups[Results[j].parent] = dist; + parentNameSeq[Results[j].parent] = Results[j].parentAligned; + } + } + } for (itDup = removeDups.begin(); itDup != removeDups.end(); itDup++) { - //Sequence* seq = getSequence(itDup->first); //makes copy so you can filter and mask and not effect template itSeq = parentNameSeq.find(itDup->first); -//cout << itDup->first << itSeq->second << endl; - Sequence* seq = new Sequence(itDup->first, itSeq->second); + Sequence seq(itDup->first, itSeq->second); - SeqDist member; + SeqCompare member; member.seq = seq; member.dist = itDup->second; - seqs.push_back(member); } //limit number of parents to explore - default 3 if (Results.size() > parents) { //sort by distance - sort(seqs.begin(), seqs.end(), compareSeqDist); + sort(seqs.begin(), seqs.end(), compareSeqCompare); //prioritize larger more similiar sequence fragments reverse(seqs.begin(), seqs.end()); - for (int k = seqs.size()-1; k > (parents-1); k--) { - delete seqs[k].seq; - seqs.pop_back(); - } + //for (int k = seqs.size()-1; k > (parents-1); k--) { + // delete seqs[k].seq; + //seqs.pop_back(); + //} } - + //put seqs into vector to send to slayer - vector seqsForSlayer; - for (int k = 0; k < seqs.size(); k++) { seqsForSlayer.push_back(seqs[k].seq); } - //mask then send to slayer... - if (seqMask != "") { - decalc->setMask(seqMask); - - //mask querys - decalc->runMask(query); - - //mask parents - for (int k = 0; k < seqsForSlayer.size(); k++) { - decalc->runMask(seqsForSlayer[k]); - } - - spotMap = decalc->getMaskMap(); +// cout << query->getAligned() << endl; + vector seqsForSlayer; + for (int k = 0; k < seqs.size(); k++) { +// cout << seqs[k].seq->getAligned() << endl; + seqsForSlayer.push_back(seqs[k].seq); +// cout << seqs[k].seq->getName() << endl; } - if (m->control_pressed) { for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; } return 0; } - + if (m->control_pressed) { return 0; } + //send to slayer - chimeraFlags = slayer.getResults(query, seqsForSlayer); + chimeraFlags = slayer.getResults(*query, seqsForSlayer); if (m->control_pressed) { return 0; } chimeraResults = slayer.getOutput(); + printResults.flag = chimeraFlags; + printResults.results = chimeraResults; + //free memory - for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; } + //for (int k = 0; k < seqs.size(); k++) { delete seqs[k].seq; } } - - //delete maligner; - //delete slayer; - + //cout << endl << endl; return 0; } catch(exception& e) { @@ -402,32 +872,66 @@ int ChimeraSlayer::getChimeras(Sequence* query) { //*************************************************************************************************************** void ChimeraSlayer::printBlock(data_struct data, string flag, ostream& out){ try { - //out << ":)\n"; - - out << querySeq->getName() << '\t'; + out << querySeq.getName() << '\t'; out << data.parentA.getName() << "\t" << data.parentB.getName() << '\t'; - //out << "Left Window: " << spotMap[data.winLStart] << " " << spotMap[data.winLEnd] << endl; - //out << "Right Window: " << spotMap[data.winRStart] << " " << spotMap[data.winREnd] << endl; - + out << data.divr_qla_qrb << '\t' << data.qla_qrb << '\t' << data.bsa << '\t'; out << data.divr_qlb_qra << '\t' << data.qlb_qra << '\t' << data.bsb << '\t'; - out << flag << '\t' << spotMap[data.winLStart] << "-" << spotMap[data.winLEnd] << '\t' << spotMap[data.winRStart] << "-" << spotMap[data.winREnd] << '\t'; + out << flag << '\t' << data.winLStart << "-" << data.winLEnd << '\t' << data.winRStart << "-" << data.winREnd << '\t'; - //out << "Similarity of parents: " << data.ab << endl; - //out << "Similarity of query to parentA: " << data.qa << endl; - //out << "Similarity of query to parentB: " << data.qb << endl; + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "printBlock"); + exit(1); + } +} +//*************************************************************************************************************** +void ChimeraSlayer::printBlock(data_results leftdata, data_results rightdata, bool leftChimeric, bool rightChimeric, string flag, ostream& out){ + try { + if ((leftChimeric) && (!rightChimeric)) { //print left + out << querySeq.getName() << '\t'; + out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName() << '\t'; + + out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t'; + out << leftdata.results[0].divr_qlb_qra << '\t' << leftdata.results[0].qlb_qra << '\t' << leftdata.results[0].bsb << '\t'; - //out << "Per_id(QL,A): " << data.qla << endl; - //out << "Per_id(QL,B): " << data.qlb << endl; - //out << "Per_id(QR,A): " << data.qra << endl; - //out << "Per_id(QR,B): " << data.qrb << endl; - + out << flag << '\t' << leftdata.results[0].winLStart << "-" << leftdata.results[0].winLEnd << '\t' << leftdata.results[0].winRStart << "-" << leftdata.results[0].winREnd << '\t'; - //out << "DeltaL: " << (data.qla - data.qlb) << endl; - //out << "DeltaR: " << (data.qra - data.qrb) << endl; - + }else if ((!leftChimeric) && (rightChimeric)) { //print right + out << querySeq.getName() << '\t'; + out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName() << '\t'; + + out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t'; + out << rightdata.results[0].divr_qlb_qra << '\t' << rightdata.results[0].qlb_qra << '\t' << rightdata.results[0].bsb << '\t'; + + out << flag << '\t' << rightdata.results[0].winLStart << "-" << rightdata.results[0].winLEnd << '\t' << rightdata.results[0].winRStart << "-" << rightdata.results[0].winREnd << '\t'; + + }else { //print both results + if (leftdata.flag == "yes") { + out << querySeq.getName() + "_LEFT" << '\t'; + out << leftdata.results[0].parentA.getName() << "\t" << leftdata.results[0].parentB.getName() << '\t'; + + out << leftdata.results[0].divr_qla_qrb << '\t' << leftdata.results[0].qla_qrb << '\t' << leftdata.results[0].bsa << '\t'; + out << leftdata.results[0].divr_qlb_qra << '\t' << leftdata.results[0].qlb_qra << '\t' << leftdata.results[0].bsb << '\t'; + + out << flag << '\t' << leftdata.results[0].winLStart << "-" << leftdata.results[0].winLEnd << '\t' << leftdata.results[0].winRStart << "-" << leftdata.results[0].winREnd << '\t'; + } + + if (rightdata.flag == "yes") { + if (leftdata.flag == "yes") { out << endl; } + + out << querySeq.getName() + "_RIGHT"<< '\t'; + out << rightdata.results[0].parentA.getName() << "\t" << rightdata.results[0].parentB.getName() << '\t'; + + out << rightdata.results[0].divr_qla_qrb << '\t' << rightdata.results[0].qla_qrb << '\t' << rightdata.results[0].bsa << '\t'; + out << rightdata.results[0].divr_qlb_qra << '\t' << rightdata.results[0].qlb_qra << '\t' << rightdata.results[0].bsb << '\t'; + + out << flag << '\t' << rightdata.results[0].winLStart << "-" << rightdata.results[0].winLEnd << '\t' << rightdata.results[0].winRStart << "-" << rightdata.results[0].winREnd << '\t'; + + } + } } catch(exception& e) { m->errorOut(e, "ChimeraSlayer", "printBlock"); @@ -435,18 +939,74 @@ void ChimeraSlayer::printBlock(data_struct data, string flag, ostream& out){ } } //*************************************************************************************************************** +string ChimeraSlayer::getBlock(data_results leftdata, data_results rightdata, bool leftChimeric, bool rightChimeric, string flag){ + try { + + string out = ""; + + if ((leftChimeric) && (!rightChimeric)) { //get left + out += querySeq.getName() + "\t"; + out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t"; + + out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t"; + out += toString(leftdata.results[0].divr_qlb_qra) + "\t" + toString(leftdata.results[0].qlb_qra) + "\t" + toString(leftdata.results[0].bsb) + "\t"; + + out += flag + "\t" + toString(leftdata.results[0].winLStart) + "-" + toString(leftdata.results[0].winLEnd) + "\t" + toString(leftdata.results[0].winRStart) + "-" + toString(leftdata.results[0].winREnd) + "\t"; + + }else if ((!leftChimeric) && (rightChimeric)) { //print right + out += querySeq.getName() + "\t"; + out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName() + "\t"; + + out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t"; + out += toString(rightdata.results[0].divr_qlb_qra) + "\t" + toString(rightdata.results[0].qlb_qra) + "\t" + toString(rightdata.results[0].bsb) + "\t"; + + out += flag + "\t" + toString(rightdata.results[0].winLStart) + "-" + toString(rightdata.results[0].winLEnd) + "\t" + toString(rightdata.results[0].winRStart) + "-" + toString(rightdata.results[0].winREnd) + "\t"; + + }else { //print both results + + if (leftdata.flag == "yes") { + out += querySeq.getName() + "_LEFT\t"; + out += leftdata.results[0].parentA.getName() + "\t" + leftdata.results[0].parentB.getName() + "\t"; + + out += toString(leftdata.results[0].divr_qla_qrb) + "\t" + toString(leftdata.results[0].qla_qrb) + "\t" + toString(leftdata.results[0].bsa) + "\t"; + out += toString(leftdata.results[0].divr_qlb_qra) + "\t" + toString(leftdata.results[0].qlb_qra) + "\t" + toString(leftdata.results[0].bsb) + "\t"; + + out += flag + "\t" + toString(leftdata.results[0].winLStart) + "-" + toString(leftdata.results[0].winLEnd) + "\t" + toString(leftdata.results[0].winRStart) + "-" + toString(leftdata.results[0].winREnd) + "\t"; + } + + if (rightdata.flag == "yes") { + if (leftdata.flag == "yes") { out += "\n"; } + out += querySeq.getName() + "_RIGHT\t"; + out += rightdata.results[0].parentA.getName() + "\t" + rightdata.results[0].parentB.getName() + "\t"; + + out += toString(rightdata.results[0].divr_qla_qrb) + "\t" + toString(rightdata.results[0].qla_qrb) + "\t" + toString(rightdata.results[0].bsa) + "\t"; + out += toString(rightdata.results[0].divr_qlb_qra) + "\t" + toString(rightdata.results[0].qlb_qra) + "\t" + toString(rightdata.results[0].bsb) + "\t"; + + out += flag + "\t" + toString(rightdata.results[0].winLStart) + "-" + toString(rightdata.results[0].winLEnd) + "\t" + toString(rightdata.results[0].winRStart) + "-" + toString(rightdata.results[0].winREnd) + "\t"; + } + } + + return out; + + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "getBlock"); + exit(1); + } +} +//*************************************************************************************************************** string ChimeraSlayer::getBlock(data_struct data, string flag){ try { string outputString = ""; - outputString += querySeq->getName() + "\t"; + outputString += querySeq.getName() + "\t"; outputString += data.parentA.getName() + "\t" + data.parentB.getName() + "\t"; outputString += toString(data.divr_qla_qrb) + "\t" + toString(data.qla_qrb) + "\t" + toString(data.bsa) + "\t"; outputString += toString(data.divr_qlb_qra) + "\t" + toString(data.qlb_qra) + "\t" + toString(data.bsb) + "\t"; - outputString += flag + "\t" + toString(spotMap[data.winLStart]) + "-" + toString(spotMap[data.winLEnd]) + "\t" + toString(spotMap[data.winRStart]) + "-" + toString(spotMap[data.winREnd]) + "\t"; + outputString += flag + "\t" + toString(data.winLStart) + "-" + toString(data.winLEnd) + "\t" + toString(data.winRStart) + "-" + toString(data.winREnd) + "\t"; return outputString; } @@ -455,5 +1015,215 @@ string ChimeraSlayer::getBlock(data_struct data, string flag){ exit(1); } } +//*************************************************************************************************************** +vector ChimeraSlayer::getRefSeqs(Sequence q, vector& thisTemplate, vector& thisFilteredTemplate){ + try { + + vector refSeqs; + + if (searchMethod == "distance") { + //find closest seqs to query in template - returns copies of seqs so trim does not destroy - remember to deallocate + Sequence* newSeq = new Sequence(q.getName(), q.getAligned()); + runFilter(newSeq); + refSeqs = decalc.findClosest(*newSeq, thisTemplate, thisFilteredTemplate, numWanted, minSim); + delete newSeq; + }else if (searchMethod == "blast") { + refSeqs = getBlastSeqs(q, thisTemplate, numWanted); //fills indexes + }else if (searchMethod == "kmer") { + refSeqs = getKmerSeqs(q, thisTemplate, numWanted); //fills indexes + }else { m->mothurOut("not valid search."); exit(1); } //should never get here + + return refSeqs; + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "getRefSeqs"); + exit(1); + } +} //***************************************************************************************************************/ +vector ChimeraSlayer::getBlastSeqs(Sequence q, vector& db, int num) { + try { + + vector refResults; + + //get parts of query + string queryUnAligned = q.getUnaligned(); + string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence + string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence +//cout << "whole length = " << queryUnAligned.length() << '\t' << "left length = " << leftQuery.length() << '\t' << "right length = "<< rightQuery.length() << endl; + Sequence* queryLeft = new Sequence(q.getName(), leftQuery); + Sequence* queryRight = new Sequence(q.getName(), rightQuery); + + vector tempIndexesLeft = databaseLeft->findClosestMegaBlast(queryLeft, num+1, minSim); + vector tempIndexesRight = databaseLeft->findClosestMegaBlast(queryRight, num+1, minSim); + + + //cout << q->getName() << '\t' << leftQuery << '\t' << "leftMatches = " << tempIndexesLeft.size() << '\t' << rightQuery << " rightMatches = " << tempIndexesRight.size() << endl; +// vector smaller; +// vector larger; +// +// if (tempIndexesRight.size() < tempIndexesLeft.size()) { smaller = tempIndexesRight; larger = tempIndexesLeft; } +// else { smaller = tempIndexesLeft; larger = tempIndexesRight; } + + //merge results + map seen; + map::iterator it; + vector mergedResults; + + int index = 0; +// for (int i = 0; i < smaller.size(); i++) { + while(index < tempIndexesLeft.size() && index < tempIndexesRight.size()){ + + if (m->control_pressed) { delete queryRight; delete queryLeft; return refResults; } + + //add left if you havent already + it = seen.find(tempIndexesLeft[index]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesLeft[index]); + seen[tempIndexesLeft[index]] = tempIndexesLeft[index]; + } + + //add right if you havent already + it = seen.find(tempIndexesRight[index]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesRight[index]); + seen[tempIndexesRight[index]] = tempIndexesRight[index]; + } + index++; + } + + + for (int i = index; i < tempIndexesLeft.size(); i++) { + if (m->control_pressed) { delete queryRight; delete queryLeft; return refResults; } + + //add right if you havent already + it = seen.find(tempIndexesLeft[i]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesLeft[i]); + seen[tempIndexesLeft[i]] = tempIndexesLeft[i]; + } + } + + for (int i = index; i < tempIndexesRight.size(); i++) { + if (m->control_pressed) { delete queryRight; delete queryLeft; return refResults; } + + //add right if you havent already + it = seen.find(tempIndexesRight[i]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesRight[i]); + seen[tempIndexesRight[i]] = tempIndexesRight[i]; + } + } + //string qname = q->getName().substr(0, q->getName().find_last_of('_')); + //cout << qname << endl; + + if (mergedResults.size() == 0) { numNoParents++; } + + for (int i = 0; i < mergedResults.size(); i++) { + //cout << q->getName() << mergedResults[i] << '\t' << db[mergedResults[i]]->getName() << endl; + if (db[mergedResults[i]]->getName() != q.getName()) { + Sequence temp(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); + refResults.push_back(temp); + } + } + //cout << endl << endl; + + delete queryRight; + delete queryLeft; + + return refResults; + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "getBlastSeqs"); + exit(1); + } +} +//*************************************************************************************************************** +vector ChimeraSlayer::getKmerSeqs(Sequence q, vector& db, int num) { + try { + vector refResults; + + //get parts of query + string queryUnAligned = q.getUnaligned(); + string leftQuery = queryUnAligned.substr(0, int(queryUnAligned.length() * 0.33)); //first 1/3 of the sequence + string rightQuery = queryUnAligned.substr(int(queryUnAligned.length() * 0.66)); //last 1/3 of the sequence + + Sequence* queryLeft = new Sequence(q.getName(), leftQuery); + Sequence* queryRight = new Sequence(q.getName(), rightQuery); + + vector tempIndexesLeft = databaseLeft->findClosestSequences(queryLeft, num); + vector tempIndexesRight = databaseRight->findClosestSequences(queryRight, num); + + //merge results + map seen; + map::iterator it; + vector mergedResults; + + int index = 0; + // for (int i = 0; i < smaller.size(); i++) { + while(index < tempIndexesLeft.size() && index < tempIndexesRight.size()){ + + if (m->control_pressed) { delete queryRight; delete queryLeft; return refResults; } + + //add left if you havent already + it = seen.find(tempIndexesLeft[index]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesLeft[index]); + seen[tempIndexesLeft[index]] = tempIndexesLeft[index]; + } + + //add right if you havent already + it = seen.find(tempIndexesRight[index]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesRight[index]); + seen[tempIndexesRight[index]] = tempIndexesRight[index]; + } + index++; + } + + + for (int i = index; i < tempIndexesLeft.size(); i++) { + if (m->control_pressed) { delete queryRight; delete queryLeft; return refResults; } + + //add right if you havent already + it = seen.find(tempIndexesLeft[i]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesLeft[i]); + seen[tempIndexesLeft[i]] = tempIndexesLeft[i]; + } + } + + for (int i = index; i < tempIndexesRight.size(); i++) { + if (m->control_pressed) { delete queryRight; delete queryLeft; return refResults; } + + //add right if you havent already + it = seen.find(tempIndexesRight[i]); + if (it == seen.end()) { + mergedResults.push_back(tempIndexesRight[i]); + seen[tempIndexesRight[i]] = tempIndexesRight[i]; + } + } + + for (int i = 0; i < mergedResults.size(); i++) { + //cout << mergedResults[i] << '\t' << db[mergedResults[i]]->getName() << endl; + if (db[mergedResults[i]]->getName() != q.getName()) { + Sequence temp(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); + refResults.push_back(temp); + + } + } + + //cout << endl; + delete queryRight; + delete queryLeft; + + return refResults; + } + catch(exception& e) { + m->errorOut(e, "ChimeraSlayer", "getKmerSeqs"); + exit(1); + } +} +//*************************************************************************************************************** +