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 = m->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 = m->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 = m->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") { m->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 = m->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 ((m->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 m->openInputFile(namefile, inNames);
249 inNames >> name; //read from first column A
250 inNames >> names; //read from second column A,B,C,D
253 //parse names into vector
254 vector<string> theseNames;
255 m->splitAtComma(names, theseNames);
257 for (int i = 0; i < theseNames.size(); i++) { nameMap[theseNames[i]] = name; }
259 if (m->control_pressed) { inNames.close(); nameMap.clear(); return 0; }
265 catch(exception& e) {
266 m->errorOut(e, "ClassifyOtuCommand", "readNamesFile");
270 //**********************************************************************************************************************
271 int ClassifyOtuCommand::readTaxonomyFile() {
275 m->openInputFile(taxfile, in);
283 //are there confidence scores, if so remove them
284 if (tax.find_first_of('(') != -1) { removeConfidences(tax); }
288 if (m->control_pressed) { in.close(); taxMap.clear(); return 0; }
294 catch(exception& e) {
295 m->errorOut(e, "ClassifyOtuCommand", "readTaxonomyFile");
299 //**********************************************************************************************************************
300 string ClassifyOtuCommand::findConsensusTaxonomy(int bin, ListVector* thisList, int& size) {
303 vector<string> names;
304 map<string, string>::iterator it;
305 map<string, string>::iterator it2;
307 //parse names into vector
308 string binnames = thisList->get(bin);
309 m->splitAtComma(binnames, names);
311 //create a tree containing sequences from this bin
312 PhyloTree* phylo = new PhyloTree();
315 for (int i = 0; i < names.size(); i++) {
317 //if namesfile include the names
318 if (namefile != "") {
319 //is this sequence in the name file - namemap maps seqName -> repSeqName
320 it2 = nameMap.find(names[i]);
322 if (it2 == nameMap.end()) { //this name is not in name file, skip it
323 m->mothurOut(names[i] + " is not in your name file. I will not include it in the consensus."); m->mothurOutEndLine();
326 //is this sequence in the taxonomy file - look for repSeqName since we are assuming the taxonomy file is unique
327 it = taxMap.find(it2->second);
329 if (it == taxMap.end()) { //this name is not in taxonomy file, skip it
330 if (names[i] != it->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(); }
331 else { m->mothurOut(names[i] + " is not in your taxonomy file. I will not include it in the consensus."); m->mothurOutEndLine(); }
335 phylo->addSeqToTree(names[i], it->second);
341 //is this sequence in the taxonomy file - look for repSeqName since we are assuming the taxonomy file is unique
342 it = taxMap.find(names[i]);
344 if (it == taxMap.end()) { //this name is not in taxonomy file, skip it
345 m->mothurOut(names[i] + " is not in your taxonomy file. I will not include it in the consensus."); m->mothurOutEndLine();
348 phylo->addSeqToTree(names[i], it->second);
354 if (m->control_pressed) { delete phylo; return conTax; }
359 phylo->assignHeirarchyIDs(0);
361 TaxNode currentNode = phylo->get(0);
364 while (currentNode.children.size() != 0) { //you still have more to explore
367 int bestChildSize = 0;
369 //go through children
370 for (map<string, int>::iterator itChild = currentNode.children.begin(); itChild != currentNode.children.end(); itChild++) {
372 TaxNode temp = phylo->get(itChild->second);
374 //select child with largest accesions - most seqs assigned to it
375 if (temp.accessions.size() > bestChildSize) {
376 bestChild = phylo->get(itChild->second);
377 bestChildSize = temp.accessions.size();
382 //is this taxonomy above cutoff
383 int consensusConfidence = ceil((bestChildSize / (float) size) * 100);
385 if (consensusConfidence >= cutoff) { //if yes, add it
387 conTax += bestChild.name + "(" + toString(consensusConfidence) + ");";
389 conTax += bestChild.name + ";";
396 currentNode = bestChild;
400 if (conTax == "") { conTax = "unclassified;"; }
407 catch(exception& e) {
408 m->errorOut(e, "ClassifyOtuCommand", "findConsensusTaxonomy");
413 //**********************************************************************************************************************
414 int ClassifyOtuCommand::process(ListVector* processList) {
420 if (outputDir == "") { outputDir += m->hasPath(listfile); }
423 string outputFile = outputDir + m->getRootName(m->getSimpleName(listfile)) + processList->getLabel() + ".cons.taxonomy";
424 m->openOutputFile(outputFile, out);
425 outputNames.push_back(outputFile);
427 //for each bin in the list vector
428 for (int i = 0; i < processList->getNumBins(); i++) {
429 conTax = findConsensusTaxonomy(i, processList, size);
431 if (m->control_pressed) { out.close(); return 0; }
433 //output to new names file
434 out << (i+1) << '\t' << size << '\t' << conTax << endl;
442 catch(exception& e) {
443 m->errorOut(e, "ClassifyOtuCommand", "process");
447 /**************************************************************************************************/
448 void ClassifyOtuCommand::removeConfidences(string& tax) {
454 while (tax.find_first_of(';') != -1) {
456 taxon = tax.substr(0,tax.find_first_of(';'));
458 int pos = taxon.find_first_of('(');
460 taxon = taxon.substr(0, pos); //rip off confidence
465 tax = tax.substr(tax.find_first_of(';')+1, tax.length());
471 catch(exception& e) {
472 m->errorOut(e, "ClassifyOtuCommand", "removeConfidences");
476 //**********************************************************************************************************************