phylodiversity.cpp

   1 /*
   2  *  phylodiversity.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 4/30/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "phylodiversity.h"
  11
  12 /**************************************************************************************************/
  13 EstOutput PhyloDiversity::getValues(Tree* t, vector<int> treeNodes) {
  14     try {
  15
  16                 map<string, float> DScore;
  17                 float totalLength = 0.0;
  18                 data.clear();
  19
  20                 //initialize Dscore
  21                 for (int i=0; i<globaldata->Groups.size(); i++) {               DScore[globaldata->Groups[i]] = 0.0;    }
  22
  23                 /********************************************************/
  24                 //calculate a D value for each group
  25                 for(int v=0;v<treeNodes.size();v++){
  26
  27                         if (m->control_pressed) { return data; }
  28
  29                         //is this node from a sequence which is in one of the users groups
  30                         if (inUsersGroups(t->tree[treeNodes[v]].getGroup(), globaldata->Groups) == true) {
  31
  32                                 //calc the branch length
  33                                 //while you aren't at root
  34                                 float sum = 0.0;
  35                                 int index = treeNodes[v];
  36
  37                                 while(t->tree[index].getParent() != -1){
  38
  39                                         //if you have a BL
  40                                         if(t->tree[index].getBranchLength() != -1){
  41                                                 sum += abs(t->tree[index].getBranchLength());
  42                                         }
  43                                         index = t->tree[index].getParent();
  44                                 }
  45
  46                                 //get last breanch length added
  47                                 if(t->tree[index].getBranchLength() != -1){
  48                                         sum += abs(t->tree[index].getBranchLength());
  49                                 }
  50
  51                                 //for each group in the groups update the total branch length accounting for the names file
  52                                 vector<string> groups = t->tree[treeNodes[v]].getGroup();
  53                                 for (int j = 0; j < groups.size(); j++) {
  54                                         int numSeqsInGroupJ = 0;
  55                                         map<string, int>::iterator it;
  56                                         it = t->tree[treeNodes[v]].pcount.find(groups[j]);
  57                                         if (it != t->tree[treeNodes[v]].pcount.end()) { //this leaf node contains seqs from group j
  58                                                 numSeqsInGroupJ = it->second;
  59                                         }
  60
  61                                         //add branch length to total for group
  62                                         DScore[groups[j]] += (numSeqsInGroupJ * sum);
  63                                 }
  64                         }
  65                 }
  66
  67
  68                 for (int i=0; i<globaldata->Groups.size(); i++) {
  69                         if (groupTotals[globaldata->Groups[i]] != 0.0) { //avoid divide by zero error
  70                                 float percent = DScore[globaldata->Groups[i]] / groupTotals[globaldata->Groups[i]];
  71                                 data.push_back(percent);
  72                         }else {  data.push_back(0.0);  }
  73                 }
  74
  75                 return data;
  76         }
  77         catch(exception& e) {
  78                 m->errorOut(e, "PhyloDiversity", "getValues");
  79                 exit(1);
  80         }
  81 }
  82 /**************************************************************************************************/
  83 void PhyloDiversity::setTotalGroupBranchLengths(Tree* t) {
  84     try {
  85
  86                 groupTotals.clear();
  87
  88                 //initialize group totals
  89                 for (int i=0; i<globaldata->Groups.size(); i++) {               groupTotals[globaldata->Groups[i]] = 0.0;       }
  90
  91
  92                 /********************************************************/
  93                 //calculate a D value for each group
  94                 for(int v=0;v<t->getNumLeaves();v++){
  95
  96                         //is this node from a sequence which is in one of the users groups
  97                         if (inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
  98
  99                                 //calc the branch length
 100                                 int index = v;
 101                                 float sum = 0.0;
 102
 103                                 while(t->tree[index].getParent() != -1){  //while you aren't at root
 104
 105                                         //if you have a BL
 106                                         if(t->tree[index].getBranchLength() != -1){
 107                                                 sum += abs(t->tree[index].getBranchLength());
 108                                         }
 109                                         index = t->tree[index].getParent();
 110                                 }
 111
 112                                 //get last breanch length added
 113                                 if(t->tree[index].getBranchLength() != -1){
 114                                         sum += abs(t->tree[index].getBranchLength());
 115                                 }
 116
 117                                 //account for the names file
 118                                 vector<string> groups = t->tree[v].getGroup();
 119                                 for (int j = 0; j < groups.size(); j++) {
 120                                         int numSeqsInGroupJ = 0;
 121                                         map<string, int>::iterator it;
 122                                         it = t->tree[v].pcount.find(groups[j]);
 123                                         if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group j
 124                                                 numSeqsInGroupJ = it->second;
 125                                         }
 126
 127                                         //add branch length to total for group
 128                                         groupTotals[groups[j]] += (numSeqsInGroupJ * sum);
 129                                 }//end for
 130                         }//end if
 131                 }//end for
 132
 133         }
 134         catch(exception& e) {
 135                 m->errorOut(e, "PhyloDiversity", "setTotalGroupBranchLengths");
 136                 exit(1);
 137         }
 138 }
 139 /**************************************************************************************************/
 140
 141