1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449
|
//
// classifytreecommand.cpp
// Mothur
//
// Created by Sarah Westcott on 2/20/12.
// Copyright (c) 2012 Schloss Lab. All rights reserved.
//
#include "classifytreecommand.h"
#include "phylotree.h"
#include "treereader.h"
//**********************************************************************************************************************
vector<string> ClassifyTreeCommand::setParameters(){
try {
CommandParameter ptree("tree", "InputTypes", "", "", "", "", "none","tree-summary",false,true,true); parameters.push_back(ptree);
CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "", "", "none","",false,true,true); parameters.push_back(ptaxonomy);
CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","",false,false,true); parameters.push_back(pname);
CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","",false,false,true); parameters.push_back(pcount);
CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","",false,false,true); parameters.push_back(pgroup);
CommandParameter pmethod("output", "Multiple", "node-taxon", "node", "", "", "","",false,false); parameters.push_back(pmethod);
CommandParameter pcutoff("cutoff", "Number", "", "51", "", "", "","",false,true); parameters.push_back(pcutoff);
CommandParameter pseed("seed", "Number", "", "0", "", "", "","",false,false); parameters.push_back(pseed);
CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
abort = false; calledHelp = false;
vector<string> tempOutNames;
outputTypes["tree"] = tempOutNames;
outputTypes["summary"] = tempOutNames;
vector<string> myArray;
for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
return myArray;
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "setParameters");
exit(1);
}
}
//**********************************************************************************************************************
string ClassifyTreeCommand::getHelpString(){
try {
string helpString = "";
helpString += "The classify.tree command reads a tree and taxonomy file and output the consensus taxonomy for each node on the tree. \n";
helpString += "If you provide a group file, the concensus for each group will also be provided. \n";
helpString += "The new tree contains labels at each internal node. The label is the node number so you can relate the tree to the summary file.\n";
helpString += "The count parameter allows you add a count file so you can have the summary totals broken up by group.\n";
helpString += "The summary file lists the concensus taxonomy for the descendants of each node.\n";
helpString += "The classify.tree command parameters are tree, group, name, count and taxonomy. The tree and taxonomy files are required.\n";
helpString += "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";
helpString += "The output parameter allows you to specify whether you want the tree node number displayed on the tree, or the taxonomy displayed. Default=node. Options are node or taxon.\n";
helpString += "The classify.tree command should be used in the following format: classify.tree(tree=test.tre, group=test.group, taxonomy=test.taxonomy)\n";
return helpString;
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "getHelpString");
exit(1);
}
}
//**********************************************************************************************************************
string ClassifyTreeCommand::getOutputPattern(string type) {
try {
string pattern = "";
if (type == "summary") { pattern = "[filename],taxonomy.summary"; } //makes file like: amazon.0.03.fasta
else if (type == "tree") { pattern = "[filename],taxonomy.tre"; }
else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->setControl_pressed(true); }
return pattern;
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "getOutputPattern");
exit(1);
}
}
//**********************************************************************************************************************
ClassifyTreeCommand::ClassifyTreeCommand(string option) : Command() {
try {
//allow user to run help
if(option == "help") { help(); abort = true; calledHelp = true; }
else if(option == "citation") { citation(); abort = true; calledHelp = true;}
else if(option == "category") { abort = true; calledHelp = true; }
else {
OptionParser parser(option, setParameters());
map<string, string> parameters = parser.getParameters();
ValidParameters validParameter;
//check for required parameters
treefile = validParameter.validFile(parameters, "tree");
if (treefile == "not open") { treefile = ""; abort = true; }
else if (treefile == "not found") { treefile = "";
treefile = current->getTreeFile();
if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter.\n"); }
else { m->mothurOut("No valid current files. You must provide a tree file.\n"); abort = true; }
}else { current->setTreeFile(treefile); }
taxonomyfile = validParameter.validFile(parameters, "taxonomy");
if (taxonomyfile == "not open") { taxonomyfile = ""; abort = true; }
else if (taxonomyfile == "not found") { taxonomyfile = "";
taxonomyfile = current->getTaxonomyFile();
if (taxonomyfile != "") { m->mothurOut("Using " + taxonomyfile + " as input file for the taxonomy parameter.\n"); }
else { m->mothurOut("No valid current files. You must provide a taxonomy file.\n"); abort = true; }
}else { current->setTaxonomyFile(taxonomyfile); }
namefile = validParameter.validFile(parameters, "name");
if (namefile == "not open") { namefile = ""; abort = true; }
else if (namefile == "not found") { namefile = ""; }
else { current->setNameFile(namefile); }
groupfile = validParameter.validFile(parameters, "group");
if (groupfile == "not open") { groupfile = ""; abort = true; }
else if (groupfile == "not found") { groupfile = ""; }
else { current->setGroupFile(groupfile); }
countfile = validParameter.validFile(parameters, "count");
if (countfile == "not open") { countfile = ""; abort = true; }
else if (countfile == "not found") { countfile = ""; }
else { current->setCountFile(countfile); }
if ((namefile != "") && (countfile != "")) {
m->mothurOut("[ERROR]: you may only use one of the following: name or count.\n"); abort = true;
}
if ((groupfile != "") && (countfile != "")) {
m->mothurOut("[ERROR]: you may only use one of the following: group or count.\n"); abort=true;
}
string temp = validParameter.valid(parameters, "cutoff"); if (temp == "not found") { temp = "51"; }
util.mothurConvert(temp, cutoff);
if ((cutoff < 51) || (cutoff > 100)) { m->mothurOut("cutoff must be above 50, and no greater than 100.\n"); abort = true; }
output = validParameter.valid(parameters, "output");
if (output == "not found") { output = "node"; }
if ((output == "node") || (output == "taxon")) {
}else { m->mothurOut("[ERROR]: " + output + "is not a valid output option. Valid output options are node or taxon.\n"); abort = true; }
}
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "ClassifyTreeCommand");
exit(1);
}
}
//**********************************************************************************************************************
int ClassifyTreeCommand::execute(){
try {
if (abort) { if (calledHelp) { return 0; } return 2; }
cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
long start = time(nullptr);
/***************************************************/
// reading tree info //
/***************************************************/
current->setTreeFile(treefile);
TreeReader* reader = new TreeReader(treefile, groupfile, namefile);
vector<Tree*> T = reader->getTrees();
CountTable* tmap = T[0]->getCountTable();
Tree* outputTree = T[0];
delete reader;
if (namefile != "") { util.readNames(namefile, nameMap, nameCount); }
if (m->getControl_pressed()) { delete tmap; delete outputTree; return 0; }
util.readTax(taxonomyfile, taxMap, true);
/***************************************************/
// get concensus taxonomies //
/***************************************************/
getClassifications(outputTree);
delete outputTree; delete tmap;
if (m->getControl_pressed()) { for (int i = 0; i < outputNames.size(); i++) { util.mothurRemove(outputNames[i]); } return 0; }
//set tree file as new current treefile
if (treefile != "") {
string currentName = "";
itTypes = outputTypes.find("tree");
if (itTypes != outputTypes.end()) {
if ((itTypes->second).size() != 0) { currentName = (itTypes->second)[0]; current->setTreeFile(currentName); }
}
}
m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(nullptr) - start) + " secs to find the concensus taxonomies.\n");
m->mothurOut("\nOutput File Names: \n");
for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i] +"\n"); } m->mothurOutEndLine();
return 0;
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "execute");
exit(1);
}
}
//**********************************************************************************************************************
//traverse tree finding concensus taxonomy at each node
//label node with a number to relate to output summary file
//report all concensus taxonomies to file
int ClassifyTreeCommand::getClassifications(Tree*& T){
try {
string thisOutputDir = outputdir;
if (outputdir == "") { thisOutputDir += util.hasPath(treefile); }
map<string, string> variables;
variables["[filename]"] = thisOutputDir + util.getRootName(util.getSimpleName(treefile));
string outputFileName = getOutputFileName("summary", variables);
outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
ofstream out;
util.openOutputFile(outputFileName, out);
out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
//print headings
out << "TreeNode\t";
if (groupfile != "") { out << "Group\t"; }
out << "NumRep\tTaxonomy" << endl;
string treeOutputDir = outputdir;
if (outputdir == "") { treeOutputDir += util.hasPath(treefile); }
variables["[filename]"] = treeOutputDir + util.getRootName(util.getSimpleName(treefile));
string outputTreeFileName = getOutputFileName("tree", variables);
//create a map from tree node index to names of descendants, save time later
map<int, map<string, set<string> > > nodeToDescendants; //node# -> (groupName -> groupMembers)
for (int i = 0; i < T->getNumNodes(); i++) {
if (m->getControl_pressed()) { return 0; }
nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants);
}
//for each node
for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
if (m->getControl_pressed()) { out.close(); return 0; }
string tax = "not classifed";
int size;
if (groupfile != "") {
for (map<string, set<string> >::iterator itGroups = nodeToDescendants[i].begin(); itGroups != nodeToDescendants[i].end(); itGroups++) {
if (itGroups->first != "AllGroups") {
tax = getTaxonomy(itGroups->second, size);
out << (i+1) << '\t' << itGroups->first << '\t' << size << '\t' << tax << endl;
}
}
}else {
string group = "AllGroups";
tax = getTaxonomy(nodeToDescendants[i][group], size);
out << (i+1) << '\t' << size << '\t' << tax << endl;
}
if (output == "node") { T->tree[i].setLabel(toString(i+1)); }
else {
string cleanedTax = tax;
util.removeConfidences(cleanedTax);
for (int j = 0; j < cleanedTax.length(); j++) {
//special chars to trees - , ) ( ; [ ] :
if ((cleanedTax[j] == ',') || (cleanedTax[j] == '(') || (cleanedTax[j] == ')') || (cleanedTax[j] == ';') || (cleanedTax[j] == ':') || (cleanedTax[j] == ']') || (cleanedTax[j] == '[')) {
cleanedTax[j] = '_'; //change any special chars to _ so the tree can be read by tree readers
}
}
T->tree[i].setLabel(cleanedTax);
}
}
out.close();
ofstream outTree;
util.openOutputFile(outputTreeFileName, outTree);
outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
T->print(outTree, "both");
outTree.close();
return 0;
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "GetConcensusTaxonomies");
exit(1);
}
}
//**********************************************************************************************************************
string ClassifyTreeCommand::getTaxonomy(set<string> names, int& size) {
try{
string conTax = "";
size = 0;
//create a tree containing sequences from this bin
PhyloTree* phylo = new PhyloTree();
for (set<string>::iterator it = names.begin(); it != names.end(); it++) {
//if namesfile include the names
if (namefile != "") {
//is this sequence in the name file - namemap maps seqName -> repSeqName
map<string, string>::iterator it2 = nameMap.find(*it);
if (it2 == nameMap.end()) { //this name is not in name file, skip it
m->mothurOut((*it) + " is not in your name file. I will not include it in the consensus.\n");
}else{
//is this sequence in the taxonomy file - look for repSeqName since we are assuming the taxonomy file is unique
map<string, string>::iterator itTax = taxMap.find((it2->second));
if (itTax == taxMap.end()) { //this name is not in taxonomy file, skip it
if ((*it) != (it2->second)) { m->mothurOut((*it) + " is represented by " + it2->second + " and is not in your taxonomy file. I will not include it in the consensus.\n"); }
else { m->mothurOut((*it) + " is not in your taxonomy file. I will not include it in the consensus.\n"); }
}else{
//add seq to tree
int num = nameCount[(*it)]; // we know its there since we found it in nameMap
for (int i = 0; i < num; i++) { phylo->addSeqToTree((*it)+toString(i), itTax->second); }
size += num;
}
}
}else{
//is this sequence in the taxonomy file - look for repSeqName since we are assuming the taxonomy file is unique
map<string, string>::iterator itTax = taxMap.find((*it));
if (itTax == taxMap.end()) { //this name is not in taxonomy file, skip it
m->mothurOut((*it) + " is not in your taxonomy file. I will not include it in the consensus.\n");
}else{
if (countfile != "") {
int numDups = ct->getNumSeqs((*it));
for (int j = 0; j < numDups; j++) { phylo->addSeqToTree((*it), itTax->second); }
size += numDups;
}else{
//add seq to tree
phylo->addSeqToTree((*it), itTax->second);
size++;
} }
}
if (m->getControl_pressed()) { delete phylo; return conTax; }
}
//build tree
phylo->assignHeirarchyIDs(0);
TaxNode currentNode = phylo->get(0);
int myLevel = 0;
//at each level
while (currentNode.children.size() != 0) { //you still have more to explore
TaxNode bestChild;
int bestChildSize = 0;
//go through children
for (map<string, int>::iterator itChild = currentNode.children.begin(); itChild != currentNode.children.end(); itChild++) {
TaxNode temp = phylo->get(itChild->second);
//select child with largest accesions - most seqs assigned to it
if (temp.accessions.size() > bestChildSize) {
bestChild = phylo->get(itChild->second);
bestChildSize = temp.accessions.size();
}
}
//is this taxonomy above cutoff
int consensusConfidence = ceil((bestChildSize / (float) size) * 100);
if (consensusConfidence >= cutoff) { //if yes, add it
conTax += bestChild.name + "(" + toString(consensusConfidence) + ");";
myLevel++;
}else{ //if no, quit
break;
}
//move down a level
currentNode = bestChild;
}
if (myLevel != phylo->getMaxLevel()) {
while (myLevel != phylo->getMaxLevel()) {
conTax += "unclassified;";
myLevel++;
}
}
if (conTax == "") { conTax = "no_consensus;"; }
delete phylo;
return conTax;
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "getTaxonomy");
exit(1);
}
}
//**********************************************************************************************************************
map<string, set<string> > ClassifyTreeCommand::getDescendantList(Tree*& T, int i, map<int, map<string, set<string> > > descendants){
try {
map<string ,set<string> > names;
map<string ,set<string> >::iterator it;
map<string ,set<string> >::iterator it2;
int lc = T->tree[i].getLChild();
int rc = T->tree[i].getRChild();
// TreeMap* tmap = T->getTreeMap();
if (lc == -1) { //you are a leaf your only descendant is yourself
vector<string> groups = T->tree[i].getGroup();
set<string> mynames; mynames.insert(T->tree[i].getName());
for (int j = 0; j < groups.size(); j++) { names[groups[j]] = mynames; } //mygroup -> me
names["AllGroups"] = mynames;
}else{ //your descedants are the combination of your childrens descendants
names = descendants[lc];
for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
it2 = names.find(it->first); //do we already have this group
if (it2 == names.end()) { //nope, so add it
names[it->first] = it->second;
}else {
for (set<string>::iterator it3 = (it->second).begin(); it3 != (it->second).end(); it3++) {
names[it->first].insert(*it3);
}
}
}
}
return names;
}
catch(exception& e) {
m->errorOut(e, "ClassifyTreeCommand", "getDescendantList");
exit(1);
}
}
/*****************************************************************/
|