5 * Created by westcott on 11/12/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "indicatorcommand.h"
11 #include "sharedutilities.h"
14 //**********************************************************************************************************************
15 vector<string> IndicatorCommand::setParameters(){
17 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
18 CommandParameter pdesign("design", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(pdesign);
19 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(pshared);
20 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(prelabund);
21 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
22 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
23 CommandParameter ptree("tree", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(ptree);
24 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
25 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
26 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
28 vector<string> myArray;
29 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
33 m->errorOut(e, "IndicatorCommand", "setParameters");
37 //**********************************************************************************************************************
38 string IndicatorCommand::getHelpString(){
40 string helpString = "";
41 helpString += "The indicator command can be run in 3 ways: with a shared or relabund file and a design file, or with a shared or relabund file and a tree file, or with a shared or relabund file, tree file and design file. \n";
42 helpString += "The indicator command outputs a .indicator.summary file and a .indicator.tre if a tree is given. \n";
43 helpString += "The new tree contains labels at each internal node. The label is the node number so you can relate the tree to the summary file.\n";
44 helpString += "The summary file lists the indicator value for each OTU for each node.\n";
45 helpString += "The indicator command parameters are tree, groups, shared, relabund, design and label. \n";
46 helpString += "The design parameter allows you to relate the tree to the shared or relabund file, if your tree contains the grouping names, or if no tree is provided to group your groups into groupings.\n";
47 helpString += "The groups parameter allows you to specify which of the groups in your shared or relabund you would like analyzed, or if you provide a design file the groups in your design file. The groups may be entered separated by dashes.\n";
48 helpString += "The label parameter indicates at what distance your tree relates to the shared or relabund.\n";
49 helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n";
50 helpString += "The iters parameter allows you to set number of randomization for the P value. The default is 1000.";
51 helpString += "The indicator command should be used in the following format: indicator(tree=test.tre, shared=test.shared, label=0.03)\n";
52 helpString += "Note: No spaces between parameter labels (i.e. tree), '=' and parameters (i.e.yourTreefile).\n";
56 m->errorOut(e, "IndicatorCommand", "getHelpString");
61 //**********************************************************************************************************************
62 IndicatorCommand::IndicatorCommand(){
64 abort = true; calledHelp = true;
66 vector<string> tempOutNames;
67 outputTypes["tree"] = tempOutNames;
68 outputTypes["summary"] = tempOutNames;
71 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
75 //**********************************************************************************************************************
76 IndicatorCommand::IndicatorCommand(string option) {
78 abort = false; calledHelp = false;
80 //allow user to run help
81 if(option == "help") { help(); abort = true; calledHelp = true; }
82 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
85 vector<string> myArray = setParameters();
87 OptionParser parser(option);
88 map<string, string> parameters = parser.getParameters();
90 ValidParameters validParameter;
91 map<string, string>::iterator it;
93 //check to make sure all parameters are valid for command
94 for (it = parameters.begin(); it != parameters.end(); it++) {
95 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
101 m->Treenames.clear();
103 vector<string> tempOutNames;
104 outputTypes["tree"] = tempOutNames;
105 outputTypes["summary"] = tempOutNames;
107 //if the user changes the input directory command factory will send this info to us in the output parameter
108 string inputDir = validParameter.validFile(parameters, "inputdir", false);
109 if (inputDir == "not found"){ inputDir = ""; }
112 it = parameters.find("tree");
113 //user has given a template file
114 if(it != parameters.end()){
115 path = m->hasPath(it->second);
116 //if the user has not given a path then, add inputdir. else leave path alone.
117 if (path == "") { parameters["tree"] = inputDir + it->second; }
120 it = parameters.find("shared");
121 //user has given a template file
122 if(it != parameters.end()){
123 path = m->hasPath(it->second);
124 //if the user has not given a path then, add inputdir. else leave path alone.
125 if (path == "") { parameters["shared"] = inputDir + it->second; }
128 it = parameters.find("relabund");
129 //user has given a template file
130 if(it != parameters.end()){
131 path = m->hasPath(it->second);
132 //if the user has not given a path then, add inputdir. else leave path alone.
133 if (path == "") { parameters["relabund"] = inputDir + it->second; }
136 it = parameters.find("design");
137 //user has given a template file
138 if(it != parameters.end()){
139 path = m->hasPath(it->second);
140 //if the user has not given a path then, add inputdir. else leave path alone.
141 if (path == "") { parameters["design"] = inputDir + it->second; }
145 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
147 //check for required parameters
148 treefile = validParameter.validFile(parameters, "tree", true);
149 if (treefile == "not open") { treefile = ""; abort = true; }
150 else if (treefile == "not found") { treefile = ""; }
151 else { m->setTreeFile(treefile); }
153 sharedfile = validParameter.validFile(parameters, "shared", true);
154 if (sharedfile == "not open") { abort = true; }
155 else if (sharedfile == "not found") { sharedfile = ""; }
156 else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
158 relabundfile = validParameter.validFile(parameters, "relabund", true);
159 if (relabundfile == "not open") { abort = true; }
160 else if (relabundfile == "not found") { relabundfile = ""; }
161 else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
163 designfile = validParameter.validFile(parameters, "design", true);
164 if (designfile == "not open") { designfile = ""; abort = true; }
165 else if (designfile == "not found") { designfile = ""; }
166 else { m->setDesignFile(designfile); }
168 groups = validParameter.validFile(parameters, "groups", false);
169 if (groups == "not found") { groups = ""; Groups.push_back("all"); }
170 else { m->splitAtDash(groups, Groups); }
171 m->setGroups(Groups);
173 label = validParameter.validFile(parameters, "label", false);
174 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=""; }
176 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
177 m->mothurConvert(temp, iters);
179 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
180 m->setProcessors(temp);
181 m->mothurConvert(temp, processors);
183 if ((relabundfile == "") && (sharedfile == "")) {
184 //is there are current file available for either of these?
185 //give priority to shared, then relabund
186 sharedfile = m->getSharedFile();
187 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
189 relabundfile = m->getRelAbundFile();
190 if (relabundfile != "") { inputFileName = relabundfile; m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
192 m->mothurOut("No valid current files. You must provide a shared or relabund."); m->mothurOutEndLine();
199 if ((designfile == "") && (treefile == "")) {
200 treefile = m->getTreeFile();
201 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
203 designfile = m->getDesignFile();
204 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
206 m->mothurOut("[ERROR]: You must provide either a tree or design file."); m->mothurOutEndLine(); abort = true;
211 if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("[ERROR]: You may not use both a shared and relabund file."); m->mothurOutEndLine(); abort = true; }
216 catch(exception& e) {
217 m->errorOut(e, "IndicatorCommand", "IndicatorCommand");
221 //**********************************************************************************************************************
223 int IndicatorCommand::execute(){
226 if (abort == true) { if (calledHelp) { return 0; } return 2; }
228 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
230 int start = time(NULL);
232 //read designfile if given and set up groups for read of sharedfiles
233 if (designfile != "") {
234 designMap = new GroupMap(designfile);
235 designMap->readDesignMap();
237 //fill Groups - checks for "all" and for any typo groups
239 vector<string> nameGroups = designMap->getNamesOfGroups();
240 util.setGroups(Groups, nameGroups);
241 designMap->setNamesOfGroups(nameGroups);
243 //loop through the Groups and fill Globaldata's Groups with the design file info
244 vector<string> namesSeqs = designMap->getNamesSeqs(Groups);
245 m->setGroups(namesSeqs);
248 /***************************************************/
249 // use smart distancing to get right sharedRabund //
250 /***************************************************/
251 if (sharedfile != "") {
253 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } return 0; }
254 if (lookup[0] == NULL) { m->mothurOut("[ERROR] reading shared file."); m->mothurOutEndLine(); return 0; }
257 if (m->control_pressed) { if (designfile != "") { delete designMap; } for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } return 0; }
258 if (lookupFloat[0] == NULL) { m->mothurOut("[ERROR] reading relabund file."); m->mothurOutEndLine(); return 0; }
261 //reset groups if needed
262 if (designfile != "") { m->setGroups(Groups); }
264 /***************************************************/
265 // reading tree info //
266 /***************************************************/
267 if (treefile != "") {
268 string groupfile = "";
269 m->setTreeFile(treefile);
270 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
271 treeMap = new TreeMap();
272 bool mismatch = false;
274 for (int i = 0; i < m->Treenames.size(); i++) {
275 //sanity check - is this a group that is not in the sharedfile?
276 if (designfile == "") {
277 if (!(m->inUsersGroups(m->Treenames[i], m->getAllGroups()))) {
278 m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
281 treeMap->addSeq(m->Treenames[i], "Group1");
283 vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
284 vector<string> myNames = designMap->getNamesSeqs(myGroups);
286 for(int k = 0; k < myNames.size(); k++) {
287 if (!(m->inUsersGroups(myNames[k], m->getAllGroups()))) {
288 m->mothurOut("[ERROR]: " + myNames[k] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
292 treeMap->addSeq(m->Treenames[i], "Group1");
296 if ((designfile != "") && (m->Treenames.size() != Groups.size())) { cout << Groups.size() << '\t' << m->Treenames.size() << endl; m->mothurOut("[ERROR]: You design file does not match your tree, aborting."); m->mothurOutEndLine(); mismatch = true; }
298 if (mismatch) { //cleanup and exit
299 if (designfile != "") { delete designMap; }
300 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
301 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
306 read = new ReadNewickTree(treefile);
307 int readOk = read->read(treeMap);
309 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
311 vector<Tree*> T = read->getTrees();
315 if (m->control_pressed) {
316 if (designfile != "") { delete designMap; }
317 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
318 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
319 for (int i = 0; i < T.size(); i++) { delete T[i]; } delete treeMap; return 0;
322 map<string, string> nameMap;
323 T[0]->assembleTree(nameMap);
325 /***************************************************/
326 // create ouptut tree - respecting pickedGroups //
327 /***************************************************/
328 Tree* outputTree = new Tree(m->getNumGroups(), treeMap);
330 outputTree->getSubTree(T[0], m->getGroups());
331 outputTree->assembleTree(nameMap);
333 //no longer need original tree, we have output tree to use and label
334 for (int i = 0; i < T.size(); i++) { delete T[i]; }
336 if (m->control_pressed) {
337 if (designfile != "") { delete designMap; }
338 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
339 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
340 delete outputTree; delete treeMap; return 0;
343 /***************************************************/
344 // get indicator species values //
345 /***************************************************/
346 GetIndicatorSpecies(outputTree);
347 delete outputTree; delete treeMap;
349 }else { //run with design file only
350 //get indicator species
351 GetIndicatorSpecies();
354 if (designfile != "") { delete designMap; }
355 if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
356 else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
358 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
360 //set tree file as new current treefile
361 if (treefile != "") {
363 itTypes = outputTypes.find("tree");
364 if (itTypes != outputTypes.end()) {
365 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTreeFile(current); }
369 m->mothurOutEndLine(); m->mothurOutEndLine(); m->mothurOut("It took " + toString(time(NULL) - start) + " secs to find the indicator species."); m->mothurOutEndLine();
370 m->mothurOutEndLine();
371 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
372 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
373 m->mothurOutEndLine();
377 catch(exception& e) {
378 m->errorOut(e, "IndicatorCommand", "execute");
382 //**********************************************************************************************************************
383 //divide shared or relabund file by groupings in the design file
384 //report all otu values to file
385 int IndicatorCommand::GetIndicatorSpecies(){
387 string thisOutputDir = outputDir;
388 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
389 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
390 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
393 m->openOutputFile(outputFileName, out);
394 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
395 m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
398 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
399 else { numBins = lookupFloat[0]->getNumBins(); }
401 if (m->control_pressed) { out.close(); return 0; }
403 /*****************************************************/
404 //create vectors containing rabund info //
405 /*****************************************************/
407 vector<float> indicatorValues; //size of numBins
408 vector<float> pValues;
409 map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
411 if (sharedfile != "") {
412 vector< vector<SharedRAbundVector*> > groupings;
413 set<string> groupsAlreadyAdded;
414 vector<SharedRAbundVector*> subset;
417 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
419 for (int k = 0; k < lookup.size(); k++) {
420 //are you from this grouping?
421 if (designMap->getGroup(lookup[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
422 subset.push_back(lookup[k]);
423 groupsAlreadyAdded.insert(lookup[k]->getGroup());
426 if (subset.size() != 0) { groupings.push_back(subset); }
430 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
432 indicatorValues = getValues(groupings, randomGroupingsMap);
434 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
436 vector< vector<SharedRAbundFloatVector*> > groupings;
437 set<string> groupsAlreadyAdded;
438 vector<SharedRAbundFloatVector*> subset;
441 for (int i = 0; i < (designMap->getNamesOfGroups()).size(); i++) {
442 for (int k = 0; k < lookupFloat.size(); k++) {
443 //are you from this grouping?
444 if (designMap->getGroup(lookupFloat[k]->getGroup()) == (designMap->getNamesOfGroups())[i]) {
445 subset.push_back(lookupFloat[k]);
446 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
449 if (subset.size() != 0) { groupings.push_back(subset); }
453 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
455 indicatorValues = getValues(groupings, randomGroupingsMap);
457 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
460 if (m->control_pressed) { out.close(); return 0; }
463 /******************************************************/
464 //output indicator values to table form //
465 /*****************************************************/
466 out << "OTU\tIndicator_Value\tpValue" << endl;
467 for (int j = 0; j < indicatorValues.size(); j++) {
469 if (m->control_pressed) { out.close(); return 0; }
471 out << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
473 if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
474 else { out << "<" << (1/(float)iters) << endl; }
476 if (pValues[j] <= 0.05) {
477 cout << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
478 string pValueString = "<" + toString((1/(float)iters));
479 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
480 else { cout << "<" << (1/(float)iters); }
481 m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
482 m->mothurOutEndLine();
490 catch(exception& e) {
491 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
495 //**********************************************************************************************************************
496 //traverse tree finding indicator species values for each otu at each node
497 //label node with otu number that has highest indicator value
498 //report all otu values to file
499 int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
502 string thisOutputDir = outputDir;
503 if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
504 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
505 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
508 m->openOutputFile(outputFileName, out);
509 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
512 if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
513 else { numBins = lookupFloat[0]->getNumBins(); }
517 for (int i = 0; i < numBins; i++) { out << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; }
520 m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
522 string treeOutputDir = outputDir;
523 if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
524 string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
527 //create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
528 map<int, set<string> > nodeToDescendants;
529 map<int, set<int> > descendantNodes;
530 for (int i = 0; i < T->getNumNodes(); i++) {
531 if (m->control_pressed) { return 0; }
533 nodeToDescendants[i] = getDescendantList(T, i, nodeToDescendants, descendantNodes);
536 //you need the distances to leaf to decide grouping below
537 //this will also set branch lengths if the tree does not include them
538 map<int, float> distToRoot = getDistToRoot(T);
541 for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
542 //cout << endl << i+1 << endl;
543 if (m->control_pressed) { out.close(); return 0; }
545 /*****************************************************/
546 //create vectors containing rabund info //
547 /*****************************************************/
549 vector<float> indicatorValues; //size of numBins
550 vector<float> pValues;
551 map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
553 if (sharedfile != "") {
554 vector< vector<SharedRAbundVector*> > groupings;
556 //get nodes that will be a valid grouping
557 //you are valid if you are not one of my descendants
558 //AND your distToRoot is >= mine
559 //AND you were not added as part of a larger grouping. Largest nodes are added first.
561 set<string> groupsAlreadyAdded;
562 //create a grouping with my grouping
563 vector<SharedRAbundVector*> subset;
565 int doneCount = nodeToDescendants[i].size();
566 for (int k = 0; k < lookup.size(); k++) {
567 //is this descendant of i
568 if ((nodeToDescendants[i].count(lookup[k]->getGroup()) != 0)) {
569 subset.push_back(lookup[k]);
570 groupsAlreadyAdded.insert(lookup[k]->getGroup());
573 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
575 if (subset.size() != 0) { groupings.push_back(subset); }
578 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
581 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
582 vector<SharedRAbundVector*> subset;
584 int doneCount = nodeToDescendants[j].size();
585 for (int k = 0; k < lookup.size(); k++) {
586 //is this descendant of j, and we didn't already add this as part of a larger grouping
587 if ((nodeToDescendants[j].count(lookup[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0)) {
588 subset.push_back(lookup[k]);
589 groupsAlreadyAdded.insert(lookup[k]->getGroup());
592 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
595 //if subset.size == 0 then the node was added as part of a larger grouping
596 if (subset.size() != 0) { groupings.push_back(subset); }
600 if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
602 indicatorValues = getValues(groupings, randomGroupingsMap);
604 pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
606 vector< vector<SharedRAbundFloatVector*> > groupings;
608 //get nodes that will be a valid grouping
609 //you are valid if you are not one of my descendants
610 //AND your distToRoot is >= mine
611 //AND you were not added as part of a larger grouping. Largest nodes are added first.
613 set<string> groupsAlreadyAdded;
614 //create a grouping with my grouping
615 vector<SharedRAbundFloatVector*> subset;
617 int doneCount = nodeToDescendants[i].size();
618 for (int k = 0; k < lookupFloat.size(); k++) {
619 //is this descendant of i
620 if ((nodeToDescendants[i].count(lookupFloat[k]->getGroup()) != 0)) {
621 subset.push_back(lookupFloat[k]);
622 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
625 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
627 if (subset.size() != 0) { groupings.push_back(subset); }
629 for (int j = (T->getNumNodes()-1); j >= 0; j--) {
630 if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
631 vector<SharedRAbundFloatVector*> subset;
633 int doneCount = nodeToDescendants[j].size();
634 for (int k = 0; k < lookupFloat.size(); k++) {
635 //is this descendant of j, and we didn't already add this as part of a larger grouping
636 if ((nodeToDescendants[j].count(lookupFloat[k]->getGroup()) != 0) && (groupsAlreadyAdded.count(lookupFloat[k]->getGroup()) == 0)) {
637 subset.push_back(lookupFloat[k]);
638 groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
641 if (count == doneCount) { break; } //quit once you get the rabunds for this grouping
644 //if subset.size == 0 then the node was added as part of a larger grouping
645 if (subset.size() != 0) { groupings.push_back(subset); }
649 if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
651 indicatorValues = getValues(groupings, randomGroupingsMap);
653 pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
656 if (m->control_pressed) { out.close(); return 0; }
659 /******************************************************/
660 //output indicator values to table form + label tree //
661 /*****************************************************/
662 out << (i+1) << '\t';
663 for (int j = 0; j < indicatorValues.size(); j++) {
665 if (m->control_pressed) { out.close(); return 0; }
667 if (pValues[j] < (1/(float)iters)) {
668 out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
670 out << indicatorValues[j] << '\t' << pValues[j] << '\t';
673 if (pValues[j] <= 0.05) {
674 cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
675 string pValueString = "<" + toString((1/(float)iters));
676 if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
677 else { cout << "<" << (1/(float)iters); }
678 m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
679 m->mothurOutEndLine();
684 T->tree[i].setLabel((i+1));
689 m->openOutputFile(outputTreeFileName, outTree);
690 outputNames.push_back(outputTreeFileName); outputTypes["tree"].push_back(outputTreeFileName);
692 T->print(outTree, "both");
697 catch(exception& e) {
698 m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
702 //**********************************************************************************************************************
703 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
705 vector<float> values;
706 map< vector<int>, vector<int> >::iterator it;
709 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
711 if (m->control_pressed) { return values; }
714 float AijDenominator = 0.0;
717 //get overall abundance of each grouping
718 for (int j = 0; j < groupings.size(); j++) {
720 float totalAbund = 0;
722 for (int k = 0; k < groupings[j].size(); k++) {
723 vector<int> temp; temp.push_back(j); temp.push_back(k);
724 it = groupingsMap.find(temp);
726 if (it == groupingsMap.end()) { //this one didnt get moved
727 totalAbund += groupings[j][k]->getAbundance(i);
728 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
730 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
731 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
737 float Aij = (totalAbund / (float) groupings[j].size());
738 terms.push_back(Aij);
740 //percentage of sites represented
741 Bij.push_back(numNotZero / (float) groupings[j].size());
743 AijDenominator += Aij;
746 float maxIndVal = 0.0;
747 for (int j = 0; j < terms.size(); j++) {
748 float thisAij = (terms[j] / AijDenominator); //relative abundance
749 float thisValue = thisAij * Bij[j] * 100.0;
752 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
755 values.push_back(maxIndVal);
760 catch(exception& e) {
761 m->errorOut(e, "IndicatorCommand", "getValues");
765 //**********************************************************************************************************************
766 //same as above, just data type difference
767 vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
769 vector<float> values;
771 /*for (int j = 0; j < groupings.size(); j++) {
772 cout << "grouping " << j << endl;
773 for (int k = 0; k < groupings[j].size(); k++) {
774 cout << groupings[j][k]->getGroup() << endl;
777 map< vector<int>, vector<int> >::iterator it;
780 for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
782 float AijDenominator = 0.0;
784 //get overall abundance of each grouping
785 for (int j = 0; j < groupings.size(); j++) {
787 int totalAbund = 0.0;
789 for (int k = 0; k < groupings[j].size(); k++) {
790 vector<int> temp; temp.push_back(j); temp.push_back(k);
791 it = groupingsMap.find(temp);
793 if (it == groupingsMap.end()) { //this one didnt get moved
794 totalAbund += groupings[j][k]->getAbundance(i);
795 if (groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
797 totalAbund += groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
798 if (groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
804 float Aij = (totalAbund / (float) groupings[j].size());
805 terms.push_back(Aij);
807 //percentage of sites represented
808 Bij.push_back(numNotZero / (float) groupings[j].size());
810 AijDenominator += Aij;
813 float maxIndVal = 0.0;
814 for (int j = 0; j < terms.size(); j++) {
815 float thisAij = (terms[j] / AijDenominator); //relative abundance
816 float thisValue = thisAij * Bij[j] * 100.0;
819 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
822 values.push_back(maxIndVal);
827 catch(exception& e) {
828 m->errorOut(e, "IndicatorCommand", "getValues");
832 //**********************************************************************************************************************
833 //you need the distances to root to decide groupings
834 //this will also set branch lengths if the tree does not include them
835 map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
837 map<int, float> dists;
839 bool hasBranchLengths = false;
840 for (int i = 0; i < T->getNumNodes(); i++) {
841 if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; }
844 //set branchlengths if needed
845 if (!hasBranchLengths) {
846 for (int i = 0; i < T->getNumNodes(); i++) {
848 int lc = T->tree[i].getLChild();
849 int rc = T->tree[i].getRChild();
851 if (lc == -1) { // you are a leaf
852 //if you are a leaf set you priliminary length to 1.0, this may adjust later
853 T->tree[i].setBranchLength(1.0);
855 }else{ // you are an internal node
856 //look at your children's length to leaf
857 float ldist = dists[lc];
858 float rdist = dists[rc];
860 float greater = ldist;
861 if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
862 else { dists[i] = rdist + 1.0; }
865 //branch length = difference + 1
866 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
867 T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
874 for (int i = 0; i < T->getNumNodes(); i++) {
879 while(T->tree[index].getParent() != -1){
880 if (T->tree[index].getBranchLength() != -1) {
881 sum += abs(T->tree[index].getBranchLength());
883 index = T->tree[index].getParent();
891 catch(exception& e) {
892 m->errorOut(e, "IndicatorCommand", "getLengthToLeaf");
896 //**********************************************************************************************************************
897 set<string> IndicatorCommand::getDescendantList(Tree*& T, int i, map<int, set<string> > descendants, map<int, set<int> >& nodes){
901 set<string>::iterator it;
903 int lc = T->tree[i].getLChild();
904 int rc = T->tree[i].getRChild();
906 if (lc == -1) { //you are a leaf your only descendant is yourself
907 set<int> temp; temp.insert(i);
910 if (designfile == "") {
911 names.insert(T->tree[i].getName());
913 vector<string> myGroup; myGroup.push_back(T->tree[i].getName());
914 vector<string> myReps = designMap->getNamesSeqs(myGroup);
915 for (int k = 0; k < myReps.size(); k++) {
916 names.insert(myReps[k]);
920 }else{ //your descedants are the combination of your childrens descendants
921 names = descendants[lc];
922 nodes[i] = nodes[lc];
923 for (it = descendants[rc].begin(); it != descendants[rc].end(); it++) {
926 for (set<int>::iterator itNum = nodes[rc].begin(); itNum != nodes[rc].end(); itNum++) {
927 nodes[i].insert(*itNum);
929 //you are your own descendant
935 catch(exception& e) {
936 m->errorOut(e, "IndicatorCommand", "getDescendantList");
940 //**********************************************************************************************************************
941 int IndicatorCommand::getShared(){
943 InputData* input = new InputData(sharedfile, "sharedfile");
944 lookup = input->getSharedRAbundVectors();
945 string lastLabel = lookup[0]->getLabel();
947 if (label == "") { label = lastLabel; delete input; return 0; }
949 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
950 set<string> labels; labels.insert(label);
951 set<string> processedLabels;
952 set<string> userLabels = labels;
954 //as long as you are not at the end of the file or done wih the lines you want
955 while((lookup[0] != NULL) && (userLabels.size() != 0)) {
956 if (m->control_pressed) { delete input; return 0; }
958 if(labels.count(lookup[0]->getLabel()) == 1){
959 processedLabels.insert(lookup[0]->getLabel());
960 userLabels.erase(lookup[0]->getLabel());
964 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
965 string saveLabel = lookup[0]->getLabel();
967 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
968 lookup = input->getSharedRAbundVectors(lastLabel);
970 processedLabels.insert(lookup[0]->getLabel());
971 userLabels.erase(lookup[0]->getLabel());
973 //restore real lastlabel to save below
974 lookup[0]->setLabel(saveLabel);
978 lastLabel = lookup[0]->getLabel();
980 //get next line to process
981 //prevent memory leak
982 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
983 lookup = input->getSharedRAbundVectors();
987 if (m->control_pressed) { delete input; return 0; }
989 //output error messages about any remaining user labels
990 set<string>::iterator it;
991 bool needToRun = false;
992 for (it = userLabels.begin(); it != userLabels.end(); it++) {
993 m->mothurOut("Your file does not include the label " + *it);
994 if (processedLabels.count(lastLabel) != 1) {
995 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
998 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1002 //run last label if you need to
1003 if (needToRun == true) {
1004 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
1005 lookup = input->getSharedRAbundVectors(lastLabel);
1011 catch(exception& e) {
1012 m->errorOut(e, "IndicatorCommand", "getShared");
1016 //**********************************************************************************************************************
1017 int IndicatorCommand::getSharedFloat(){
1019 InputData* input = new InputData(relabundfile, "relabund");
1020 lookupFloat = input->getSharedRAbundFloatVectors();
1021 string lastLabel = lookupFloat[0]->getLabel();
1023 if (label == "") { label = lastLabel; delete input; return 0; }
1025 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1026 set<string> labels; labels.insert(label);
1027 set<string> processedLabels;
1028 set<string> userLabels = labels;
1030 //as long as you are not at the end of the file or done wih the lines you want
1031 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
1033 if (m->control_pressed) { delete input; return 0; }
1035 if(labels.count(lookupFloat[0]->getLabel()) == 1){
1036 processedLabels.insert(lookupFloat[0]->getLabel());
1037 userLabels.erase(lookupFloat[0]->getLabel());
1041 if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1042 string saveLabel = lookupFloat[0]->getLabel();
1044 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1045 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1047 processedLabels.insert(lookupFloat[0]->getLabel());
1048 userLabels.erase(lookupFloat[0]->getLabel());
1050 //restore real lastlabel to save below
1051 lookupFloat[0]->setLabel(saveLabel);
1055 lastLabel = lookupFloat[0]->getLabel();
1057 //get next line to process
1058 //prevent memory leak
1059 for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; }
1060 lookupFloat = input->getSharedRAbundFloatVectors();
1064 if (m->control_pressed) { delete input; return 0; }
1066 //output error messages about any remaining user labels
1067 set<string>::iterator it;
1068 bool needToRun = false;
1069 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1070 m->mothurOut("Your file does not include the label " + *it);
1071 if (processedLabels.count(lastLabel) != 1) {
1072 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1075 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1079 //run last label if you need to
1080 if (needToRun == true) {
1081 for (int i = 0; i < lookupFloat.size(); i++) { if (lookupFloat[i] != NULL) { delete lookupFloat[i]; } }
1082 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
1088 catch(exception& e) {
1089 m->errorOut(e, "IndicatorCommand", "getShared");
1093 //**********************************************************************************************************************
1094 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1096 vector<float> pvalues;
1097 pvalues.resize(indicatorValues.size(), 0);
1099 for(int i=0;i<numIters;i++){
1100 if (m->control_pressed) { break; }
1101 groupingsMap = randomizeGroupings(groupings, num);
1102 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1104 for (int j = 0; j < indicatorValues.size(); j++) {
1105 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1111 }catch(exception& e) {
1112 m->errorOut(e, "IndicatorCommand", "driver");
1116 //**********************************************************************************************************************
1117 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1119 vector<float> pvalues;
1121 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1122 if(processors == 1){
1123 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1124 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1127 //divide iters between processors
1128 int numItersPerProcessor = iters / processors;
1130 vector<int> processIDS;
1134 //loop through and create all the processes you want
1135 while (process != processors) {
1139 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1141 }else if (pid == 0){
1142 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1144 //pass pvalues to parent
1146 string tempFile = toString(getpid()) + ".pvalues.temp";
1147 m->openOutputFile(tempFile, out);
1150 for (int i = 0; i < pvalues.size(); i++) {
1151 out << pvalues[i] << '\t';
1159 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1160 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1166 //special case for last processor in case it doesn't divide evenly
1167 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1168 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1170 //force parent to wait until all the processes are done
1171 for (int i=0;i<processIDS.size();i++) {
1172 int temp = processIDS[i];
1177 for (int i = 0; i < processIDS.size(); i++) {
1179 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1180 m->openInputFile(tempFile, in);
1182 ////// to do ///////////
1183 int numTemp; numTemp = 0;
1184 for (int j = 0; j < pvalues.size(); j++) {
1185 in >> numTemp; m->gobble(in);
1186 pvalues[j] += numTemp;
1188 in.close(); m->mothurRemove(tempFile);
1190 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1193 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1194 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1199 catch(exception& e) {
1200 m->errorOut(e, "IndicatorCommand", "getPValues");
1205 //**********************************************************************************************************************
1206 //same as above, just data type difference
1207 vector<float> IndicatorCommand::driver(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues, int numIters){
1209 vector<float> pvalues;
1210 pvalues.resize(indicatorValues.size(), 0);
1212 for(int i=0;i<numIters;i++){
1213 if (m->control_pressed) { break; }
1214 groupingsMap = randomizeGroupings(groupings, num);
1215 vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
1217 for (int j = 0; j < indicatorValues.size(); j++) {
1218 if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
1224 }catch(exception& e) {
1225 m->errorOut(e, "IndicatorCommand", "driver");
1229 //**********************************************************************************************************************
1230 //same as above, just data type difference
1231 vector<float> IndicatorCommand::getPValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap, int num, vector<float> indicatorValues){
1233 vector<float> pvalues;
1235 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1236 if(processors == 1){
1237 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1238 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1241 //divide iters between processors
1242 int numItersPerProcessor = iters / processors;
1244 vector<int> processIDS;
1247 //loop through and create all the processes you want
1248 while (process != processors) {
1252 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1254 }else if (pid == 0){
1255 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1257 //pass pvalues to parent
1259 string tempFile = toString(getpid()) + ".pvalues.temp";
1260 m->openOutputFile(tempFile, out);
1263 for (int i = 0; i < pvalues.size(); i++) {
1264 out << pvalues[i] << '\t';
1272 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1273 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1279 //special case for last processor in case it doesn't divide evenly
1280 numItersPerProcessor = iters - ((processors-1) * numItersPerProcessor);
1281 pvalues = driver(groupings, groupingsMap, num, indicatorValues, numItersPerProcessor);
1283 //force parent to wait until all the processes are done
1284 for (int i=0;i<processIDS.size();i++) {
1285 int temp = processIDS[i];
1290 for (int i = 0; i < processIDS.size(); i++) {
1292 string tempFile = toString(processIDS[i]) + ".pvalues.temp";
1293 m->openInputFile(tempFile, in);
1295 ////// to do ///////////
1296 int numTemp; numTemp = 0;
1297 for (int j = 0; j < pvalues.size(); j++) {
1298 in >> numTemp; m->gobble(in);
1299 pvalues[j] += numTemp;
1301 in.close(); m->mothurRemove(tempFile);
1303 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1306 pvalues = driver(groupings, groupingsMap, num, indicatorValues, iters);
1307 for (int i = 0; i < pvalues.size(); i++) { pvalues[i] /= (double)iters; }
1312 catch(exception& e) {
1313 m->errorOut(e, "IndicatorCommand", "getPValues");
1317 //**********************************************************************************************************************
1318 //swap groups between groupings, in essence randomizing the second column of the design file
1319 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundVector*> >& groupings, int numLookupGroups){
1322 map< vector<int>, vector<int> > randomGroupings;
1324 for (int i = 0; i < numLookupGroups; i++) {
1325 if (m->control_pressed) {break;}
1327 //get random groups to swap to switch with
1328 //generate random int between 0 and groupings.size()-1
1329 int z = m->getRandomIndex(groupings.size()-1);
1330 int x = m->getRandomIndex(groupings.size()-1);
1331 int a = m->getRandomIndex(groupings[z].size()-1);
1332 int b = m->getRandomIndex(groupings[x].size()-1);
1333 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1334 //if ((z < 0) || (z > 1) || x<0 || x>1 || a <0 || a>groupings[z].size()-1 || b<0 || b>groupings[x].size()-1) { cout << "probelm" << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl; }
1339 from.push_back(z); from.push_back(a);
1340 to.push_back(x); to.push_back(b);
1342 randomGroupings[from] = to;
1344 //cout << "done" << endl;
1345 return randomGroupings;
1347 catch(exception& e) {
1348 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1352 //**********************************************************************************************************************
1353 //swap groups between groupings, in essence randomizing the second column of the design file
1354 map< vector<int>, vector<int> > IndicatorCommand::randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >& groupings, int numLookupGroups){
1357 map< vector<int>, vector<int> > randomGroupings;
1359 for (int i = 0; i < numLookupGroups; i++) {
1361 //get random groups to swap to switch with
1362 //generate random int between 0 and groupings.size()-1
1363 int z = m->getRandomIndex(groupings.size()-1);
1364 int x = m->getRandomIndex(groupings.size()-1);
1365 int a = m->getRandomIndex(groupings[z].size()-1);
1366 int b = m->getRandomIndex(groupings[x].size()-1);
1367 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
1372 from.push_back(z); from.push_back(a);
1373 to.push_back(x); to.push_back(b);
1375 randomGroupings[from] = to;
1378 return randomGroupings;
1380 catch(exception& e) {
1381 m->errorOut(e, "IndicatorCommand", "randomizeGroupings");
1386 /*****************************************************************/