indicatorcommand.cpp

   1 /*
   2  *  indicatorcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 11/12/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "indicatorcommand.h"
  11 #include "sharedutilities.h"
  12
  13
  14 //**********************************************************************************************************************
  15 vector<string> IndicatorCommand::setParameters(){
  16         try {
  17                 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
  18                 CommandParameter pdesign("design", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(pdesign);
  19                 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(pshared);
  20                 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(prelabund);
  21                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  22                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
  23                 CommandParameter ptree("tree", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(ptree);
  24                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  25                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  26                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
  27
  28                 vector<string> myArray;
  29                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  30                 return myArray;
  31         }
  32         catch(exception& e) {
  33                 m->errorOut(e, "IndicatorCommand", "setParameters");
  34                 exit(1);
  35         }
  36 }
  37 //**********************************************************************************************************************
  38 string IndicatorCommand::getHelpString(){
  39         try {
  40                 string helpString = "";
  41                 helpString += "The indicator command can be run in 3 ways: with a shared or relabund file and a design file, or with a shared or relabund file and a tree file, or with a shared or relabund file, tree file and design file. \n";
  42                 helpString += "The indicator command outputs a .indicator.summary file and a .indicator.tre if a tree is given. \n";
  43                 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";
  44                 helpString += "The summary file lists the indicator value for each OTU for each node.\n";
  45                 helpString += "The indicator command parameters are tree, groups, shared, relabund, design and label. \n";
  46                 helpString += "The design parameter allows you to relate the tree to the shared or relabund file, if your tree contains the grouping names, or if no tree is provided to group your groups into groupings.\n";
  47                 helpString += "The groups parameter allows you to specify which of the groups in your shared or relabund you would like analyzed, or if you provide a design file the groups in your design file.  The groups may be entered separated by dashes.\n";
  48                 helpString += "The label parameter indicates at what distance your tree relates to the shared or relabund.\n";
  49                 helpString += "The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n";
  50                 helpString += "The iters parameter allows you to set number of randomization for the P value.  The default is 1000.";
  51                 helpString += "The indicator command should be used in the following format: indicator(tree=test.tre, shared=test.shared, label=0.03)\n";
  52                 helpString += "Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreefile).\n";
  53                 return helpString;
  54         }
  55         catch(exception& e) {
  56                 m->errorOut(e, "IndicatorCommand", "getHelpString");
  57                 exit(1);
  58         }
  59 }
  60
  61 //**********************************************************************************************************************
  62 IndicatorCommand::IndicatorCommand(){
  63         try {
  64                 abort = true; calledHelp = true;
  65                 setParameters();
  66                 vector<string> tempOutNames;
  67                 outputTypes["tree"] = tempOutNames;
  68                 outputTypes["summary"] = tempOutNames;
  69         }
  70         catch(exception& e) {
  71                 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
  72                 exit(1);
  73         }
  74 }
  75 //**********************************************************************************************************************
  76 IndicatorCommand::IndicatorCommand(string option)  {
  77         try {
  78                 abort = false; calledHelp = false;
  79
  80                 //allow user to run help
  81                 if(option == "help") { help(); abort = true; calledHelp = true; }
  82                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  83
  84                 else {
  85                         vector<string> myArray = setParameters();
  86
  87                         OptionParser parser(option);
  88                         map<string, string> parameters = parser.getParameters();
  89
  90                         ValidParameters validParameter;
  91                         map<string, string>::iterator it;
  92
  93                         //check to make sure all parameters are valid for command
  94                         for (it = parameters.begin(); it != parameters.end(); it++) {
  95                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  96                         }
  97
  98                         m->runParse = true;
  99                         m->clearGroups();
 100                         m->clearAllGroups();
 101                         m->Treenames.clear();
 102                         m->names.clear();
 103
 104                         vector<string> tempOutNames;
 105                         outputTypes["tree"] = tempOutNames;
 106                         outputTypes["summary"] = tempOutNames;
 107
 108                         //if the user changes the input directory command factory will send this info to us in the output parameter
 109                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 110                         if (inputDir == "not found"){   inputDir = "";          }
 111                         else {
 112                                 string path;
 113                                 it = parameters.find("tree");
 114                                 //user has given a template file
 115                                 if(it != parameters.end()){
 116                                         path = m->hasPath(it->second);
 117                                         //if the user has not given a path then, add inputdir. else leave path alone.
 118                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
 119                                 }
 120
 121                                 it = parameters.find("shared");
 122                                 //user has given a template file
 123                                 if(it != parameters.end()){
 124                                         path = m->hasPath(it->second);
 125                                         //if the user has not given a path then, add inputdir. else leave path alone.
 126                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 127                                 }
 128
 129                                 it = parameters.find("relabund");
 130                                 //user has given a template file
 131                                 if(it != parameters.end()){
 132                                         path = m->hasPath(it->second);
 133                                         //if the user has not given a path then, add inputdir. else leave path alone.
 134                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
 135                                 }
 136
 137                                 it = parameters.find("design");
 138                                 //user has given a template file
 139                                 if(it != parameters.end()){
 140                                         path = m->hasPath(it->second);
 141                                         //if the user has not given a path then, add inputdir. else leave path alone.
 142                                         if (path == "") {       parameters["design"] = inputDir + it->second;           }
 143                                 }
 144                         }
 145
 146                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 147
 148                         //check for required parameters
 149                         treefile = validParameter.validFile(parameters, "tree", true);
 150                         if (treefile == "not open") { treefile = ""; abort = true; }
 151                         else if (treefile == "not found") { treefile = "";      }
 152                         else { m->setTreeFile(treefile); }
 153
 154                         sharedfile = validParameter.validFile(parameters, "shared", true);
 155                         if (sharedfile == "not open") { abort = true; }
 156                         else if (sharedfile == "not found") { sharedfile = ""; }
 157                         else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
 158
 159                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 160                         if (relabundfile == "not open") { abort = true; }
 161                         else if (relabundfile == "not found") { relabundfile = ""; }
 162                         else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
 163
 164                         designfile = validParameter.validFile(parameters, "design", true);
 165                         if (designfile == "not open") { designfile = ""; abort = true; }
 166                         else if (designfile == "not found") { designfile = ""; }
 167                         else { m->setDesignFile(designfile); }
 168
 169                         groups = validParameter.validFile(parameters, "groups", false);
 170                         if (groups == "not found") { groups = "";  Groups.push_back("all"); }
 171                         else { m->splitAtDash(groups, Groups);  }
 172                         m->setGroups(Groups);
 173
 174                         label = validParameter.validFile(parameters, "label", false);
 175                         if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }
 176
 177                         string temp = validParameter.validFile(parameters, "iters", false);             if (temp == "not found") { temp = "1000"; }
 178                         m->mothurConvert(temp, iters);
 179
 180                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 181                         m->setProcessors(temp);
 182                         m->mothurConvert(temp, processors);
 183
 184                         if ((relabundfile == "") && (sharedfile == "")) {
 185                                 //is there are current file available for either of these?
 186                                 //give priority to shared, then relabund
 187                                 sharedfile = m->getSharedFile();
 188                                 if (sharedfile != "") {  inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 189                                 else {
 190                                         relabundfile = m->getRelAbundFile();
 191                                         if (relabundfile != "") {  inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
 192                                         else {
 193                                                 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
 194                                                 abort = true;
 195                                         }
 196                                 }
 197                         }
 198
 199
 200                         if ((designfile == "") && (treefile == "")) {
 201                                 treefile = m->getTreeFile();
 202                                 if (treefile != "") {  m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
 203                                 else {
 204                                         designfile = m->getDesignFile();
 205                                         if (designfile != "") {  m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
 206                                         else {
 207                                                 m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true;
 208                                         }
 209                                 }
 210                         }
 211
 212                         if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true;  }
 213
 214
 215                 }
 216         }
 217         catch(exception& e) {
 218                 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
 219                 exit(1);
 220         }
 221 }
 222 //**********************************************************************************************************************
 223
 224 int IndicatorCommand::execute(){
 225         try {
 226
 227                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 228
 229                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 230
 231                 int start = time(NULL);
 232
 233                 //read designfile if given and set up groups for read of sharedfiles
 234                 if (designfile != "") {
 235                         designMap = new GroupMap(designfile);
 236                         designMap->readDesignMap();
 237
 238                         //fill Groups - checks for "all" and for any typo groups
 239                         SharedUtil* util = new SharedUtil();
 240                         vector<string> nameGroups = designMap->getNamesOfGroups();
 241                         util->setGroups(Groups, nameGroups);
 242                         designMap->setNamesOfGroups(nameGroups);
 243                         delete util;
 244
 245                         //loop through the Groups and fill Globaldata's Groups with the design file info
 246                         vector<string> namesSeqs = designMap->getNamesSeqs(Groups);
 247                         m->setGroups(namesSeqs);
 248                 }
 249
 250                 /***************************************************/
 251                 // use smart distancing to get right sharedRabund  //
 252                 /***************************************************/
 253                 if (sharedfile != "") {
 254                         getShared();
 255                         if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } return 0; }
 256                         if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
 257                 }else {
 258                         getSharedFloat();
 259                         if (m->control_pressed) {  if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 260                         if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
 261                 }
 262
 263                 //reset groups if needed
 264                 if (designfile != "") { m->setGroups(Groups); }
 265
 266                 /***************************************************/
 267                 //    reading tree info                                                    //
 268                 /***************************************************/
 269                 if (treefile != "") {
 270                         string groupfile = "";
 271                         m->setTreeFile(treefile);
 272                         Tree* tree = new Tree(treefile); delete tree;  //extracts names from tree to make faked out groupmap
 273                         treeMap = new TreeMap();
 274                         bool mismatch = false;
 275
 276                         for (int i = 0; i < m->Treenames.size(); i++) {
 277                                 //sanity check - is this a group that is not in the sharedfile?
 278                                 if (designfile == "") {
 279                                         if (!(m->inUsersGroups(m->Treenames[i], m->getAllGroups()))) {
 280                                                 m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
 281                                                 mismatch = true;
 282                                         }
 283                                         treeMap->addSeq(m->Treenames[i], "Group1");
 284                                 }else{
 285                                         vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
 286                                         vector<string> myNames = designMap->getNamesSeqs(myGroups);
 287
 288                                         for(int k = 0; k < myNames.size(); k++) {
 289                                                 if (!(m->inUsersGroups(myNames[k], m->getAllGroups()))) {
 290                                                         m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
 291                                                         mismatch = true;
 292                                                 }
 293                                         }
 294                                         treeMap->addSeq(m->Treenames[i], "Group1");
 295                                 }
 296                         }
 297
 298                         if ((designfile != "") && (m->Treenames.size() != Groups.size())) { cout << Groups.size() << '\t' << m->Treenames.size() << endl; m->mothurOut("[ERROR]: You design file does not match your tree, aborting."); m->mothurOutEndLine(); mismatch = true; }
 299
 300                         if (mismatch) { //cleanup and exit
 301                                 if (designfile != "") { delete designMap; }
 302                                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 303                                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 304                                 delete treeMap;
 305                                 return 0;
 306                         }
 307
 308                         read = new ReadNewickTree(treefile);
 309                         int readOk = read->read(treeMap);
 310
 311                         if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
 312
 313                         vector<Tree*> T = read->getTrees();
 314
 315                         delete read;
 316
 317                         if (m->control_pressed) {
 318                                 if (designfile != "") { delete designMap; }
 319                                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 320                                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 321                                 for (int i = 0; i < T.size(); i++) {  delete T[i];  }  delete treeMap; return 0;
 322                         }
 323
 324                         T[0]->assembleTree();
 325
 326                         /***************************************************/
 327                         //    create ouptut tree - respecting pickedGroups //
 328                         /***************************************************/
 329                         Tree* outputTree = new Tree(m->getNumGroups(), treeMap);
 330
 331                         outputTree->getSubTree(T[0], m->getGroups());
 332                         outputTree->assembleTree();
 333
 334                         //no longer need original tree, we have output tree to use and label
 335                         for (int i = 0; i < T.size(); i++) {  delete T[i];  }
 336
 337
 338                         if (m->control_pressed) {
 339                                 if (designfile != "") { delete designMap; }
 340                                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 341                                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 342                                 delete outputTree; delete treeMap;  return 0;
 343                         }
 344
 345                         /***************************************************/
 346                         //              get indicator species values                       //
 347                         /***************************************************/
 348                         GetIndicatorSpecies(outputTree);
 349                         delete outputTree; delete treeMap;
 350
 351                 }else { //run with design file only
 352                         //get indicator species
 353                         GetIndicatorSpecies();
 354                 }
 355
 356                 if (designfile != "") { delete designMap; }
 357                 if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 358                 else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 359
 360                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]);        } return 0; }
 361
 362                 //set tree file as new current treefile
 363                 if (treefile != "") {
 364                         string current = "";
 365                         itTypes = outputTypes.find("tree");
 366                         if (itTypes != outputTypes.end()) {
 367                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
 368                         }
 369                 }
 370
 371                 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
 372                 m->mothurOutEndLine();
 373                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 374                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 375                 m->mothurOutEndLine();
 376
 377                 return 0;
 378         }
 379         catch(exception& e) {
 380                 m->errorOut(e, "IndicatorCommand", "execute");
 381                 exit(1);
 382         }
 383 }
 384 //**********************************************************************************************************************
 385 //divide shared or relabund file by groupings in the design file
 386 //report all otu values to file
 387 int IndicatorCommand::GetIndicatorSpecies(){
 388         try {
 389                 string thisOutputDir = outputDir;
 390                 if (outputDir == "") {  thisOutputDir += m->hasPath(inputFileName);  }
 391                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
 392                 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
 393
 394                 ofstream out;
 395                 m->openOutputFile(outputFileName, out);
 396                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 397                 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
 398
 399                 int numBins = 0;
 400                 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
 401                 else { numBins = lookupFloat[0]->getNumBins(); }
 402
 403                 if (m->control_pressed) { out.close(); return 0; }
 404
 405                 /*****************************************************/
 406                 //create vectors containing rabund info                          //
 407                 /*****************************************************/
 408
 409                 vector<float> indicatorValues; //size of numBins
 410                 vector<float> pValues;
 411                 map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
 412
 413                 if (sharedfile != "") {
 414                         vector< vector<SharedRAbundVector*> > groupings;
 415                         set<string> groupsAlreadyAdded;
 416                         vector<SharedRAbundVector*> subset;
 417
 418                         //for each grouping
 419                         for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
 420
 421                                 for (int k = 0; k < lookup.size(); k++) {
 422                                         //are you from this grouping?
 423                                         if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
 424                                                 subset.push_back(lookup[k]);
 425                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 426                                         }
 427                                 }
 428                                 if (subset.size() != 0) { groupings.push_back(subset); }
 429                                 subset.clear();
 430                         }
 431
 432                         if (groupsAlreadyAdded.size() != lookup.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 433
 434                         indicatorValues = getValues(groupings, randomGroupingsMap);
 435
 436                         pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
 437                 }else {
 438                         vector< vector<SharedRAbundFloatVector*> > groupings;
 439                         set<string> groupsAlreadyAdded;
 440                         vector<SharedRAbundFloatVector*> subset;
 441
 442                         //for each grouping
 443                         for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
 444                                 for (int k = 0; k < lookupFloat.size(); k++) {
 445                                         //are you from this grouping?
 446                                         if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
 447                                                 subset.push_back(lookupFloat[k]);
 448                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 449                                         }
 450                                 }
 451                                 if (subset.size() != 0) { groupings.push_back(subset); }
 452                                 subset.clear();
 453                         }
 454
 455                         if (groupsAlreadyAdded.size() != lookupFloat.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 456
 457                         indicatorValues = getValues(groupings, randomGroupingsMap);
 458
 459                         pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
 460                 }
 461
 462                 if (m->control_pressed) { out.close(); return 0; }
 463
 464
 465                 /******************************************************/
 466                 //output indicator values to table form               //
 467                 /*****************************************************/
 468                 out << "OTU\tIndicator_Value\tpValue" << endl;
 469                 for (int j = 0; j < indicatorValues.size(); j++) {
 470
 471                         if (m->control_pressed) { out.close(); return 0; }
 472
 473                         out << (j+1) << '\t' << indicatorValues[j] << '\t';
 474
 475                         if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
 476                         else { out << "<" << (1/(float)iters) << endl; }
 477
 478                         if (pValues[j] <= 0.05) {
 479                                 cout << "OTU" << j+1 << '\t' << indicatorValues[j]  << '\t';
 480                                 string pValueString = "<" + toString((1/(float)iters));
 481                                 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
 482                                 else { cout << "<" << (1/(float)iters); }
 483                                 m->mothurOutJustToLog("OTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
 484                                 m->mothurOutEndLine();
 485                         }
 486                 }
 487
 488                 out.close();
 489
 490                 return 0;
 491         }
 492         catch(exception& e) {
 493                 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
 494                 exit(1);
 495         }
 496 }
 497 //**********************************************************************************************************************
 498 //traverse tree finding indicator species values for each otu at each node
 499 //label node with otu number that has highest indicator value
 500 //report all otu values to file
 501 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
 502         try {
 503
 504                 string thisOutputDir = outputDir;
 505                 if (outputDir == "") {  thisOutputDir += m->hasPath(inputFileName);  }
 506                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
 507                 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
 508
 509                 ofstream out;
 510                 m->openOutputFile(outputFileName, out);
 511                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 512
 513                 int numBins = 0;
 514                 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
 515                 else { numBins = lookupFloat[0]->getNumBins(); }
 516
 517                 //print headings
 518                 out << "TreeNode\t";
 519                 for (int i = 0; i < numBins; i++) { out << "OTU" << (i+1) << "_IndValue" << '\t' << "pValue" << '\t'; }
 520                 out << endl;
 521
 522                 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
 523
 524                 string treeOutputDir = outputDir;
 525                 if (outputDir == "") {  treeOutputDir += m->hasPath(treefile);  }
 526                 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
 527
 528
 529                 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
 530                 map<int, set<string> > nodeToDescendants;
 531                 map<int, set<int> > descendantNodes;
 532                 for (int i = 0; i < T->getNumNodes(); i++) {
 533                         if (m->control_pressed) { return 0; }
 534
 535                         nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
 536                 }
 537
 538                 //you need the distances to leaf to decide grouping below
 539                 //this will also set branch lengths if the tree does not include them
 540                 map<int, float> distToRoot = getDistToRoot(T);
 541
 542                 //for each node
 543                 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
 544                         //cout << endl << i+1 << endl;
 545                         if (m->control_pressed) { out.close(); return 0; }
 546
 547                         /*****************************************************/
 548                         //create vectors containing rabund info                          //
 549                         /*****************************************************/
 550
 551                         vector<float> indicatorValues; //size of numBins
 552                         vector<float> pValues;
 553                         map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
 554
 555                         if (sharedfile != "") {
 556                                 vector< vector<SharedRAbundVector*> > groupings;
 557
 558                                 //get nodes that will be a valid grouping
 559                                 //you are valid if you are not one of my descendants
 560                                 //AND your distToRoot is >= mine
 561                                 //AND you were not added as part of a larger grouping. Largest nodes are added first.
 562
 563                                 set<string> groupsAlreadyAdded;
 564                                 //create a grouping with my grouping
 565                                 vector<SharedRAbundVector*> subset;
 566                                 int count = 0;
 567                                 int doneCount = nodeToDescendants[i].size();
 568                                 for (int k = 0; k < lookup.size(); k++) {
 569                                         //is this descendant of i
 570                                         if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
 571                                                 subset.push_back(lookup[k]);
 572                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 573                                                 count++;
 574                                         }
 575                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 576                                 }
 577                                 if (subset.size() != 0) { groupings.push_back(subset); }
 578
 579
 580                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 581
 582
 583                                         if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
 584                                                 vector<SharedRAbundVector*> subset;
 585                                                 int count = 0;
 586                                                 int doneCount = nodeToDescendants[j].size();
 587                                                 for (int k = 0; k < lookup.size(); k++) {
 588                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 589                                                         if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
 590                                                                 subset.push_back(lookup[k]);
 591                                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 592                                                                 count++;
 593                                                         }
 594                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 595                                                 }
 596
 597                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 598                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 599                                         }
 600                                 }
 601
 602                                 if (groupsAlreadyAdded.size() != lookup.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 603
 604                                 indicatorValues = getValues(groupings, randomGroupingsMap);
 605
 606                                 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
 607                         }else {
 608                                 vector< vector<SharedRAbundFloatVector*> > groupings;
 609
 610                                 //get nodes that will be a valid grouping
 611                                 //you are valid if you are not one of my descendants
 612                                 //AND your distToRoot is >= mine
 613                                 //AND you were not added as part of a larger grouping. Largest nodes are added first.
 614
 615                                 set<string> groupsAlreadyAdded;
 616                                 //create a grouping with my grouping
 617                                 vector<SharedRAbundFloatVector*> subset;
 618                                 int count = 0;
 619                                 int doneCount = nodeToDescendants[i].size();
 620                                 for (int k = 0; k < lookupFloat.size(); k++) {
 621                                         //is this descendant of i
 622                                         if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
 623                                                 subset.push_back(lookupFloat[k]);
 624                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 625                                                 count++;
 626                                         }
 627                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 628                                 }
 629                                 if (subset.size() != 0) { groupings.push_back(subset); }
 630
 631                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 632                                         if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
 633                                                 vector<SharedRAbundFloatVector*> subset;
 634                                                 int count = 0;
 635                                                 int doneCount = nodeToDescendants[j].size();
 636                                                 for (int k = 0; k < lookupFloat.size(); k++) {
 637                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 638                                                         if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
 639                                                                 subset.push_back(lookupFloat[k]);
 640                                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 641                                                                 count++;
 642                                                         }
 643                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 644                                                 }
 645
 646                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 647                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 648                                         }
 649                                 }
 650
 651                                 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 652
 653                                 indicatorValues = getValues(groupings, randomGroupingsMap);
 654
 655                                 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
 656                         }
 657
 658                         if (m->control_pressed) { out.close(); return 0; }
 659
 660
 661                         /******************************************************/
 662                         //output indicator values to table form + label tree  //
 663                         /*****************************************************/
 664                         out << (i+1) << '\t';
 665                         for (int j = 0; j < indicatorValues.size(); j++) {
 666
 667                                 if (m->control_pressed) { out.close(); return 0; }
 668
 669                                 if (pValues[j] < (1/(float)iters)) {
 670                                         out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
 671                                 }else {
 672                                         out << indicatorValues[j] << '\t' << pValues[j] << '\t';
 673                                 }
 674
 675                                 if (pValues[j] <= 0.05) {
 676                                         cout << i+1 << "\tOTU" << j+1 << '\t' << indicatorValues[j]  << '\t';
 677                                         string pValueString = "<" + toString((1/(float)iters));
 678                                         if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
 679                                         else { cout << "<" << (1/(float)iters); }
 680                                         m->mothurOutJustToLog(toString(i) + "\tOTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
 681                                         m->mothurOutEndLine();
 682                                 }
 683                         }
 684                         out << endl;
 685
 686                         T->tree[i].setLabel((i+1));
 687                 }
 688                 out.close();
 689
 690                 ofstream outTree;
 691                 m->openOutputFile(outputTreeFileName, outTree);
 692                 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
 693
 694                 T->print(outTree, "both");
 695                 outTree.close();
 696
 697                 return 0;
 698         }
 699         catch(exception& e) {
 700                 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
 701                 exit(1);
 702         }
 703 }
 704 //**********************************************************************************************************************
 705 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
 706         try {
 707                 vector<float> values;
 708                 map< vector<int>, vector<int> >::iterator it;
 709
 710                 //for each otu
 711                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 712
 713                         if (m->control_pressed) { return values; }
 714
 715                         vector<float> terms;
 716                         float AijDenominator = 0.0;
 717                         vector<float> Bij;
 718
 719                         //get overall abundance of each grouping
 720                         for (int j = 0; j < groupings.size(); j++) {
 721
 722                                 float totalAbund = 0;
 723                                 int numNotZero = 0;
 724                                 for (int k = 0; k < groupings[j].size(); k++) {
 725                                         vector<int> temp; temp.push_back(j); temp.push_back(k);
 726                                         it = groupingsMap.find(temp);
 727
 728                                         if (it == groupingsMap.end()) { //this one didnt get moved
 729                                                 totalAbund += groupings[j][k]->getAbundance(i);
 730                                                 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
 731                                         }else {
 732                                                 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
 733                                                 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
 734                                         }
 735
 736                                 }
 737
 738                                 //mean abundance
 739                                 float Aij = (totalAbund / (float) groupings[j].size());
 740                                 terms.push_back(Aij);
 741
 742                                 //percentage of sites represented
 743                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 744
 745                                 AijDenominator += Aij;
 746                         }
 747
 748                         float maxIndVal = 0.0;
 749                         for (int j = 0; j < terms.size(); j++) {
 750                                 float thisAij = (terms[j] / AijDenominator); //relative abundance
 751                                 float thisValue = thisAij * Bij[j] * 100.0;
 752
 753                                 //save largest
 754                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 755                         }
 756
 757                         values.push_back(maxIndVal);
 758                 }
 759
 760                 return values;
 761         }
 762         catch(exception& e) {
 763                 m->errorOut(e, "IndicatorCommand", "getValues");
 764                 exit(1);
 765         }
 766 }
 767 //**********************************************************************************************************************
 768 //same as above, just data type difference
 769 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
 770         try {
 771                 vector<float> values;
 772
 773         /*for (int j = 0; j < groupings.size(); j++) {
 774                 cout << "grouping " << j << endl;
 775                 for (int k = 0; k < groupings[j].size(); k++) {
 776                         cout << groupings[j][k]->getGroup() << endl;
 777                 }
 778         }*/
 779                 map< vector<int>, vector<int> >::iterator it;
 780
 781                 //for each otu
 782                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 783                         vector<float> terms;
 784                         float AijDenominator = 0.0;
 785                         vector<float> Bij;
 786                         //get overall abundance of each grouping
 787                         for (int j = 0; j < groupings.size(); j++) {
 788
 789                                 int totalAbund = 0.0;
 790                                 int numNotZero = 0;
 791                                 for (int k = 0; k < groupings[j].size(); k++) {
 792                                         vector<int> temp; temp.push_back(j); temp.push_back(k);
 793                                         it = groupingsMap.find(temp);
 794
 795                                         if (it == groupingsMap.end()) { //this one didnt get moved
 796                                                 totalAbund += groupings[j][k]->getAbundance(i);
 797                                                 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
 798                                         }else {
 799                                                 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
 800                                                 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
 801                                         }
 802
 803                                 }
 804
 805                                 //mean abundance
 806                                 float Aij = (totalAbund / (float) groupings[j].size());
 807                                 terms.push_back(Aij);
 808
 809                                 //percentage of sites represented
 810                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 811
 812                                 AijDenominator += Aij;
 813                         }
 814
 815                         float maxIndVal = 0.0;
 816                         for (int j = 0; j < terms.size(); j++) {
 817                                 float thisAij = (terms[j] / AijDenominator); //relative abundance
 818                                 float thisValue = thisAij * Bij[j] * 100.0;
 819
 820                                 //save largest
 821                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 822                         }
 823
 824                         values.push_back(maxIndVal);
 825                 }
 826
 827                 return values;
 828         }
 829         catch(exception& e) {
 830                 m->errorOut(e, "IndicatorCommand", "getValues");
 831                 exit(1);
 832         }
 833 }
 834 //**********************************************************************************************************************
 835 //you need the distances to root to decide groupings
 836 //this will also set branch lengths if the tree does not include them
 837 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
 838         try {
 839                 map<int, float> dists;
 840
 841                 bool hasBranchLengths = false;
 842                 for (int i = 0; i < T->getNumNodes(); i++) {
 843                         if (T->tree[i].getBranchLength() > 0.0) {  hasBranchLengths = true; break; }
 844                 }
 845
 846                 //set branchlengths if needed
 847                 if (!hasBranchLengths) {
 848                         for (int i = 0; i < T->getNumNodes(); i++) {
 849
 850                                 int lc = T->tree[i].getLChild();
 851                                 int rc = T->tree[i].getRChild();
 852
 853                                 if (lc == -1) { // you are a leaf
 854                                         //if you are a leaf set you priliminary length to 1.0, this may adjust later
 855                                         T->tree[i].setBranchLength(1.0);
 856                                         dists[i] = 1.0;
 857                                 }else{ // you are an internal node
 858                                         //look at your children's length to leaf
 859                                         float ldist = dists[lc];
 860                                         float rdist = dists[rc];
 861
 862                                         float greater = ldist;
 863                                         if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
 864                                         else { dists[i] = rdist + 1.0; }
 865
 866
 867                                         //branch length = difference + 1
 868                                         T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
 869                                         T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
 870                                 }
 871                         }
 872                 }
 873
 874                 dists.clear();
 875
 876                 for (int i = 0; i < T->getNumNodes(); i++) {
 877
 878                         double sum = 0.0;
 879                         int index = i;
 880
 881                         while(T->tree[index].getParent() != -1){
 882                                 if (T->tree[index].getBranchLength() != -1) {
 883                                         sum += abs(T->tree[index].getBranchLength());
 884                                 }
 885                                 index = T->tree[index].getParent();
 886                         }
 887
 888                         dists[i] = sum;
 889                 }
 890
 891                 return dists;
 892         }
 893         catch(exception& e) {
 894                 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
 895                 exit(1);
 896         }
 897 }
 898 //**********************************************************************************************************************
 899 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
 900         try {
 901                 set<string> names;
 902
 903                 set<string>::iterator it;
 904
 905                 int lc = T->tree[i].getLChild();
 906                 int rc = T->tree[i].getRChild();
 907
 908                 if (lc == -1) { //you are a leaf your only descendant is yourself
 909                         set<int> temp; temp.insert(i);
 910                         nodes[i] = temp;
 911
 912                         if (designfile == "") {
 913                                 names.insert(T->tree[i].getName());
 914                         }else {
 915                                 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
 916                                 vector<string> myReps = designMap->getNamesSeqs(myGroup);
 917                                 for (int k = 0; k < myReps.size(); k++) {
 918                                         names.insert(myReps[k]);
 919                                 }
 920                         }
 921
 922                 }else{ //your descedants are the combination of your childrens descendants
 923                         names = descendants[lc];
 924                         nodes[i] = nodes[lc];
 925                         for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
 926                                 names.insert(*it);
 927                         }
 928                         for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
 929                                 nodes[i].insert(*itNum);
 930                         }
 931                         //you are your own descendant
 932                         nodes[i].insert(i);
 933                 }
 934
 935                 return names;
 936         }
 937         catch(exception& e) {
 938                 m->errorOut(e, "IndicatorCommand", "getDescendantList");
 939                 exit(1);
 940         }
 941 }
 942 //**********************************************************************************************************************
 943 int IndicatorCommand::getShared(){
 944         try {
 945                 InputData* input = new InputData(sharedfile, "sharedfile");
 946                 lookup = input->getSharedRAbundVectors();
 947                 string lastLabel = lookup[0]->getLabel();
 948
 949                 if (label == "") { label = lastLabel; delete input; return 0; }
 950
 951                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 952                 set<string> labels; labels.insert(label);
 953                 set<string> processedLabels;
 954                 set<string> userLabels = labels;
 955
 956                 //as long as you are not at the end of the file or done wih the lines you want
 957                 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
 958                         if (m->control_pressed) {  delete input; return 0;  }
 959
 960                         if(labels.count(lookup[0]->getLabel()) == 1){
 961                                 processedLabels.insert(lookup[0]->getLabel());
 962                                 userLabels.erase(lookup[0]->getLabel());
 963                                 break;
 964                         }
 965
 966                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 967                                 string saveLabel = lookup[0]->getLabel();
 968
 969                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 970                                 lookup = input->getSharedRAbundVectors(lastLabel);
 971
 972                                 processedLabels.insert(lookup[0]->getLabel());
 973                                 userLabels.erase(lookup[0]->getLabel());
 974
 975                                 //restore real lastlabel to save below
 976                                 lookup[0]->setLabel(saveLabel);
 977                                 break;
 978                         }
 979
 980                         lastLabel = lookup[0]->getLabel();
 981
 982                         //get next line to process
 983                         //prevent memory leak
 984                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 985                         lookup = input->getSharedRAbundVectors();
 986                 }
 987
 988
 989                 if (m->control_pressed) { delete input; return 0;  }
 990
 991                 //output error messages about any remaining user labels
 992                 set<string>::iterator it;
 993                 bool needToRun = false;
 994                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 995                         m->mothurOut("Your file does not include the label " + *it);
 996                         if (processedLabels.count(lastLabel) != 1) {
 997                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 998                                 needToRun = true;
 999                         }else {
1000                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1001                         }
1002                 }
1003
1004                 //run last label if you need to
1005                 if (needToRun == true)  {
1006                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
1007                         lookup = input->getSharedRAbundVectors(lastLabel);
1008                 }
1009
1010                 delete input;
1011                 return 0;
1012         }
1013         catch(exception& e) {
1014                 m->errorOut(e, "IndicatorCommand", "getShared");
1015                 exit(1);
1016         }
1017 }
1018 //**********************************************************************************************************************
1019 int IndicatorCommand::getSharedFloat(){
1020         try {
1021                 InputData* input = new InputData(relabundfile, "relabund");
1022                 lookupFloat = input->getSharedRAbundFloatVectors();
1023                 string lastLabel = lookupFloat[0]->getLabel();
1024
1025                 if (label == "") { label = lastLabel; delete input; return 0; }
1026
1027                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1028                 set<string> labels; labels.insert(label);
1029                 set<string> processedLabels;
1030                 set<string> userLabels = labels;
1031
1032                 //as long as you are not at the end of the file or done wih the lines you want
1033                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1034
1035                         if (m->control_pressed) {  delete input; return 0;  }
1036
1037                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
1038                                 processedLabels.insert(lookupFloat[0]->getLabel());
1039                                 userLabels.erase(lookupFloat[0]->getLabel());
1040                                 break;
1041                         }
1042
1043                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1044                                 string saveLabel = lookupFloat[0]->getLabel();
1045
1046                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
1047                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1048
1049                                 processedLabels.insert(lookupFloat[0]->getLabel());
1050                                 userLabels.erase(lookupFloat[0]->getLabel());
1051
1052                                 //restore real lastlabel to save below
1053                                 lookupFloat[0]->setLabel(saveLabel);
1054                                 break;
1055                         }
1056
1057                         lastLabel = lookupFloat[0]->getLabel();
1058
1059                         //get next line to process
1060                         //prevent memory leak
1061                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
1062                         lookupFloat = input->getSharedRAbundFloatVectors();
1063                 }
1064
1065
1066                 if (m->control_pressed) { delete input; return 0;  }
1067
1068                 //output error messages about any remaining user labels
1069                 set<string>::iterator it;
1070                 bool needToRun = false;
1071                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1072                         m->mothurOut("Your file does not include the label " + *it);
1073                         if (processedLabels.count(lastLabel) != 1) {
1074                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1075                                 needToRun = true;
1076                         }else {
1077                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1078                         }
1079                 }
1080
1081                 //run last label if you need to
1082                 if (needToRun == true)  {
1083                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }
1084                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1085                 }
1086
1087                 delete input;
1088                 return 0;
1089         }
1090         catch(exception& e) {
1091                 m->errorOut(e, "IndicatorCommand", "getShared");
1092         exit(1);
1093         }
1094 }
1095 //**********************************************************************************************************************
1096 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1097         try {
1098                 vector<float> pvalues;
1099                 pvalues.resize(indicatorValues.size(), 0);
1100
1101                 for(int i=0;i<numIters;i++){
1102                         if (m->control_pressed) { break; }
1103                         groupingsMap = randomizeGroupings(groupings, num);
1104                         vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1105
1106                         for (int j = 0; j < indicatorValues.size(); j++) {
1107                                 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1108                         }
1109                 }
1110
1111                 return pvalues;
1112
1113         }catch(exception& e) {
1114                 m->errorOut(e, "IndicatorCommand", "driver");
1115                 exit(1);
1116         }
1117 }
1118 //**********************************************************************************************************************
1119 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1120         try {
1121                 vector<float> pvalues;
1122
1123 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1124                 if(processors == 1){
1125                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1126                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1127                 }else{
1128
1129                         //divide iters between processors
1130                         int numItersPerProcessor = iters / processors;
1131
1132                         vector<int> processIDS;
1133                         int process = 1;
1134                         int num = 0;
1135
1136                         //loop through and create all the processes you want
1137                         while (process != processors) {
1138                                 int pid = fork();
1139
1140                                 if (pid > 0) {
1141                                         processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1142                                         process++;
1143                                 }else if (pid == 0){
1144                                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1145
1146                                         //pass pvalues to parent
1147                                         ofstream out;
1148                                         string tempFile = toString(getpid()) + ".pvalues.temp";
1149                                         m->openOutputFile(tempFile, out);
1150
1151                                         //pass values
1152                                         for (int i = 0; i < pvalues.size(); i++) {
1153                                                 out << pvalues[i] << '\t';
1154                                         }
1155                                         out << endl;
1156
1157                                         out.close();
1158
1159                                         exit(0);
1160                                 }else {
1161                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1162                                         for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1163                                         exit(0);
1164                                 }
1165                         }
1166
1167                         //do my part
1168                         //special case for last processor in case it doesn't divide evenly
1169                         numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1170                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1171
1172                         //force parent to wait until all the processes are done
1173                         for (int i=0;i<processIDS.size();i++) {
1174                                 int temp = processIDS[i];
1175                                 wait(&temp);
1176                         }
1177
1178                         //combine results
1179                         for (int i = 0; i < processIDS.size(); i++) {
1180                                 ifstream in;
1181                                 string tempFile =  toString(processIDS[i]) + ".pvalues.temp";
1182                                 m->openInputFile(tempFile, in);
1183
1184                                 ////// to do ///////////
1185                                 int numTemp; numTemp = 0;
1186                                 for (int j = 0; j < pvalues.size(); j++) {
1187                                         in >> numTemp; m->gobble(in);
1188                                         pvalues[j] += numTemp;
1189                                 }
1190                                 in.close(); m->mothurRemove(tempFile);
1191                         }
1192                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1193                 }
1194 #else
1195                 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1196                 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1197 #endif
1198
1199                 return pvalues;
1200         }
1201         catch(exception& e) {
1202                 m->errorOut(e, "IndicatorCommand", "getPValues");
1203                 exit(1);
1204         }
1205 }
1206
1207 //**********************************************************************************************************************
1208 //same as above, just data type difference
1209 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1210         try {
1211                 vector<float> pvalues;
1212                 pvalues.resize(indicatorValues.size(), 0);
1213
1214                 for(int i=0;i<numIters;i++){
1215                         if (m->control_pressed) { break; }
1216                         groupingsMap = randomizeGroupings(groupings, num);
1217                         vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1218
1219                         for (int j = 0; j < indicatorValues.size(); j++) {
1220                                 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1221                         }
1222                 }
1223
1224                 return pvalues;
1225
1226         }catch(exception& e) {
1227                 m->errorOut(e, "IndicatorCommand", "driver");
1228                 exit(1);
1229         }
1230 }
1231 //**********************************************************************************************************************
1232 //same as above, just data type difference
1233 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1234         try {
1235                 vector<float> pvalues;
1236
1237 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1238                 if(processors == 1){
1239                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1240                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1241                 }else{
1242
1243                         //divide iters between processors
1244                         int numItersPerProcessor = iters / processors;
1245
1246                         vector<int> processIDS;
1247                         int process = 1;
1248
1249                         //loop through and create all the processes you want
1250                         while (process != processors) {
1251                                 int pid = fork();
1252
1253                                 if (pid > 0) {
1254                                         processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1255                                         process++;
1256                                 }else if (pid == 0){
1257                                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1258
1259                                         //pass pvalues to parent
1260                                         ofstream out;
1261                                         string tempFile = toString(getpid()) + ".pvalues.temp";
1262                                         m->openOutputFile(tempFile, out);
1263
1264                                         //pass values
1265                                         for (int i = 0; i < pvalues.size(); i++) {
1266                                                 out << pvalues[i] << '\t';
1267                                         }
1268                                         out << endl;
1269
1270                                         out.close();
1271
1272                                         exit(0);
1273                                 }else {
1274                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1275                                         for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1276                                         exit(0);
1277                                 }
1278                         }
1279
1280                         //do my part
1281                         //special case for last processor in case it doesn't divide evenly
1282                         numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1283                         pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1284
1285                         //force parent to wait until all the processes are done
1286                         for (int i=0;i<processIDS.size();i++) {
1287                                 int temp = processIDS[i];
1288                                 wait(&temp);
1289                         }
1290
1291                         //combine results
1292                         for (int i = 0; i < processIDS.size(); i++) {
1293                                 ifstream in;
1294                                 string tempFile =  toString(processIDS[i]) + ".pvalues.temp";
1295                                 m->openInputFile(tempFile, in);
1296
1297                                 ////// to do ///////////
1298                                 int numTemp; numTemp = 0;
1299                                 for (int j = 0; j < pvalues.size(); j++) {
1300                                         in >> numTemp; m->gobble(in);
1301                                         pvalues[j] += numTemp;
1302                                 }
1303                                 in.close(); m->mothurRemove(tempFile);
1304                         }
1305                         for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1306                 }
1307 #else
1308                 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1309                 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1310 #endif
1311
1312                 return pvalues;
1313         }
1314         catch(exception& e) {
1315                 m->errorOut(e, "IndicatorCommand", "getPValues");
1316                 exit(1);
1317         }
1318 }
1319 //**********************************************************************************************************************
1320 //swap groups between groupings, in essence randomizing the second column of the design file
1321 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1322         try {
1323
1324                 map< vector<int>, vector<int> > randomGroupings;
1325
1326                 for (int i = 0; i < numLookupGroups; i++) {
1327                         if (m->control_pressed) {break;}
1328
1329                         //get random groups to swap to switch with
1330                         //generate random int between 0 and groupings.size()-1
1331                         int z = m->getRandomIndex(groupings.size()-1);
1332                         int x = m->getRandomIndex(groupings.size()-1);
1333                         int a = m->getRandomIndex(groupings[z].size()-1);
1334                         int b = m->getRandomIndex(groupings[x].size()-1);
1335                         //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1336                         //if ((z < 0) || (z > 1) || x<0 || x>1 || a <0 || a>groupings[z].size()-1 || b<0 || b>groupings[x].size()-1) { cout << "probelm" << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;      }
1337
1338                         vector<int> from;
1339                         vector<int> to;
1340
1341                         from.push_back(z); from.push_back(a);
1342                         to.push_back(x); to.push_back(b);
1343
1344                         randomGroupings[from] = to;
1345                 }
1346                 //cout << "done" << endl;
1347                 return randomGroupings;
1348         }
1349         catch(exception& e) {
1350                 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1351                 exit(1);
1352         }
1353 }
1354 //**********************************************************************************************************************
1355 //swap groups between groupings, in essence randomizing the second column of the design file
1356 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1357         try {
1358
1359                 map< vector<int>, vector<int> > randomGroupings;
1360
1361                 for (int i = 0; i < numLookupGroups; i++) {
1362
1363                         //get random groups to swap to switch with
1364                         //generate random int between 0 and groupings.size()-1
1365                         int z = m->getRandomIndex(groupings.size()-1);
1366                         int x = m->getRandomIndex(groupings.size()-1);
1367                         int a = m->getRandomIndex(groupings[z].size()-1);
1368                         int b = m->getRandomIndex(groupings[x].size()-1);
1369                         //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1370
1371                         vector<int> from;
1372                         vector<int> to;
1373
1374                         from.push_back(z); from.push_back(a);
1375                         to.push_back(x); to.push_back(b);
1376
1377                         randomGroupings[from] = to;
1378                 }
1379
1380                 return randomGroupings;
1381         }
1382         catch(exception& e) {
1383                 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1384                 exit(1);
1385         }
1386 }
1387
1388 /*****************************************************************/
1389
1390