1 #ifndef PHYLODIVERSITYCOMMAND_H
2 #define PHYLODIVERSITYCOMMAND_H
5 * phylodiversitycommand.h
8 * Created by westcott on 4/30/10.
9 * Copyright 2010 Schloss Lab. All rights reserved.
13 #include "command.hpp"
14 #include "counttable.h"
15 #include "sharedutilities.h"
18 class PhyloDiversityCommand : public Command {
21 PhyloDiversityCommand(string);
22 PhyloDiversityCommand();
23 ~PhyloDiversityCommand(){}
25 vector<string> setParameters();
26 string getCommandName() { return "phylo.diversity"; }
27 string getCommandCategory() { return "Hypothesis Testing"; }
29 string getHelpString();
30 string getOutputPattern(string);
31 string getCitation() { return "Faith DP (1994). Phylogenetic pattern and the quantification of organismal biodiversity. Philos Trans R Soc Lond B Biol Sci 345: 45-58. \nhttp://www.mothur.org/wiki/Phylo.diversity"; }
32 string getDescription() { return "phylo.diversity"; }
35 void help() { m->mothurOut(getHelpString()); }
39 int iters, processors, numUniquesInName;
40 bool abort, rarefy, summary, collect, scale;
41 string groups, outputDir, treefile, groupfile, namefile, countfile;
42 vector<string> Groups, outputNames; //holds groups to be used, and outputFile names
44 map<string, int> getRootForGroups(Tree* t);
46 void printData(set<int>&, map< string, vector<float> >&, ofstream&, int);
47 void printSumData(map< string, vector<float> >&, ofstream&, int);
48 vector<float> calcBranchLength(Tree*, int, vector< map<string, bool> >&, map<string, int>);
49 int driver(Tree*, map< string, vector<float> >&, map<string, vector<float> >&, int, int, vector<int>&, set<int>&, ofstream&, ofstream&, bool);
50 int createProcesses(vector<int>&, Tree*, map< string, vector<float> >&, map<string, vector<float> >&, int, int, vector<int>&, set<int>&, ofstream&, ofstream&);
54 /***********************************************************************/
58 map< string, vector<float> > div;
59 map<string, vector<float> > sumDiv;
60 map<string, int> rootForGroup;
61 vector<int> randomLeaf;
62 set<int> numSampledList;
70 phylodivData(MothurOut* mout, int ni, map< string, vector<float> > cd, map< string, vector<float> > csd, Tree* tree, CountTable* count, int incre, vector<int> crl, set<int> nsl, map<string, int> rfg) {
84 /**************************************************************************************************/
85 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
87 static DWORD WINAPI MyPhyloDivThreadFunction(LPVOID lpParam){
88 phylodivData* pDataArray;
89 pDataArray = (phylodivData*)lpParam;
91 int numLeafNodes = pDataArray->randomLeaf.size();
92 vector<string> mGroups = pDataArray->m->getGroups();
94 for (int l = 0; l < pDataArray->numIters; l++) {
95 random_shuffle(pDataArray->randomLeaf.begin(), pDataArray->randomLeaf.end());
98 map<string, int> counts;
99 vector< map<string, bool> > countedBranch;
100 for (int i = 0; i < pDataArray->t->getNumNodes(); i++) {
101 map<string, bool> temp;
102 for (int j = 0; j < mGroups.size(); j++) { temp[mGroups[j]] = false; }
103 countedBranch.push_back(temp);
106 for (int j = 0; j < mGroups.size(); j++) { counts[mGroups[j]] = 0; }
108 for(int k = 0; k < numLeafNodes; k++){
110 if (pDataArray->m->control_pressed) { return 0; }
112 //calc branch length of randomLeaf k
113 //vector<float> br = calcBranchLength(t, randomLeaf[k], countedBranch, rootForGroup);
114 //(Tree* t, int leaf, vector< map<string, bool> >& counted, map<string, int> roots
115 /////////////////////////////////////////////////////////////////////////////////////
117 int index = pDataArray->randomLeaf[k];
119 vector<string> groups = pDataArray->t->tree[pDataArray->randomLeaf[k]].getGroup();
120 br.resize(groups.size(), 0.0);
123 if(pDataArray->t->tree[index].getBranchLength() != -1){
124 for (int k = 0; k < groups.size(); k++) {
125 br[k] += abs(pDataArray->t->tree[index].getBranchLength());
129 index = pDataArray->t->tree[index].getParent();
131 //while you aren't at root
132 while(pDataArray->t->tree[index].getParent() != -1){
134 if (pDataArray->m->control_pressed) { return 0; }
136 for (int k = 0; k < groups.size(); k++) {
138 if (index >= pDataArray->rootForGroup[groups[k]]) { countedBranch[index][groups[k]] = true; } //if you are at this groups "root", then say we are done
140 if (!countedBranch[index][groups[k]]){ //if counted[index][groups[k] is true this groups has already added all br from here to root, so quit early
141 if (pDataArray->t->tree[index].getBranchLength() != -1) {
142 br[k] += abs(pDataArray->t->tree[index].getBranchLength());
144 countedBranch[index][groups[k]] = true;
147 index = pDataArray->t->tree[index].getParent();
149 /////////////////////////////////////////////////////////////////////////////////////
151 //for each group in the groups update the total branch length accounting for the names file
152 groups = pDataArray->t->tree[pDataArray->randomLeaf[k]].getGroup();
154 for (int j = 0; j < groups.size(); j++) {
156 if (pDataArray->m->inUsersGroups(groups[j], mGroups)) {
157 int numSeqsInGroupJ = 0;
158 map<string, int>::iterator it;
159 it = pDataArray->t->tree[pDataArray->randomLeaf[k]].pcount.find(groups[j]);
160 if (it != pDataArray->t->tree[pDataArray->randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j
161 numSeqsInGroupJ = it->second;
164 if (numSeqsInGroupJ != 0) { pDataArray->div[groups[j]][(counts[groups[j]]+1)] = pDataArray->div[groups[j]][counts[groups[j]]] + br[j]; }
166 for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) {
167 pDataArray->div[groups[j]][s] = pDataArray->div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths
169 counts[groups[j]] += numSeqsInGroupJ;
175 //add this diversity to the sum
176 for (int j = 0; j < mGroups.size(); j++) {
177 for (int g = 0; g < pDataArray->div[mGroups[j]].size(); g++) {
178 pDataArray->sumDiv[mGroups[j]][g] += pDataArray->div[mGroups[j]][g];
186 catch(exception& e) {
187 pDataArray->m->errorOut(e, "PhyloDiversityCommand", "MyPhyloDivThreadFunction");