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, 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");
 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                         //parse names into vector
 254                         vector<string> theseNames;
 255                         splitAtComma(names, theseNames);
 256
 257                         for (int i = 0; i < theseNames.size(); i++) {  nameMap[theseNames[i]] = name;  }
 258
 259                         if (m->control_pressed) { inNames.close(); nameMap.clear(); return 0; }
 260                 }
 261                 inNames.close();
 262
 263                 return 0;
 264         }
 265         catch(exception& e) {
 266                 m->errorOut(e, "ClassifyOtuCommand", "readNamesFile");
 267                 exit(1);
 268         }
 269 }
 270 //**********************************************************************************************************************
 271 int ClassifyOtuCommand::readTaxonomyFile() {
 272         try {
 273
 274                 ifstream in;
 275                 openInputFile(taxfile, in);
 276
 277                 string name, tax;
 278
 279                 while(!in.eof()){
 280                         in >> name >> tax;
 281                         gobble(in);
 282
 283                         //are there confidence scores, if so remove them
 284                         if (tax.find_first_of('(') != -1) {  removeConfidences(tax);    }
 285
 286                         taxMap[name] = tax;
 287
 288                         if (m->control_pressed) { in.close(); taxMap.clear(); return 0; }
 289                 }
 290                 in.close();
 291
 292                 return 0;
 293         }
 294         catch(exception& e) {
 295                 m->errorOut(e, "ClassifyOtuCommand", "readTaxonomyFile");
 296                 exit(1);
 297         }
 298 }
 299 //**********************************************************************************************************************
 300 string ClassifyOtuCommand::findConsensusTaxonomy(int bin, ListVector* thisList, int& size) {
 301         try{
 302                 string conTax = "";
 303                 vector<string> names;
 304                 map<string, string>::iterator it;
 305                 map<string, string>::iterator it2;
 306
 307                 //parse names into vector
 308                 string binnames = thisList->get(bin);
 309                 splitAtComma(binnames, names);
 310
 311                 //create a tree containing sequences from this bin
 312                 PhyloTree* phylo = new PhyloTree();
 313
 314                 size = 0;
 315                 for (int i = 0; i < names.size(); i++) {
 316
 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]);
 321
 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();
 324                                 }else{
 325
 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);
 328
 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(); }
 332                                         }else{
 333
 334                                                 //add seq to tree
 335                                                 phylo->addSeqToTree(names[i], it->second);
 336                                                 size++;
 337                                         }
 338                                 }
 339
 340                         }else{
 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]);
 343
 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();
 346                                 }else{
 347                                         //add seq to tree
 348                                         phylo->addSeqToTree(names[i], it->second);
 349                                         size++;
 350                                 }
 351                         }
 352
 353
 354                         if (m->control_pressed) { delete phylo; return conTax; }
 355
 356                 }
 357
 358                 //build tree
 359                 phylo->assignHeirarchyIDs(0);
 360
 361                 TaxNode currentNode = phylo->get(0);
 362
 363                 //at each level
 364                 while (currentNode.children.size() != 0) { //you still have more to explore
 365
 366                         TaxNode bestChild;
 367                         int bestChildSize = 0;
 368
 369                         //go through children
 370                         for (map<string, int>::iterator itChild = currentNode.children.begin(); itChild != currentNode.children.end(); itChild++) {
 371
 372                                 TaxNode temp = phylo->get(itChild->second);
 373
 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();
 378                                 }
 379
 380                         }
 381
 382                         //is this taxonomy above cutoff
 383                         int consensusConfidence = ceil((bestChildSize / (float) size) * 100);
 384
 385                         if (consensusConfidence >= cutoff) { //if yes, add it
 386                                 if (probs) {
 387                                         conTax += bestChild.name + "(" + toString(consensusConfidence) + ");";
 388                                 }else{
 389                                         conTax += bestChild.name + ";";
 390                                 }
 391                         }else{ //if no, quit
 392                                 break;
 393                         }
 394
 395                         //move down a level
 396                         currentNode = bestChild;
 397                 }
 398
 399
 400                 if (conTax == "") {  conTax = "unclassified;";  }
 401
 402                 delete phylo;
 403
 404                 return conTax;
 405
 406         }
 407         catch(exception& e) {
 408                 m->errorOut(e, "ClassifyOtuCommand", "findConsensusTaxonomy");
 409                 exit(1);
 410         }
 411 }
 412
 413 //**********************************************************************************************************************
 414 int ClassifyOtuCommand::process(ListVector* processList) {
 415         try{
 416                 string conTax;
 417                 int size;
 418
 419                 //create output file
 420                 if (outputDir == "") { outputDir += hasPath(listfile); }
 421
 422                 ofstream out;
 423                 string outputFile = outputDir + getRootName(getSimpleName(listfile)) + processList->getLabel() + ".cons.taxonomy";
 424                 openOutputFile(outputFile, out);
 425                 outputNames.push_back(outputFile);
 426
 427                 //for each bin in the list vector
 428                 for (int i = 0; i < processList->getNumBins(); i++) {
 429                         conTax  = findConsensusTaxonomy(i, processList, size);
 430
 431                         if (m->control_pressed) { out.close();  return 0; }
 432
 433                         //output to new names file
 434                         out << (i+1) << '\t' << size << '\t' << conTax << endl;
 435                 }
 436
 437                 out.close();
 438
 439                 return 0;
 440
 441         }
 442         catch(exception& e) {
 443                 m->errorOut(e, "ClassifyOtuCommand", "process");
 444                 exit(1);
 445         }
 446 }
 447 /**************************************************************************************************/
 448 void ClassifyOtuCommand::removeConfidences(string& tax) {
 449         try {
 450
 451                 string taxon;
 452                 string newTax = "";
 453
 454                 while (tax.find_first_of(';') != -1) {
 455                         //get taxon
 456                         taxon = tax.substr(0,tax.find_first_of(';'));
 457
 458                         int pos = taxon.find_first_of('(');
 459                         if (pos != -1) {
 460                                 taxon = taxon.substr(0, pos); //rip off confidence
 461                         }
 462
 463                         taxon += ";";
 464
 465                         tax = tax.substr(tax.find_first_of(';')+1, tax.length());
 466                         newTax += taxon;
 467                 }
 468
 469                 tax = newTax;
 470         }
 471         catch(exception& e) {
 472                 m->errorOut(e, "ClassifyOtuCommand", "removeConfidences");
 473                 exit(1);
 474         }
 475 }
 476 //**********************************************************************************************************************
 477
 478