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