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(SharedRAbundVector* sharedA, SharedRAbundVector* sharedB){
  14         try {
  15                 data.resize(1,0);
  16
  17                 int f11, f1A, f2A, f1B, f2B, S12, tempA, tempB;
  18                 f11 = 0; f1A = 0; f2A = 0; f1B = 0; f2B = 0; S12 = 0;
  19                 float ChaoAB, part1, part2, part3;
  20
  21         /*      f11 = number of shared OTUs with one observed individual in A and B
  22                 f1A, f2A = number of shared OTUs with one or two individuals observed in A
  23                 f1B, f2B = number of shared OTUs with one or two individuals observed in B
  24                 S12 = number of shared OTUs in A and B  */
  25
  26                 //loop through vectors calculating the f11, f1A, f2A, f1B, f2B, S12 values
  27                 for (int i = 0; i< sharedA->size(); i++) {
  28                         tempA = sharedA->getAbundance(i); //store in temps to avoid calling getAbundance multiple times
  29                         tempB = sharedB->getAbundance(i);
  30                         if ((tempA != 0) && (tempB != 0)) {//they are shared
  31                                 S12++; //they are shared
  32                                 //do both A and B have one
  33                                 if ((tempA == 1) && (tempB == 1)) { f11++; }
  34
  35                                 //does A have one or two
  36                                 if (tempA == 1)                 { f1A++; }
  37                                 else if (tempA == 2)    { f2A++; }
  38
  39                                 //does B have one or two
  40                                 if (tempB == 1)                 { f1B++; }
  41                                 else if (tempB == 2)    { f2B++; }
  42
  43                         }
  44                 }
  45
  46                 //checks for divide by zero error
  47                 if ((f2A == 0) || (f2B == 0)) {
  48                         part1 = ((float)(f1A*f1B)/(float)(4*(f2A+1)*(f2B+1)));
  49                         part2 = ((float)(f1A*(f1A-1))/(float)(2*f2A+2));
  50                         part3 = ((float)(f1B*(f1B-1))/(float)(2*f2B+2));
  51                 }else {
  52                         part1 = ((float)(f1A*f1B)/(float)(4*f2A*f2B));
  53                         part2 = ((float)(f1A*f1A)/(float)(2*f2A));
  54                         part3 = ((float)(f1B*f1B)/(float)(2*f2B));
  55                 }
  56
  57                 ChaoAB = (float)S12 + (float)(f11*part1) + part2 + part3;
  58                 data[0] = ChaoAB;
  59
  60                 return data;
  61         }
  62         catch(exception& e) {
  63                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  64                 exit(1);
  65         }
  66         catch(...) {
  67                 cout << "An unknown error has occurred in the SharedChao1 class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  68                 exit(1);
  69         }
  70
  71 }
  72
  73 /***********************************************************************/