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
 305                 string treeOutputDir = outputDir;
 306                 if (outputDir == "") {  treeOutputDir += m->hasPath(treefile);  }
 307                 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
 308
 309
 310                 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
 311                 map<int, set<string> > nodeToDescendants;
 312                 map<int, set<int> > descendantNodes;
 313                 for (int i = 0; i < T->getNumNodes(); i++) {
 314                         if (m->control_pressed) { return 0; }
 315
 316                         nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
 317                 }
 318
 319                 //you need the distances to leaf to decide grouping below
 320                 //this will also set branch lengths if the tree does not include them
 321                 map<int, float> distToRoot = getDistToRoot(T);
 322
 323                 //for each node
 324                 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
 325
 326                         if (m->control_pressed) { out.close(); return 0; }
 327
 328                         /*****************************************************/
 329                         //create vectors containing rabund info                          //
 330                         /*****************************************************/
 331
 332                         vector<float> indicatorValues; //size of numBins
 333
 334                         if (sharedfile != "") {
 335                                 vector< vector<SharedRAbundVector*> > groupings;
 336
 337                                 //get nodes that will be a valid grouping
 338                                 //you are valid if you are not one of my descendants
 339                                 //AND your distToRoot is >= mine
 340                                 //AND you were not added as part of a larger grouping. Largest nodes are added first.
 341
 342                                 set<string> groupsAlreadyAdded;
 343                                 //create a grouping with my grouping
 344                                 vector<SharedRAbundVector*> subset;
 345                                 int count = 0;
 346                                 int doneCount = nodeToDescendants[i].size();
 347                                 for (int k = 0; k < lookup.size(); k++) {
 348                                         //is this descendant of i
 349                                         if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
 350                                                 subset.push_back(lookup[k]);
 351                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 352                                                 count++;
 353                                         }
 354                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 355                                 }
 356                                 if (subset.size() != 0) { groupings.push_back(subset); }
 357
 358
 359                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 360
 361
 362                                         if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
 363                                                 vector<SharedRAbundVector*> subset;
 364                                                 int count = 0;
 365                                                 int doneCount = nodeToDescendants[j].size();
 366                                                 for (int k = 0; k < lookup.size(); k++) {
 367                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 368                                                         if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
 369                                                                 subset.push_back(lookup[k]);
 370                                                                 groupsAlreadyAdded.insert(lookup[k]->getGroup());
 371                                                                 count++;
 372                                                         }
 373                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 374                                                 }
 375
 376                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 377                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 378                                         }
 379                                 }
 380
 381                                 if (groupsAlreadyAdded.size() != lookup.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 382                                 //for (int k = 0; k < lookup.size(); k++) {
 383                                 //      if (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0) { cout << lookup[k]->getGroup() << endl; }
 384                                 //}
 385
 386                                 indicatorValues = getValues(groupings);
 387
 388                         }else {
 389                                 vector< vector<SharedRAbundFloatVector*> > groupings;
 390
 391                                 //get nodes that will be a valid grouping
 392                                 //you are valid if you are not one of my descendants
 393                                 //AND your distToRoot is >= mine
 394                                 //AND you were not added as part of a larger grouping. Largest nodes are added first.
 395
 396                                 set<string> groupsAlreadyAdded;
 397                                 //create a grouping with my grouping
 398                                 vector<SharedRAbundFloatVector*> subset;
 399                                 int count = 0;
 400                                 int doneCount = nodeToDescendants[i].size();
 401                                 for (int k = 0; k < lookupFloat.size(); k++) {
 402                                         //is this descendant of i
 403                                         if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
 404                                                 subset.push_back(lookupFloat[k]);
 405                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 406                                                 count++;
 407                                         }
 408                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 409                                 }
 410                                 if (subset.size() != 0) { groupings.push_back(subset); }
 411
 412                                 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
 413                                         if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
 414                                                 vector<SharedRAbundFloatVector*> subset;
 415                                                 int count = 0;
 416                                                 int doneCount = nodeToDescendants[j].size();
 417                                                 for (int k = 0; k < lookupFloat.size(); k++) {
 418                                                         //is this descendant of j, and we didn't already add this as part of a larger grouping
 419                                                         if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
 420                                                                 subset.push_back(lookupFloat[k]);
 421                                                                 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
 422                                                                 count++;
 423                                                         }
 424                                                         if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
 425                                                 }
 426
 427                                                 //if subset.size == 0 then the node was added as part of a larger grouping
 428                                                 if (subset.size() != 0) { groupings.push_back(subset); }
 429                                         }
 430                                 }
 431
 432                                 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
 433
 434                                 indicatorValues = getValues(groupings);
 435                         }
 436
 437                         if (m->control_pressed) { out.close(); return 0; }
 438
 439
 440                         /******************************************************/
 441                         //output indicator values to table form + label tree  //
 442                         /*****************************************************/
 443                         out << (i+1) << '\t';
 444                         for (int j = 0; j < indicatorValues.size(); j++) {
 445
 446                                 if (m->control_pressed) { out.close(); return 0; }
 447
 448                                 out << indicatorValues[j] << '\t';
 449                         }
 450                         out << endl;
 451
 452                         T->tree[i].setLabel((i+1));
 453
 454                 }
 455                 out.close();
 456
 457                 ofstream outTree;
 458                 m->openOutputFile(outputTreeFileName, outTree);
 459                 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
 460
 461                 T->print(outTree, "both");
 462                 outTree.close();
 463
 464                 return 0;
 465         }
 466         catch(exception& e) {
 467                 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
 468                 exit(1);
 469         }
 470 }
 471 //**********************************************************************************************************************
 472 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings){
 473         try {
 474                 vector<float> values;
 475
 476                 //for each otu
 477                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 478
 479                         if (m->control_pressed) { return values; }
 480
 481                         vector<float> terms;
 482                         float AijDenominator = 0.0;
 483                         vector<float> Bij;
 484                         //get overall abundance of each grouping
 485                         for (int j = 0; j < groupings.size(); j++) {
 486
 487                                 float totalAbund = 0;
 488                                 int numNotZero = 0;
 489                                 for (int k = 0; k < groupings[j].size(); k++) {
 490                                         totalAbund += groupings[j][k]->getAbundance(i);
 491                                         if (groupings[j][k]->getAbundance(i) != 0) { numNotZero++; }
 492                                 }
 493
 494                                 float Aij = (totalAbund / (float) groupings[j].size());
 495                                 terms.push_back(Aij);
 496
 497                                 //percentage of sites represented
 498                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 499
 500                                 AijDenominator += Aij;
 501                         }
 502
 503                         float maxIndVal = 0.0;
 504                         for (int j = 0; j < terms.size(); j++) {
 505                                 float thisAij = (terms[j] / AijDenominator);
 506                                 float thisValue = thisAij * Bij[j] * 100.0;
 507
 508                                 //save largest
 509                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 510                         }
 511
 512                         values.push_back(maxIndVal);
 513                 }
 514
 515                 return values;
 516         }
 517         catch(exception& e) {
 518                 m->errorOut(e, "IndicatorCommand", "getValues");
 519                 exit(1);
 520         }
 521 }
 522 //**********************************************************************************************************************
 523 //same as above, just data type difference
 524 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings){
 525         try {
 526                 vector<float> values;
 527
 528         /*for (int j = 0; j < groupings.size(); j++) {
 529                 cout << "grouping " << j << endl;
 530                 for (int k = 0; k < groupings[j].size(); k++) {
 531                         cout << groupings[j][k]->getGroup() << endl;
 532                 }
 533         }*/
 534                 //for each otu
 535                 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
 536                         vector<float> terms;
 537                         float AijDenominator = 0.0;
 538                         vector<float> Bij;
 539                         //get overall abundance of each grouping
 540                         for (int j = 0; j < groupings.size(); j++) {
 541
 542                                 int totalAbund = 0.0;
 543                                 int numNotZero = 0;
 544                                 for (int k = 0; k < groupings[j].size(); k++) {
 545                                         totalAbund += groupings[j][k]->getAbundance(i);
 546                                         if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
 547
 548                                 }
 549
 550                                 float Aij = (totalAbund / (float) groupings[j].size());
 551                                 terms.push_back(Aij);
 552
 553                                 //percentage of sites represented
 554                                 Bij.push_back(numNotZero / (float) groupings[j].size());
 555
 556                                 AijDenominator += Aij;
 557                         }
 558
 559                         float maxIndVal = 0.0;
 560                         for (int j = 0; j < terms.size(); j++) {
 561                                 float thisAij = (terms[j] / AijDenominator);
 562                                 float thisValue = thisAij * Bij[j] * 100.0;
 563
 564                                 //save largest
 565                                 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 566                         }
 567
 568                         values.push_back(maxIndVal);
 569                 }
 570
 571                 return values;
 572         }
 573         catch(exception& e) {
 574                 m->errorOut(e, "IndicatorCommand", "getValues");
 575                 exit(1);
 576         }
 577 }
 578 //**********************************************************************************************************************
 579 //you need the distances to root to decide groupings
 580 //this will also set branch lengths if the tree does not include them
 581 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
 582         try {
 583                 map<int, float> dists;
 584
 585                 bool hasBranchLengths = false;
 586                 for (int i = 0; i < T->getNumNodes(); i++) {
 587                         if (T->tree[i].getBranchLength() > 0.0) {  hasBranchLengths = true; break; }
 588                 }
 589
 590                 //set branchlengths if needed
 591                 if (!hasBranchLengths) {
 592                         for (int i = 0; i < T->getNumNodes(); i++) {
 593
 594                                 int lc = T->tree[i].getLChild();
 595                                 int rc = T->tree[i].getRChild();
 596
 597                                 if (lc == -1) { // you are a leaf
 598                                         //if you are a leaf set you priliminary length to 1.0, this may adjust later
 599                                         T->tree[i].setBranchLength(1.0);
 600                                         dists[i] = 1.0;
 601                                 }else{ // you are an internal node
 602                                         //look at your children's length to leaf
 603                                         float ldist = dists[lc];
 604                                         float rdist = dists[rc];
 605
 606                                         float greater = ldist;
 607                                         if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
 608                                         else { dists[i] = rdist + 1.0; }
 609
 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                 }
 617
 618                 dists.clear();
 619
 620                 for (int i = 0; i < T->getNumNodes(); i++) {
 621
 622                         double sum = 0.0;
 623                         int index = i;
 624
 625                         while(T->tree[index].getParent() != -1){
 626                                 if (T->tree[index].getBranchLength() != -1) {
 627                                         sum += abs(T->tree[index].getBranchLength());
 628                                 }
 629                                 index = T->tree[index].getParent();
 630                         }
 631
 632                         dists[i] = sum;
 633                 }
 634
 635                 return dists;
 636         }
 637         catch(exception& e) {
 638                 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
 639                 exit(1);
 640         }
 641 }
 642 //**********************************************************************************************************************
 643 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
 644         try {
 645                 set<string> names;
 646
 647                 set<string>::iterator it;
 648
 649                 int lc = T->tree[i].getLChild();
 650                 int rc = T->tree[i].getRChild();
 651
 652                 if (lc == -1) { //you are a leaf your only descendant is yourself
 653                         set<int> temp; temp.insert(i);
 654                         names.insert(T->tree[i].getName());
 655                         nodes[i] = temp;
 656                 }else{ //your descedants are the combination of your childrens descendants
 657                         names = descendants[lc];
 658                         nodes[i] = nodes[lc];
 659                         for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
 660                                 names.insert(*it);
 661                         }
 662                         for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
 663                                 nodes[i].insert(*itNum);
 664                         }
 665                         //you are your own descendant
 666                         nodes[i].insert(i);
 667                 }
 668
 669                 return names;
 670         }
 671         catch(exception& e) {
 672                 m->errorOut(e, "IndicatorCommand", "getDescendantList");
 673                 exit(1);
 674         }
 675 }
 676 //**********************************************************************************************************************
 677 int IndicatorCommand::getShared(){
 678         try {
 679                 InputData* input = new InputData(sharedfile, "sharedfile");
 680                 lookup = input->getSharedRAbundVectors();
 681                 string lastLabel = lookup[0]->getLabel();
 682
 683                 if (label == "") { label = lastLabel; delete input; return 0; }
 684
 685                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 686                 set<string> labels; labels.insert(label);
 687                 set<string> processedLabels;
 688                 set<string> userLabels = labels;
 689
 690                 //as long as you are not at the end of the file or done wih the lines you want
 691                 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
 692                         if (m->control_pressed) {  delete input; return 0;  }
 693
 694                         if(labels.count(lookup[0]->getLabel()) == 1){
 695                                 processedLabels.insert(lookup[0]->getLabel());
 696                                 userLabels.erase(lookup[0]->getLabel());
 697                                 break;
 698                         }
 699
 700                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 701                                 string saveLabel = lookup[0]->getLabel();
 702
 703                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 704                                 lookup = input->getSharedRAbundVectors(lastLabel);
 705
 706                                 processedLabels.insert(lookup[0]->getLabel());
 707                                 userLabels.erase(lookup[0]->getLabel());
 708
 709                                 //restore real lastlabel to save below
 710                                 lookup[0]->setLabel(saveLabel);
 711                                 break;
 712                         }
 713
 714                         lastLabel = lookup[0]->getLabel();
 715
 716                         //get next line to process
 717                         //prevent memory leak
 718                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 719                         lookup = input->getSharedRAbundVectors();
 720                 }
 721
 722
 723                 if (m->control_pressed) { delete input; return 0;  }
 724
 725                 //output error messages about any remaining user labels
 726                 set<string>::iterator it;
 727                 bool needToRun = false;
 728                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 729                         m->mothurOut("Your file does not include the label " + *it);
 730                         if (processedLabels.count(lastLabel) != 1) {
 731                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 732                                 needToRun = true;
 733                         }else {
 734                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 735                         }
 736                 }
 737
 738                 //run last label if you need to
 739                 if (needToRun == true)  {
 740                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } }
 741                         lookup = input->getSharedRAbundVectors(lastLabel);
 742                 }
 743
 744                 delete input;
 745                 return 0;
 746         }
 747         catch(exception& e) {
 748                 m->errorOut(e, "IndicatorCommand", "getShared");
 749                 exit(1);
 750         }
 751 }
 752 //**********************************************************************************************************************
 753 int IndicatorCommand::getSharedFloat(){
 754         try {
 755                 InputData* input = new InputData(relabundfile, "relabund");
 756                 lookupFloat = input->getSharedRAbundFloatVectors();
 757                 string lastLabel = lookupFloat[0]->getLabel();
 758
 759                 if (label == "") { label = lastLabel; delete input; return 0; }
 760
 761                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 762                 set<string> labels; labels.insert(label);
 763                 set<string> processedLabels;
 764                 set<string> userLabels = labels;
 765
 766                 //as long as you are not at the end of the file or done wih the lines you want
 767                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
 768
 769                         if (m->control_pressed) {  delete input; return 0;  }
 770
 771                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
 772                                 processedLabels.insert(lookupFloat[0]->getLabel());
 773                                 userLabels.erase(lookupFloat[0]->getLabel());
 774                                 break;
 775                         }
 776
 777                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 778                                 string saveLabel = lookupFloat[0]->getLabel();
 779
 780                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 781                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 782
 783                                 processedLabels.insert(lookupFloat[0]->getLabel());
 784                                 userLabels.erase(lookupFloat[0]->getLabel());
 785
 786                                 //restore real lastlabel to save below
 787                                 lookupFloat[0]->setLabel(saveLabel);
 788                                 break;
 789                         }
 790
 791                         lastLabel = lookupFloat[0]->getLabel();
 792
 793                         //get next line to process
 794                         //prevent memory leak
 795                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }
 796                         lookupFloat = input->getSharedRAbundFloatVectors();
 797                 }
 798
 799
 800                 if (m->control_pressed) { delete input; return 0;  }
 801
 802                 //output error messages about any remaining user labels
 803                 set<string>::iterator it;
 804                 bool needToRun = false;
 805                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 806                         m->mothurOut("Your file does not include the label " + *it);
 807                         if (processedLabels.count(lastLabel) != 1) {
 808                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 809                                 needToRun = true;
 810                         }else {
 811                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 812                         }
 813                 }
 814
 815                 //run last label if you need to
 816                 if (needToRun == true)  {
 817                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }
 818                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 819                 }
 820
 821                 delete input;
 822                 return 0;
 823         }
 824         catch(exception& e) {
 825                 m->errorOut(e, "IndicatorCommand", "getShared");
 826         exit(1);
 827         }
 828 }
 829 /*****************************************************************/
 830
 831