sharedthetayc.cpp

   1 /*
   2  *  sharedthetayc.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 "sharedthetayc.h"
  11
  12 /***********************************************************************/
  13 EstOutput ThetaYC::getValues(vector<SharedRAbundVector*> shared) {
  14         try {
  15                 data.resize(3,0.0000);
  16
  17                 double Atotal = 0;
  18                 double Btotal = 0;
  19                 double thetaYC = 0;
  20                 double pi = 0;
  21                 double qi = 0;
  22                 double a = 0;
  23                 double b = 0;
  24                 double d = 0;
  25
  26                 double sumPcubed = 0;
  27                 double sumQcubed = 0;
  28                 double sumPQsq = 0;
  29                 double sumPsqQ = 0;
  30
  31                 //get the total values we need to calculate the theta denominator sums
  32                 for (int i = 0; i < shared[0]->getNumBins(); i++) {
  33                         //store in temps to avoid multiple repetitive function calls
  34                         Atotal += (double)shared[0]->getAbundance(i);
  35                         Btotal += (double)shared[1]->getAbundance(i);
  36                 }
  37
  38                 //calculate the theta denominator sums
  39                 for (int j = 0; j < shared[0]->getNumBins(); j++) {
  40                         //store in temps to avoid multiple repetitive function calls
  41                         pi = shared[0]->getAbundance(j) / Atotal;
  42                         qi = shared[1]->getAbundance(j) / Btotal;
  43
  44                         a += pi * pi;
  45                         b += qi * qi;
  46                         d += pi * qi;
  47
  48                         sumPcubed += pi * pi * pi;
  49                         sumQcubed += qi * qi * qi;
  50                         sumPQsq += pi * qi * qi;
  51                         sumPsqQ += pi * pi * qi;
  52                 }
  53
  54                 thetaYC = d / (a + b - d);
  55
  56                 if (isnan(thetaYC) || isinf(thetaYC)) { thetaYC = 0; }
  57
  58                 double varA = 4 / Atotal * (sumPcubed - a * a);
  59                 double varB = 4 / Btotal * (sumQcubed - b * b);
  60                 double varD = sumPQsq / Atotal + sumPsqQ / Btotal - d * d * (1/Atotal + 1/Btotal);
  61                 double covAD = 2 / Atotal * (sumPsqQ - a * d);
  62                 double covBD = 2 / Btotal * (sumPQsq - b* d);
  63
  64                 double varT = d * d * (varA + varB) / pow(a + b - d, (double)4.0) + pow(a+b, (double)2.0) * varD / pow(a+b-d, (double)4.0)
  65                                                 - 2.0 * (a + b) * d / pow(a + b - d, (double)4.0) * (covAD + covBD);
  66
  67                 double ci = 1.95 * sqrt(varT);
  68
  69                 data[0] = thetaYC;
  70                 data[1] = thetaYC - ci;
  71                 data[2] = thetaYC + ci;
  72
  73                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
  74                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
  75                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
  76
  77                 data[0] = 1.0 - data[0];
  78         double hold = data[1];
  79         data[1] = 1.0 - data[2];
  80         data[2] = 1.0 - hold;
  81
  82                 return data;
  83         }
  84         catch(exception& e) {
  85                 m->errorOut(e, "ThetaYC", "getValues");
  86                 exit(1);
  87         }
  88 }
  89
  90 /***********************************************************************/