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