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