weighted.cpp

   1 /*
   2  *  weighted.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 2/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "weighted.h"
  11
  12 /**************************************************************************************************/
  13
  14 EstOutput Weighted::getValues(Tree* t) {
  15     try {
  16                 globaldata = GlobalData::getInstance();
  17                 int numGroups;
  18                 vector<double> D;
  19
  20                 //if the user has not entered specific groups to analyze then do them all
  21                 if (globaldata->Groups.size() == 0) {
  22                         numGroups = tmap->getNumGroups();
  23                 }else {
  24                         numGroups = globaldata->Groups.size();
  25                 }
  26
  27                 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
  28                 int n = 1;
  29                 for (int i=1; i<numGroups; i++) {
  30                         for (int l = n; l < numGroups; l++) {
  31                                 //initialize weighted scores
  32                                 if (globaldata->Groups.size() == 0) {
  33                                         WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] = 0.0;
  34                                         D.push_back(0.0000);  //initialize a spot in D for each combination
  35                                 }else {
  36                                         WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] = 0.0;
  37                                         D.push_back(0.0000); //initialize a spot in D for each combination
  38                                 }
  39                         }
  40                 }
  41
  42                 data.clear(); //clear out old values
  43
  44                 for(int i=0;i<t->getNumNodes();i++){
  45                         //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC.
  46                         n = 1;
  47                         for (int b=1; b<numGroups; b++) {
  48                                 for (int l = n; l < numGroups; l++) {
  49
  50                                 /********************************************************/
  51                                 //calculate a D value for each group combo
  52                                         for(int v=0;v<t->getNumLeaves();v++){
  53                                                 int index = v;
  54                                                 double sum = 0.0000;
  55
  56                                                 //while you aren't at root
  57                                                 while(t->tree[index].getParent() != -1){
  58
  59                                                         //if you have a BL
  60                                                         if(t->tree[index].getBranchLength() != -1){
  61                                                                 sum += t->tree[index].getBranchLength();
  62                                                         }
  63
  64                                                         index = t->tree[index].getParent();
  65                                                 }
  66
  67                                                 //get last breanch length added
  68                                                 if(t->tree[index].getBranchLength() != -1){
  69                                                                 sum += t->tree[index].getBranchLength();
  70                                                 }
  71
  72                                                 if (globaldata->Groups.size() == 0) {
  73                                                         //is this sum from a sequence which is in one of the users groups
  74                                                         if (inUsersGroups(t->tree[v].getGroup(), tmap->namesOfGroups) == true) {
  75                                                                 //is this sum from a sequence which is in this groupCombo
  76                                                                 if ((t->tree[v].getGroup() == tmap->namesOfGroups[b-1]) || (t->tree[v].getGroup() == tmap->namesOfGroups[l])) {
  77                                                                         sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
  78                                                                         D[n-1] += sum;
  79                                                                 }
  80                                                         }
  81                                                 }else {
  82                                                         //is this sum from a sequence which is in one of the users groups
  83                                                         if (inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
  84                                                                 //is this sum from a sequence which is in this groupCombo
  85                                                                 if ((t->tree[v].getGroup() == globaldata->Groups[b-1]) || (t->tree[v].getGroup() == globaldata->Groups[l])) {
  86                                                                         sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
  87                                                                         D[n-1] += sum;
  88                                                                 }
  89                                                         }
  90                                                 }
  91                                         }
  92                                 /*********************************************************/
  93                                 //calculate a u value for each combo
  94                                         double u;
  95                                         //the user has not entered specific groups
  96                                         if (globaldata->Groups.size() == 0) {
  97                                                 //does this node have descendants from group b-1
  98                                                 it = t->tree[i].pcount.find(tmap->namesOfGroups[b-1]);
  99                                                 //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 100                                                 if (it != t->tree[i].pcount.end()) {
 101                                                         u = (double) t->tree[i].pcount[tmap->namesOfGroups[b-1]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[b-1]];
 102                                                 }else { u = 0.00; }
 103
 104                                                 //does this node have descendants from group l
 105                                                 it = t->tree[i].pcount.find(tmap->namesOfGroups[l]);
 106                                                 //if it does subtract their percentage from u
 107                                                 if (it != t->tree[i].pcount.end()) {
 108                                                         u -= (double) t->tree[i].pcount[tmap->namesOfGroups[l]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[l]];
 109                                                 }
 110
 111                                                 u = abs(u) * t->tree[i].getBranchLength();
 112
 113                                                 //save groupcombs u value
 114                                                 WScore[tmap->namesOfGroups[b-1]+tmap->namesOfGroups[l]] += u;
 115
 116                                         //the user has entered specific groups
 117                                         }else {
 118                                                 //does this node have descendants from group b-1
 119                                                 it = t->tree[i].pcount.find(globaldata->Groups[b-1]);
 120                                                 //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 121                                                 if (it != t->tree[i].pcount.end()) {
 122                                                         u = (double) t->tree[i].pcount[globaldata->Groups[b-1]] / (double) tmap->seqsPerGroup[globaldata->Groups[b-1]];
 123                                                 }else { u = 0.00; }
 124
 125                                                 //does this node have descendants from group l
 126                                                 it = t->tree[i].pcount.find(globaldata->Groups[l]);
 127                                                 //if it does subtract their percentage from u
 128                                                 if (it != t->tree[i].pcount.end()) {
 129                                                         u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]];
 130                                                 }
 131
 132                                                 u = abs(u) * t->tree[i].getBranchLength();
 133
 134                                                 //save groupcombs u value
 135                                                 WScore[globaldata->Groups[b-1]+globaldata->Groups[l]] += u;
 136                                         }
 137                                 /*********************************************************/
 138                                 }
 139                                 n++;
 140                         }
 141                 }
 142
 143                 //calculate weighted score for each group combination
 144                 double UN;
 145                 n = 1;
 146                 for (int i=1; i<numGroups; i++) {
 147                         for (int l = n; l < numGroups; l++) {
 148                                 //the user has not entered specific groups
 149                                 if (globaldata->Groups.size() == 0) {
 150                                         UN = (WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] / D[n-1]);
 151                                 }else {//they have entered specific groups
 152                                         UN = (WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] / D[n-1]);
 153                                 }
 154                                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 155                                 data.push_back(UN);
 156                         }
 157                 }
 158
 159                 return data;
 160         }
 161         catch(exception& e) {
 162                 cout << "Standard Error: " << e.what() << " has occurred in the Weighted class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 163                 exit(1);
 164         }
 165         catch(...) {
 166                 cout << "An unknown error has occurred in the Weighted class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 167                 exit(1);
 168         }
 169
 170 }
 171
 172 /**************************************************************************************************/