X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=subsamplecommand.cpp;h=3e954759b078ff9371f5c28b64b3de86ef3c89cc;hp=cfc7b1d9283779fdf9cc4a21a60c3f6d49d5177b;hb=cf9987b67aa49777a4c91c2d21f96e58bf17aa82;hpb=1a5c2356c1b955c6ec024b2baf9f46377ee7c72e diff --git a/subsamplecommand.cpp b/subsamplecommand.cpp index cfc7b1d..3e95475 100644 --- a/subsamplecommand.cpp +++ b/subsamplecommand.cpp @@ -10,6 +10,7 @@ #include "subsamplecommand.h" #include "sharedutilities.h" #include "deconvolutecommand.h" +#include "getseqscommand.h" #include "subsample.h" //********************************************************************************************************************** @@ -17,6 +18,7 @@ vector SubSampleCommand::setParameters(){ try { CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "FLSSR", "none","fasta",false,false,true); parameters.push_back(pfasta); CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","name",false,false,true); parameters.push_back(pname); + CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "none", "none", "none","taxonomy",false,false,true); parameters.push_back(ptaxonomy); CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","count",false,false,true); parameters.push_back(pcount); CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","group",false,false,true); parameters.push_back(pgroup); CommandParameter plist("list", "InputTypes", "", "", "none", "FLSSR", "none","list",false,false,true); parameters.push_back(plist); @@ -44,7 +46,7 @@ string SubSampleCommand::getHelpString(){ try { string helpString = ""; helpString += "The sub.sample command is designed to be used as a way to normalize your data, or create a smaller set from your original set.\n"; - helpString += "The sub.sample command parameters are fasta, name, list, group, count, rabund, sabund, shared, groups, size, persample and label. You must provide a fasta, list, sabund, rabund or shared file as an input file.\n"; + helpString += "The sub.sample command parameters are fasta, name, list, group, count, rabund, sabund, shared, taxonomy, groups, size, persample and label. You must provide a fasta, list, sabund, rabund or shared file as an input file.\n"; helpString += "The namefile is only used with the fasta file, not with the listfile, because the list file should contain all sequences.\n"; helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like included. The group names are separated by dashes.\n"; helpString += "The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n"; @@ -70,11 +72,12 @@ string SubSampleCommand::getOutputPattern(string type) { string pattern = ""; if (type == "fasta") { pattern = "[filename],subsample,[extension]"; } - else if (type == "sabund") { pattern = "[filename],subsample,[extension]"; } + else if (type == "sabund") { pattern = "[filename],subsample,[extension]"; } else if (type == "name") { pattern = "[filename],subsample,[extension]"; } else if (type == "group") { pattern = "[filename],subsample,[extension]"; } else if (type == "count") { pattern = "[filename],subsample,[extension]"; } - else if (type == "list") { pattern = "[filename],subsample,[extension]"; } + else if (type == "list") { pattern = "[filename],[distance],subsample,[extension]"; } + else if (type == "taxonomy") { pattern = "[filename],subsample,[extension]"; } else if (type == "shared") { pattern = "[filename],[distance],subsample,[extension]"; } else if (type == "rabund") { pattern = "[filename],subsample,[extension]"; } else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } @@ -100,6 +103,7 @@ SubSampleCommand::SubSampleCommand(){ outputTypes["name"] = tempOutNames; outputTypes["group"] = tempOutNames; outputTypes["count"] = tempOutNames; + outputTypes["taxonomy"] = tempOutNames; } catch(exception& e) { m->errorOut(e, "SubSampleCommand", "GetRelAbundCommand"); @@ -140,6 +144,7 @@ SubSampleCommand::SubSampleCommand(string option) { outputTypes["name"] = tempOutNames; outputTypes["group"] = tempOutNames; outputTypes["count"] = tempOutNames; + outputTypes["taxonomy"] = tempOutNames; //if the user changes the output directory command factory will send this info to us in the output parameter outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; } @@ -212,6 +217,14 @@ SubSampleCommand::SubSampleCommand(string option) { //if the user has not given a path then, add inputdir. else leave path alone. if (path == "") { parameters["count"] = inputDir + it->second; } } + + it = parameters.find("taxonomy"); + //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["taxonomy"] = inputDir + it->second; } + } } //check for required parameters @@ -249,13 +262,18 @@ SubSampleCommand::SubSampleCommand(string option) { if (groupfile == "not open") { groupfile = ""; abort = true; } else if (groupfile == "not found") { groupfile = ""; } else { m->setGroupFile(groupfile); } + + taxonomyfile = validParameter.validFile(parameters, "taxonomy", true); + if (taxonomyfile == "not open") { taxonomyfile = ""; abort = true; } + else if (taxonomyfile == "not found") { taxonomyfile = ""; } + else { m->setTaxonomyFile(taxonomyfile); } countfile = validParameter.validFile(parameters, "count", true); if (countfile == "not open") { countfile = ""; abort = true; } else if (countfile == "not found") { countfile = ""; } else { m->setCountTableFile(countfile); - ct.readTable(countfile); + ct.readTable(countfile, true, false); } if ((namefile != "") && (countfile != "")) { @@ -297,7 +315,10 @@ SubSampleCommand::SubSampleCommand(string option) { } } - if ((namefile != "") && (fastafile == "")) { m->mothurOut("You may only use a namefile with a fastafile."); m->mothurOutEndLine(); abort = true; } + if ((namefile != "") && ((fastafile == "") && (taxonomyfile == ""))) { m->mothurOut("You may only use a name file with a fasta file or taxonomy file."); m->mothurOutEndLine(); abort = true; } + + if ((taxonomyfile != "") && ((fastafile == "") && (listfile == ""))) { m->mothurOut("You may only use a taxonomyfile with a fastafile or listfile."); m->mothurOutEndLine(); abort = true; } + if ((fastafile == "") && (listfile == "") && (sabundfile == "") && (rabundfile == "") && (sharedfile == "")) { m->mothurOut("You must provide a fasta, list, sabund, rabund or shared file as an input file."); m->mothurOutEndLine(); abort = true; } @@ -389,6 +410,10 @@ int SubSampleCommand::execute(){ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); } } + itTypes = outputTypes.find("taxonomy"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); } + } m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); @@ -527,7 +552,7 @@ int SubSampleCommand::getSubSampleFasta() { else{ itGroupCounts = groupCounts.find(group); if (itGroupCounts != groupCounts.end()) { - if (groupCounts[group] < size) { subset.insert(names[j]); groupCounts[group]++; } + if (itGroupCounts->second < size) { subset.insert(names[j]); (itGroupCounts->second)++; } } } } @@ -632,6 +657,7 @@ int SubSampleCommand::getSubSampleFasta() { } in.close(); out.close(); + if (count != subset.size()) { m->mothurOut("[ERROR]: The subset selected contained " + toString(subset.size()) + " sequences, but I only found " + toString(count) + " of those in the fastafile."); m->mothurOutEndLine(); @@ -665,7 +691,33 @@ int SubSampleCommand::getSubSampleFasta() { m->mothurOut("/******************************************/"); m->mothurOutEndLine(); m->mothurOut("Done."); m->mothurOutEndLine(); - } + + if (taxonomyfile != "") { + set tempSubset; + //get new unique names from fasta file + //read through fasta file outputting only the names on the subsample list after deconvolute + ifstream in2; + m->openInputFile(outputFileName, in2); + + while (!in2.eof()) { + Sequence seq(in2); m->gobble(in2); + if (seq.getName() != "") { + tempSubset.insert(seq.getName()); + } + } + in2.close(); + + //send that list to getTax + int tcount = getTax(tempSubset); + if (tcount != tempSubset.size()) { m->mothurOut("[ERROR]: subsampled fasta file contains " + toString(tempSubset.size()) + " sequences, but I only found " + toString(tcount) + " in your taxonomy file, please correct."); m->mothurOutEndLine(); } + } + }else { + if (taxonomyfile != "") { + int tcount = getTax(subset); + if (tcount != subset.size()) { m->mothurOut("[ERROR]: subsampled fasta file contains " + toString(subset.size()) + " sequences, but I only found " + toString(tcount) + " in your taxonomy file, please correct."); m->mothurOutEndLine(); } + + } //should only contain uniques. + } outputTypes["fasta"].push_back(outputFileName); outputNames.push_back(outputFileName); @@ -895,7 +947,7 @@ int SubSampleCommand::processShared(vector& thislookup) { try { //save mothurOut's binLabels to restore for next label - vector saveBinLabels = m->currentBinLabels; + vector saveBinLabels = m->currentSharedBinLabels; string thisOutputDir = outputDir; if (outputDir == "") { thisOutputDir += m->hasPath(sharedfile); } @@ -913,7 +965,7 @@ int SubSampleCommand::processShared(vector& thislookup) { m->openOutputFile(outputFileName, out); outputTypes["shared"].push_back(outputFileName); outputNames.push_back(outputFileName); - m->currentBinLabels = subsampledLabels; + m->currentSharedBinLabels = subsampledLabels; thislookup[0]->printHeaders(out); @@ -925,7 +977,7 @@ int SubSampleCommand::processShared(vector& thislookup) { //save mothurOut's binLabels to restore for next label - m->currentBinLabels = saveBinLabels; + m->currentSharedBinLabels = saveBinLabels; return 0; @@ -938,22 +990,14 @@ int SubSampleCommand::processShared(vector& thislookup) { //********************************************************************************************************************** int SubSampleCommand::getSubSampleList() { try { - - string thisOutputDir = outputDir; - if (outputDir == "") { thisOutputDir += m->hasPath(listfile); } - map variables; - variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(listfile)); - variables["[extension]"] = m->getExtension(listfile); - string outputFileName = getOutputFileName("list", variables); - ofstream out; - m->openOutputFile(outputFileName, out); - outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName); - + + if (namefile != "") { m->readNames(namefile, nameMap); } + InputData* input = new InputData(listfile, "list"); ListVector* list = input->getListVector(); string lastLabel = list->getLabel(); - //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label. + //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label. set processedLabels; set userLabels = labels; @@ -975,7 +1019,7 @@ int SubSampleCommand::getSubSampleList() { //file mismatch quit if (list->getNumSeqs() != groupMap.getNumSeqs()) { m->mothurOut("[ERROR]: your list file contains " + toString(list->getNumSeqs()) + " sequences, and your groupfile contains " + toString(groupMap.getNumSeqs()) + ", please correct."); - m->mothurOutEndLine(); delete list; delete input; out.close(); outGroup.close(); return 0; + m->mothurOutEndLine(); delete list; delete input; outGroup.close(); return 0; } }else if (countfile != "") { if (ct.hasGroupInfo()) { @@ -1139,13 +1183,13 @@ int SubSampleCommand::getSubSampleList() { //as long as you are not at the end of the file or done wih the lines you want while((list != NULL) && ((allLines == 1) || (userLabels.size() != 0))) { - if (m->control_pressed) { delete list; delete input; out.close(); return 0; } + if (m->control_pressed) { delete list; delete input; return 0; } if(allLines == 1 || labels.count(list->getLabel()) == 1){ m->mothurOut(list->getLabel()); m->mothurOutEndLine(); - processList(list, out, subset); + processList(list, subset); processedLabels.insert(list->getLabel()); userLabels.erase(list->getLabel()); @@ -1159,7 +1203,7 @@ int SubSampleCommand::getSubSampleList() { list = input->getListVector(lastLabel); m->mothurOut(list->getLabel()); m->mothurOutEndLine(); - processList(list, out, subset); + processList(list, subset); processedLabels.insert(list->getLabel()); userLabels.erase(list->getLabel()); @@ -1177,7 +1221,7 @@ int SubSampleCommand::getSubSampleList() { } - if (m->control_pressed) { if (list != NULL) { delete list; } delete input; out.close(); return 0; } + if (m->control_pressed) { if (list != NULL) { delete list; } delete input; return 0; } //output error messages about any remaining user labels set::iterator it; @@ -1200,14 +1244,77 @@ int SubSampleCommand::getSubSampleList() { m->mothurOut(list->getLabel()); m->mothurOutEndLine(); - processList(list, out, subset); + processList(list, subset); delete list; list = NULL; } - out.close(); if (list != NULL) { delete list; } delete input; + + if (taxonomyfile != "") { + if (namefile == "") { + InputData input(listfile, "list"); + ListVector* list = input.getListVector(); + string lastLabel = list->getLabel(); + + for (int i = 0; i < list->getNumBins(); i++) { + vector temp; + string bin = list->get(i); + m->splitAtComma(bin, temp); + for (int j = 0; j < temp.size(); j++) { vector tempFakeOut; tempFakeOut.push_back(temp[j]); nameMap[temp[j]] = tempFakeOut; } + } + delete list; + + int tcount = getTax(subset); + if (tcount != subset.size()) { m->mothurOut("[ERROR]: subsampled list file contains " + toString(subset.size()) + " sequences, but I only found " + toString(tcount) + " in your taxonomy file, did you forget a name file? Please correct."); m->mothurOutEndLine(); } + }else { + string tempAccnos = "temp.accnos"; + ofstream outAccnos; + m->openOutputFile(tempAccnos, outAccnos); + for (set::iterator it = subset.begin(); it != subset.end(); it++) { outAccnos << *it << endl; } + outAccnos.close(); + + m->mothurOut("Sampling taxonomy and name file... "); m->mothurOutEndLine(); + string thisNameOutputDir = outputDir; + if (outputDir == "") { thisNameOutputDir += m->hasPath(namefile); } + map variables; + variables["[filename]"] = thisNameOutputDir + m->getRootName(m->getSimpleName(namefile)); + variables["[extension]"] = m->getExtension(namefile); + string outputNameFileName = getOutputFileName("name", variables); + + string thisTaxOutputDir = outputDir; + if (outputDir == "") { thisTaxOutputDir += m->hasPath(taxonomyfile); } + variables["[filename]"] = thisTaxOutputDir + m->getRootName(m->getSimpleName(taxonomyfile)); + variables["[extension]"] = m->getExtension(taxonomyfile); + string outputTaxFileName = getOutputFileName("taxonomy", variables); + + + //use unique.seqs to create new name and fastafile + string inputString = "dups=f, name=" + namefile + ", taxonomy=" + taxonomyfile + ", accnos=" + tempAccnos; + m->mothurOut("/******************************************/"); m->mothurOutEndLine(); + m->mothurOut("Running command: get.seqs(" + inputString + ")"); m->mothurOutEndLine(); + m->mothurCalling = true; + + Command* getCommand = new GetSeqsCommand(inputString); + getCommand->execute(); + + map > filenames = getCommand->getOutputFiles(); + + delete getCommand; + m->mothurCalling = false; + + m->renameFile(filenames["name"][0], outputNameFileName); + m->renameFile(filenames["taxonomy"][0], outputTaxFileName); + + outputTypes["name"].push_back(outputNameFileName); outputNames.push_back(outputNameFileName); + outputNames.push_back(outputTaxFileName); outputTypes["taxonomy"].push_back(outputTaxFileName); + + m->mothurOut("/******************************************/"); m->mothurOutEndLine(); + + m->mothurOut("Done."); m->mothurOutEndLine(); + } + } return 0; @@ -1218,14 +1325,26 @@ int SubSampleCommand::getSubSampleList() { } } //********************************************************************************************************************** -int SubSampleCommand::processList(ListVector*& list, ofstream& out, set& subset) { +int SubSampleCommand::processList(ListVector*& list, set& subset) { try { - + string thisOutputDir = outputDir; + if (outputDir == "") { thisOutputDir += m->hasPath(listfile); } + map variables; + variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(listfile)); + variables["[extension]"] = m->getExtension(listfile); + variables["[distance]"] = list->getLabel(); + string outputFileName = getOutputFileName("list", variables); + ofstream out; + m->openOutputFile(outputFileName, out); + outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName); + int numBins = list->getNumBins(); ListVector* temp = new ListVector(); temp->setLabel(list->getLabel()); + vector binLabels = list->getLabels(); + vector newLabels; for (int i = 0; i < numBins; i++) { if (m->control_pressed) { break; } @@ -1241,15 +1360,19 @@ int SubSampleCommand::processList(ListVector*& list, ofstream& out, set& if (newNames != "") { newNames = newNames.substr(0, newNames.length()-1); //rip off extra comma temp->push_back(newNames); + newLabels.push_back(binLabels[i]); } } + temp->setLabels(newLabels); delete list; list = temp; - if (m->control_pressed) { return 0; } + if (m->control_pressed) { out.close(); return 0; } + list->printHeaders(out); list->print(out); + out.close(); return 0; @@ -1579,7 +1702,63 @@ int SubSampleCommand::processSabund(SAbundVector*& sabund, ofstream& out) { m->errorOut(e, "SubSampleCommand", "processSabund"); exit(1); } -} +} + +//********************************************************************************************************************** +int SubSampleCommand::getTax(set& subset) { + try { + + string thisTaxOutputDir = outputDir; + if (outputDir == "") { thisTaxOutputDir += m->hasPath(taxonomyfile); } + map variables; + variables["[filename]"] = thisTaxOutputDir + m->getRootName(m->getSimpleName(taxonomyfile)); + variables["[extension]"] = m->getExtension(taxonomyfile); + string outputTaxFileName = getOutputFileName("taxonomy", variables); + ofstream outTax; + m->openOutputFile(outputTaxFileName, outTax); + outputNames.push_back(outputTaxFileName); outputTypes["taxonomy"].push_back(outputTaxFileName); + + //read through fasta file outputting only the names on the subsample list + ifstream inTax; + m->openInputFile(taxonomyfile, inTax); + + string tname, tax; + int tcount = 0; + map >::iterator itNameMap; + + while(!inTax.eof()){ + + if (m->control_pressed) { inTax.close(); outTax.close(); return 0; } + + inTax >> tname; m->gobble(inTax); //read from first column + inTax >> tax; m->gobble(inTax); //read from second column + + //does the subset contain a sequence that this sequence represents + itNameMap = nameMap.find(tname); + if (itNameMap != nameMap.end()) { + vector nameRepresents = itNameMap->second; + + + for (int i = 0; i < nameRepresents.size(); i++){ + if (subset.count(nameRepresents[i]) != 0) { + outTax << nameRepresents[i] << '\t' << tax << endl; + tcount++; + + } + } + }else{ m->mothurOut("[ERROR]: " + tname + " is missing, please correct."); m->mothurOutEndLine(); } + } + inTax.close(); + outTax.close(); + + return tcount; + } + catch(exception& e) { + m->errorOut(e, "SubSampleCommand", "getTax"); + exit(1); + } +} + //**********************************************************************************************************************