classifyotucommand.cpp

   1 /*
   2  *  classifyotucommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 6/1/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "classifyotucommand.h"
  11 #include "phylotree.h"
  12
  13
  14 //**********************************************************************************************************************
  15 ClassifyOtuCommand::ClassifyOtuCommand(string option)  {
  16         try{
  17                 abort = false;
  18                 allLines = 1;
  19                 labels.clear();
  20
  21                 //allow user to run help
  22                 if (option == "help") {
  23                         help(); abort = true;
  24                 } else {
  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)));
  28
  29                         OptionParser parser(option);
  30                         map<string, string> parameters = parser.getParameters();
  31
  32                         ValidParameters validParameter;
  33                         map<string, string>::iterator it;
  34
  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;  }
  38                         }
  39
  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 = "";          }
  43                         else {
  44                                 string path;
  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;             }
  51                                 }
  52
  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;             }
  59                                 }
  60
  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;         }
  67                                 }
  68                         }
  69
  70
  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 = "";         }
  73
  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; }
  78
  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; }
  82
  83                         namefile = validParameter.validFile(parameters, "name", true);
  84                         if (namefile == "not open") { abort = true; }
  85                         else if (namefile == "not found") { namefile = ""; }
  86
  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;  }
  91                         else {
  92                                 if(label != "all") {  splitAtDash(label, labels);  allLines = 0;  }
  93                                 else { allLines = 1;  }
  94                         }
  95
  96                         string temp = validParameter.validFile(parameters, "cutoff", false);                    if (temp == "not found") { temp = "51"; }
  97                         convert(temp, cutoff);
  98
  99                         temp = validParameter.validFile(parameters, "probs", false);                                    if (temp == "not found"){       temp = "true";                  }
 100                         probs = isTrue(temp);
 101
 102
 103                         if ((cutoff < 51) || (cutoff > 100)) { m->mothurOut("cutoff must be above 50, and no greater than 100."); m->mothurOutEndLine(); abort = true;  }
 104
 105                 }
 106         }
 107         catch(exception& e) {
 108                 m->errorOut(e, "ClassifyOtuCommand", "ClassifyOtuCommand");
 109                 exit(1);
 110         }
 111 }
 112
 113 //**********************************************************************************************************************
 114
 115 void ClassifyOtuCommand::help(){
 116         try {
 117                 m->mothurOut("The classify.otu command parameters are list, taxonomy, name, cutoff, and label.  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 concensus 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 concensus 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");
 126         }
 127         catch(exception& e) {
 128                 m->errorOut(e, "ClassifyOtuCommand", "help");
 129                 exit(1);
 130         }
 131 }
 132
 133 //**********************************************************************************************************************
 134
 135 ClassifyOtuCommand::~ClassifyOtuCommand(){}
 136
 137 //**********************************************************************************************************************
 138
 139 int ClassifyOtuCommand::execute(){
 140         try {
 141
 142                 if (abort == true) { return 0; }
 143
 144                 //if user gave a namesfile then use it
 145                 if (namefile != "") {   readNamesFile();        }
 146
 147                 //read taxonomy file and save in map for easy access in building bin trees
 148                 readTaxonomyFile();
 149
 150                 if (m->control_pressed) { return 0; }
 151
 152                 input = new InputData(listfile, "list");
 153                 list = input->getListVector();
 154                 string lastLabel = list->getLabel();
 155
 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;
 159
 160                 if (m->control_pressed) { delete input; delete list; for (int i = 0; i < outputNames.size(); i++) {     remove(outputNames[i].c_str());  }  return 0; }
 161
 162                 while((list != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 163
 164                         if (allLines == 1 || labels.count(list->getLabel()) == 1){
 165
 166                                         m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
 167                                         process(list);
 168                                         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } delete input; delete list; return 0; }
 169
 170                                         processedLabels.insert(list->getLabel());
 171                                         userLabels.erase(list->getLabel());
 172                         }
 173
 174                         if ((anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 175                                         string saveLabel = list->getLabel();
 176
 177                                         delete list;
 178                                         list = input->getListVector(lastLabel);
 179                                         m->mothurOut(list->getLabel() + "\t" + toString(list->size())); m->mothurOutEndLine();
 180                                         process(list);
 181
 182
 183                                         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } delete input; delete list; return 0; }
 184
 185                                         processedLabels.insert(list->getLabel());
 186                                         userLabels.erase(list->getLabel());
 187
 188                                         //restore real lastlabel to save below
 189                                         list->setLabel(saveLabel);
 190                         }
 191
 192                         lastLabel = list->getLabel();
 193
 194                         delete list;
 195                         list = input->getListVector();
 196                 }
 197
 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();
 204                                 needToRun = true;
 205                         }else {
 206                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 207                         }
 208                 }
 209
 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();
 215
 216                         process(list);
 217                         delete list;
 218
 219                         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } delete input; delete list; return 0; }
 220                 }
 221
 222                 delete input;
 223
 224                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str());  } return 0; }
 225
 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();
 230
 231                 return 0;
 232         }
 233         catch(exception& e) {
 234                 m->errorOut(e, "ClassifyOtuCommand", "execute");
 235                 exit(1);
 236         }
 237 }
 238
 239 //**********************************************************************************************************************
 240 int ClassifyOtuCommand::readNamesFile() {
 241         try {
 242
 243                 ifstream inNames;
 244                 openInputFile(namefile, inNames);
 245
 246                 string name, names;
 247
 248                 while(inNames){
 249                         inNames >> name;                        //read from first column  A
 250                         inNames >> names;               //read from second column  A,B,C,D
 251                         gobble(inNames);
 252
 253                         nameMap[name] = names;
 254
 255                         if (m->control_pressed) { inNames.close(); nameMap.clear(); return 0; }
 256                 }
 257                 inNames.close();
 258
 259                 return 0;
 260         }
 261         catch(exception& e) {
 262                 m->errorOut(e, "ClassifyOtuCommand", "readNamesFile");
 263                 exit(1);
 264         }
 265 }
 266 //**********************************************************************************************************************
 267 int ClassifyOtuCommand::readTaxonomyFile() {
 268         try {
 269
 270                 ifstream in;
 271                 openInputFile(taxfile, in);
 272
 273                 string name, tax;
 274
 275                 while(!in.eof()){
 276                         in >> name >> tax;
 277                         gobble(in);
 278
 279                         //are there confidence scores, if so remove them
 280                         if (tax.find_first_of('(') != -1) {  removeConfidences(tax);    }
 281
 282                         taxMap[name] = tax;
 283
 284                         if (m->control_pressed) { in.close(); taxMap.clear(); return 0; }
 285                 }
 286                 in.close();
 287
 288                 return 0;
 289         }
 290         catch(exception& e) {
 291                 m->errorOut(e, "ClassifyOtuCommand", "readTaxonomyFile");
 292                 exit(1);
 293         }
 294 }
 295 //**********************************************************************************************************************
 296 string ClassifyOtuCommand::findConcensusTaxonomy(int bin, ListVector* thisList, int& size) {
 297         try{
 298                 string conTax = "";
 299                 vector<string> names;
 300                 map<string, string>::iterator it;
 301                 map<string, string>::iterator it2;
 302
 303                 //parse names into vector
 304                 string binnames = thisList->get(bin);
 305                 splitAtComma(binnames, names);
 306
 307                 //create a tree containing sequences from this bin
 308                 PhyloTree* phylo = new PhyloTree();
 309
 310                 size = 0;
 311                 for (int i = 0; i < names.size(); i++) {
 312
 313                         if (m->control_pressed) { delete phylo; return conTax; }
 314
 315                         //is this sequence in the taxonomy file
 316                         it = taxMap.find(names[i]);
 317
 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 concensus."); m->mothurOutEndLine();
 320                         }else{
 321
 322                                 //if namesfile include the names
 323                                 if (namefile != "") {
 324                                         //is this sequence in the name file
 325                                         it2 = nameMap.find(names[i]);
 326
 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 concensus."); m->mothurOutEndLine();
 329                                         }else{
 330
 331                                                 vector<string> nameFileNames;
 332                                                 splitAtComma(it2->second, nameFileNames);
 333
 334                                                 for (int j = 0; j < nameFileNames.size(); j++) {
 335                                                         //add seq to tree
 336                                                         phylo->addSeqToTree(nameFileNames[j], it->second);
 337                                                         size++;
 338                                                 }
 339                                         }
 340
 341                                 }else{
 342                                         //add seq to tree
 343                                         phylo->addSeqToTree(names[i], it->second);
 344                                         size++;
 345                                 }
 346                         }
 347                 }
 348
 349                 //build tree
 350                 phylo->assignHeirarchyIDs(0);
 351
 352                 TaxNode currentNode = phylo->get(0);
 353
 354                 //at each level
 355                 while (currentNode.children.size() != 0) { //you still have more to explore
 356
 357                         TaxNode bestChild;
 358                         int bestChildSize = 0;
 359
 360                         //go through children
 361                         for (map<string, int>::iterator itChild = currentNode.children.begin(); itChild != currentNode.children.end(); itChild++) {
 362
 363                                 TaxNode temp = phylo->get(itChild->second);
 364
 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();
 369                                 }
 370
 371                         }
 372
 373                         //is this taxonomy above cutoff
 374                         int concensusConfidence = ceil((bestChildSize / (float) size) * 100);
 375
 376                         if (concensusConfidence >= cutoff) { //if yes, add it
 377                                 if (probs) {
 378                                         conTax += bestChild.name + "(" + toString(concensusConfidence) + ");";
 379                                 }else{
 380                                         conTax += bestChild.name + ";";
 381                                 }
 382                         }else{ //if no, quit
 383                                 break;
 384                         }
 385
 386                         //move down a level
 387                         currentNode = bestChild;
 388                 }
 389
 390
 391                 if (conTax == "") {  conTax = "unclassified;";  }
 392
 393                 delete phylo;
 394
 395                 return conTax;
 396
 397         }
 398         catch(exception& e) {
 399                 m->errorOut(e, "ClassifyOtuCommand", "findConcensusTaxonomy");
 400                 exit(1);
 401         }
 402 }
 403
 404 //**********************************************************************************************************************
 405 int ClassifyOtuCommand::process(ListVector* processList) {
 406         try{
 407                 string conTax;
 408                 int size;
 409
 410                 //create output file
 411                 if (outputDir == "") { outputDir += hasPath(listfile); }
 412
 413                 ofstream out;
 414                 string outputFile = outputDir + getRootName(getSimpleName(listfile)) + processList->getLabel() + ".cons.taxonomy";
 415                 openOutputFile(outputFile, out);
 416                 outputNames.push_back(outputFile);
 417
 418                 //for each bin in the list vector
 419                 for (int i = 0; i < processList->getNumBins(); i++) {
 420                         conTax  = findConcensusTaxonomy(i, processList, size);
 421
 422                         if (m->control_pressed) { out.close();  return 0; }
 423
 424                         //output to new names file
 425                         out << (i+1) << '\t' << size << '\t' << conTax << endl;
 426                 }
 427
 428                 out.close();
 429
 430                 return 0;
 431
 432         }
 433         catch(exception& e) {
 434                 m->errorOut(e, "ClassifyOtuCommand", "process");
 435                 exit(1);
 436         }
 437 }
 438 /**************************************************************************************************/
 439 void ClassifyOtuCommand::removeConfidences(string& tax) {
 440         try {
 441
 442                 string taxon;
 443                 string newTax = "";
 444
 445                 while (tax.find_first_of(';') != -1) {
 446                         //get taxon
 447                         taxon = tax.substr(0,tax.find_first_of(';'));
 448
 449                         int pos = taxon.find_first_of('(');
 450                         if (pos != -1) {
 451                                 taxon = taxon.substr(0, pos); //rip off confidence
 452                         }
 453
 454                         taxon += ";";
 455
 456                         tax = tax.substr(tax.find_first_of(';')+1, tax.length());
 457                         newTax += taxon;
 458                 }
 459
 460                 tax = newTax;
 461         }
 462         catch(exception& e) {
 463                 m->errorOut(e, "ClassifyOtuCommand", "removeConfidences");
 464                 exit(1);
 465         }
 466 }
 467 //**********************************************************************************************************************
 468
 469