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