From: westcott Date: Thu, 23 Dec 2010 19:05:17 +0000 (+0000) Subject: fixed memory leak in decalc findClosest function used by chimera slayer X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=1f7448f6ec056b30ed146a5f779dfd6029788316 fixed memory leak in decalc findClosest function used by chimera slayer --- diff --git a/commandfactory.cpp b/commandfactory.cpp index b7b4992..f1a7105 100644 --- a/commandfactory.cpp +++ b/commandfactory.cpp @@ -103,6 +103,7 @@ #include "removeotuscommand.h" #include "indicatorcommand.h" #include "consensusseqscommand.h" +#include "corraxescommand.h" /*******************************************************/ @@ -210,6 +211,7 @@ CommandFactory::CommandFactory(){ commands["remove.otus"] = "remove.otus"; commands["indicator"] = "indicator"; commands["consensus.seqs"] = "consensus.seqs"; + commands["corr.axes"] = "corr.axes"; commands["pairwise.seqs"] = "MPIEnabled"; commands["pipeline.pds"] = "MPIEnabled"; commands["classify.seqs"] = "MPIEnabled"; @@ -364,6 +366,7 @@ Command* CommandFactory::getCommand(string commandName, string optionString){ else if(commandName == "sub.sample") { command = new SubSampleCommand(optionString); } else if(commandName == "indicator") { command = new IndicatorCommand(optionString); } else if(commandName == "consensus.seqs") { command = new ConsensusSeqsCommand(optionString); } + else if(commandName == "corr.axes") { command = new CorrAxesCommand(optionString); } else { command = new NoCommand(optionString); } return command; @@ -485,6 +488,7 @@ Command* CommandFactory::getCommand(string commandName, string optionString, str else if(commandName == "sub.sample") { pipecommand = new SubSampleCommand(optionString); } else if(commandName == "indicator") { pipecommand = new IndicatorCommand(optionString); } else if(commandName == "consensus.seqs") { pipecommand = new ConsensusSeqsCommand(optionString); } + else if(commandName == "corr.axes") { pipecommand = new CorrAxesCommand(optionString); } else { pipecommand = new NoCommand(optionString); } return pipecommand; @@ -594,6 +598,7 @@ Command* CommandFactory::getCommand(string commandName){ else if(commandName == "sub.sample") { shellcommand = new SubSampleCommand(); } else if(commandName == "indicator") { shellcommand = new IndicatorCommand(); } else if(commandName == "consensus.seqs") { shellcommand = new ConsensusSeqsCommand(); } + else if(commandName == "corr.axes") { shellcommand = new CorrAxesCommand(); } else { shellcommand = new NoCommand(); } return shellcommand; diff --git a/corraxescommand.cpp b/corraxescommand.cpp index 322b7f9..a36f95d 100644 --- a/corraxescommand.cpp +++ b/corraxescommand.cpp @@ -12,7 +12,7 @@ //********************************************************************************************************************** vector CorrAxesCommand::getValidParameters(){ try { - string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metastats","outputdir","inputdir"}; + string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"}; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); return myArray; } @@ -69,7 +69,7 @@ CorrAxesCommand::CorrAxesCommand(string option) { else { //valid paramters for this command - string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metastats","outputdir","inputdir"}; + string Array[] = {"axes","shared","relabund","numaxes","label","groups","method","metadata","outputdir","inputdir"}; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); OptionParser parser(option); @@ -114,6 +114,14 @@ CorrAxesCommand::CorrAxesCommand(string option) { //if the user has not given a path then, add inputdir. else leave path alone. if (path == "") { parameters["relabund"] = inputDir + it->second; } } + + it = parameters.find("metadata"); + //user has given a template file + if(it != parameters.end()){ + path = m->hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["metadata"] = inputDir + it->second; } + } } @@ -132,6 +140,10 @@ CorrAxesCommand::CorrAxesCommand(string option) { else if (relabundfile == "not found") { relabundfile = ""; } else { inputFileName = relabundfile; } + metadatafile = validParameter.validFile(parameters, "metadata", true); + if (metadatafile == "not open") { abort = true; } + else if (metadatafile == "not found") { metadatafile = ""; } + groups = validParameter.validFile(parameters, "groups", false); if (groups == "not found") { groups = ""; pickedGroups = false; } else { @@ -233,13 +245,23 @@ int CorrAxesCommand::execute(){ if (m->control_pressed) { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; } + /*************************************************************************************/ + // calc the r values // + /************************************************************************************/ + string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + "corr.axes"; outputNames.push_back(outputFileName); outputTypes["corr.axes"].push_back(outputFileName); ofstream out; m->openOutputFile(outputFileName, out); + out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); + + //output headings + out << "OTU\t"; + for (int i = 0; i < numaxes; i++) { out << "axis" << (i+1) << '\t'; } + out << endl; if (method == "pearson") { calcPearson(axes, out); } - //else if (method == "spearman") { calcSpearman(axes, out); } + else if (method == "spearman") { calcSpearman(axes, out); } //else if (method == "kendall") { calcKendal(axes, out); } //else { m->mothurOut("[ERROR]: Invalid method."); m->mothurOutEndLine(); } @@ -278,6 +300,8 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou //for each otu for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) { + out << i+1 << '\t'; + //find the averages this otu - Y float sumOtu = 0.0; for (int j = 0; j < lookupFloat.size(); j++) { @@ -285,7 +309,30 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou } float Ybar = sumOtu / (float) lookupFloat.size(); - //find the average + //find r value for each axis + for (int k = 0; k < averageAxes.size(); k++) { + + double r = 0.0; + double numerator = 0.0; + double denomTerm1 = 0.0; + double denomTerm2 = 0.0; + for (int j = 0; j < lookupFloat.size(); j++) { + float Yi = lookupFloat[j]->getAbundance(i); + float Xi = axes[lookupFloat[j]->getGroup()][k]; + + numerator += ((Xi - averageAxes[k]) * (Yi - Ybar)); + denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k])); + denomTerm2 += ((Yi - Ybar) * (Yi - Ybar)); + } + + double denom = (sqrt(denomTerm1) * sqrt(denomTerm2)); + + r = numerator / denom; + + out << r << '\t'; + } + + out << endl; } return 0; @@ -296,6 +343,66 @@ int CorrAxesCommand::calcPearson(map >& axes, ofstream& ou } } //********************************************************************************************************************** +int CorrAxesCommand::calcSpearman(map >& axes, ofstream& out) { + try { + + //find average of each axis - X + vector averageAxes; averageAxes.resize(numaxes, 0.0); + for (map >::iterator it = axes.begin(); it != axes.end(); it++) { + vector temp = it->second; + for (int i = 0; i < temp.size(); i++) { + averageAxes[i] += temp[i]; + } + } + + for (int i = 0; i < averageAxes.size(); i++) { averageAxes[i] = averageAxes[i] / (float) axes.size(); } + + //for each otu + for (int i = 0; i < lookupFloat[0]->getNumBins(); i++) { + + out << i+1 << '\t'; + + //find the averages this otu - Y + float sumOtu = 0.0; + for (int j = 0; j < lookupFloat.size(); j++) { + sumOtu += lookupFloat[j]->getAbundance(i); + } + float Ybar = sumOtu / (float) lookupFloat.size(); + + //find r value for each axis + for (int k = 0; k < averageAxes.size(); k++) { + + double r = 0.0; + double numerator = 0.0; + double denomTerm1 = 0.0; + double denomTerm2 = 0.0; + for (int j = 0; j < lookupFloat.size(); j++) { + float Yi = lookupFloat[j]->getAbundance(i); + float Xi = axes[lookupFloat[j]->getGroup()][k]; + + numerator += ((Xi - averageAxes[k]) * (Yi - Ybar)); + denomTerm1 += ((Xi - averageAxes[k]) * (Xi - averageAxes[k])); + denomTerm2 += ((Yi - Ybar) * (Yi - Ybar)); + } + + double denom = (sqrt(denomTerm1 * denomTerm2)); + + r = numerator / denom; + + out << r << '\t'; + } + + out << endl; + } + + return 0; + } + catch(exception& e) { + m->errorOut(e, "CorrAxesCommand", "calcSpearman"); + exit(1); + } +} +//********************************************************************************************************************** int CorrAxesCommand::getShared(){ try { InputData* input = new InputData(sharedfile, "sharedfile"); @@ -547,7 +654,6 @@ map > CorrAxesCommand::readAxes(){ headerLine = headerLine.substr(pos+4); }else { done = true; } } - cout << "count = " << count << endl; while (!in.eof()) { diff --git a/corraxescommand.h b/corraxescommand.h index 507399c..01c8d2b 100644 --- a/corraxescommand.h +++ b/corraxescommand.h @@ -29,7 +29,7 @@ public: private: GlobalData* globaldata; - string axesfile, sharedfile, relabundfile, groups, label, inputFileName, outputDir, method; + string axesfile, sharedfile, relabundfile, metadatafile, groups, label, inputFileName, outputDir, method; bool abort, pickedGroups; int numaxes; set names; @@ -45,6 +45,7 @@ private: int eliminateZeroOTUS(vector&); map > readAxes(); int calcPearson(map >&, ofstream&); + int calcSpearman(map >&, ofstream&); }; diff --git a/decalc.cpp b/decalc.cpp index 9c84429..ff68ca5 100644 --- a/decalc.cpp +++ b/decalc.cpp @@ -688,6 +688,7 @@ vector DeCalculator::findClosest(Sequence* querySeq, vector seqsMatches; + vector distsLeft; vector distsRight; @@ -759,14 +760,14 @@ vector DeCalculator::findClosest(Sequence* querySeq, vectorgetDist(); SeqDist subjectLeft; - subjectLeft.seq = db[j]; + subjectLeft.seq = NULL; subjectLeft.dist = distLeft; subjectLeft.index = j; distsLeft.push_back(subjectLeft); SeqDist subjectRight; - subjectRight.seq = db[j]; + subjectRight.seq = NULL; subjectRight.dist = distRight; subjectRight.index = j; @@ -790,18 +791,18 @@ vector DeCalculator::findClosest(Sequence* querySeq, vectorgetName()); + it = seen.find(db[distsLeft[i].index]->getName()); if (it == seen.end()) { dists.push_back(distsLeft[i]); - seen[distsLeft[i].seq->getName()] = distsLeft[i].seq->getName(); + seen[db[distsLeft[i].index]->getName()] = db[distsLeft[i].index]->getName(); lastLeft = distsLeft[i].dist; } //add right if you havent already - it = seen.find(distsRight[i].seq->getName()); + it = seen.find(db[distsRight[i].index]->getName()); if (it == seen.end()) { dists.push_back(distsRight[i]); - seen[distsRight[i].seq->getName()] = distsRight[i].seq->getName(); + seen[db[distsRight[i].index]->getName()] = db[distsRight[i].index]->getName(); lastRight = distsRight[i].dist; } @@ -829,11 +830,13 @@ vector DeCalculator::findClosest(Sequence* querySeq, vectorgetName() << '\t' << dists[i].dist << endl; - Sequence* temp = new Sequence(dists[i].seq->getName(), dists[i].seq->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. + Sequence* temp = new Sequence(db[dists[i].index]->getName(), db[dists[i].index]->getAligned()); //have to make a copy so you can trim and filter without stepping on eachother. seqsMatches.push_back(temp); indexes.push_back(dists[i].index); } + dists.clear(); distsLeft.clear(); distsRight.clear(); + return seqsMatches; } catch(exception& e) { @@ -951,7 +954,7 @@ map DeCalculator::trimSeqs(Sequence* query, vector topMatch map trimmedPos; //check to make sure that is not whole seq if ((rearPos - frontPos - 1) <= 0) { - m->mothurOut("[ERROR]: when I trim your sequences, the entire sequence is trimmed."); m->mothurOutEndLine(); + m->mothurOut("[ERROR]: when I trim " + query->getName() + ", the entire sequence is trimmed. Skipping."); m->mothurOutEndLine(); query->setAligned(""); //trim topMatches for (int i = 0; i < topMatches.size(); i++) { diff --git a/mothur b/mothur index 3ee5a5a..6347591 100755 Binary files a/mothur and b/mothur differ