2 * phylodiversitycommand.cpp
5 * Created by westcott on 4/30/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "phylodiversitycommand.h"
12 //**********************************************************************************************************************
13 vector<string> PhyloDiversityCommand::setParameters(){
16 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
17 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup);
18 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
19 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
20 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
21 CommandParameter pfreq("freq", "Number", "", "100", "", "", "",false,false); parameters.push_back(pfreq);
22 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
23 CommandParameter prarefy("rarefy", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prarefy);
24 CommandParameter psummary("summary", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(psummary);
25 CommandParameter pcollect("collect", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcollect);
26 CommandParameter pscale("scale", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pscale);
27 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
28 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
30 vector<string> myArray;
31 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
35 m->errorOut(e, "PhyloDiversityCommand", "setParameters");
39 //**********************************************************************************************************************
40 string PhyloDiversityCommand::getHelpString(){
42 string helpString = "";
43 helpString += "The phylo.diversity command parameters are tree, group, name, groups, iters, freq, processors, scale, rarefy, collect and summary. tree and group are required, unless you have valid current files.\n";
44 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed. The group names are separated by dashes. By default all groups are used.\n";
45 helpString += "The iters parameter allows you to specify the number of randomizations to preform, by default iters=1000, if you set rarefy to true.\n";
46 helpString += "The freq parameter is used indicate when to output your data, by default it is set to 100. But you can set it to a percentage of the number of sequence. For example freq=0.10, means 10%. \n";
47 helpString += "The scale parameter is used indicate that you want your output scaled to the number of sequences sampled, default = false. \n";
48 helpString += "The rarefy parameter allows you to create a rarefaction curve. The default is false.\n";
49 helpString += "The collect parameter allows you to create a collectors curve. The default is false.\n";
50 helpString += "The summary parameter allows you to create a .summary file. The default is true.\n";
51 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
52 helpString += "The phylo.diversity command should be in the following format: phylo.diversity(groups=yourGroups, rarefy=yourRarefy, iters=yourIters).\n";
53 helpString += "Example phylo.diversity(groups=A-B-C, rarefy=T, iters=500).\n";
54 helpString += "The phylo.diversity command output two files: .phylo.diversity and if rarefy=T, .rarefaction.\n";
55 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
59 m->errorOut(e, "PhyloDiversityCommand", "getHelpString");
65 //**********************************************************************************************************************
66 PhyloDiversityCommand::PhyloDiversityCommand(){
68 abort = true; calledHelp = true;
70 vector<string> tempOutNames;
71 outputTypes["phylodiv"] = tempOutNames;
72 outputTypes["rarefy"] = tempOutNames;
73 outputTypes["summary"] = tempOutNames;
76 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
80 //**********************************************************************************************************************
81 PhyloDiversityCommand::PhyloDiversityCommand(string option) {
83 abort = false; calledHelp = false;
85 //allow user to run help
86 if(option == "help") { help(); abort = true; calledHelp = true; }
87 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
90 vector<string> myArray = setParameters();;
92 OptionParser parser(option);
93 map<string,string> parameters = parser.getParameters();
94 map<string,string>::iterator it;
96 ValidParameters validParameter;
98 //check to make sure all parameters are valid for command
99 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
100 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
103 //initialize outputTypes
104 vector<string> tempOutNames;
105 outputTypes["phylodiv"] = tempOutNames;
106 outputTypes["rarefy"] = tempOutNames;
107 outputTypes["summary"] = tempOutNames;
109 //if the user changes the input directory command factory will send this info to us in the output parameter
110 string inputDir = validParameter.validFile(parameters, "inputdir", false);
111 if (inputDir == "not found"){ inputDir = ""; }
114 it = parameters.find("tree");
115 //user has given a template file
116 if(it != parameters.end()){
117 path = m->hasPath(it->second);
118 //if the user has not given a path then, add inputdir. else leave path alone.
119 if (path == "") { parameters["tree"] = inputDir + it->second; }
122 it = parameters.find("group");
123 //user has given a template file
124 if(it != parameters.end()){
125 path = m->hasPath(it->second);
126 //if the user has not given a path then, add inputdir. else leave path alone.
127 if (path == "") { parameters["group"] = inputDir + it->second; }
130 it = parameters.find("name");
131 //user has given a template file
132 if(it != parameters.end()){
133 path = m->hasPath(it->second);
134 //if the user has not given a path then, add inputdir. else leave path alone.
135 if (path == "") { parameters["name"] = inputDir + it->second; }
142 m->Treenames.clear();
145 //check for required parameters
146 treefile = validParameter.validFile(parameters, "tree", true);
147 if (treefile == "not open") { abort = true; }
148 else if (treefile == "not found") {
149 //if there is a current design file, use it
150 treefile = m->getTreeFile();
151 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
152 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
153 }else { m->setTreeFile(treefile); }
155 //check for required parameters
156 groupfile = validParameter.validFile(parameters, "group", true);
157 if (groupfile == "not open") { groupfile = ""; abort = true; }
158 else if (groupfile == "not found") { groupfile = ""; }
159 else { m->setGroupFile(groupfile); }
161 namefile = validParameter.validFile(parameters, "name", true);
162 if (namefile == "not open") { abort = true; }
163 else if (namefile == "not found") { namefile = ""; }
164 else { m->setNameFile(namefile); }
166 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(treefile); }
169 temp = validParameter.validFile(parameters, "freq", false); if (temp == "not found") { temp = "100"; }
172 temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
173 convert(temp, iters);
175 temp = validParameter.validFile(parameters, "rarefy", false); if (temp == "not found") { temp = "F"; }
176 rarefy = m->isTrue(temp);
177 if (!rarefy) { iters = 1; }
179 temp = validParameter.validFile(parameters, "summary", false); if (temp == "not found") { temp = "T"; }
180 summary = m->isTrue(temp);
182 temp = validParameter.validFile(parameters, "scale", false); if (temp == "not found") { temp = "F"; }
183 scale = m->isTrue(temp);
185 temp = validParameter.validFile(parameters, "collect", false); if (temp == "not found") { temp = "F"; }
186 collect = m->isTrue(temp);
188 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
189 m->setProcessors(temp);
190 convert(temp, processors);
192 groups = validParameter.validFile(parameters, "groups", false);
193 if (groups == "not found") { groups = ""; }
195 m->splitAtDash(groups, Groups);
196 m->setGroups(Groups);
199 if ((!collect) && (!rarefy) && (!summary)) { m->mothurOut("No outputs selected. You must set either collect, rarefy or summary to true, summary=T by default."); m->mothurOutEndLine(); abort=true; }
203 catch(exception& e) {
204 m->errorOut(e, "PhyloDiversityCommand", "PhyloDiversityCommand");
208 //**********************************************************************************************************************
210 int PhyloDiversityCommand::execute(){
213 if (abort == true) { if (calledHelp) { return 0; } return 2; }
215 m->setTreeFile(treefile);
217 if (groupfile != "") {
218 //read in group map info.
219 tmap = new TreeMap(groupfile);
221 }else{ //fake out by putting everyone in one group
222 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
223 tmap = new TreeMap();
225 for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
228 if (namefile != "") { readNamesFile(); }
230 read = new ReadNewickTree(treefile);
231 int readOk = read->read(tmap);
233 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
235 read->AssembleTrees();
236 vector<Tree*> trees = read->getTrees();
239 //make sure all files match
240 //if you provide a namefile we will use the numNames in the namefile as long as the number of unique match the tree names size.
242 if (namefile != "") {
243 if (numUniquesInName == m->Treenames.size()) { numNamesInTree = nameMap.size(); }
244 else { numNamesInTree = m->Treenames.size(); }
245 }else { numNamesInTree = m->Treenames.size(); }
248 //output any names that are in group file but not in tree
249 if (numNamesInTree < tmap->getNumSeqs()) {
250 for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
251 //is that name in the tree?
253 for (int j = 0; j < m->Treenames.size(); j++) {
254 if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
258 if (m->control_pressed) {
259 delete tmap; for (int i = 0; i < trees.size(); i++) { delete trees[i]; }
260 for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } outputTypes.clear();
265 //then you did not find it so report it
266 if (count == m->Treenames.size()) {
267 //if it is in your namefile then don't remove
268 map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
270 if (it == nameMap.end()) {
271 m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
272 tmap->removeSeq(tmap->namesOfSeqs[i]);
273 i--; //need this because removeSeq removes name from namesOfSeqs
279 SharedUtil* util = new SharedUtil();
280 vector<string> mGroups = m->getGroups();
281 vector<string> tGroups = tmap->getNamesOfGroups();
282 util->setGroups(mGroups, tGroups, "phylo.diversity"); //sets the groups the user wants to analyze
285 //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed
286 for (int i = 0; i < mGroups.size(); i++) { if (mGroups[i] == "xxx") { mGroups.erase(mGroups.begin()+i); break; } }
287 m->setGroups(mGroups);
289 vector<string> outputNames;
291 //for each of the users trees
292 for(int i = 0; i < trees.size(); i++) {
294 if (m->control_pressed) { delete tmap; for (int j = 0; j < trees.size(); j++) { delete trees[j]; } for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; }
296 ofstream outSum, outRare, outCollect;
297 string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.summary";
298 string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.rarefaction";
299 string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv";
301 if (summary) { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); outputTypes["summary"].push_back(outSumFile); }
302 if (rarefy) { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); outputTypes["rarefy"].push_back(outRareFile); }
303 if (collect) { m->openOutputFile(outCollectFile, outCollect); outputNames.push_back(outCollectFile); outputTypes["phylodiv"].push_back(outCollectFile); }
305 int numLeafNodes = trees[i]->getNumLeaves();
307 //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator
308 vector<int> randomLeaf;
309 for (int j = 0; j < numLeafNodes; j++) {
310 if (m->inUsersGroups(trees[i]->tree[j].getGroup(), mGroups) == true) { //is this a node from the group the user selected.
311 randomLeaf.push_back(j);
315 numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using
317 //each group, each sampling, if no rarefy iters = 1;
318 map<string, vector<float> > diversity;
320 //each group, each sampling, if no rarefy iters = 1;
321 map<string, vector<float> > sumDiversity;
323 //find largest group total
324 int largestGroup = 0;
325 for (int j = 0; j < mGroups.size(); j++) {
326 if (tmap->seqsPerGroup[mGroups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[mGroups[j]]; }
328 //initialize diversity
329 diversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0); //numSampled
332 //initialize sumDiversity
333 sumDiversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0);
336 //convert freq percentage to number
338 if (freq < 1.0) { increment = largestGroup * freq;
339 }else { increment = freq; }
341 //initialize sampling spots
342 set<int> numSampledList;
343 for(int k = 1; k <= largestGroup; k++){ if((k == 1) || (k % increment == 0)){ numSampledList.insert(k); } }
344 if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); }
346 //add other groups ending points
347 for (int j = 0; j < mGroups.size(); j++) {
348 if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) { numSampledList.insert(diversity[mGroups[j]].size()-1); }
351 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
353 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
356 vector<int> procIters;
358 int numItersPerProcessor = iters / processors;
360 //divide iters between processes
361 for (int h = 0; h < processors; h++) {
362 if(h == processors - 1){
363 numItersPerProcessor = iters - h * numItersPerProcessor;
365 procIters.push_back(numItersPerProcessor);
368 createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum);
370 }else{ //no need to paralellize if you dont want to rarefy
371 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
376 driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true);
379 if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); }
383 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
385 m->mothurOutEndLine();
386 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
387 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
388 m->mothurOutEndLine();
393 catch(exception& e) {
394 m->errorOut(e, "PhyloDiversityCommand", "execute");
398 //**********************************************************************************************************************
399 int PhyloDiversityCommand::createProcesses(vector<int>& procIters, Tree* t, map< string, vector<float> >& div, map<string, vector<float> >& sumDiv, int numIters, int increment, vector<int>& randomLeaf, set<int>& numSampledList, ofstream& outCollect, ofstream& outSum){
401 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
404 vector<int> processIDS;
405 map< string, vector<float> >::iterator itSum;
407 //loop through and create all the processes you want
408 while (process != processors) {
412 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
415 driver(t, div, sumDiv, procIters[process], increment, randomLeaf, numSampledList, outCollect, outSum, false);
417 string outTemp = outputDir + toString(getpid()) + ".sumDiv.temp";
419 m->openOutputFile(outTemp, out);
421 //output the sumDIversity
422 for (itSum = sumDiv.begin(); itSum != sumDiv.end(); itSum++) {
423 out << itSum->first << '\t' << (itSum->second).size() << '\t';
424 for (int k = 0; k < (itSum->second).size(); k++) {
425 out << (itSum->second)[k] << '\t';
434 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
435 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
440 driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true);
442 //force parent to wait until all the processes are done
443 for (int i=0;i<(processors-1);i++) {
444 int temp = processIDS[i];
448 //get data created by processes
449 for (int i=0;i<(processors-1);i++) {
451 //input the sumDIversity
452 string inTemp = outputDir + toString(processIDS[i]) + ".sumDiv.temp";
454 m->openInputFile(inTemp, in);
456 //output the sumDIversity
457 for (int j = 0; j < sumDiv.size(); j++) {
461 in >> group >> size; m->gobble(in);
463 for (int k = 0; k < size; k++) {
467 sumDiv[group][k] += tempVal;
473 m->mothurRemove(inTemp);
481 catch(exception& e) {
482 m->errorOut(e, "PhyloDiversityCommand", "createProcesses");
486 //**********************************************************************************************************************
487 int PhyloDiversityCommand::driver(Tree* t, map< string, vector<float> >& div, map<string, vector<float> >& sumDiv, int numIters, int increment, vector<int>& randomLeaf, set<int>& numSampledList, ofstream& outCollect, ofstream& outSum, bool doSumCollect){
489 int numLeafNodes = randomLeaf.size();
490 vector<string> mGroups = m->getGroups();
492 for (int l = 0; l < numIters; l++) {
493 random_shuffle(randomLeaf.begin(), randomLeaf.end());
496 map<string, int> counts;
497 map< string, set<int> > countedBranch;
498 for (int j = 0; j < mGroups.size(); j++) { counts[mGroups[j]] = 0; countedBranch[mGroups[j]].insert(-2); } //add dummy index to initialize countedBranch sets
500 for(int k = 0; k < numLeafNodes; k++){
502 if (m->control_pressed) { return 0; }
504 //calc branch length of randomLeaf k
505 vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch);
507 //for each group in the groups update the total branch length accounting for the names file
508 vector<string> groups = t->tree[randomLeaf[k]].getGroup();
510 for (int j = 0; j < groups.size(); j++) {
511 int numSeqsInGroupJ = 0;
512 map<string, int>::iterator it;
513 it = t->tree[randomLeaf[k]].pcount.find(groups[j]);
514 if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
515 numSeqsInGroupJ = it->second;
518 if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; }
520 for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
521 div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
523 counts[groups[j]] += numSeqsInGroupJ;
528 //add this diversity to the sum
529 for (int j = 0; j < mGroups.size(); j++) {
530 for (int g = 0; g < div[mGroups[j]].size(); g++) {
531 sumDiv[mGroups[j]][g] += div[mGroups[j]][g];
536 if ((collect) && (l == 0) && doSumCollect) { printData(numSampledList, div, outCollect, 1); }
537 if ((summary) && (l == 0) && doSumCollect) { printSumData(div, outSum, 1); }
543 catch(exception& e) {
544 m->errorOut(e, "PhyloDiversityCommand", "driver");
549 //**********************************************************************************************************************
551 void PhyloDiversityCommand::printSumData(map< string, vector<float> >& div, ofstream& out, int numIters){
554 out << "Groups\tnumSampled\tphyloDiversity" << endl;
556 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
558 vector<string> mGroups = m->getGroups();
559 for (int j = 0; j < mGroups.size(); j++) {
560 int numSampled = (div[mGroups[j]].size()-1);
561 out << mGroups[j] << '\t' << numSampled << '\t';
565 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
566 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
568 out << setprecision(4) << score << endl;
574 catch(exception& e) {
575 m->errorOut(e, "PhyloDiversityCommand", "printSumData");
579 //**********************************************************************************************************************
581 void PhyloDiversityCommand::printData(set<int>& num, map< string, vector<float> >& div, ofstream& out, int numIters){
584 out << "numSampled\t";
585 vector<string> mGroups = m->getGroups();
586 for (int i = 0; i < mGroups.size(); i++) { out << mGroups[i] << '\t'; }
589 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
591 for (set<int>::iterator it = num.begin(); it != num.end(); it++) {
592 int numSampled = *it;
594 out << numSampled << '\t';
596 for (int j = 0; j < mGroups.size(); j++) {
597 if (numSampled < div[mGroups[j]].size()) {
599 if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; }
600 else { score = div[mGroups[j]][numSampled] / (float)numIters; }
602 out << setprecision(4) << score << '\t';
603 }else { out << "NA" << '\t'; }
611 catch(exception& e) {
612 m->errorOut(e, "PhyloDiversityCommand", "printData");
616 //**********************************************************************************************************************
617 //need a vector of floats one branch length for every group the node represents.
618 vector<float> PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, map< string, set<int> >& counted){
621 //calc the branch length
622 //while you aren't at root
626 vector<string> groups = t->tree[leaf].getGroup();
627 sums.resize(groups.size(), 0.0);
629 map<string, map<int, double> > tempTotals; //maps node to total Branch Length
630 map< string, set<int> > tempCounted;
631 set<int>::iterator it;
634 if(t->tree[index].getBranchLength() != -1){
635 for (int k = 0; k < groups.size(); k++) {
636 sums[k] += abs(t->tree[index].getBranchLength());
637 counted[groups[k]].insert(index);
641 for (int k = 0; k < groups.size(); k++) {
642 tempTotals[groups[k]][index] = 0.0;
645 index = t->tree[index].getParent();
647 //while you aren't at root
648 while(t->tree[index].getParent() != -1){
650 if (m->control_pressed) { return sums; }
653 for (int k = 0; k < groups.size(); k++) {
654 map<string, int>::iterator itGroup = t->tree[index].pcount.find(groups[k]);
655 if (itGroup != t->tree[index].pcount.end()) { pcountSize++; }
657 //do both your chidren have have descendants from the users groups?
658 int lc = t->tree[index].getLChild();
659 int rc = t->tree[index].getRChild();
662 itGroup = t->tree[lc].pcount.find(groups[k]);
663 if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; }
666 itGroup = t->tree[rc].pcount.find(groups[k]);
667 if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; }
669 //if yes, add your childrens tempTotals
670 if ((LpcountSize != 0) && (RpcountSize != 0)) {
671 sums[k] += tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
673 for (it = tempCounted[groups[k]].begin(); it != tempCounted[groups[k]].end(); it++) { counted[groups[k]].insert(*it); }
675 //cout << "added to total " << tempTotals[lc] << '\t' << tempTotals[rc] << endl;
676 if (t->tree[index].getBranchLength() != -1) {
677 if (counted[groups[k]].count(index) == 0) {
678 tempTotals[groups[k]][index] = abs(t->tree[index].getBranchLength());
679 tempCounted[groups[k]].insert(index);
681 tempTotals[groups[k]][index] = 0.0;
684 tempTotals[groups[k]][index] = 0.0;
686 }else { //if no, your tempTotal is your childrens temp totals + your branch length
687 tempTotals[groups[k]][index] = tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc];
689 if (counted[groups[k]].count(index) == 0) {
690 tempTotals[groups[k]][index] += abs(t->tree[index].getBranchLength());
691 tempCounted[groups[k]].insert(index);
695 //cout << "temptotal = "<< tempTotals[i] << endl;
698 index = t->tree[index].getParent();
704 catch(exception& e) {
705 m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength");
709 /*****************************************************************/
710 int PhyloDiversityCommand::readNamesFile() {
713 numUniquesInName = 0;
716 m->openInputFile(namefile, in);
718 string first, second;
719 map<string, string>::iterator itNames;
722 in >> first >> second; m->gobble(in);
726 itNames = m->names.find(first);
727 if (itNames == m->names.end()) {
728 m->names[first] = second;
730 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
731 vector<string> dupNames;
732 m->splitAtComma(second, dupNames);
734 for (int i = 0; i < dupNames.size(); i++) {
735 nameMap[dupNames[i]] = dupNames[i];
736 if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
738 }else { m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
744 catch(exception& e) {
745 m->errorOut(e, "PhyloDiversityCommand", "readNamesFile");
750 //**********************************************************************************************************************