2 * classifyotucommand.cpp
5 * Created by westcott on 6/1/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "classifyotucommand.h"
11 #include "phylotree.h"
13 //**********************************************************************************************************************
14 vector<string> ClassifyOtuCommand::getValidParameters(){
16 string AlignArray[] = {"list","label","name","taxonomy","cutoff","probs","outputdir","inputdir"};
17 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
21 m->errorOut(e, "ClassifyOtuCommand", "getValidParameters");
25 //**********************************************************************************************************************
26 ClassifyOtuCommand::ClassifyOtuCommand(){
29 //initialize outputTypes
30 vector<string> tempOutNames;
31 outputTypes["constaxonomy"] = tempOutNames;
34 m->errorOut(e, "ClassifyOtuCommand", "ClassifyOtuCommand");
38 //**********************************************************************************************************************
39 vector<string> ClassifyOtuCommand::getRequiredParameters(){
41 string Array[] = {"list","taxonomy"};
42 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
46 m->errorOut(e, "ClassifyOtuCommand", "getRequiredParameters");
50 //**********************************************************************************************************************
51 vector<string> ClassifyOtuCommand::getRequiredFiles(){
53 vector<string> myArray;
57 m->errorOut(e, "ClassifyOtuCommand", "getRequiredFiles");
61 //**********************************************************************************************************************
62 ClassifyOtuCommand::ClassifyOtuCommand(string option) {
68 //allow user to run help
69 if (option == "help") {
72 //valid paramters for this command
73 string Array[] = {"list","label","name","taxonomy","cutoff","probs","outputdir","inputdir"};
74 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
76 OptionParser parser(option);
77 map<string, string> parameters = parser.getParameters();
79 ValidParameters validParameter;
80 map<string, string>::iterator it;
82 //check to make sure all parameters are valid for command
83 for (it = parameters.begin(); it != parameters.end(); it++) {
84 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
87 //initialize outputTypes
88 vector<string> tempOutNames;
89 outputTypes["constaxonomy"] = tempOutNames;
91 //if the user changes the input directory command factory will send this info to us in the output parameter
92 string inputDir = validParameter.validFile(parameters, "inputdir", false);
93 if (inputDir == "not found"){ inputDir = ""; }
96 it = parameters.find("list");
97 //user has given a template file
98 if(it != parameters.end()){
99 path = m->hasPath(it->second);
100 //if the user has not given a path then, add inputdir. else leave path alone.
101 if (path == "") { parameters["list"] = inputDir + it->second; }
104 it = parameters.find("name");
105 //user has given a template file
106 if(it != parameters.end()){
107 path = m->hasPath(it->second);
108 //if the user has not given a path then, add inputdir. else leave path alone.
109 if (path == "") { parameters["name"] = inputDir + it->second; }
112 it = parameters.find("taxonomy");
113 //user has given a template file
114 if(it != parameters.end()){
115 path = m->hasPath(it->second);
116 //if the user has not given a path then, add inputdir. else leave path alone.
117 if (path == "") { parameters["taxonomy"] = inputDir + it->second; }
122 //if the user changes the output directory command factory will send this info to us in the output parameter
123 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
125 //check for required parameters
126 listfile = validParameter.validFile(parameters, "list", true);
127 if (listfile == "not found") { m->mothurOut("list is a required parameter for the classify.otu command."); m->mothurOutEndLine(); abort = true; }
128 else if (listfile == "not open") { abort = true; }
130 taxfile = validParameter.validFile(parameters, "taxonomy", true);
131 if (taxfile == "not found") { m->mothurOut("taxonomy is a required parameter for the classify.otu command."); m->mothurOutEndLine(); abort = true; }
132 else if (taxfile == "not open") { abort = true; }
134 namefile = validParameter.validFile(parameters, "name", true);
135 if (namefile == "not open") { abort = true; }
136 else if (namefile == "not found") { namefile = ""; }
138 //check for optional parameter and set defaults
139 // ...at some point should added some additional type checking...
140 label = validParameter.validFile(parameters, "label", false);
141 if (label == "not found") { label = ""; allLines = 1; }
143 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
144 else { allLines = 1; }
147 string temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found") { temp = "51"; }
148 convert(temp, cutoff);
150 temp = validParameter.validFile(parameters, "probs", false); if (temp == "not found"){ temp = "true"; }
151 probs = m->isTrue(temp);
154 if ((cutoff < 51) || (cutoff > 100)) { m->mothurOut("cutoff must be above 50, and no greater than 100."); m->mothurOutEndLine(); abort = true; }
158 catch(exception& e) {
159 m->errorOut(e, "ClassifyOtuCommand", "ClassifyOtuCommand");
164 //**********************************************************************************************************************
166 void ClassifyOtuCommand::help(){
168 m->mothurOut("The classify.otu command parameters are list, taxonomy, name, cutoff, label and probs. The taxonomy and list parameters are required.\n");
169 m->mothurOut("The name parameter allows you add a names file with your taxonomy file.\n");
170 m->mothurOut("The label parameter allows you to select what distance levels you would like a output files created for, and is separated by dashes.\n");
171 m->mothurOut("The default value for label is all labels in your inputfile.\n");
172 m->mothurOut("The cutoff parameter allows you to specify a consensus confidence threshold for your taxonomy. The default is 51, meaning 51%. Cutoff cannot be below 51.\n");
173 m->mothurOut("The probs parameter shuts off the outputting of the consensus confidence results. The default is true, meaning you want the confidence to be shown.\n");
174 m->mothurOut("The classify.otu command should be in the following format: classify.otu(taxonomy=yourTaxonomyFile, list=yourListFile, name=yourNamesFile, label=yourLabels).\n");
175 m->mothurOut("Example classify.otu(taxonomy=abrecovery.silva.full.taxonomy, list=abrecovery.fn.list, label=0.10).\n");
176 m->mothurOut("Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListFile).\n\n");
178 catch(exception& e) {
179 m->errorOut(e, "ClassifyOtuCommand", "help");
184 //**********************************************************************************************************************
186 ClassifyOtuCommand::~ClassifyOtuCommand(){}
188 //**********************************************************************************************************************
190 int ClassifyOtuCommand::execute(){
193 if (abort == true) { return 0; }
195 //if user gave a namesfile then use it
196 if (namefile != "") { readNamesFile(); }
198 //read taxonomy file and save in map for easy access in building bin trees
201 if (m->control_pressed) { return 0; }
203 input = new InputData(listfile, "list");
204 list = input->getListVector();
205 string lastLabel = list->getLabel();
207 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
208 set<string> processedLabels;
209 set<string> userLabels = labels;
211 if (m->control_pressed) { outputTypes.clear(); delete input; delete list; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
213 while((list != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
215 if (allLines == 1 || labels.count(list->getLabel()) == 1){
217 m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
219 if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; delete list; return 0; }
221 processedLabels.insert(list->getLabel());
222 userLabels.erase(list->getLabel());
225 if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
226 string saveLabel = list->getLabel();
229 list = input->getListVector(lastLabel);
230 m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
234 if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; delete list; return 0; }
236 processedLabels.insert(list->getLabel());
237 userLabels.erase(list->getLabel());
239 //restore real lastlabel to save below
240 list->setLabel(saveLabel);
243 lastLabel = list->getLabel();
246 list = input->getListVector();
249 //output error messages about any remaining user labels
250 bool needToRun = false;
251 for (set<string>::iterator it = userLabels.begin(); it != userLabels.end(); it++) {
252 m->mothurOut("Your file does not include the label " + (*it));
253 if (processedLabels.count(lastLabel) != 1) {
254 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
257 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
261 //run last label if you need to
262 if (needToRun == true) {
263 if (list != NULL) { delete list; }
264 list = input->getListVector(lastLabel);
265 m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
270 if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; delete list; return 0; }
275 if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
277 m->mothurOutEndLine();
278 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
279 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
280 m->mothurOutEndLine();
284 catch(exception& e) {
285 m->errorOut(e, "ClassifyOtuCommand", "execute");
290 //**********************************************************************************************************************
291 int ClassifyOtuCommand::readNamesFile() {
295 m->openInputFile(namefile, inNames);
299 while(!inNames.eof()){
300 inNames >> name; //read from first column A
301 inNames >> names; //read from second column A,B,C,D
304 //parse names into vector
305 vector<string> theseNames;
306 m->splitAtComma(names, theseNames);
308 for (int i = 0; i < theseNames.size(); i++) { nameMap[theseNames[i]] = name; }
310 if (m->control_pressed) { inNames.close(); nameMap.clear(); return 0; }
316 catch(exception& e) {
317 m->errorOut(e, "ClassifyOtuCommand", "readNamesFile");
321 //**********************************************************************************************************************
322 int ClassifyOtuCommand::readTaxonomyFile() {
326 m->openInputFile(taxfile, in);
334 //are there confidence scores, if so remove them
335 if (tax.find_first_of('(') != -1) { removeConfidences(tax); }
339 if (m->control_pressed) { in.close(); taxMap.clear(); return 0; }
345 catch(exception& e) {
346 m->errorOut(e, "ClassifyOtuCommand", "readTaxonomyFile");
350 //**********************************************************************************************************************
351 string ClassifyOtuCommand::findConsensusTaxonomy(int bin, ListVector* thisList, int& size) {
354 vector<string> names;
355 map<string, string>::iterator it;
356 map<string, string>::iterator it2;
358 //parse names into vector
359 string binnames = thisList->get(bin);
360 m->splitAtComma(binnames, names);
362 //create a tree containing sequences from this bin
363 PhyloTree* phylo = new PhyloTree();
366 for (int i = 0; i < names.size(); i++) {
368 //if namesfile include the names
369 if (namefile != "") {
371 //is this sequence in the name file - namemap maps seqName -> repSeqName
372 it2 = nameMap.find(names[i]);
374 if (it2 == nameMap.end()) { //this name is not in name file, skip it
375 m->mothurOut(names[i] + " is not in your name file. I will not include it in the consensus."); m->mothurOutEndLine();
378 //is this sequence in the taxonomy file - look for repSeqName since we are assuming the taxonomy file is unique
379 it = taxMap.find(it2->second);
381 if (it == taxMap.end()) { //this name is not in taxonomy file, skip it
383 if (names[i] != it2->second) { m->mothurOut(names[i] + " is represented by " + it2->second + " and is not in your taxonomy file. I will not include it in the consensus."); m->mothurOutEndLine(); }
384 else { m->mothurOut(names[i] + " is not in your taxonomy file. I will not include it in the consensus."); m->mothurOutEndLine(); }
388 phylo->addSeqToTree(names[i], it->second);
394 //is this sequence in the taxonomy file - look for repSeqName since we are assuming the taxonomy file is unique
395 it = taxMap.find(names[i]);
397 if (it == taxMap.end()) { //this name is not in taxonomy file, skip it
398 m->mothurOut(names[i] + " is not in your taxonomy file. I will not include it in the consensus."); m->mothurOutEndLine();
401 phylo->addSeqToTree(names[i], it->second);
407 if (m->control_pressed) { delete phylo; return conTax; }
412 phylo->assignHeirarchyIDs(0);
414 TaxNode currentNode = phylo->get(0);
417 while (currentNode.children.size() != 0) { //you still have more to explore
420 int bestChildSize = 0;
422 //go through children
423 for (map<string, int>::iterator itChild = currentNode.children.begin(); itChild != currentNode.children.end(); itChild++) {
425 TaxNode temp = phylo->get(itChild->second);
427 //select child with largest accesions - most seqs assigned to it
428 if (temp.accessions.size() > bestChildSize) {
429 bestChild = phylo->get(itChild->second);
430 bestChildSize = temp.accessions.size();
435 //is this taxonomy above cutoff
436 int consensusConfidence = ceil((bestChildSize / (float) size) * 100);
438 if (consensusConfidence >= cutoff) { //if yes, add it
440 conTax += bestChild.name + "(" + toString(consensusConfidence) + ");";
442 conTax += bestChild.name + ";";
449 currentNode = bestChild;
453 if (conTax == "") { conTax = "no_consensus;"; }
460 catch(exception& e) {
461 m->errorOut(e, "ClassifyOtuCommand", "findConsensusTaxonomy");
466 //**********************************************************************************************************************
467 int ClassifyOtuCommand::process(ListVector* processList) {
473 if (outputDir == "") { outputDir += m->hasPath(listfile); }
476 string outputFile = outputDir + m->getRootName(m->getSimpleName(listfile)) + processList->getLabel() + ".cons.taxonomy";
477 m->openOutputFile(outputFile, out);
478 outputNames.push_back(outputFile); outputTypes["constaxonomy"].push_back(outputFile);
480 out << "OTU\tSize\tTaxonomy" << endl;
482 //for each bin in the list vector
483 for (int i = 0; i < processList->getNumBins(); i++) {
485 conTax = findConsensusTaxonomy(i, processList, size);
487 if (m->control_pressed) { out.close(); return 0; }
489 //output to new names file
490 out << (i+1) << '\t' << size << '\t' << conTax << endl;
498 catch(exception& e) {
499 m->errorOut(e, "ClassifyOtuCommand", "process");
503 /**************************************************************************************************/
504 void ClassifyOtuCommand::removeConfidences(string& tax) {
510 while (tax.find_first_of(';') != -1) {
512 taxon = tax.substr(0,tax.find_first_of(';'));
514 int pos = taxon.find_first_of('(');
516 taxon = taxon.substr(0, pos); //rip off confidence
521 tax = tax.substr(tax.find_first_of(';')+1, tax.length());
527 catch(exception& e) {
528 m->errorOut(e, "ClassifyOtuCommand", "removeConfidences");
532 //**********************************************************************************************************************