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"
14 //**********************************************************************************************************************
15 ClassifyOtuCommand::ClassifyOtuCommand(string option) {
21 //allow user to run help
22 if (option == "help") {
25 //valid paramters for this command
26 string Array[] = {"list","label","name","taxonomy","cutoff","probs","outputdir","inputdir"};
27 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
29 OptionParser parser(option);
30 map<string, string> parameters = parser.getParameters();
32 ValidParameters validParameter;
33 map<string, string>::iterator it;
35 //check to make sure all parameters are valid for command
36 for (it = parameters.begin(); it != parameters.end(); it++) {
37 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
40 //if the user changes the input directory command factory will send this info to us in the output parameter
41 string inputDir = validParameter.validFile(parameters, "inputdir", false);
42 if (inputDir == "not found"){ inputDir = ""; }
45 it = parameters.find("list");
46 //user has given a template file
47 if(it != parameters.end()){
48 path = hasPath(it->second);
49 //if the user has not given a path then, add inputdir. else leave path alone.
50 if (path == "") { parameters["list"] = inputDir + it->second; }
53 it = parameters.find("name");
54 //user has given a template file
55 if(it != parameters.end()){
56 path = hasPath(it->second);
57 //if the user has not given a path then, add inputdir. else leave path alone.
58 if (path == "") { parameters["name"] = inputDir + it->second; }
61 it = parameters.find("taxonomy");
62 //user has given a template file
63 if(it != parameters.end()){
64 path = hasPath(it->second);
65 //if the user has not given a path then, add inputdir. else leave path alone.
66 if (path == "") { parameters["taxonomy"] = inputDir + it->second; }
71 //if the user changes the output directory command factory will send this info to us in the output parameter
72 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
74 //check for required parameters
75 listfile = validParameter.validFile(parameters, "list", true);
76 if (listfile == "not found") { m->mothurOut("list is a required parameter for the classify.otu command."); m->mothurOutEndLine(); abort = true; }
77 else if (listfile == "not open") { abort = true; }
79 taxfile = validParameter.validFile(parameters, "taxonomy", true);
80 if (taxfile == "not found") { m->mothurOut("taxonomy is a required parameter for the classify.otu command."); m->mothurOutEndLine(); abort = true; }
81 else if (taxfile == "not open") { abort = true; }
83 namefile = validParameter.validFile(parameters, "name", true);
84 if (namefile == "not open") { abort = true; }
85 else if (namefile == "not found") { namefile = ""; }
87 //check for optional parameter and set defaults
88 // ...at some point should added some additional type checking...
89 label = validParameter.validFile(parameters, "label", false);
90 if (label == "not found") { label = ""; allLines = 1; }
92 if(label != "all") { splitAtDash(label, labels); allLines = 0; }
93 else { allLines = 1; }
96 string temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found") { temp = "51"; }
97 convert(temp, cutoff);
99 temp = validParameter.validFile(parameters, "probs", false); if (temp == "not found"){ temp = "true"; }
100 probs = isTrue(temp);
103 if ((cutoff < 51) || (cutoff > 100)) { m->mothurOut("cutoff must be above 50, and no greater than 100."); m->mothurOutEndLine(); abort = true; }
107 catch(exception& e) {
108 m->errorOut(e, "ClassifyOtuCommand", "ClassifyOtuCommand");
113 //**********************************************************************************************************************
115 void ClassifyOtuCommand::help(){
117 m->mothurOut("The classify.otu command parameters are list, taxonomy, name, cutoff, label and probs. The taxonomy and list parameters are required.\n");
118 m->mothurOut("The name parameter allows you add a names file with your taxonomy file.\n");
119 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");
120 m->mothurOut("The default value for label is all labels in your inputfile.\n");
121 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");
122 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");
123 m->mothurOut("The classify.otu command should be in the following format: classify.otu(taxonomy=yourTaxonomyFile, list=yourListFile, name=yourNamesFile, label=yourLabels).\n");
124 m->mothurOut("Example classify.otu(taxonomy=abrecovery.silva.full.taxonomy, list=abrecovery.fn.list, label=0.10).\n");
125 m->mothurOut("Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListFile).\n\n");
127 catch(exception& e) {
128 m->errorOut(e, "ClassifyOtuCommand", "help");
133 //**********************************************************************************************************************
135 ClassifyOtuCommand::~ClassifyOtuCommand(){}
137 //**********************************************************************************************************************
139 int ClassifyOtuCommand::execute(){
142 if (abort == true) { return 0; }
144 //if user gave a namesfile then use it
145 if (namefile != "") { readNamesFile(); }
147 //read taxonomy file and save in map for easy access in building bin trees
150 if (m->control_pressed) { return 0; }
152 input = new InputData(listfile, "list");
153 list = input->getListVector();
154 string lastLabel = list->getLabel();
156 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
157 set<string> processedLabels;
158 set<string> userLabels = labels;
160 if (m->control_pressed) { delete input; delete list; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
162 while((list != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
164 if (allLines == 1 || labels.count(list->getLabel()) == 1){
166 m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
168 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; delete list; return 0; }
170 processedLabels.insert(list->getLabel());
171 userLabels.erase(list->getLabel());
174 if ((anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
175 string saveLabel = list->getLabel();
178 list = input->getListVector(lastLabel);
179 m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
183 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; delete list; return 0; }
185 processedLabels.insert(list->getLabel());
186 userLabels.erase(list->getLabel());
188 //restore real lastlabel to save below
189 list->setLabel(saveLabel);
192 lastLabel = list->getLabel();
195 list = input->getListVector();
198 //output error messages about any remaining user labels
199 bool needToRun = false;
200 for (set<string>::iterator it = userLabels.begin(); it != userLabels.end(); it++) {
201 m->mothurOut("Your file does not include the label " + (*it));
202 if (processedLabels.count(lastLabel) != 1) {
203 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
206 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
210 //run last label if you need to
211 if (needToRun == true) {
212 if (list != NULL) { delete list; }
213 list = input->getListVector(lastLabel);
214 m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
219 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } delete input; delete list; return 0; }
224 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
226 m->mothurOutEndLine();
227 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
228 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
229 m->mothurOutEndLine();
233 catch(exception& e) {
234 m->errorOut(e, "ClassifyOtuCommand", "execute");
239 //**********************************************************************************************************************
240 int ClassifyOtuCommand::readNamesFile() {
244 openInputFile(namefile, inNames);
249 inNames >> name; //read from first column A
250 inNames >> names; //read from second column A,B,C,D
253 nameMap[name] = names;
255 if (m->control_pressed) { inNames.close(); nameMap.clear(); return 0; }
261 catch(exception& e) {
262 m->errorOut(e, "ClassifyOtuCommand", "readNamesFile");
266 //**********************************************************************************************************************
267 int ClassifyOtuCommand::readTaxonomyFile() {
271 openInputFile(taxfile, in);
279 //are there confidence scores, if so remove them
280 if (tax.find_first_of('(') != -1) { removeConfidences(tax); }
284 if (m->control_pressed) { in.close(); taxMap.clear(); return 0; }
290 catch(exception& e) {
291 m->errorOut(e, "ClassifyOtuCommand", "readTaxonomyFile");
295 //**********************************************************************************************************************
296 string ClassifyOtuCommand::findConsensusTaxonomy(int bin, ListVector* thisList, int& size) {
299 vector<string> names;
300 map<string, string>::iterator it;
301 map<string, string>::iterator it2;
303 //parse names into vector
304 string binnames = thisList->get(bin);
305 splitAtComma(binnames, names);
307 //create a tree containing sequences from this bin
308 PhyloTree* phylo = new PhyloTree();
311 for (int i = 0; i < names.size(); i++) {
313 if (m->control_pressed) { delete phylo; return conTax; }
315 //is this sequence in the taxonomy file
316 it = taxMap.find(names[i]);
318 if (it == taxMap.end()) { //this name is not in taxonomy file, skip it
319 m->mothurOut(names[i] + " is not in your taxonomy file. I will not include it in the consensus."); m->mothurOutEndLine();
322 //if namesfile include the names
323 if (namefile != "") {
324 //is this sequence in the name file
325 it2 = nameMap.find(names[i]);
327 if (it2 == nameMap.end()) { //this name is not in name file, skip it
328 m->mothurOut(names[i] + " is not in your name file. I will not include it in the consensus."); m->mothurOutEndLine();
331 vector<string> nameFileNames;
332 splitAtComma(it2->second, nameFileNames);
334 for (int j = 0; j < nameFileNames.size(); j++) {
336 phylo->addSeqToTree(nameFileNames[j], it->second);
343 phylo->addSeqToTree(names[i], it->second);
350 phylo->assignHeirarchyIDs(0);
352 TaxNode currentNode = phylo->get(0);
355 while (currentNode.children.size() != 0) { //you still have more to explore
358 int bestChildSize = 0;
360 //go through children
361 for (map<string, int>::iterator itChild = currentNode.children.begin(); itChild != currentNode.children.end(); itChild++) {
363 TaxNode temp = phylo->get(itChild->second);
365 //select child with largest accesions - most seqs assigned to it
366 if (temp.accessions.size() > bestChildSize) {
367 bestChild = phylo->get(itChild->second);
368 bestChildSize = temp.accessions.size();
373 //is this taxonomy above cutoff
374 int consensusConfidence = ceil((bestChildSize / (float) size) * 100);
376 if (consensusConfidence >= cutoff) { //if yes, add it
378 conTax += bestChild.name + "(" + toString(consensusConfidence) + ");";
380 conTax += bestChild.name + ";";
387 currentNode = bestChild;
391 if (conTax == "") { conTax = "unclassified;"; }
398 catch(exception& e) {
399 m->errorOut(e, "ClassifyOtuCommand", "findConsensusTaxonomy");
404 //**********************************************************************************************************************
405 int ClassifyOtuCommand::process(ListVector* processList) {
411 if (outputDir == "") { outputDir += hasPath(listfile); }
414 string outputFile = outputDir + getRootName(getSimpleName(listfile)) + processList->getLabel() + ".cons.taxonomy";
415 openOutputFile(outputFile, out);
416 outputNames.push_back(outputFile);
418 //for each bin in the list vector
419 for (int i = 0; i < processList->getNumBins(); i++) {
420 conTax = findConsensusTaxonomy(i, processList, size);
422 if (m->control_pressed) { out.close(); return 0; }
424 //output to new names file
425 out << (i+1) << '\t' << size << '\t' << conTax << endl;
433 catch(exception& e) {
434 m->errorOut(e, "ClassifyOtuCommand", "process");
438 /**************************************************************************************************/
439 void ClassifyOtuCommand::removeConfidences(string& tax) {
445 while (tax.find_first_of(';') != -1) {
447 taxon = tax.substr(0,tax.find_first_of(';'));
449 int pos = taxon.find_first_of('(');
451 taxon = taxon.substr(0, pos); //rip off confidence
456 tax = tax.substr(tax.find_first_of(';')+1, tax.length());
462 catch(exception& e) {
463 m->errorOut(e, "ClassifyOtuCommand", "removeConfidences");
467 //**********************************************************************************************************************