phylodiversitycommand.h

   1 #ifndef PHYLODIVERSITYCOMMAND_H
   2 #define PHYLODIVERSITYCOMMAND_H
   3
   4 /*
   5  *  phylodiversitycommand.h
   6  *  Mothur
   7  *
   8  *  Created by westcott on 4/30/10.
   9  *  Copyright 2010 Schloss Lab. All rights reserved.
  10  *
  11  */
  12
  13 #include "command.hpp"
  14 #include "counttable.h"
  15 #include "sharedutilities.h"
  16 #include "tree.h"
  17
  18 class PhyloDiversityCommand : public Command {
  19
  20         public:
  21                 PhyloDiversityCommand(string);
  22                 PhyloDiversityCommand();
  23                 ~PhyloDiversityCommand(){}
  24
  25                 vector<string> setParameters();
  26                 string getCommandName()                 { return "phylo.diversity";                     }
  27                 string getCommandCategory()             { return "Hypothesis Testing";          }
  28
  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"; }
  33
  34                 int execute();
  35                 void help() { m->mothurOut(getHelpString()); }
  36 private:
  37                 CountTable* ct;
  38                 float freq;
  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
  43
  44         map<string, int> getRootForGroups(Tree* t);
  45                 int readNamesFile();
  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&);
  51
  52 };
  53
  54 /***********************************************************************/
  55 struct phylodivData {
  56     int numIters;
  57         MothurOut* m;
  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;
  63     int increment;
  64     Tree* t;
  65     CountTable* ct;
  66     bool includeRoot;
  67
  68
  69         phylodivData(){}
  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) {
  71         m = mout;
  72         t = tree;
  73         ct = count;
  74         div = cd;
  75         numIters = ni;
  76         sumDiv = csd;
  77         increment = incre;
  78         randomLeaf = crl;
  79         numSampledList = nsl;
  80         rootForGroup = rfg;
  81         }
  82 };
  83
  84 /**************************************************************************************************/
  85 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
  86 #else
  87 static DWORD WINAPI MyPhyloDivThreadFunction(LPVOID lpParam){
  88         phylodivData* pDataArray;
  89         pDataArray = (phylodivData*)lpParam;
  90         try {
  91         int numLeafNodes = pDataArray->randomLeaf.size();
  92                 vector<string> mGroups = pDataArray->m->getGroups();
  93
  94                 for (int l = 0; l < pDataArray->numIters; l++) {
  95             random_shuffle(pDataArray->randomLeaf.begin(), pDataArray->randomLeaf.end());
  96
  97             //initialize counts
  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);
 104             }
 105
 106             for (int j = 0; j < mGroups.size(); j++) {  counts[mGroups[j]] = 0;   }
 107
 108             for(int k = 0; k < numLeafNodes; k++){
 109
 110                 if (pDataArray->m->control_pressed) { return 0; }
 111
 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                 /////////////////////////////////////////////////////////////////////////////////////
 116                 vector<float> br;
 117                 int index = pDataArray->randomLeaf[k];
 118
 119                 vector<string> groups = pDataArray->t->tree[pDataArray->randomLeaf[k]].getGroup();
 120                 br.resize(groups.size(), 0.0);
 121
 122                 //you are a leaf
 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());
 126                     }
 127                 }
 128
 129                 index = pDataArray->t->tree[index].getParent();
 130
 131                 //while you aren't at root
 132                 while(pDataArray->t->tree[index].getParent() != -1){
 133
 134                     if (pDataArray->m->control_pressed) {  return 0; }
 135
 136                     for (int k = 0; k < groups.size(); k++) {
 137
 138                         if (index >= pDataArray->rootForGroup[groups[k]]) { countedBranch[index][groups[k]] = true; } //if you are at this groups "root", then say we are done
 139
 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());
 143                             }
 144                             countedBranch[index][groups[k]] = true;
 145                         }
 146                     }
 147                     index = pDataArray->t->tree[index].getParent();
 148                 }
 149                 /////////////////////////////////////////////////////////////////////////////////////
 150
 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();
 153
 154                 for (int j = 0; j < groups.size(); j++) {
 155
 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;
 162                         }
 163
 164                         if (numSeqsInGroupJ != 0) {     pDataArray->div[groups[j]][(counts[groups[j]]+1)] = pDataArray->div[groups[j]][counts[groups[j]]] + br[j];  }
 165
 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
 168                         }
 169                         counts[groups[j]] += numSeqsInGroupJ;
 170                     }
 171                 }
 172             }
 173
 174
 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];
 179                 }
 180             }
 181
 182         }
 183
 184         return 0;
 185     }
 186         catch(exception& e) {
 187                 pDataArray->m->errorOut(e, "PhyloDiversityCommand", "MyPhyloDivThreadFunction");
 188                 exit(1);
 189         }
 190 }
 191 #endif
 192
 193 #endif
 194