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                 numGroups = globaldata->Groups.size();
  21
  22                 vector<double> sums = getBranchLengthSums(t);
  23
  24                 if (m->control_pressed) { return data; }
  25
  26                 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
  27                 int count = 0;
  28                 for (int i=0; i<numGroups; i++) {
  29                         for (int l = i+1; l < numGroups; l++) {
  30                                 //initialize weighted scores
  31                                 WScore[globaldata->Groups[i]+globaldata->Groups[l]] = 0.0;
  32
  33                                 vector<string> groups; groups.push_back(globaldata->Groups[i]); groups.push_back(globaldata->Groups[l]);
  34
  35                                 D.push_back(0.0000); //initialize a spot in D for each combination
  36
  37                                 //adding the wieghted sums from group i
  38                                 for (int j = 0; j < t->groupNodeInfo[groups[0]].size(); j++) { //the leaf nodes that have seqs from group i
  39                                         map<string, int>::iterator it = t->tree[t->groupNodeInfo[groups[0]][j]].pcount.find(groups[0]);
  40                                         int numSeqsInGroupI = it->second;
  41
  42                                         double weightedSum = ((numSeqsInGroupI * sums[t->groupNodeInfo[groups[0]][j]]) / (double)tmap->seqsPerGroup[groups[0]]);
  43
  44                                         D[count] += weightedSum;
  45                                 }
  46
  47                                 //adding the wieghted sums from group l
  48                                 for (int j = 0; j < t->groupNodeInfo[groups[1]].size(); j++) { //the leaf nodes that have seqs from group l
  49                                         map<string, int>::iterator it = t->tree[t->groupNodeInfo[groups[1]][j]].pcount.find(groups[1]);
  50                                         int numSeqsInGroupL = it->second;
  51
  52                                         double weightedSum = ((numSeqsInGroupL * sums[t->groupNodeInfo[groups[1]][j]]) / (double)tmap->seqsPerGroup[groups[1]]);
  53
  54                                         D[count] += weightedSum;
  55                                 }
  56
  57                                 /********************************************************
  58                                 //calculate a D value for each group combo
  59                                 for(int v=0;v<t->getNumLeaves();v++){
  60
  61                                         if (m->control_pressed) { return data; }
  62
  63                                         //is this sum from a sequence which is in one of the users groups
  64                                         if (m->inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
  65                                                 //is this sum from a sequence which is in this groupCombo
  66                                                 if (m->inUsersGroups(t->tree[v].getGroup(), groups)) {
  67                                                         int numSeqsInGroupI, numSeqsInGroupL;
  68
  69                                                         map<string, int>::iterator it;
  70                                                         it = t->tree[v].pcount.find(groups[0]);
  71                                                         if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i
  72                                                                 numSeqsInGroupI = it->second;
  73                                                         }else{ numSeqsInGroupI = 0; }
  74
  75                                                         it = t->tree[v].pcount.find(groups[1]);
  76                                                         if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l
  77                                                                 numSeqsInGroupL = it->second;
  78                                                         }else{ numSeqsInGroupL = 0; }
  79
  80                                                         double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]);
  81
  82                                                         //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
  83
  84                                                         D[count] += weightedSum;
  85                                                 }
  86                                         }
  87                                 }
  88                                 /*********************************************************/
  89                                 count++;
  90                         }
  91                 }
  92
  93                 data.clear(); //clear out old values
  94
  95                 for(int i=0;i<t->getNumNodes();i++){
  96                         //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC.
  97                         for (int b=0; b<numGroups; b++) {
  98                                 for (int l = b+1; l < numGroups; l++) {
  99
 100                                         if (m->control_pressed) { return data; }
 101
 102                                         //calculate a u value for each combo
 103                                         double u;
 104                                         //does this node have descendants from group b-1
 105                                         it = t->tree[i].pcount.find(globaldata->Groups[b]);
 106                                         //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 107                                         if (it != t->tree[i].pcount.end()) {
 108                                                 u = (double) t->tree[i].pcount[globaldata->Groups[b]] / (double) tmap->seqsPerGroup[globaldata->Groups[b]];
 109                                         }else { u = 0.00; }
 110
 111                                         //does this node have descendants from group l
 112                                         it = t->tree[i].pcount.find(globaldata->Groups[l]);
 113                                         //if it does subtract their percentage from u
 114                                         if (it != t->tree[i].pcount.end()) {
 115                                                 u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]];
 116                                         }
 117
 118                                         u = abs(u * t->tree[i].getBranchLength());
 119
 120                                         //save groupcombs u value
 121                                         WScore[globaldata->Groups[b]+globaldata->Groups[l]] += u;
 122                                 /*********************************************************/
 123                                 }
 124                         }
 125                 }
 126
 127                 //calculate weighted score for each group combination
 128                 double UN;
 129                 count = 0;
 130                 for (int i=0; i<numGroups; i++) {
 131                         for (int l = i+1; l < numGroups; l++) {
 132                                 UN = (WScore[globaldata->Groups[i]+globaldata->Groups[l]] / D[count]);
 133
 134                                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 135                                 data.push_back(UN);
 136                                 count++;
 137                         }
 138                 }
 139                 return data;
 140         }
 141         catch(exception& e) {
 142                 m->errorOut(e, "Weighted", "getValues");
 143                 exit(1);
 144         }
 145 }
 146
 147 /**************************************************************************************************/
 148 EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
 149  try {
 150                 globaldata = GlobalData::getInstance();
 151
 152                 data.clear(); //clear out old values
 153
 154                 vector<double> sums = getBranchLengthSums(t);
 155
 156                 if (m->control_pressed) { return data; }
 157
 158                 //initialize weighted score
 159                 WScore[(groupA+groupB)] = 0.0;
 160                 double D = 0.0;
 161
 162                 vector<string> groups; groups.push_back(groupA); groups.push_back(groupB);
 163
 164                 //adding the wieghted sums from groupA
 165                 for (int j = 0; j < t->groupNodeInfo[groupA].size(); j++) { //the leaf nodes that have seqs from group i
 166                         map<string, int>::iterator it = t->tree[j].pcount.find(groupA);
 167                         int numSeqsInGroupA = it->second;
 168
 169                         double weightedSum = ((numSeqsInGroupA * sums[j]) / (double)tmap->seqsPerGroup[groupA]);
 170
 171                         D += weightedSum;
 172                 }
 173
 174                 //adding the wieghted sums from groupB
 175                 for (int j = 0; j < t->groupNodeInfo[groupB].size(); j++) { //the leaf nodes that have seqs from group l
 176                         map<string, int>::iterator it = t->tree[j].pcount.find(groupB);
 177                         int numSeqsInGroupB = it->second;
 178
 179                         double weightedSum = ((numSeqsInGroupB * sums[j]) / (double)tmap->seqsPerGroup[groupB]);
 180
 181                         D += weightedSum;
 182                 }
 183
 184
 185                 /********************************************************
 186                 //calculate a D value for the group combo
 187                 for(int v=0;v<t->getNumLeaves();v++){
 188                         if (m->control_pressed) { return data; }
 189
 190                         int index = v;
 191                         double sum = 0.0000;
 192
 193                         //while you aren't at root
 194                         while(t->tree[index].getParent() != -1){
 195
 196                                 //if you have a BL
 197                                 if(t->tree[index].getBranchLength() != -1){
 198                                         sum += abs(t->tree[index].getBranchLength());
 199                                 }
 200                                 index = t->tree[index].getParent();
 201                         }
 202
 203                         //get last breanch length added
 204                         if(t->tree[index].getBranchLength() != -1){
 205                                 sum += abs(t->tree[index].getBranchLength());
 206                         }
 207
 208                         if (m->inUsersGroups(t->tree[v].getGroup(), groups)) {
 209                                 int numSeqsInGroupI, numSeqsInGroupL;
 210
 211                                 map<string, int>::iterator it;
 212                                 it = t->tree[v].pcount.find(groups[0]);
 213                                 if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i
 214                                         numSeqsInGroupI = it->second;
 215                                 }else{ numSeqsInGroupI = 0; }
 216
 217                                 it = t->tree[v].pcount.find(groups[1]);
 218                                 if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l
 219                                         numSeqsInGroupL = it->second;
 220                                 }else{ numSeqsInGroupL = 0; }
 221
 222                                 double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]);
 223
 224                                 //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
 225
 226                                 D += weightedSum;
 227                         }
 228                 }
 229                 /********************************************************/
 230
 231                 //calculate u for the group comb
 232                 for(int i=0;i<t->getNumNodes();i++){
 233
 234                         if (m->control_pressed) { return data; }
 235
 236                         double u;
 237                         //does this node have descendants from groupA
 238                         it = t->tree[i].pcount.find(groupA);
 239                         //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 240                         if (it != t->tree[i].pcount.end()) {
 241                                 u = (double) t->tree[i].pcount[groupA] / (double) tmap->seqsPerGroup[groupA];
 242                         }else { u = 0.00; }
 243
 244
 245                         //does this node have descendants from group l
 246                         it = t->tree[i].pcount.find(groupB);
 247                         //if it does subtract their percentage from u
 248                         if (it != t->tree[i].pcount.end()) {
 249                                 u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB];
 250                         }
 251
 252                         u = abs(u * t->tree[i].getBranchLength());
 253
 254                         //save groupcombs u value
 255                         WScore[(groupA+groupB)] += u;
 256                 }
 257
 258                 /********************************************************/
 259
 260                 //calculate weighted score for the group combination
 261                 double UN;
 262                 UN = (WScore[(groupA+groupB)] / D);
 263
 264                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 265                 data.push_back(UN);
 266
 267                 return data;
 268         }
 269         catch(exception& e) {
 270                 m->errorOut(e, "Weighted", "getValues");
 271                 exit(1);
 272         }
 273 }
 274 /**************************************************************************************************/
 275 vector<double> Weighted::getBranchLengthSums(Tree* t) {
 276         try {
 277
 278                 vector<double> sums;
 279
 280                 for(int v=0;v<t->getNumLeaves();v++){
 281
 282                         if (m->control_pressed) { return sums; }
 283
 284                         int index = v;
 285                         double sum = 0.0000;
 286
 287                         //while you aren't at root
 288                         while(t->tree[index].getParent() != -1){
 289
 290                                 //if you have a BL
 291                                 if(t->tree[index].getBranchLength() != -1){
 292                                         sum += abs(t->tree[index].getBranchLength());
 293                                 }
 294                                 index = t->tree[index].getParent();
 295                         }
 296
 297                         //get last breanch length added
 298                         if(t->tree[index].getBranchLength() != -1){
 299                                 sum += abs(t->tree[index].getBranchLength());
 300                         }
 301
 302                         sums.push_back(sum);
 303                 }
 304
 305                 return sums;
 306         }
 307         catch(exception& e) {
 308                 m->errorOut(e, "Weighted", "getBranchLengthSums");
 309                 exit(1);
 310         }
 311 }
 312 /**************************************************************************************************/
 313
 314
 315