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 //**********************************************************************************************************************
  12 vector<string> IndicatorCommand::getValidParameters(){
  13         try {
  14                 string Array[] =  {"tree","shared","relabund","label","groups","outputdir","inputdir"};
  15                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  16                 return myArray;
  17         }
  18         catch(exception& e) {
  19                 m->errorOut(e, "IndicatorCommand", "getValidParameters");
  20                 exit(1);
  21         }
  22 }
  23 //**********************************************************************************************************************
  24 vector<string> IndicatorCommand::getRequiredParameters(){
  25         try {
  26                 string Array[] =  {"tree"};
  27                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  28                 return myArray;
  29         }
  30         catch(exception& e) {
  31                 m->errorOut(e, "IndicatorCommand", "getRequiredParameters");
  32                 exit(1);
  33         }
  34 }
  35 //**********************************************************************************************************************
  36 IndicatorCommand::IndicatorCommand(){
  37         try {
  38                 abort = true;
  39                 //initialize outputTypes
  40                 vector<string> tempOutNames;
  41                 outputTypes["tree"] = tempOutNames;
  42                 outputTypes["summary"] = tempOutNames;
  43         }
  44         catch(exception& e) {
  45                 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
  46                 exit(1);
  47         }
  48 }
  49
  50 //**********************************************************************************************************************
  51 vector<string> IndicatorCommand::getRequiredFiles(){
  52         try {
  53                 vector<string> myArray;
  54                 return myArray;
  55         }
  56         catch(exception& e) {
  57                 m->errorOut(e, "IndicatorCommand", "getRequiredFiles");
  58                 exit(1);
  59         }
  60 }
  61 //**********************************************************************************************************************
  62 IndicatorCommand::IndicatorCommand(string option)  {
  63         try {
  64                 globaldata = GlobalData::getInstance();
  65                 abort = false;
  66
  67                 //allow user to run help
  68                 if(option == "help") { help(); abort = true; }
  69
  70                 else {
  71                         //valid paramters for this command
  72                         string Array[] =  {"tree","shared","relabund","groups","label","outputdir","inputdir"};
  73                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  74
  75                         OptionParser parser(option);
  76                         map<string, string> parameters = parser.getParameters();
  77
  78                         ValidParameters validParameter;
  79                         map<string, string>::iterator it;
  80
  81                         //check to make sure all parameters are valid for command
  82                         for (it = parameters.begin(); it != parameters.end(); it++) {
  83                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  84                         }
  85
  86                         globaldata->newRead();
  87
  88                         vector<string> tempOutNames;
  89                         outputTypes["tree"] = tempOutNames;
  90                         outputTypes["summary"] = tempOutNames;
  91
  92                         //if the user changes the input directory command factory will send this info to us in the output parameter
  93                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  94                         if (inputDir == "not found"){   inputDir = "";          }
  95                         else {
  96                                 string path;
  97                                 it = parameters.find("tree");
  98                                 //user has given a template file
  99                                 if(it != parameters.end()){
 100                                         path = m->hasPath(it->second);
 101                                         //if the user has not given a path then, add inputdir. else leave path alone.
 102                                         if (path == "") {       parameters["tree"] = inputDir + it->second;             }
 103                                 }
 104
 105                                 it = parameters.find("shared");
 106                                 //user has given a template file
 107                                 if(it != parameters.end()){
 108                                         path = m->hasPath(it->second);
 109                                         //if the user has not given a path then, add inputdir. else leave path alone.
 110                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 111                                 }
 112
 113                                 it = parameters.find("relabund");
 114                                 //user has given a template file
 115                                 if(it != parameters.end()){
 116                                         path = m->hasPath(it->second);
 117                                         //if the user has not given a path then, add inputdir. else leave path alone.
 118                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
 119                                 }
 120
 121                         }
 122
 123                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 124
 125                         //check for required parameters
 126                         treefile = validParameter.validFile(parameters, "tree", true);
 127                         if (treefile == "not open") { abort = true; }
 128                         else if (treefile == "not found") { treefile = ""; m->mothurOut("tree is a required parameter for the indicator command."); m->mothurOutEndLine(); abort = true;  }
 129                         else {  globaldata->setTreeFile(treefile);  globaldata->setFormat("tree");      }
 130
 131                         sharedfile = validParameter.validFile(parameters, "shared", true);
 132                         if (sharedfile == "not open") { abort = true; }
 133                         else if (sharedfile == "not found") { sharedfile = ""; }
 134                         else { inputFileName = sharedfile; }
 135
 136                         relabundfile = validParameter.validFile(parameters, "relabund", true);
 137                         if (relabundfile == "not open") { abort = true; }
 138                         else if (relabundfile == "not found") { relabundfile = ""; }
 139                         else { inputFileName = relabundfile; }
 140
 141                         groups = validParameter.validFile(parameters, "groups", false);
 142                         if (groups == "not found") { groups = "";  Groups.push_back("all"); }
 143                         else { m->splitAtDash(groups, Groups);  }
 144                         globaldata->Groups = Groups;
 145
 146                         label = validParameter.validFile(parameters, "label", false);
 147                         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=""; }
 148
 149                         if ((relabundfile == "") && (sharedfile == "")) { m->mothurOut("You must provide either a shared or relabund file."); m->mothurOutEndLine(); abort = true;  }
 150
 151                         if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true;  }
 152
 153                 }
 154         }
 155         catch(exception& e) {
 156                 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
 157                 exit(1);
 158         }
 159 }
 160 //**********************************************************************************************************************
 161
 162 void IndicatorCommand::help(){
 163         try {
 164                 m->mothurOut("The indicator command reads a shared or relabund file and a tree file, and outputs a .indicator.tre and .indicator.summary file. \n");
 165                 m->mothurOut("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");
 166                 m->mothurOut("The summary file lists the indicator value for each OTU for each node.\n");
 167                 m->mothurOut("The indicator command parameters are tree, groups, shared, relabund and label. The tree parameter is required as well as either shared or relabund.\n");
 168                 m->mothurOut("The groups parameter allows you to specify which of the groups in your shared or relabund you would like analyzed.  The groups may be entered separated by dashes.\n");
 169                 m->mothurOut("The label parameter indicates at what distance your tree relates to the shared or relabund.\n");
 170                 m->mothurOut("The indicator command should be used in the following format: indicator(tree=test.tre, shared=test.shared, label=0.03)\n");
 171                 m->mothurOut("Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreefile).\n\n");
 172         }
 173         catch(exception& e) {
 174                 m->errorOut(e, "IndicatorCommand", "help");
 175                 exit(1);
 176         }
 177 }
 178
 179 //**********************************************************************************************************************
 180
 181 IndicatorCommand::~IndicatorCommand(){}
 182
 183 //**********************************************************************************************************************
 184
 185 int IndicatorCommand::execute(){
 186         try {
 187
 188                 if (abort == true) { return 0; }
 189
 190                 /***************************************************/
 191                 // use smart distancing to get right sharedRabund  //
 192                 /***************************************************/
 193                 if (sharedfile != "") {
 194                         getShared();
 195                         if (m->control_pressed) {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } return 0; }
 196                         if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
 197                 }else {
 198                         getSharedFloat();
 199                         if (m->control_pressed) {  for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } return 0; }
 200                         if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
 201                 }
 202
 203                 /***************************************************/
 204                 //    reading tree info                                                    //
 205                 /***************************************************/
 206                 string groupfile = "";
 207                 Tree* tree = new Tree(treefile); delete tree;  //extracts names from tree to make faked out groupmap
 208
 209                 globaldata->setGroupFile(groupfile);
 210                 treeMap = new TreeMap();
 211                 bool mismatch = false;
 212                 for (int i = 0; i < globaldata->Treenames.size(); i++) {
 213                         //sanity check - is this a group that is not in the sharedfile?
 214                         if (!(m->inUsersGroups(globaldata->Treenames[i], globaldata->gGroupmap->namesOfGroups))) {
 215                                 m->mothurOut("[ERROR]: " + globaldata->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
 216                                 mismatch = true;
 217                         }
 218                         treeMap->addSeq(globaldata->Treenames[i], "Group1");
 219                 }
 220
 221                 if (mismatch) { //cleanup and exit
 222                         if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 223                         else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 224                         delete treeMap;
 225                         return 0;
 226                 }
 227
 228                 globaldata->gTreemap = treeMap;
 229
 230                 read = new ReadNewickTree(treefile);
 231                 int readOk = read->read();
 232
 233                 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); globaldata->gTree.clear(); delete globaldata->gTreemap; delete read; return 0; }
 234
 235                 vector<Tree*> T = globaldata->gTree;
 236
 237                 delete read;
 238
 239                 if (m->control_pressed) {
 240                         if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 241                         else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 242                         for (int i = 0; i < T.size(); i++) {  delete T[i];  } globaldata->gTree.clear(); delete globaldata->gTreemap; return 0;
 243                 }
 244
 245                 T[0]->assembleTree();
 246
 247                 /***************************************************/
 248                 //    create ouptut tree - respecting pickedGroups //
 249                 /***************************************************/
 250                 Tree* outputTree = new Tree(globaldata->Groups.size());
 251
 252                 outputTree->getSubTree(T[0], globaldata->Groups);
 253                 outputTree->assembleTree();
 254
 255                 //no longer need original tree, we have output tree to use and label
 256                 for (int i = 0; i < T.size(); i++) {  delete T[i];  } globaldata->gTree.clear();
 257
 258
 259                 if (m->control_pressed) {
 260                         if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 261                         else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 262                         delete outputTree; delete globaldata->gTreemap;  return 0;
 263                 }
 264
 265                 /***************************************************/
 266                 //              get indicator species values                       //
 267                 /***************************************************/
 268                 GetIndicatorSpecies(outputTree);
 269
 270                 if (m->control_pressed) {
 271                         if (sharedfile != "") {  for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } }
 272                         else { for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } }
 273                         for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); }
 274                         delete outputTree; delete globaldata->gTreemap;  return 0;
 275                 }
 276
 277                 m->mothurOutEndLine();
 278                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 279                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 280                 m->mothurOutEndLine();
 281
 282                 return 0;
 283         }
 284         catch(exception& e) {
 285                 m->errorOut(e, "IndicatorCommand", "execute");
 286                 exit(1);
 287         }
 288 }
 289 //**********************************************************************************************************************
 290 //traverse tree finding indicator species values for each otu at each node
 291 //label node with otu number that has highest indicator value
 292 //report all otu values to file
 293 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
 294         try {
 295
 296                 string thisOutputDir = outputDir;
 297                 if (outputDir == "") {  thisOutputDir += m->hasPath(inputFileName);  }
 298                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
 299                 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
 300
 301                 ofstream out;
 302                 m->openOutputFile(outputFileName, out);
 303                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 304                 out << "Node\tOTU#\tIndVal" << endl;
 305
 306                 string treeOutputDir = outputDir;
 307                 if (outputDir == "") {  treeOutputDir += m->hasPath(treefile);  }
 308                 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
 309
 310
 311                 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
 312                 map<int, set<string> > nodeToDescendants;
 313                 map<int, set<int> > descendantNodes;
 314                 for (int i = 0; i < T->getNumNodes(); i++) {
 315                         if (m->control_pressed) { return 0; }
 316
 317                         nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
 318                 }
 319
 320                 //you need the distances to leaf to decide grouping below
 321                 //this will also set branch lengths if the tree does not include them
 322                 map<int, float> distToLeaf = getLengthToLeaf(T);
 323
 324                 //for each node
 325                 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
 326
 327                         if (m->control_pressed) { out.close(); return 0; }
 328
 329                         /*****************************************************/
 330                         //create vectors containing rabund info                          //
 331                         /*****************************************************/
 332
 333                         vector<float> indicatorValues; //size of numBins
 334
 335                         if (sharedfile != "") {
 336                                 vector< vector<SharedRAbundVector*> > groupings;
 337
 338                                 /*groupings.resize(1);
 339                                 groupings[0].push_back(lookup[0]);
 340                                 groupings[0].push_back(lookup[1]);
 341                                 groupings[0].push_back(lookup[2]);
 342                                 groupings[0].push_back(lookup[3]);
 343                                 groupings[0].push_back(lookup[4]);*/
 344
 345                                 //get nodes that will be a valid grouping
 346                                 //you are valid if you are not one of my descendants
 347                                 //AND your distToLeaf is <= mine
 348                                 //AND your distToLeaf is >= my smallest childs
 349                                 //AND you were not added as part of a larger groupings
 350
 351                                 set<string> groupsAlreadyAdded;
 352                                 //create a grouping with my grouping
 353                                 vector<SharedRAbundVector*> subset;
 354                                 int count = 0;
 355                                 int doneCount = nodeToDescendants[i].size();
 356                                 for (int k = 0; k < lookup.size(); k++) {
 357                                         //is this descendant of i
 358                                         if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
 359                                                 subset.push_back(lookup[k]);
 360                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 361                                                 count++;
 362                                         }
 363                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 364                                 }
 365                                 if (subset.size() != 0) { groupings.push_back(subset); }
 366
 367
 368                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 369                                         if ((descendantNodes[i].count(j) == 0) && (distToLeaf[j] <= distToLeaf[i]) && ((distToLeaf[j] >= distToLeaf[T->tree[i].getLChild()]) || (distToLeaf[j] >= distToLeaf[T->tree[i].getRChild()]))) {
 370                                                 vector<SharedRAbundVector*> subset;
 371                                                 int count = 0;
 372                                                 int doneCount = nodeToDescendants[j].size();
 373                                                 for (int k = 0; k < lookup.size(); k++) {
 374                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 375                                                         if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
 376                                                                 subset.push_back(lookup[k]);
 377                                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 378                                                                 count++;
 379                                                         }
 380                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 381                                                 }
 382
 383                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 384                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 385                                         }
 386                                 }
 387
 388                                 if (groupsAlreadyAdded.size() != lookup.size()) { cout << i << '\t' << groupsAlreadyAdded.size() << '\t' << lookup.size() << endl; m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 389                                 for (int k = 0; k < lookup.size(); k++) {
 390                                         if (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0) { cout << lookup[k]->getGroup() << endl; }
 391                                 }
 392
 393                                 indicatorValues = getValues(groupings);
 394
 395                         }else {
 396                                 vector< vector<SharedRAbundFloatVector*> > groupings;
 397
 398                                 //get nodes that will be a valid grouping
 399                                 //you are valid if you are not one of my descendants
 400                                 //AND your distToLeaf is <= mine
 401                                 //AND your distToLeaf is >= my smallest childs
 402                                 //AND you were not added as part of a larger grouping
 403
 404                                 set<string> groupsAlreadyAdded;
 405                                 //create a grouping with my grouping
 406                                 vector<SharedRAbundFloatVector*> subset;
 407                                 int count = 0;
 408                                 int doneCount = nodeToDescendants[i].size();
 409                                 for (int k = 0; k < lookupFloat.size(); k++) {
 410                                         //is this descendant of i
 411                                         if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
 412                                                 subset.push_back(lookupFloat[k]);
 413                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 414                                                 count++;
 415                                         }
 416                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 417                                 }
 418                                 if (subset.size() != 0) { groupings.push_back(subset); }
 419
 420                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 421                                         if ((descendantNodes[i].count(j) == 0) && (distToLeaf[j] <= distToLeaf[i]) && ((distToLeaf[j] >= distToLeaf[T->tree[i].getLChild()]) || (distToLeaf[j] >= distToLeaf[T->tree[i].getRChild()]))) {
 422                                                 vector<SharedRAbundFloatVector*> subset;
 423                                                 int count = 0;
 424                                                 int doneCount = nodeToDescendants[j].size();
 425                                                 for (int k = 0; k < lookupFloat.size(); k++) {
 426                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 427                                                         if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
 428                                                                 subset.push_back(lookupFloat[k]);
 429                                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 430                                                                 count++;
 431                                                         }
 432                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 433                                                 }
 434
 435                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 436                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 437                                         }
 438                                 }
 439
 440                                 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 441
 442                                 indicatorValues = getValues(groupings);
 443                         }
 444
 445                         if (m->control_pressed) { out.close(); return 0; }
 446
 447                         /******************************************************/
 448                         //output indicator values to table form + label tree  //
 449                         /*****************************************************/
 450                         for (int j = 0; j < indicatorValues.size(); j++) {
 451
 452                                 if (m->control_pressed) { out.close(); return 0; }
 453
 454                                 out << (i+1) << '\t' << (j+1) << '\t' << indicatorValues[j] << endl;
 455                         }
 456
 457                         T->tree[i].setLabel((i+1));
 458
 459                 }
 460                 out.close();
 461
 462                 ofstream outTree;
 463                 m->openOutputFile(outputTreeFileName, outTree);
 464                 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
 465
 466                 T->print(outTree, "both");
 467                 outTree.close();
 468
 469                 return 0;
 470         }
 471         catch(exception& e) {
 472                 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
 473                 exit(1);
 474         }
 475 }
 476 //**********************************************************************************************************************
 477 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings){
 478         try {
 479                 vector<float> values;
 480
 481                 //for each otu
 482                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 483
 484                         if (m->control_pressed) { return values; }
 485
 486                         vector<float> terms;
 487                         float AijDenominator = 0.0;
 488                         vector<float> Bij;
 489                         //get overall abundance of each grouping
 490                         for (int j = 0; j < groupings.size(); j++) {
 491
 492                                 float totalAbund = 0;
 493                                 int numNotZero = 0;
 494                                 for (int k = 0; k < groupings[j].size(); k++) {
 495                                         totalAbund += groupings[j][k]->getAbundance(i);
 496                                         if (groupings[j][k]->getAbundance(i) != 0) { numNotZero++; }
 497                                 }
 498
 499                                 float Aij = (totalAbund / (float) groupings[j].size());
 500                                 terms.push_back(Aij);
 501
 502                                 //percentage of sites represented
 503                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 504
 505                                 AijDenominator += Aij;
 506                         }
 507
 508                         float maxIndVal = 0.0;
 509                         for (int j = 0; j < terms.size(); j++) {
 510                                 float thisAij = (terms[j] / AijDenominator);
 511                                 float thisValue = thisAij * Bij[j] * 100.0;
 512
 513                                 //save largest
 514                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 515                         }
 516
 517                         values.push_back(maxIndVal);
 518                 }
 519
 520                 return values;
 521         }
 522         catch(exception& e) {
 523                 m->errorOut(e, "IndicatorCommand", "getValues");
 524                 exit(1);
 525         }
 526 }
 527 //**********************************************************************************************************************
 528 //same as above, just data type difference
 529 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings){
 530         try {
 531                 vector<float> values;
 532
 533         /*for (int j = 0; j < groupings.size(); j++) {
 534                 cout << "grouping " << j << endl;
 535                 for (int k = 0; k < groupings[j].size(); k++) {
 536                         cout << groupings[j][k]->getGroup() << endl;
 537                 }
 538         }*/
 539                 //for each otu
 540                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 541                         vector<float> terms;
 542                         float AijDenominator = 0.0;
 543                         vector<float> Bij;
 544                         //get overall abundance of each grouping
 545                         for (int j = 0; j < groupings.size(); j++) {
 546
 547                                 int totalAbund = 0.0;
 548                                 int numNotZero = 0;
 549                                 for (int k = 0; k < groupings[j].size(); k++) {
 550                                         totalAbund += groupings[j][k]->getAbundance(i);
 551                                         if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
 552
 553                                 }
 554
 555                                 float Aij = (totalAbund / (float) groupings[j].size());
 556                                 terms.push_back(Aij);
 557
 558                                 //percentage of sites represented
 559                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 560
 561                                 AijDenominator += Aij;
 562                         }
 563
 564                         float maxIndVal = 0.0;
 565                         for (int j = 0; j < terms.size(); j++) {
 566                                 float thisAij = (terms[j] / AijDenominator);
 567                                 float thisValue = thisAij * Bij[j] * 100.0;
 568
 569                                 //save largest
 570                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 571                         }
 572
 573                         values.push_back(maxIndVal);
 574                 }
 575
 576                 return values;
 577         }
 578         catch(exception& e) {
 579                 m->errorOut(e, "IndicatorCommand", "getValues");
 580                 exit(1);
 581         }
 582 }
 583 //**********************************************************************************************************************
 584 //you need the distances to leaf to decide groupings
 585 //this will also set branch lengths if the tree does not include them
 586 map<int, float> IndicatorCommand::getLengthToLeaf(Tree*& T){
 587         try {
 588                 map<int, float> dists;
 589
 590                 for (int i = 0; i < T->getNumNodes(); i++) {
 591
 592                         int lc = T->tree[i].getLChild();
 593                         int rc = T->tree[i].getRChild();
 594
 595                         //if you have no branch length, set it then calc
 596                         if (T->tree[i].getBranchLength() <= 0.0) {
 597
 598                                 if (lc == -1) { // you are a leaf
 599                                         //if you are a leaf set you priliminary length to 1.0, this may adjust later
 600                                         T->tree[i].setBranchLength(1.0);
 601                                         dists[i] = 1.0;
 602                                 }else{ // you are an internal node
 603                                         //look at your children's length to leaf
 604                                         float ldist = dists[lc];
 605                                         float rdist = dists[rc];
 606
 607                                         float greater = ldist;
 608                                         if (rdist > greater) { greater = rdist; dists[i] = ldist; }
 609                                         else { dists[i] = rdist; }
 610
 611                                         //branch length = difference + 1
 612                                         T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
 613                                         T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
 614                                 }
 615
 616                         }else{
 617                                 if (lc == -1) { dists[i] = T->tree[i].getBranchLength(); }
 618                                 else { //smaller of my two children distances plus my branch length
 619                                         //look at your children's length to leaf
 620                                         float ldist = dists[lc];
 621                                         float rdist = dists[rc];
 622
 623                                         float smaller = ldist;
 624                                         if (rdist < smaller) { smaller = rdist;  }
 625
 626                                         dists[i] = smaller + T->tree[i].getBranchLength();
 627                                 }
 628                         }
 629
 630                 }
 631
 632                 return dists;
 633         }
 634         catch(exception& e) {
 635                 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
 636                 exit(1);
 637         }
 638 }
 639 //**********************************************************************************************************************
 640 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
 641         try {
 642                 set<string> names;
 643
 644                 set<string>::iterator it;
 645
 646                 int lc = T->tree[i].getLChild();
 647                 int rc = T->tree[i].getRChild();
 648
 649                 if (lc == -1) { //you are a leaf your only descendant is yourself
 650                         set<int> temp; temp.insert(i);
 651                         names.insert(T->tree[i].getName());
 652                         nodes[i] = temp;
 653                 }else{ //your descedants are the combination of your childrens descendants
 654                         names = descendants[lc];
 655                         nodes[i] = nodes[lc];
 656                         for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
 657                                 names.insert(*it);
 658                         }
 659                         for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
 660                                 nodes[i].insert(*itNum);
 661                         }
 662                         //you are your own descendant
 663                         nodes[i].insert(i);
 664                 }
 665
 666                 return names;
 667         }
 668         catch(exception& e) {
 669                 m->errorOut(e, "IndicatorCommand", "getDescendantList");
 670                 exit(1);
 671         }
 672 }
 673 //**********************************************************************************************************************
 674 int IndicatorCommand::getShared(){
 675         try {
 676                 InputData* input = new InputData(sharedfile, "sharedfile");
 677                 lookup = input->getSharedRAbundVectors();
 678                 string lastLabel = lookup[0]->getLabel();
 679
 680                 if (label == "") { label = lastLabel; delete input; return 0; }
 681
 682                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 683                 set<string> labels; labels.insert(label);
 684                 set<string> processedLabels;
 685                 set<string> userLabels = labels;
 686
 687                 //as long as you are not at the end of the file or done wih the lines you want
 688                 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
 689                         if (m->control_pressed) {  delete input; return 0;  }
 690
 691                         if(labels.count(lookup[0]->getLabel()) == 1){
 692                                 processedLabels.insert(lookup[0]->getLabel());
 693                                 userLabels.erase(lookup[0]->getLabel());
 694                                 break;
 695                         }
 696
 697                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 698                                 string saveLabel = lookup[0]->getLabel();
 699
 700                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 701                                 lookup = input->getSharedRAbundVectors(lastLabel);
 702
 703                                 processedLabels.insert(lookup[0]->getLabel());
 704                                 userLabels.erase(lookup[0]->getLabel());
 705
 706                                 //restore real lastlabel to save below
 707                                 lookup[0]->setLabel(saveLabel);
 708                                 break;
 709                         }
 710
 711                         lastLabel = lookup[0]->getLabel();
 712
 713                         //get next line to process
 714                         //prevent memory leak
 715                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 716                         lookup = input->getSharedRAbundVectors();
 717                 }
 718
 719
 720                 if (m->control_pressed) { delete input; return 0;  }
 721
 722                 //output error messages about any remaining user labels
 723                 set<string>::iterator it;
 724                 bool needToRun = false;
 725                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 726                         m->mothurOut("Your file does not include the label " + *it);
 727                         if (processedLabels.count(lastLabel) != 1) {
 728                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 729                                 needToRun = true;
 730                         }else {
 731                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 732                         }
 733                 }
 734
 735                 //run last label if you need to
 736                 if (needToRun == true)  {
 737                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
 738                         lookup = input->getSharedRAbundVectors(lastLabel);
 739                 }
 740
 741                 delete input;
 742                 return 0;
 743         }
 744         catch(exception& e) {
 745                 m->errorOut(e, "IndicatorCommand", "getShared");
 746                 exit(1);
 747         }
 748 }
 749 //**********************************************************************************************************************
 750 int IndicatorCommand::getSharedFloat(){
 751         try {
 752                 InputData* input = new InputData(relabundfile, "relabund");
 753                 lookupFloat = input->getSharedRAbundFloatVectors();
 754                 string lastLabel = lookupFloat[0]->getLabel();
 755
 756                 if (label == "") { label = lastLabel; delete input; return 0; }
 757
 758                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 759                 set<string> labels; labels.insert(label);
 760                 set<string> processedLabels;
 761                 set<string> userLabels = labels;
 762
 763                 //as long as you are not at the end of the file or done wih the lines you want
 764                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
 765
 766                         if (m->control_pressed) {  delete input; return 0;  }
 767
 768                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
 769                                 processedLabels.insert(lookupFloat[0]->getLabel());
 770                                 userLabels.erase(lookupFloat[0]->getLabel());
 771                                 break;
 772                         }
 773
 774                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 775                                 string saveLabel = lookupFloat[0]->getLabel();
 776
 777                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 778                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 779
 780                                 processedLabels.insert(lookupFloat[0]->getLabel());
 781                                 userLabels.erase(lookupFloat[0]->getLabel());
 782
 783                                 //restore real lastlabel to save below
 784                                 lookupFloat[0]->setLabel(saveLabel);
 785                                 break;
 786                         }
 787
 788                         lastLabel = lookupFloat[0]->getLabel();
 789
 790                         //get next line to process
 791                         //prevent memory leak
 792                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 793                         lookupFloat = input->getSharedRAbundFloatVectors();
 794                 }
 795
 796
 797                 if (m->control_pressed) { delete input; return 0;  }
 798
 799                 //output error messages about any remaining user labels
 800                 set<string>::iterator it;
 801                 bool needToRun = false;
 802                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 803                         m->mothurOut("Your file does not include the label " + *it);
 804                         if (processedLabels.count(lastLabel) != 1) {
 805                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 806                                 needToRun = true;
 807                         }else {
 808                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 809                         }
 810                 }
 811
 812                 //run last label if you need to
 813                 if (needToRun == true)  {
 814                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }
 815                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 816                 }
 817
 818                 delete input;
 819                 return 0;
 820         }
 821         catch(exception& e) {
 822                 m->errorOut(e, "IndicatorCommand", "getShared");
 823         exit(1);
 824         }
 825 }
 826 /*****************************************************************/
 827
 828