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                         float percent = DScore[globaldata->Groups[i]];
  72                         data.push_back(percent);
  73                         //}else {  data.push_back(0.0);  }
  74                 }
  75
  76                 return data;
  77         }
  78         catch(exception& e) {
  79                 m->errorOut(e, "PhyloDiversity", "getValues");
  80                 exit(1);
  81         }
  82 }
  83 /**************************************************************************************************/
  84 void PhyloDiversity::setTotalGroupBranchLengths(Tree* t) {
  85     try {
  86
  87                 groupTotals.clear();
  88
  89                 //initialize group totals
  90                 for (int i=0; i<globaldata->Groups.size(); i++) {               groupTotals[globaldata->Groups[i]] = 0.0;       }
  91
  92
  93                 /********************************************************/
  94                 //calculate a D value for each group
  95                 for(int v=0;v<t->getNumLeaves();v++){
  96
  97                         //is this node from a sequence which is in one of the users groups
  98                         if (inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
  99
 100                                 //calc the branch length
 101                                 int index = v;
 102                                 float sum = 0.0;
 103
 104                                 while(t->tree[index].getParent() != -1){  //while you aren't at root
 105
 106                                         //if you have a BL
 107                                         if(t->tree[index].getBranchLength() != -1){
 108                                                 sum += abs(t->tree[index].getBranchLength());
 109                                         }
 110                                         index = t->tree[index].getParent();
 111                                 }
 112
 113                                 //get last breanch length added
 114                                 if(t->tree[index].getBranchLength() != -1){
 115                                         sum += abs(t->tree[index].getBranchLength());
 116                                 }
 117
 118                                 //account for the names file
 119                                 vector<string> groups = t->tree[v].getGroup();
 120                                 for (int j = 0; j < groups.size(); j++) {
 121                                         int numSeqsInGroupJ = 0;
 122                                         map<string, int>::iterator it;
 123                                         it = t->tree[v].pcount.find(groups[j]);
 124                                         if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group j
 125                                                 numSeqsInGroupJ = it->second;
 126                                         }
 127
 128                                         //add branch length to total for group
 129                                         groupTotals[groups[j]] += (numSeqsInGroupJ * sum);
 130                                 }//end for
 131                         }//end if
 132                 }//end for
 133
 134         }
 135         catch(exception& e) {
 136                 m->errorOut(e, "PhyloDiversity", "setTotalGroupBranchLengths");
 137                 exit(1);
 138         }
 139 }
 140 /**************************************************************************************************/
 141
 142