deuniquetreecommand.cpp

   1 /*
   2  *  deuniquetreecommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 5/27/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "deuniquetreecommand.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> DeuniqueTreeCommand::setParameters(){
  14         try {
  15                 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
  16                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
  17                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  18                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  19
  20                 vector<string> myArray;
  21                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  22                 return myArray;
  23         }
  24         catch(exception& e) {
  25                 m->errorOut(e, "DeuniqueTreeCommand", "setParameters");
  26                 exit(1);
  27         }
  28 }
  29 //**********************************************************************************************************************
  30 string DeuniqueTreeCommand::getHelpString(){
  31         try {
  32                 string helpString = "";
  33                 helpString += "The deunique.tree command parameters are tree and name.  Both parameters are required unless you have valid current files.\n";
  34                 helpString += "The deunique.tree command should be in the following format: deunique.tree(tree=yourTreeFile, name=yourNameFile).\n";
  35                 helpString += "Example deunique.tree(tree=abrecovery.tree, name=abrecovery.names).\n";
  36                 helpString += "Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreeFile).\n";
  37                 return helpString;
  38         }
  39         catch(exception& e) {
  40                 m->errorOut(e, "DeuniqueTreeCommand", "getHelpString");
  41                 exit(1);
  42         }
  43 }
  44 //**********************************************************************************************************************
  45 DeuniqueTreeCommand::DeuniqueTreeCommand(){
  46         try {
  47                 abort = true; calledHelp = true;
  48                 setParameters();
  49                 vector<string> tempOutNames;
  50                 outputTypes["tree"] = tempOutNames;
  51         }
  52         catch(exception& e) {
  53                 m->errorOut(e, "DeuniqueTreeCommand", "DeuniqueTreeCommand");
  54                 exit(1);
  55         }
  56 }
  57 /***********************************************************/
  58 DeuniqueTreeCommand::DeuniqueTreeCommand(string option)  {
  59         try {
  60                 abort = false; calledHelp = false;
  61
  62                 //allow user to run help
  63                 if(option == "help") { help(); abort = true; calledHelp = true; }
  64                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  65
  66                 else {
  67                         vector<string> myArray = setParameters();
  68
  69                         OptionParser parser(option);
  70                         map<string,string> parameters = parser.getParameters();
  71                         map<string,string>::iterator it;
  72
  73                         ValidParameters validParameter;
  74
  75                         //check to make sure all parameters are valid for command
  76                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  77                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  78                         }
  79
  80                         //initialize outputTypes
  81                         vector<string> tempOutNames;
  82                         outputTypes["tree"] = tempOutNames;
  83
  84                         //if the user changes the input directory command factory will send this info to us in the output parameter
  85                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  86                         if (inputDir == "not found"){   inputDir = "";          }
  87                         else {
  88                                 string path;
  89                                 it = parameters.find("tree");
  90                                 //user has given a template file
  91                                 if(it != parameters.end()){
  92                                         path = m->hasPath(it->second);
  93                                         //if the user has not given a path then, add inputdir. else leave path alone.
  94                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
  95                                 }
  96
  97                                 it = parameters.find("name");
  98                                 //user has given a template file
  99                                 if(it != parameters.end()){
 100                                         path = m->hasPath(it->second);
 101                                         //if the user has not given a path then, add inputdir. else leave path alone.
 102                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 103                                 }
 104                         }
 105
 106                         m->runParse = true;
 107                         m->clearGroups();
 108                         m->clearAllGroups();
 109                         m->Treenames.clear();
 110                         m->names.clear();
 111
 112                         //check for required parameters
 113                         treefile = validParameter.validFile(parameters, "tree", true);
 114                         if (treefile == "not open") { abort = true; }
 115                         else if (treefile == "not found") {                             //if there is a current design file, use it
 116                                 treefile = m->getTreeFile();
 117                                 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
 118                                 else {  m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
 119                         }else { m->setTreeFile(treefile); }
 120
 121                         namefile = validParameter.validFile(parameters, "name", true);
 122                         if (namefile == "not open") { abort = true; }
 123                         else if (namefile == "not found") {                             //if there is a current design file, use it
 124                                 namefile = m->getNameFile();
 125                                 if (namefile != "") { m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
 126                                 else {  m->mothurOut("You have no current name file and the name parameter is required."); m->mothurOutEndLine(); abort = true; }
 127                         }else { m->setNameFile(namefile); }
 128
 129                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(treefile);       }
 130                 }
 131
 132         }
 133         catch(exception& e) {
 134                 m->errorOut(e, "DeuniqueTreeCommand", "DeuniqueTreeCommand");
 135                 exit(1);
 136         }
 137 }
 138
 139 /***********************************************************/
 140 int DeuniqueTreeCommand::execute() {
 141         try {
 142
 143                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 144
 145                 m->setTreeFile(treefile);
 146
 147                 //extracts names from tree to make faked out groupmap
 148                 Tree* tree = new Tree(treefile); delete tree;
 149                 tmap = new TreeMap();
 150                 for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
 151
 152                 if (m->control_pressed) {  delete tmap;  return 0; }
 153
 154                 readNamesFile();
 155
 156                 if (m->control_pressed) {  delete tmap;  return 0; }
 157
 158                 ReadTree* read = new ReadNewickTree(treefile);
 159                 int readOk = read->read(tmap);
 160                 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
 161
 162                 read->AssembleTrees();
 163                 vector<Tree*> T = read->getTrees();
 164                 delete read;
 165
 166                 //make sure all files match
 167                 //if you provide a namefile we will use the numNames in the namefile as long as the number of unique match the tree names size.
 168                 int numNamesInTree;
 169                 if (numUniquesInName == m->Treenames.size()) {  numNamesInTree = nameMap.size();  }
 170                 else {   numNamesInTree = m->Treenames.size();  }
 171
 172                 //output any names that are in group file but not in tree
 173                 if (numNamesInTree < tmap->getNumSeqs()) {
 174                         for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
 175                                 //is that name in the tree?
 176                                 int count = 0;
 177                                 for (int j = 0; j < m->Treenames.size(); j++) {
 178                                         if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
 179                                         count++;
 180                                 }
 181
 182                                 if (m->control_pressed) {
 183                                         delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 184                                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } outputTypes.clear();
 185                                         m->clearGroups();
 186                                         return 0;
 187                                 }
 188
 189                                 //then you did not find it so report it
 190                                 if (count == m->Treenames.size()) {
 191                                         //if it is in your namefile then don't remove
 192                                         map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
 193
 194                                         if (it == nameMap.end()) {
 195                                                 m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
 196                                                 tmap->removeSeq(tmap->namesOfSeqs[i]);
 197                                                 i--; //need this because removeSeq removes name from namesOfSeqs
 198                                         }
 199                                 }
 200                         }
 201                 }
 202
 203
 204                 //print new Tree
 205                 string outputFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + "deunique.tre";
 206                 outputNames.push_back(outputFile); outputTypes["tree"].push_back(outputFile);
 207                 ofstream out;
 208                 m->openOutputFile(outputFile, out);
 209                 T[0]->print(out, "deunique");
 210                 out.close();
 211
 212                 delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
 213
 214                 //set phylip file as new current phylipfile
 215                 string current = "";
 216                 itTypes = outputTypes.find("tree");
 217                 if (itTypes != outputTypes.end()) {
 218                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
 219                 }
 220
 221                 m->mothurOutEndLine();
 222                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 223                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 224                 m->mothurOutEndLine();
 225
 226                 return 0;
 227
 228         }
 229         catch(exception& e) {
 230                 m->errorOut(e, "DeuniqueTreeCommand", "execute");
 231                 exit(1);
 232         }
 233 }
 234 /*****************************************************************/
 235 int DeuniqueTreeCommand::readNamesFile() {
 236         try {
 237                 m->names.clear();
 238                 numUniquesInName = 0;
 239
 240                 ifstream in;
 241                 m->openInputFile(namefile, in);
 242
 243                 string first, second;
 244                 map<string, string>::iterator itNames;
 245
 246                 while(!in.eof()) {
 247                         in >> first >> second; m->gobble(in);
 248
 249                         numUniquesInName++;
 250
 251                         itNames = m->names.find(first);
 252                         if (itNames == m->names.end()) {
 253                                 m->names[first] = second;
 254
 255                                 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
 256                                 vector<string> dupNames;
 257                                 m->splitAtComma(second, dupNames);
 258
 259                                 for (int i = 0; i < dupNames.size(); i++) {
 260                                         nameMap[dupNames[i]] = dupNames[i];
 261                                         if (i != 0) { tmap->addSeq(dupNames[i], "Group1"); }
 262                                 }
 263                         }else {  m->mothurOut(first + " has already been seen in namefile, aborting."); m->mothurOutEndLine(); in.close(); m->names.clear(); m->control_pressed = true; return 1; }
 264                 }
 265                 in.close();
 266
 267                 return 0;
 268         }
 269         catch(exception& e) {
 270                 m->errorOut(e, "DeuniqueTreeCommand", "readNamesFile");
 271                 exit(1);
 272         }
 273 }
 274 /***********************************************************/
 275
 276
 277
 278