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