sharedchao1.cpp

   1 /*
   2  *  sharedchao1.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 1/8/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "sharedchao1.h"
  11
  12 /***********************************************************************/
  13 EstOutput SharedChao1::getValues(vector<SharedRAbundVector*> shared){
  14         try {
  15                 data.resize(1,0);
  16
  17                 vector<int> temp;
  18                 int numGroups = shared.size();
  19                 float Chao = 0.0; float leftvalue, rightvalue;
  20
  21                 // IntNode is defined in mothur.h
  22                 // The tree used here is a binary tree used to represent the f1+++, f+1++, f++1+, f+++1, f11++, f1+1+...
  23                 // combinations required to solve the chao estimator equation for any number of groups.  Conceptually, think
  24                 // of each node as having a 1 and a + value, or for f2 values a 2 and a + value, and 2 pointers to intnodes, and 2 coeffient values.
  25                 // The coeffient value is how many times you chose branch 1 to get to that fvalue.
  26                 // If you choose left you are selecting the 1 or 2 value and right means the + value.  For instance, to find
  27                 // the number of bins that have f1+1+ you would start at the root, go left, right, left, and select the rightvalue.
  28                 // the coeffient is 2.  Note: we only set the coeffient in f2 values.
  29
  30                 //create and initialize trees to 0.
  31                 initialTree(numGroups);
  32
  33                 //loop through vectors calculating the f11, f1A, f2A, f1B, f2B, S12 values
  34                 for (int i = 0; i < shared[0]->size(); i++) {
  35                         //get bin values and calc shared
  36                         bool sharedByAll = true;
  37                         temp.clear();
  38                         for (int j = 0; j < numGroups; j++) {
  39                                 temp.push_back(shared[j]->getAbundance(i));
  40                                 if (temp[j] == 0) { sharedByAll = false; }
  41                         }
  42
  43                         //they are shared
  44                         if (sharedByAll == true) {
  45                                 // cout << "temp = ";
  46                                 // for (int h = 0; h < temp.size(); h++) { cout << temp[h] << " "; } cout << endl;
  47                                 //find f1 and f2values
  48                                 updateTree(temp);
  49                         }
  50                 }
  51
  52
  53                 //cout << "Entering " << endl;
  54                 //calculate chao1, (numleaves-1) because numleaves contains the ++ values.
  55                 bool bias;
  56                 for(int i=0;i<numLeaves;i++){
  57                         if (f2leaves[i]->lvalue == 0) { bias = true;}// break;}
  58                 }
  59
  60                 if(bias){
  61                         for (int i = 0; i < numLeaves; i++) {
  62
  63                                 leftvalue = (float)(f1leaves[i]->lvalue * (f1leaves[i]->lvalue - 1)) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * (f2leaves[i]->lvalue + 1));
  64                                 if (i != (numLeaves-1)) {
  65                                         rightvalue = (float)(f1leaves[i]->rvalue * (f1leaves[i]->rvalue - 1)) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * (f2leaves[i]->rvalue + 1));
  66                                 }else{
  67                                         rightvalue = (float)(f1leaves[i]->rvalue);
  68                                 }
  69                                 Chao += leftvalue + rightvalue;
  70                         }
  71                 }
  72                 else{
  73
  74                         for (int i = 0; i < numLeaves; i++) {
  75
  76                                 leftvalue = (float)(f1leaves[i]->lvalue * f1leaves[i]->lvalue) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * f2leaves[i]->lvalue);
  77                                 if (i != (numLeaves-1)) {
  78                                         rightvalue = (float)(f1leaves[i]->rvalue * f1leaves[i]->rvalue) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * f2leaves[i]->rvalue);
  79                                 }else{
  80                                         rightvalue = (float)(f1leaves[i]->rvalue);
  81                                 }
  82                                 Chao += leftvalue + rightvalue;
  83                         }
  84                 }
  85
  86                 for (int i = 0; i < numNodes; i++) {
  87                         delete f1leaves[i];
  88                         delete f2leaves[i];
  89                 }
  90
  91         //      cout << "exiting " << endl;
  92                 data[0] = Chao;
  93                 return data;
  94         }
  95         catch(exception& e) {
  96                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  97                 exit(1);
  98         }
  99         catch(...) {
 100                 cout << "An unknown error has occurred in the SharedChao1 class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 101                 exit(1);
 102         }
 103 }
 104
 105 /***********************************************************************/
 106 //builds trees structure with n leaf nodes initialized to 0.
 107 void SharedChao1::initialTree(int n) {
 108         try {
 109                 // (2^n) / 2. Divide by 2 because each leaf node contains 2 values. One for + and one for 1 or 2.
 110                 numLeaves = pow(2, (float)n) / 2;
 111                 numNodes = 2*numLeaves - 1;
 112                 int countleft = 0;
 113                 int countright = 1;
 114
 115                 f1leaves.resize(numNodes);
 116                 f2leaves.resize(numNodes);
 117
 118                 //initialize leaf values
 119                 for (int i = 0; i < numLeaves; i++) {
 120                         f1leaves[i] = new IntNode;
 121                         f1leaves[i]->lvalue = 0;
 122                         f1leaves[i]->rvalue = 0;
 123                         f1leaves[i]->left = NULL;
 124                         f1leaves[i]->right = NULL;
 125
 126                         f2leaves[i] = new IntNode;
 127                         f2leaves[i]->lvalue = 0;
 128                         f2leaves[i]->rvalue = 0;
 129                         f2leaves[i]->left = NULL;
 130                         f2leaves[i]->right = NULL;
 131                 }
 132
 133                 //set pointers to children
 134                 for (int j = numLeaves; j < numNodes; j++) {
 135                         f1leaves[j] = new IntNode;
 136                         f1leaves[j]->left = f1leaves[countleft];
 137                         f1leaves[j]->right = f1leaves[countright];
 138
 139                         f2leaves[j] = new IntNode;
 140                         f2leaves[j]->left = f2leaves[countleft];
 141                         f2leaves[j]->right =f2leaves[countright];
 142
 143                         countleft = countleft + 2;
 144                         countright = countright + 2;
 145                 }
 146
 147                 //point to root
 148                 f1root = f1leaves[numNodes-1];
 149
 150                 //point to root
 151                 f2root = f2leaves[numNodes-1];
 152
 153                 //set coeffients
 154                 setCoef(f2root, 0);
 155         }
 156         catch(exception& e) {
 157                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 158                 exit(1);
 159         }
 160         catch(...) {
 161                 cout << "An unknown error has occurred in the SharedChao1 class function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 162                 exit(1);
 163         }
 164 }
 165
 166 /***********************************************************************/
 167 //take vector containing the abundance info. for a bin and updates trees.
 168 void SharedChao1::updateTree(vector<int> bin) {
 169         try {
 170                 updateBranchf1(f1root, bin, 0);
 171                 updateBranchf2(f2root, bin, 0);
 172         }
 173         catch(exception& e) {
 174                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 175                 exit(1);
 176         }
 177         catch(...) {
 178                 cout << "An unknown error has occurred in the SharedChao1 class function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 179                 exit(1);
 180         }
 181 }
 182
 183 /***********************************************************************/
 184 void SharedChao1::updateBranchf1(IntNode* node, vector<int> bin, int index) {
 185         try {
 186                 //if you have more than one group
 187                 if (index == (bin.size()-1)) {
 188                         if (bin[index] == 1) { node->lvalue++; node->rvalue++; }
 189                         else { node->rvalue++;  }
 190                 }else {
 191                         if (bin[index] == 1) {
 192                                 //follow path as if you are 1
 193                                 updateBranchf1(node->left, bin, index+1);
 194                         }
 195                         //follow path as if you are +
 196                         updateBranchf1(node->right, bin, index+1);
 197                 }
 198         }
 199         catch(exception& e) {
 200                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 201                 exit(1);
 202         }
 203         catch(...) {
 204                 cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 205                 exit(1);
 206         }
 207 }
 208
 209 /***********************************************************************/
 210 void SharedChao1::updateBranchf2(IntNode* node, vector<int> bin, int index) {
 211         try {
 212                 //if you have more than one group
 213                 if (index == (bin.size()-1)) {
 214                         if (bin[index] == 2) { node->lvalue++; node->rvalue++; }
 215                         else { node->rvalue++;  }
 216                 }else {
 217                         if (bin[index] == 2) {
 218                                 //follow path as if you are 1
 219                                 updateBranchf2(node->left, bin, index+1);
 220                         }
 221                         //follow path as if you are +
 222                         updateBranchf2(node->right, bin, index+1);
 223                 }
 224         }
 225         catch(exception& e) {
 226                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 227                 exit(1);
 228         }
 229         catch(...) {
 230                 cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 231                 exit(1);
 232         }
 233 }
 234
 235 /***********************************************************************/
 236 void SharedChao1::setCoef(IntNode* node, int coef) {
 237         try {
 238                 if (node->left != NULL) {
 239                         setCoef(node->left, coef+1);
 240                         setCoef(node->right, coef);
 241                 }else {
 242                         node->lcoef = coef+1;
 243                         node->rcoef = coef;
 244                 }
 245         }
 246         catch(exception& e) {
 247                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 248                 exit(1);
 249         }
 250         catch(...) {
 251                 cout << "An unknown error has occurred in the SharedChao1 class function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 252                 exit(1);
 253         }
 254 }
 255
 256 /***********************************************************************/
 257 //for debugging purposes
 258 void SharedChao1::printTree() {
 259
 260         cout << "F1 leaves" << endl;
 261         printBranch(f1root);
 262
 263         cout << "F2 leaves" << endl;
 264         printBranch(f2root);
 265
 266
 267 }
 268 /*****************************************************************/
 269 void SharedChao1::printBranch(IntNode* node) {
 270         try {
 271
 272                 // you are not a leaf
 273                 if (node->left != NULL) {
 274                         printBranch(node->left);
 275                         printBranch(node->right);
 276                 }else { //you are a leaf
 277                         cout << node->lvalue << endl;
 278                         cout << node->rvalue << endl;
 279                 }
 280
 281         }
 282         catch(exception& e) {
 283                 cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 284                 exit(1);
 285         }
 286         catch(...) {
 287                 cout << "An unknown error has occurred in the Tree class function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 288                 exit(1);
 289         }
 290 }
 291
 292 /*****************************************************************/
 293
 294
 295
 296