From 2e28b53eb15d5dc11653247ee12aed2c7e5aac43 Mon Sep 17 00:00:00 2001 From: westcott Date: Tue, 31 May 2011 15:36:25 +0000 Subject: [PATCH] working on removing pointers from chimera.slayer to eliminate pesky memory leaks --- ccode.cpp | 8 +- ccode.h | 4 +- chimera.h | 20 +++-- chimeracheckrdp.cpp | 8 +- chimeracheckrdp.h | 4 +- chimeraslayer.cpp | 174 +++++++++++++++++++-------------------- chimeraslayer.h | 22 ++--- chimeraslayercommand.cpp | 18 ++-- decalc.cpp | 40 ++++----- decalc.h | 4 +- maligner.cpp | 106 ++++++++++++------------ maligner.h | 24 +++--- pintail.cpp | 4 +- pintail.h | 4 +- slayer.cpp | 152 +++++++++++++++++++++------------- slayer.h | 6 +- 16 files changed, 321 insertions(+), 277 deletions(-) diff --git a/ccode.cpp b/ccode.cpp index e820184..ee88ba1 100644 --- a/ccode.cpp +++ b/ccode.cpp @@ -74,7 +74,7 @@ Ccode::~Ccode() { #endif } //*************************************************************************************************************** -Sequence* Ccode::print(ostream& out, ostream& outAcc) { +Sequence Ccode::print(ostream& out, ostream& outAcc) { try { ofstream out2; @@ -155,7 +155,7 @@ Sequence* Ccode::print(ostream& out, ostream& outAcc) { //free memory for (int i = 0; i < closest.size(); i++) { delete closest[i].seq; } - return NULL; + return *querySeq; } catch(exception& e) { m->errorOut(e, "Ccode", "print"); @@ -164,7 +164,7 @@ Sequence* Ccode::print(ostream& out, ostream& outAcc) { } #ifdef USE_MPI //*************************************************************************************************************** -Sequence* Ccode::print(MPI_File& out, MPI_File& outAcc) { +Sequence Ccode::print(MPI_File& out, MPI_File& outAcc) { try { string outMapString = ""; @@ -263,7 +263,7 @@ Sequence* Ccode::print(MPI_File& out, MPI_File& outAcc) { //free memory for (int i = 0; i < closest.size(); i++) { delete closest[i].seq; } - return NULL; + return *querySeq; } catch(exception& e) { m->errorOut(e, "Ccode", "print"); diff --git a/ccode.h b/ccode.h index 60feaf4..456b735 100644 --- a/ccode.h +++ b/ccode.h @@ -28,10 +28,10 @@ class Ccode : public Chimera { ~Ccode(); int getChimeras(Sequence* query); - Sequence* print(ostream&, ostream&); + Sequence print(ostream&, ostream&); #ifdef USE_MPI - Sequence* print(MPI_File&, MPI_File&); + Sequence print(MPI_File&, MPI_File&); #endif private: diff --git a/chimera.h b/chimera.h index 3f4eb4e..9c5a219 100644 --- a/chimera.h +++ b/chimera.h @@ -103,6 +103,12 @@ struct SeqDist { float dist; int index; }; +/***********************************************************************/ +struct SeqCompare { + Sequence seq; + float dist; + int index; +}; //******************************************************************************************************************** //sorts lowest to highest inline bool compareRegionStart(results left, results right){ @@ -114,7 +120,11 @@ inline bool compareSeqDist(SeqDist left, SeqDist right){ return (left.dist < right.dist); } //******************************************************************************************************************** - +//sorts lowest to highest +inline bool compareSeqCompare(SeqCompare left, SeqCompare right){ + return (left.dist < right.dist); +} +//******************************************************************************************************************** struct sim { string leftParent; string rightParent; @@ -147,14 +157,14 @@ class Chimera { virtual void printHeader(ostream&){}; virtual int getChimeras(Sequence*){ return 0; } virtual int getChimeras(){ return 0; } - virtual Sequence* print(ostream&, ostream&){ return NULL; } - virtual Sequence* print(ostream&, ostream&, data_results, data_results) { return NULL; } + virtual Sequence print(ostream&, ostream&){ Sequence temp; return temp; } + virtual Sequence print(ostream&, ostream&, data_results, data_results) { Sequence temp; return temp; } virtual int print(ostream&, ostream&, string){ return 0; } virtual data_results getResults() { data_results results; return results; } #ifdef USE_MPI - virtual Sequence* print(MPI_File&, MPI_File&){ return 0; } - virtual Sequence* print(MPI_File&, MPI_File&, data_results, data_results){ return NULL; } + virtual Sequence print(MPI_File&, MPI_File&){ Sequence temp; return temp; } + virtual Sequence print(MPI_File&, MPI_File&, data_results, data_results){ Sequence temp; return temp; } virtual int print(MPI_File&, MPI_File&, string){ return 0; } #endif diff --git a/chimeracheckrdp.cpp b/chimeracheckrdp.cpp index 16d5f92..be59315 100644 --- a/chimeracheckrdp.cpp +++ b/chimeracheckrdp.cpp @@ -47,7 +47,7 @@ ChimeraCheckRDP::~ChimeraCheckRDP() { } } //*************************************************************************************************************** -Sequence* ChimeraCheckRDP::print(ostream& out, ostream& outAcc) { +Sequence ChimeraCheckRDP::print(ostream& out, ostream& outAcc) { try { m->mothurOut("Processing: " + querySeq->getName()); m->mothurOutEndLine(); @@ -72,7 +72,7 @@ Sequence* ChimeraCheckRDP::print(ostream& out, ostream& outAcc) { } } - return NULL; + return *querySeq; } catch(exception& e) { m->errorOut(e, "ChimeraCheckRDP", "print"); @@ -81,7 +81,7 @@ Sequence* ChimeraCheckRDP::print(ostream& out, ostream& outAcc) { } #ifdef USE_MPI //*************************************************************************************************************** -Sequence* ChimeraCheckRDP::print(MPI_File& out, MPI_File& outAcc) { +Sequence ChimeraCheckRDP::print(MPI_File& out, MPI_File& outAcc) { try { cout << "Processing: " << querySeq->getName() << endl; @@ -115,7 +115,7 @@ Sequence* ChimeraCheckRDP::print(MPI_File& out, MPI_File& outAcc) { } } - return NULL; + return *querySeq; } catch(exception& e) { m->errorOut(e, "ChimeraCheckRDP", "print"); diff --git a/chimeracheckrdp.h b/chimeracheckrdp.h index 4805ba1..cc23c2e 100644 --- a/chimeracheckrdp.h +++ b/chimeracheckrdp.h @@ -29,10 +29,10 @@ class ChimeraCheckRDP : public Chimera { ~ChimeraCheckRDP(); int getChimeras(Sequence*); - Sequence* print(ostream&, ostream&); + Sequence print(ostream&, ostream&); #ifdef USE_MPI - Sequence* print(MPI_File&, MPI_File&); + Sequence print(MPI_File&, MPI_File&); #endif private: diff --git a/chimeraslayer.cpp b/chimeraslayer.cpp index 3f45e17..b00cedf 100644 --- a/chimeraslayer.cpp +++ b/chimeraslayer.cpp @@ -35,8 +35,6 @@ int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int num realign = r; trimChimera = trim; - decalc = new DeCalculator(); - doPrep(); } catch(exception& e) { @@ -68,8 +66,6 @@ ChimeraSlayer::ChimeraSlayer(string file, string temp, bool trim, map ChimeraSlayer::getTemplate(Sequence* q, vector& userTemplateFiltered) { +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()]; + int myAbund = priority[q.getName()]; for (int i = 0; i < templateSeqs.size(); i++) { @@ -346,7 +342,6 @@ vector ChimeraSlayer::getTemplate(Sequence* q, vector& use //*************************************************************************************************************** ChimeraSlayer::~ChimeraSlayer() { - delete decalc; if (templateFileName != "self") { if (searchMethod == "kmer") { delete databaseRight; delete databaseLeft; } else if (searchMethod == "blast") { delete databaseLeft; } @@ -361,10 +356,10 @@ void ChimeraSlayer::printHeader(ostream& out) { out << "Name\tLeftParent\tRightParent\tDivQLAQRB\tPerIDQLAQRB\tBootStrapA\tDivQLBQRA\tPerIDQLBQRA\tBootStrapB\tFlag\tLeftWindow\tRightWindow\n"; } //*************************************************************************************************************** -Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc) { +Sequence ChimeraSlayer::print(ostream& out, ostream& outAcc) { try { - Sequence* trim = NULL; - if (trimChimera) { trim = new Sequence(trimQuery.getName(), trimQuery.getAligned()); } + Sequence trim; + if (trimChimera) { trim.setName(trimQuery.getName()); trim.setAligned(trimQuery.getAligned()); } if (chimeraFlags == "yes") { string chimeraFlag = "no"; @@ -375,23 +370,23 @@ Sequence* 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 (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(); + 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); + trim.setAligned(newAligned); } } } @@ -399,7 +394,7 @@ Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc) { printBlock(chimeraResults[0], chimeraFlag, out); out << endl; }else { - out << querySeq->getName() << "\tno" << endl; + out << querySeq.getName() << "\tno" << endl; } return trim; @@ -411,13 +406,13 @@ Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc) { } } //*************************************************************************************************************** -Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftPiece, data_results rightPiece) { +Sequence ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftPiece, data_results rightPiece) { try { - Sequence* trim = NULL; + Sequence trim; if (trimChimera) { string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned(); - trim = new Sequence(leftPiece.trimQuery.getName(), aligned); + trim.setName(leftPiece.trimQuery.getName()); trim.setAligned(aligned); } if ((leftPiece.flag == "yes") || (rightPiece.flag == "yes")) { @@ -445,13 +440,13 @@ Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftP 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; + m->mothurOut(querySeq.getName() + "\tyes"); m->mothurOutEndLine(); + outAcc << querySeq.getName() << endl; - if (templateFileName == "self") { chimericSeqs.insert(querySeq->getName()); } + if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); } if (trimChimera) { - string newAligned = trim->getAligned(); + string newAligned = trim.getAligned(); //right side is fine so keep that if ((leftChimeric) && (!rightChimeric)) { @@ -493,7 +488,7 @@ Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftP } } - trim->setAligned(newAligned); + trim.setAligned(newAligned); } } @@ -502,7 +497,7 @@ Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftP printBlock(leftPiece, rightPiece, leftChimeric, rightChimeric, chimeraFlag, out); out << endl; }else { - out << querySeq->getName() << "\tno" << endl; + out << querySeq.getName() << "\tno" << endl; } return trim; @@ -516,18 +511,18 @@ Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc, data_results leftP #ifdef USE_MPI //*************************************************************************************************************** -Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results leftPiece, data_results rightPiece) { +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 = NULL; + Sequence trim; if (trimChimera) { string aligned = leftPiece.trimQuery.getAligned() + rightPiece.trimQuery.getAligned(); - trim = new Sequence(leftPiece.trimQuery.getName(), aligned); + trim.setName(leftPiece.trimQuery.getName()); trim.setAligned(aligned); } @@ -558,11 +553,11 @@ Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results lef 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"; + cout << querySeq.getName() << "\tyes" << endl; + outAccString += querySeq.getName() + "\n"; results = true; - if (templateFileName == "self") { chimericSeqs.insert(querySeq->getName()); } + if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); } //write to accnos file int length = outAccString.length(); @@ -573,7 +568,7 @@ Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results lef delete buf2; if (trimChimera) { - string newAligned = trim->getAligned(); + string newAligned = trim.getAligned(); //right side is fine so keep that if ((leftChimeric) && (!rightChimeric)) { @@ -615,7 +610,7 @@ Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results lef } } - trim->setAligned(newAligned); + trim.setAligned(newAligned); } } @@ -633,7 +628,7 @@ Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results lef delete buf; }else { - outputString += querySeq->getName() + "\tno\n"; + outputString += querySeq.getName() + "\tno\n"; //write to output file int length = outputString.length(); @@ -653,15 +648,15 @@ Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc, data_results lef } } //*************************************************************************************************************** -Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { +Sequence ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { try { MPI_Status status; bool results = false; string outAccString = ""; string outputString = ""; - Sequence* trim = NULL; - if (trimChimera) { trim = new Sequence(trimQuery.getName(), trimQuery.getAligned()); } + Sequence trim; + if (trimChimera) { trim.setName(trimQuery.getName()); trim.setAligned(trimQuery.getAligned()); } if (chimeraFlags == "yes") { string chimeraFlag = "no"; @@ -672,11 +667,11 @@ Sequence* 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()); } + if (templateFileName == "self") { chimericSeqs.insert(querySeq.getName()); } //write to accnos file int length = outAccString.length(); @@ -690,13 +685,13 @@ Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { int lengthLeft = chimeraResults[0].winLEnd - chimeraResults[0].winLStart; int lengthRight = chimeraResults[0].winREnd - chimeraResults[0].winRStart; - string newAligned = trim->getAligned(); + 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); + trim.setAligned(newAligned); } } } @@ -713,7 +708,7 @@ Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) { delete buf; }else { - outputString += querySeq->getName() + "\tno\n"; + outputString += querySeq.getName() + "\tno\n"; //write to output file int length = outputString.length(); @@ -743,20 +738,20 @@ int ChimeraSlayer::getChimeras(Sequence* query) { chimeraFlags = "no"; 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 + else { thisTemplate = getTemplate(*query, thisFilteredTemplate); } //fills this template and creates the databases if (m->control_pressed) { return 0; } if (thisTemplate.size() == 0) { return 0; } //not chimeric //moved this out of maligner - 4/29/11 - vector refSeqs = getRefSeqs(query, thisTemplate, thisFilteredTemplate); + vector refSeqs = getRefSeqs(*query, thisTemplate, thisFilteredTemplate); Maligner maligner(refSeqs, match, misMatch, divR, minSim, minCov); Slayer slayer(window, increment, minSim, divR, iters, minSNP, minBS); @@ -768,13 +763,13 @@ int ChimeraSlayer::getChimeras(Sequence* query) { if (m->control_pressed) { return 0; } - string chimeraFlag = maligner.getResults(query, decalc); + 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]; } + //for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } if (chimeraFlag == "yes") { @@ -791,7 +786,7 @@ int ChimeraSlayer::getChimeras(Sequence* query) { // cout << query->getAligned() << endl; //get sequence that were given from maligner results - vector seqs; + vector seqs; map removeDups; map::iterator itDup; map parentNameSeq; @@ -820,9 +815,9 @@ int ChimeraSlayer::getChimeras(Sequence* query) { for (itDup = removeDups.begin(); itDup != removeDups.end(); itDup++) { itSeq = parentNameSeq.find(itDup->first); - 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); @@ -831,30 +826,30 @@ int ChimeraSlayer::getChimeras(Sequence* query) { //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 // cout << query->getAligned() << endl; - vector seqsForSlayer; + 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(); @@ -862,7 +857,7 @@ int ChimeraSlayer::getChimeras(Sequence* query) { 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; } } //cout << endl << endl; return 0; @@ -875,7 +870,7 @@ int ChimeraSlayer::getChimeras(Sequence* query) { //*************************************************************************************************************** void ChimeraSlayer::printBlock(data_struct data, string flag, ostream& out){ try { - out << querySeq->getName() << '\t'; + out << querySeq.getName() << '\t'; out << data.parentA.getName() << "\t" << data.parentB.getName() << '\t'; out << data.divr_qla_qrb << '\t' << data.qla_qrb << '\t' << data.bsa << '\t'; @@ -894,7 +889,7 @@ void ChimeraSlayer::printBlock(data_results leftdata, data_results rightdata, bo try { if ((leftChimeric) && (!rightChimeric)) { //print left - out << querySeq->getName() << '\t'; + 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'; @@ -903,7 +898,7 @@ void ChimeraSlayer::printBlock(data_results leftdata, data_results rightdata, bo out << flag << '\t' << leftdata.results[0].winLStart << "-" << leftdata.results[0].winLEnd << '\t' << leftdata.results[0].winRStart << "-" << leftdata.results[0].winREnd << '\t'; }else if ((!leftChimeric) && (rightChimeric)) { //print right - out << querySeq->getName() << '\t'; + 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'; @@ -913,7 +908,7 @@ void ChimeraSlayer::printBlock(data_results leftdata, data_results rightdata, bo }else { //print both results if (leftdata.flag == "yes") { - out << querySeq->getName() + "_LEFT" << '\t'; + 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'; @@ -925,7 +920,7 @@ void ChimeraSlayer::printBlock(data_results leftdata, data_results rightdata, bo if (rightdata.flag == "yes") { if (leftdata.flag == "yes") { out << endl; } - out << querySeq->getName() + "_RIGHT"<< '\t'; + 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'; @@ -948,7 +943,7 @@ string ChimeraSlayer::getBlock(data_results leftdata, data_results rightdata, bo string out = ""; if ((leftChimeric) && (!rightChimeric)) { //get left - out += querySeq->getName() + "\t"; + 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"; @@ -957,7 +952,7 @@ string ChimeraSlayer::getBlock(data_results leftdata, data_results rightdata, bo 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 += 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"; @@ -968,7 +963,7 @@ string ChimeraSlayer::getBlock(data_results leftdata, data_results rightdata, bo }else { //print both results if (leftdata.flag == "yes") { - out += querySeq->getName() + "_LEFT\t"; + 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"; @@ -979,7 +974,7 @@ string ChimeraSlayer::getBlock(data_results leftdata, data_results rightdata, bo if (rightdata.flag == "yes") { if (leftdata.flag == "yes") { out += "\n"; } - out += querySeq->getName() + "_RIGHT\t"; + 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"; @@ -1003,7 +998,7 @@ string ChimeraSlayer::getBlock(data_struct data, string flag){ 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"; @@ -1019,16 +1014,16 @@ string ChimeraSlayer::getBlock(data_struct data, string flag){ } } //*************************************************************************************************************** -vector ChimeraSlayer::getRefSeqs(Sequence* q, vector& thisTemplate, vector& thisFilteredTemplate){ +vector ChimeraSlayer::getRefSeqs(Sequence q, vector& thisTemplate, vector& thisFilteredTemplate){ try { - vector refSeqs; + 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()); + Sequence* newSeq = new Sequence(q.getName(), q.getAligned()); runFilter(newSeq); - refSeqs = decalc->findClosest(newSeq, thisTemplate, thisFilteredTemplate, numWanted, minSim); + refSeqs = decalc.findClosest(*newSeq, thisTemplate, thisFilteredTemplate, numWanted, minSim); delete newSeq; }else if (searchMethod == "blast") { refSeqs = getBlastSeqs(q, thisTemplate, numWanted); //fills indexes @@ -1044,18 +1039,18 @@ vector ChimeraSlayer::getRefSeqs(Sequence* q, vector& this } } //***************************************************************************************************************/ -vector ChimeraSlayer::getBlastSeqs(Sequence* q, vector& db, int num) { +vector ChimeraSlayer::getBlastSeqs(Sequence q, vector& db, int num) { try { - vector refResults; + vector refResults; //get parts of query - string queryUnAligned = q->getUnaligned(); + 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); + 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); @@ -1122,10 +1117,9 @@ vector ChimeraSlayer::getBlastSeqs(Sequence* q, vector& db 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 = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); + if (db[mergedResults[i]]->getName() != q.getName()) { + Sequence temp(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); refResults.push_back(temp); - } } //cout << endl << endl; @@ -1133,7 +1127,7 @@ vector ChimeraSlayer::getBlastSeqs(Sequence* q, vector& db delete queryRight; delete queryLeft; - if (refResults.size() == 0) { m->mothurOut("[WARNING]: mothur found 0 potential parents, so we are not able to check " + q->getName() + ". This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); } + if (refResults.size() == 0) { m->mothurOut("[WARNING]: megablast found 0 potential parents, so we are not able to check " + q.getName() + ". This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); } return refResults; } @@ -1143,17 +1137,17 @@ vector ChimeraSlayer::getBlastSeqs(Sequence* q, vector& db } } //*************************************************************************************************************** -vector ChimeraSlayer::getKmerSeqs(Sequence* q, vector& db, int num) { +vector ChimeraSlayer::getKmerSeqs(Sequence q, vector& db, int num) { try { - vector refResults; + vector refResults; //get parts of query - string queryUnAligned = q->getUnaligned(); + 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); + 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); @@ -1210,8 +1204,8 @@ vector ChimeraSlayer::getKmerSeqs(Sequence* q, vector& db, 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 = new Sequence(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); + if (db[mergedResults[i]]->getName() != q.getName()) { + Sequence temp(db[mergedResults[i]]->getName(), db[mergedResults[i]]->getAligned()); refResults.push_back(temp); } diff --git a/chimeraslayer.h b/chimeraslayer.h index ded0f89..304ea10 100644 --- a/chimeraslayer.h +++ b/chimeraslayer.h @@ -15,6 +15,8 @@ #include "maligner.h" #include "slayer.h" + + //***********************************************************************/ //This class was modeled after the chimeraSlayer written by the Broad Institute /***********************************************************************/ @@ -28,21 +30,21 @@ class ChimeraSlayer : public Chimera { ~ChimeraSlayer(); int getChimeras(Sequence*); - Sequence* print(ostream&, ostream&); - Sequence* print(ostream&, ostream&, data_results, data_results); + Sequence print(ostream&, ostream&); + Sequence print(ostream&, ostream&, data_results, data_results); void printHeader(ostream&); int doPrep(); data_results getResults() { return printResults; } #ifdef USE_MPI - Sequence* print(MPI_File&, MPI_File&); - Sequence* print(MPI_File&, MPI_File&, data_results, data_results); + Sequence print(MPI_File&, MPI_File&); + Sequence print(MPI_File&, MPI_File&, data_results, data_results); #endif private: - Sequence* querySeq; + Sequence querySeq; Sequence trimQuery; - DeCalculator* decalc; + DeCalculator decalc; Database* databaseRight; Database* databaseLeft; map priority; //for template=self, seqname, seqAligned, abundance @@ -60,10 +62,10 @@ class ChimeraSlayer : public Chimera { string getBlock(data_struct, string); string getBlock(data_results, data_results, bool, bool, string); //int readNameFile(string); - vector getTemplate(Sequence*, vector&); - vector getRefSeqs(Sequence*, vector&, vector&); - vector getBlastSeqs(Sequence*, vector&, int); - vector getKmerSeqs(Sequence*, vector&, int); + vector getTemplate(Sequence, vector&); + vector getRefSeqs(Sequence, vector&, vector&); + vector getBlastSeqs(Sequence, vector&, int); + vector getKmerSeqs(Sequence, vector&, int); }; diff --git a/chimeraslayercommand.cpp b/chimeraslayercommand.cpp index a422c88..0e7e19b 100644 --- a/chimeraslayercommand.cpp +++ b/chimeraslayercommand.cpp @@ -716,15 +716,15 @@ int ChimeraSlayerCommand::driver(linePair* filePos, string outputFName, string f data_results rightResults = chimera->getResults(); //if either piece is chimeric then report - Sequence* trimmed = chimera->print(out, out2, leftResults, rightResults); - if (trim) { trimmed->printSequence(out3); delete trimmed; } + Sequence trimmed = chimera->print(out, out2, leftResults, rightResults); + if (trim) { trimmed.printSequence(out3); } delete left; delete right; }else { //already chimeric //print results - Sequence* trimmed = chimera->print(out, out2); - if (trim) { trimmed->printSequence(out3); delete trimmed; } + Sequence trimmed = chimera->print(out, out2); + if (trim) { trimmed.printSequence(out3); } } @@ -832,10 +832,9 @@ int ChimeraSlayerCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_Fil data_results rightResults = chimera->getResults(); //if either piece is chimeric then report - Sequence* trimmed = chimera->print(outMPI, outAccMPI, leftResults, rightResults); + Sequence trimmed = chimera->print(outMPI, outAccMPI, leftResults, rightResults); if (trim) { - string outputString = ">" + trimmed->getName() + "\n" + trimmed->getAligned() + "\n"; - delete trimmed; + string outputString = ">" + trimmed.getName() + "\n" + trimmed.getAligned() + "\n"; //write to accnos file int length = outputString.length(); @@ -850,11 +849,10 @@ int ChimeraSlayerCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_Fil }else { //print results - Sequence* trimmed = chimera->print(outMPI, outAccMPI); + Sequence trimmed = chimera->print(outMPI, outAccMPI); if (trim) { - string outputString = ">" + trimmed->getName() + "\n" + trimmed->getAligned() + "\n"; - delete trimmed; + string outputString = ">" + trimmed.getName() + "\n" + trimmed.getAligned() + "\n"; //write to accnos file int length = outputString.length(); diff --git a/decalc.cpp b/decalc.cpp index 3a8c9d5..973340c 100644 --- a/decalc.cpp +++ b/decalc.cpp @@ -683,23 +683,23 @@ float DeCalculator::getCoef(vector obs, vector qav) { } //*************************************************************************************************************** //gets closest matches to each end, since chimeras will most likely have different parents on each end -vector DeCalculator::findClosest(Sequence* querySeq, vector& thisTemplate, vector& thisFilteredTemplate, int numWanted, int minSim) { +vector DeCalculator::findClosest(Sequence querySeq, vector& thisTemplate, vector& thisFilteredTemplate, int numWanted, int minSim) { try { //indexes.clear(); - vector seqsMatches; + vector seqsMatches; vector distsLeft; vector distsRight; Dist* distcalculator = new eachGapDist(); - string queryUnAligned = querySeq->getUnaligned(); + string queryUnAligned = querySeq.getUnaligned(); int numBases = int(queryUnAligned.length() * 0.33); string leftQuery = ""; //first 1/3 of the sequence string rightQuery = ""; //last 1/3 of the sequence - string queryAligned = querySeq->getAligned(); + string queryAligned = querySeq.getAligned(); //left side bool foundFirstBase = false; @@ -740,8 +740,8 @@ vector DeCalculator::findClosest(Sequence* querySeq, vectorgetName(), leftQuery); - Sequence queryRight(querySeq->getName(), rightQuery); + Sequence queryLeft(querySeq.getName(), leftQuery); + Sequence queryRight(querySeq.getName(), rightQuery); //cout << querySeq->getName() << '\t' << leftSpot << '\t' << rightSpot << '\t' << firstBaseSpot << '\t' << lastBaseSpot << endl; //cout << queryUnAligned.length() << '\t' << queryLeft.getUnaligned().length() << '\t' << queryRight.getUnaligned().length() << endl; @@ -841,8 +841,8 @@ vector DeCalculator::findClosest(Sequence* querySeq, vectorgetName() << '\t' << dists[i].dist << endl; - if ((thisTemplate[dists[i].index]->getName() != querySeq->getName()) && (((1.0-dists[i].dist)*100) >= minSim)) { - Sequence* temp = new Sequence(thisTemplate[dists[i].index]->getName(), thisTemplate[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. + if ((thisTemplate[dists[i].index]->getName() != querySeq.getName()) && (((1.0-dists[i].dist)*100) >= minSim)) { + Sequence temp(thisTemplate[dists[i].index]->getName(), thisTemplate[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. //cout << querySeq->getName() << '\t' << thisTemplate[dists[i].index]->getName() << '\t' << dists[i].dist << endl; seqsMatches.push_back(temp); } @@ -889,17 +889,17 @@ Sequence* DeCalculator::findClosest(Sequence* querySeq, vector db) { } } /***************************************************************************************************************/ -map DeCalculator::trimSeqs(Sequence* query, vector topMatches) { +map DeCalculator::trimSeqs(Sequence& query, vector& topMatches) { try { int frontPos = 0; //should contain first position in all seqs that is not a gap character - int rearPos = query->getAligned().length(); + int rearPos = query.getAligned().length(); //********find first position in topMatches that is a non gap character***********// //find first position all query seqs that is a non gap character for (int i = 0; i < topMatches.size(); i++) { - string aligned = topMatches[i]->getAligned(); + string aligned = topMatches[i].getAligned(); int pos = 0; //find first spot in this seq @@ -915,7 +915,7 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch } - string aligned = query->getAligned(); + string aligned = query.getAligned(); int pos = 0; //find first position in query that is a non gap character @@ -933,7 +933,7 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch //********find last position in topMatches that is a non gap character***********// for (int i = 0; i < topMatches.size(); i++) { - string aligned = topMatches[i]->getAligned(); + string aligned = topMatches[i].getAligned(); int pos = aligned.length(); //find first spot in this seq @@ -949,7 +949,7 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch } - aligned = query->getAligned(); + aligned = query.getAligned(); pos = aligned.length(); //find last position in query that is a non gap character @@ -966,24 +966,24 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch map trimmedPos; //check to make sure that is not whole seq if ((rearPos - frontPos - 1) <= 0) { - query->setAligned(""); + query.setAligned(""); //trim topMatches for (int i = 0; i < topMatches.size(); i++) { - topMatches[i]->setAligned(""); + topMatches[i].setAligned(""); } }else { //trim query - string newAligned = query->getAligned(); + string newAligned = query.getAligned(); newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1)); - query->setAligned(newAligned); + query.setAligned(newAligned); //trim topMatches for (int i = 0; i < topMatches.size(); i++) { - newAligned = topMatches[i]->getAligned(); + newAligned = topMatches[i].getAligned(); newAligned = newAligned.substr(frontPos, (rearPos-frontPos+1)); - topMatches[i]->setAligned(newAligned); + topMatches[i].setAligned(newAligned); } for (int i = 0; i < newAligned.length(); i++) { diff --git a/decalc.h b/decalc.h index 40b3a10..d6cca18 100644 --- a/decalc.h +++ b/decalc.h @@ -39,14 +39,14 @@ class DeCalculator { DeCalculator() { m = MothurOut::getInstance(); } ~DeCalculator() {}; - vector findClosest(Sequence*, vector&, vector&, int, int); //takes querySeq, a reference db, filteredRefDB, numWanted, minSim + vector findClosest(Sequence, vector&, vector&, int, int); //takes querySeq, a reference db, filteredRefDB, numWanted, minSim Sequence* findClosest(Sequence*, vector); set getPos() { return h; } void setMask(string); void setAlignmentLength(int l) { alignLength = l; } void runMask(Sequence*); void trimSeqs(Sequence*, Sequence*, map&); - map trimSeqs(Sequence*, vector); + map trimSeqs(Sequence&, vector&); void removeObviousOutliers(vector< vector >&, int); vector calcFreq(vector, string); vector findWindows(Sequence*, int, int, int&, int); diff --git a/maligner.cpp b/maligner.cpp index f4c4a98..f31e482 100644 --- a/maligner.cpp +++ b/maligner.cpp @@ -10,33 +10,33 @@ #include "maligner.h" /***********************************************************************/ //int num, int match, int misMatch, , string mode, Database* dataLeft, Database* dataRight -Maligner::Maligner(vector temp, int match, int misMatch, float div, int ms, int minCov) : db(temp), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov) { +Maligner::Maligner(vector temp, int match, int misMatch, float div, int ms, int minCov) : db(temp), matchScore(match), misMatchPenalty(misMatch), minDivR(div), minSimilarity(ms), minCoverage(minCov) { //numWanted(num), , searchMethod(mode), databaseLeft(dataLeft), databaseRight(dataRight) m = MothurOut::getInstance(); } /***********************************************************************/ -string Maligner::getResults(Sequence* q, DeCalculator* decalc) { +string Maligner::getResults(Sequence q, DeCalculator decalc) { try { outputResults.clear(); //make copy so trimming doesn't destroy query from calling class - remember to deallocate - query = new Sequence(q->getName(), q->getAligned()); + query.setName(q.getName()); query.setAligned(q.getAligned()); string chimera; //copy refSeqs so that filter does not effect original for(int i = 0; i < db.size(); i++) { - Sequence* newSeq = new Sequence(db[i]->getName(), db[i]->getAligned()); + Sequence newSeq(db[i].getName(), db[i].getAligned()); refSeqs.push_back(newSeq); } refSeqs = minCoverageFilter(refSeqs); if (refSeqs.size() < 2) { - for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } + //for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } percentIdenticalQueryChimera = 0.0; return "unknown"; } @@ -49,9 +49,9 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) { if (m->control_pressed) { return chimera; } //free memory - delete query; + //delete query; - for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } + //for (int i = 0; i < refSeqs.size(); i++) { delete refSeqs[i]; } return chimera; } @@ -61,30 +61,26 @@ string Maligner::getResults(Sequence* q, DeCalculator* decalc) { } } /***********************************************************************/ -string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { +string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator decalc) { try { string chimera; - //trims seqs to first non gap char in all seqs and last non gap char in all seqs - spotMap = decalc->trimSeqs(query, refSeqs); - + spotMap = decalc.trimSeqs(query, refSeqs); + //you trimmed the whole sequence, skip - if (query->getAligned() == "") { return "no"; } + if (query.getAligned() == "") { return "no"; } - vector temp = refSeqs; - -// for(int i=0;igetName() << endl; -// } - + vector temp = refSeqs; temp.push_back(query); - verticalFilter(temp); + temp = verticalFilter(temp); + query = temp[temp.size()-1]; + for (int i = 0; i < temp.size()-1; i++) { refSeqs[i] = temp[i]; } //for (int i = 0; i < refSeqs.size(); i++) { cout << refSeqs[i]->getName() << endl ; }//<< refSeqs[i]->getAligned() << endl - vector< vector > matrix = buildScoreMatrix(query->getAligned().length(), refSeqs.size()); //builds and initializes + vector< vector > matrix = buildScoreMatrix(query.getAligned().length(), refSeqs.size()); //builds and initializes if (m->control_pressed) { return chimera; } @@ -94,14 +90,14 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { if (m->control_pressed) { return chimera; } - vector trace = mapTraceRegionsToAlignment(path, refSeqs); + vector trace = mapTraceRegionsToAlignment(path); if (trace.size() > 1) { chimera = "yes"; } else { chimera = "no"; return chimera; } int traceStart = path[0].col; int traceEnd = path[path.size()-1].col; - string queryInRange = query->getAligned(); + string queryInRange = query.getAligned(); queryInRange = queryInRange.substr(traceStart, (traceEnd-traceStart+1)); // cout << queryInRange << endl; string chimeraSeq = constructChimericSeq(trace, refSeqs); @@ -140,23 +136,23 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { // cout << regionStart << '\t' << regionEnd << '\t' << seqIndex << endl; results temp; - temp.parent = refSeqs[seqIndex]->getName(); - temp.parentAligned = db[seqIndex]->getAligned(); + temp.parent = refSeqs[seqIndex].getName(); + temp.parentAligned = db[seqIndex].getAligned(); temp.nastRegionStart = spotMap[regionStart]; temp.nastRegionEnd = spotMap[regionEnd]; temp.regionStart = unalignedMap[regionStart]; temp.regionEnd = unalignedMap[regionEnd]; - string parentInRange = refSeqs[seqIndex]->getAligned(); + string parentInRange = refSeqs[seqIndex].getAligned(); parentInRange = parentInRange.substr(traceStart, (traceEnd-traceStart+1)); temp.queryToParent = computePercentID(queryInRange, parentInRange); temp.divR = (percentIdenticalQueryChimera / temp.queryToParent); - string queryInRegion = query->getAligned(); + string queryInRegion = query.getAligned(); queryInRegion = queryInRegion.substr(regionStart, (regionEnd-regionStart+1)); - string parentInRegion = refSeqs[seqIndex]->getAligned(); + string parentInRegion = refSeqs[seqIndex].getAligned(); parentInRegion = parentInRegion.substr(regionStart, (regionEnd-regionStart+1)); temp.queryToParentLocal = computePercentID(queryInRegion, parentInRegion); @@ -175,15 +171,15 @@ string Maligner::chimeraMaligner(int chimeraPenalty, DeCalculator* decalc) { } /***********************************************************************/ //removes top matches that do not have minimum coverage with query. -vector Maligner::minCoverageFilter(vector ref){ +vector Maligner::minCoverageFilter(vector ref){ try { - vector newRefs; + vector newRefs; - string queryAligned = query->getAligned(); + string queryAligned = query.getAligned(); for (int i = 0; i < ref.size(); i++) { - string refAligned = ref[i]->getAligned(); + string refAligned = ref[i].getAligned(); int numBases = 0; int numCovered = 0; @@ -205,9 +201,9 @@ vector Maligner::minCoverageFilter(vector ref){ //if coverage above minimum if (coverage > minCoverage) { newRefs.push_back(ref[i]); - }else { - delete ref[i]; - } + }//else { + //delete ref[i]; + //} } return newRefs; @@ -222,7 +218,7 @@ vector Maligner::minCoverageFilter(vector ref){ int Maligner::computeChimeraPenalty() { try { - int numAllowable = ((1.0 - (1.0/minDivR)) * query->getNumBases()); + int numAllowable = ((1.0 - (1.0/minDivR)) * query.getNumBases()); // if(numAllowable < 1){ numAllowable = 1; } @@ -238,16 +234,16 @@ int Maligner::computeChimeraPenalty() { } /***********************************************************************/ //this is a vertical filter -void Maligner::verticalFilter(vector seqs) { +vector Maligner::verticalFilter(vector seqs) { try { - vector gaps; gaps.resize(query->getAligned().length(), 0); + vector gaps; gaps.resize(query.getAligned().length(), 0); - string filterString = (string(query->getAligned().length(), '1')); + string filterString = (string(query.getAligned().length(), '1')); //for each sequence for (int i = 0; i < seqs.size(); i++) { - string seqAligned = seqs[i]->getAligned(); + string seqAligned = seqs[i].getAligned(); for (int j = 0; j < seqAligned.length(); j++) { //if this spot is a gap @@ -257,7 +253,7 @@ void Maligner::verticalFilter(vector seqs) { //zero out spot where all sequences have blanks int numColRemoved = 0; - for(int i = 0; i < seqs[0]->getAligned().length(); i++){ + for(int i = 0; i < seqs[0].getAligned().length(); i++){ if(gaps[i] == seqs.size()) { filterString[i] = '0'; numColRemoved++; } } @@ -265,7 +261,7 @@ void Maligner::verticalFilter(vector seqs) { //for each sequence for (int i = 0; i < seqs.size(); i++) { - string seqAligned = seqs[i]->getAligned(); + string seqAligned = seqs[i].getAligned(); string newAligned = ""; int count = 0; @@ -278,10 +274,10 @@ void Maligner::verticalFilter(vector seqs) { } } - seqs[i]->setAligned(newAligned); + seqs[i].setAligned(newAligned); } - string query = seqs[seqs.size()-1]->getAligned(); + string query = seqs[seqs.size()-1].getAligned(); int queryLength = query.length(); unalignedMap.resize(queryLength, 0); @@ -297,6 +293,8 @@ void Maligner::verticalFilter(vector seqs) { } spotMap = newMap; + + return seqs; } catch(exception& e) { m->errorOut(e, "Maligner", "verticalFilter"); @@ -333,19 +331,19 @@ vector< vector > Maligner::buildScoreMatrix(int cols, int rows) { //*************************************************************************************************************** -void Maligner::fillScoreMatrix(vector >& ms, vector seqs, int penalty) { +void Maligner::fillScoreMatrix(vector >& ms, vector seqs, int penalty) { try{ //get matrix dimensions - int numCols = query->getAligned().length(); + int numCols = query.getAligned().length(); int numRows = seqs.size(); // cout << numRows << endl; //initialize first col - string queryAligned = query->getAligned(); + string queryAligned = query.getAligned(); for (int i = 0; i < numRows; i++) { - string subjectAligned = seqs[i]->getAligned(); + string subjectAligned = seqs[i].getAligned(); //are you both gaps? if ((!isalpha(queryAligned[0])) && (!isalpha(subjectAligned[0]))) { @@ -366,7 +364,7 @@ void Maligner::fillScoreMatrix(vector >& ms, vectorgetAligned(); + string subjectAligned = seqs[i].getAligned(); int matchMisMatchScore = 0; //are you both gaps? @@ -434,7 +432,7 @@ vector Maligner::extractHighestPath(vector > try { //get matrix dimensions - int numCols = query->getAligned().length(); + int numCols = query.getAligned().length(); int numRows = ms.size(); @@ -480,7 +478,7 @@ vector Maligner::extractHighestPath(vector > } } //*************************************************************************************************************** -vector Maligner::mapTraceRegionsToAlignment(vector path, vector seqs) { +vector Maligner::mapTraceRegionsToAlignment(vector path) { try { vector trace; @@ -668,14 +666,14 @@ vector Maligner::mapTraceRegionsToAlignment(vector p */ //*************************************************************************************************************** -string Maligner::constructChimericSeq(vector trace, vector seqs) { +string Maligner::constructChimericSeq(vector trace, vector seqs) { try { string chimera = ""; for (int i = 0; i < trace.size(); i++) { // cout << i << '\t' << trace[i].row << '\t' << trace[i].col << '\t' << trace[i].oldCol << endl; - string seqAlign = seqs[trace[i].row]->getAligned(); + string seqAlign = seqs[trace[i].row].getAligned(); seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1)); chimera += seqAlign; } @@ -692,7 +690,7 @@ string Maligner::constructChimericSeq(vector trace, vector trace, vector seqs) { +string Maligner::constructAntiChimericSeq(vector trace, vector seqs) { try { string antiChimera = ""; @@ -701,7 +699,7 @@ string Maligner::constructAntiChimericSeq(vector trace, vectorgetAligned(); + string seqAlign = seqs[trace[oppositeIndex].row].getAligned(); seqAlign = seqAlign.substr(trace[i].col, (trace[i].oldCol-trace[i].col+1)); antiChimera += seqAlign; } diff --git a/maligner.h b/maligner.h index 4f5c949..fc00e99 100644 --- a/maligner.h +++ b/maligner.h @@ -20,35 +20,35 @@ class Maligner { public: - Maligner(vector, int, int, float, int, int); //int, int, int, , string, Database*, Database* + Maligner(vector, int, int, float, int, int); //int, int, int, , string, Database*, Database* ~Maligner() {}; - string getResults(Sequence*, DeCalculator*); + string getResults(Sequence, DeCalculator); float getPercentID() { return percentIdenticalQueryChimera; } vector getOutput() { return outputResults; } private: - Sequence* query; - vector refSeqs; - vector db; + Sequence query; + vector refSeqs; + vector db; int minCoverage, minSimilarity, matchScore, misMatchPenalty; float minDivR, percentIdenticalQueryChimera; vector outputResults; map spotMap; vector unalignedMap; - vector minCoverageFilter(vector); //removes top matches that do not have minimum coverage with query. + vector minCoverageFilter(vector); //removes top matches that do not have minimum coverage with query. int computeChimeraPenalty(); - void verticalFilter(vector); + vector verticalFilter(vector); vector< vector > buildScoreMatrix(int, int); - void fillScoreMatrix(vector >&, vector, int); + void fillScoreMatrix(vector >&, vector, int); vector extractHighestPath(vector >); - vector mapTraceRegionsToAlignment(vector, vector); - string constructChimericSeq(vector, vector); - string constructAntiChimericSeq(vector, vector); + vector mapTraceRegionsToAlignment(vector); + string constructChimericSeq(vector, vector); + string constructAntiChimericSeq(vector, vector); float computePercentID(string, string); - string chimeraMaligner(int, DeCalculator*); + string chimeraMaligner(int, DeCalculator); MothurOut* m; }; diff --git a/pintail.cpp b/pintail.cpp index b04efea..10f931f 100644 --- a/pintail.cpp +++ b/pintail.cpp @@ -249,7 +249,7 @@ int Pintail::doPrep() { } } //*************************************************************************************************************** -Sequence* Pintail::print(ostream& out, ostream& outAcc) { +Sequence Pintail::print(ostream& out, ostream& outAcc) { try { int index = ceil(deviation); @@ -280,7 +280,7 @@ Sequence* Pintail::print(ostream& out, ostream& outAcc) { for (int m = 0; m < expectedDistance.size(); m++) { out << expectedDistance[m] << '\t'; } out << endl; - return NULL; + return *querySeq; } catch(exception& e) { diff --git a/pintail.h b/pintail.h index 6bbc1a3..92c3998 100644 --- a/pintail.h +++ b/pintail.h @@ -28,13 +28,13 @@ class Pintail : public Chimera { ~Pintail(); int getChimeras(Sequence*); - Sequence* print(ostream&, ostream&); + Sequence print(ostream&, ostream&); void setCons(string c) { consfile = c; } void setQuantiles(string q) { quanfile = q; } #ifdef USE_MPI - Sequence* print(MPI_File&, MPI_File&); + Sequence print(MPI_File&, MPI_File&); #endif private: diff --git a/slayer.cpp b/slayer.cpp index 4987d41..a787f27 100644 --- a/slayer.cpp +++ b/slayer.cpp @@ -13,10 +13,10 @@ Slayer::Slayer(int win, int increment, int parentThreshold, float div, int i, int snp, int mi) : minBS(mi), windowSize(win), windowStep(increment), parentFragmentThreshold(parentThreshold), divRThreshold(div), iters(i), percentSNPSample(snp){ m = MothurOut::getInstance(); } /***********************************************************************/ -string Slayer::getResults(Sequence* query, vector refSeqs) { +string Slayer::getResults(Sequence query, vector refSeqs) { try { vector all; all.clear(); - myQuery = *query; + myQuery = query; for (int i = 0; i < refSeqs.size(); i++) { @@ -25,19 +25,21 @@ string Slayer::getResults(Sequence* query, vector refSeqs) { if (m->control_pressed) { return "no"; } //make copies of query and each parent because runBellerophon removes gaps and messes them up - Sequence* q = new Sequence(query->getName(), query->getAligned()); - Sequence* leftParent = new Sequence(refSeqs[i]->getName(), refSeqs[i]->getAligned()); - Sequence* rightParent = new Sequence(refSeqs[j]->getName(), refSeqs[j]->getAligned()); + Sequence q(query.getName(), query.getAligned()); + Sequence leftParent(refSeqs[i].getName(), refSeqs[i].getAligned()); + Sequence rightParent(refSeqs[j].getName(), refSeqs[j].getAligned()); //cout << q->getName() << endl << q->getAligned() << endl << endl; - //cout << leftParent->getName() << endl << leftParent->getAligned() << endl << endl; + //cout << leftParent.getName() << '\t' << leftParent.getAligned().length() << endl << endl; + //cout << rightParent.getName() << '\t' << rightParent.getAligned().length() << endl << endl; + //cout << q.getName() << '\t' << q.getAligned().length() << endl << endl; //cout << rightParent->getName() << endl << rightParent->getAligned() << endl << endl; //cout << " length = " << rightParent->getAligned().length() << endl; map spots; //map from spot in original sequence to spot in filtered sequence for query and both parents vector divs = runBellerophon(q, leftParent, rightParent, spots); - if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; } + if (m->control_pressed) { return "no"; } // cout << "examining:\t" << refSeqs[i]->getName() << '\t' << refSeqs[j]->getName() << endl; vector selectedDivs; for (int k = 0; k < divs.size(); k++) { @@ -45,7 +47,7 @@ string Slayer::getResults(Sequence* query, vector refSeqs) { vector snpsLeft = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winLStart, divs[k].winLEnd); vector snpsRight = getSNPS(divs[k].parentA.getAligned(), divs[k].querySeq.getAligned(), divs[k].parentB.getAligned(), divs[k].winRStart, divs[k].winREnd); - if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; } + if (m->control_pressed) { return "no"; } int numSNPSLeft = snpsLeft.size(); int numSNPSRight = snpsRight.size(); @@ -57,7 +59,7 @@ string Slayer::getResults(Sequence* query, vector refSeqs) { float BS_A, BS_B; bootstrapSNPS(snpsLeft, snpsRight, BS_A, BS_B, iters); - if (m->control_pressed) { delete q; delete leftParent; delete rightParent; return "no"; } + if (m->control_pressed) { return "no"; } divs[k].bsa = BS_A; divs[k].bsb = BS_B; @@ -89,10 +91,6 @@ string Slayer::getResults(Sequence* query, vector refSeqs) { //save selected for (int mi = 0; mi < selectedDivs.size(); mi++) { all.push_back(selectedDivs[mi]); } - - delete q; - delete leftParent; - delete rightParent; } } @@ -116,26 +114,27 @@ string Slayer::getResults(Sequence* query, vector refSeqs) { } } /***********************************************************************/ -vector Slayer::runBellerophon(Sequence* q, Sequence* pA, Sequence* pB, map& spots) { +vector Slayer::runBellerophon(Sequence q, Sequence pA, Sequence pB, map& spots) { try{ vector data; //vertical filter - vector temp; - temp.push_back(q); temp.push_back(pA); temp.push_back(pB); + //cout << q.getName() << endl << q.getAligned() << endl << endl; + //cout << pA.getName() << endl << pA.getUnaligned() << endl << endl; + //cout << pB.getName() << endl << pB.getUnaligned() << endl << endl; //maps spot in new alignment to spot in alignment before filter - spots = verticalFilter(temp); //fills baseSpots + spots = verticalFilter(q, pA, pB); //fills baseSpots //get these to avoid numerous function calls - string query = q->getAligned(); - string parentA = pA->getAligned(); - string parentB = pB->getAligned(); + string query = q.getAligned(); + string parentA = pA.getAligned(); + string parentB = pB.getAligned(); int length = query.length(); -//cout << q->getName() << endl << q->getAligned() << endl << endl; -//cout << pA->getName() << endl << pA->getUnaligned() << endl << endl; -//cout << pB->getName() << endl << pB->getUnaligned() << endl << endl; +//cout << q.getName() << endl << q.getAligned() << endl << endl; +//cout << pA.getName() << endl << pA.getUnaligned() << endl << endl; +//cout << pB.getName() << endl << pB.getUnaligned() << endl << endl; //cout << " length = " << length << endl; //check window size @@ -209,9 +208,9 @@ vector Slayer::runBellerophon(Sequence* q, Sequence* pA, Sequence* member.winLEnd = breakpoint; member.winRStart = breakpoint+1; member.winREnd = length-1; - member.querySeq = *(q); - member.parentA = *(pA); - member.parentB = *(pB); + member.querySeq = q; + member.parentA = pA; + member.parentB = pB; member.bsa = 0; member.bsb = 0; member.bsMax = 0; @@ -478,33 +477,41 @@ float Slayer::computePercentID(string queryAlign, string chimera, int left, int } /***********************************************************************/ //remove columns that contain any gaps -map Slayer::verticalFilter(vector seqs) { +map Slayer::verticalFilter(Sequence& q, Sequence& pA, Sequence& pB) { try { //find baseSpots baseSpots.clear(); baseSpots.resize(3); //query, parentA, parentB - vector gaps; gaps.resize(seqs[0]->getAligned().length(), 0); + vector gaps; gaps.resize(q.getAligned().length(), 0); - string filterString = (string(seqs[0]->getAligned().length(), '1')); + string filterString = (string(q.getAligned().length(), '1')); - //for each sequence - for (int i = 0; i < seqs.size(); i++) { + string seqAligned = q.getAligned(); + for (int j = 0; j < seqAligned.length(); j++) { + //if this spot is a gap + if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; } + } - string seqAligned = seqs[i]->getAligned(); - - for (int j = 0; j < seqAligned.length(); j++) { - //if this spot is a gap - if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; } - } + seqAligned = pA.getAligned(); + for (int j = 0; j < seqAligned.length(); j++) { + //if this spot is a gap + if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; } } + seqAligned = pB.getAligned(); + for (int j = 0; j < seqAligned.length(); j++) { + //if this spot is a gap + if ((seqAligned[j] == '-') || (seqAligned[j] == '.') || (toupper(seqAligned[j]) == 'N')) { gaps[j]++; } + } + + //zero out spot where any sequences have blanks int numColRemoved = 0; int count = 0; map maskMap; maskMap.clear(); - for(int i = 0; i < seqs[0]->getAligned().length(); i++){ + for(int i = 0; i < q.getAligned().length(); i++){ if(gaps[i] != 0) { filterString[i] = '0'; numColRemoved++; } else { maskMap[count] = i; @@ -512,29 +519,64 @@ map Slayer::verticalFilter(vector seqs) { } } - //for each sequence - for (int i = 0; i < seqs.size(); i++) { - - string seqAligned = seqs[i]->getAligned(); - string newAligned = ""; + seqAligned = q.getAligned(); + string newAligned = ""; - int baseCount = 0; - int count = 0; - for (int j = 0; j < seqAligned.length(); j++) { - //are you a base - if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; } + int baseCount = 0; + count = 0; + for (int j = 0; j < seqAligned.length(); j++) { + //are you a base + if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; } - //if this spot is not a gap - if (filterString[j] == '1') { - newAligned += seqAligned[j]; - baseSpots[i][count] = baseCount; - count++; - } + //if this spot is not a gap + if (filterString[j] == '1') { + newAligned += seqAligned[j]; + baseSpots[0][count] = baseCount; + count++; } + } - seqs[i]->setAligned(newAligned); + q.setAligned(newAligned); + + seqAligned = pA.getAligned(); + newAligned = ""; + + baseCount = 0; + count = 0; + for (int j = 0; j < seqAligned.length(); j++) { + //are you a base + if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; } + + //if this spot is not a gap + if (filterString[j] == '1') { + newAligned += seqAligned[j]; + baseSpots[1][count] = baseCount; + count++; + } + } + + pA.setAligned(newAligned); + + seqAligned = pB.getAligned(); + newAligned = ""; + + baseCount = 0; + count = 0; + for (int j = 0; j < seqAligned.length(); j++) { + //are you a base + if ((seqAligned[j] != '-') && (seqAligned[j] != '.') && (toupper(seqAligned[j]) != 'N')) { baseCount++; } + + //if this spot is not a gap + if (filterString[j] == '1') { + newAligned += seqAligned[j]; + baseSpots[2][count] = baseCount; + count++; + } } + pB.setAligned(newAligned); + + return maskMap; } catch(exception& e) { diff --git a/slayer.h b/slayer.h index 107b027..eeefc62 100644 --- a/slayer.h +++ b/slayer.h @@ -32,7 +32,7 @@ class Slayer { Slayer(int, int, int, float, int, int, int); ~Slayer() {}; - string getResults(Sequence*, vector); + string getResults(Sequence, vector); vector getOutput() { return outputResults; } @@ -44,10 +44,10 @@ class Slayer { vector< map > baseSpots; Sequence myQuery; - map verticalFilter(vector); + map verticalFilter(Sequence&, Sequence&, Sequence&); float computePercentID(string, string, int, int); - vector runBellerophon(Sequence*, Sequence*, Sequence*, map&); + vector runBellerophon(Sequence, Sequence, Sequence, map&); vector getSNPS(string, string, string, int, int); int bootstrapSNPS(vector, vector, float&, float&, int); float snpQA(vector); -- 2.39.2