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