uvest.cpp

   1 /*
   2  *  uvest.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 "uvest.h"
  11
  12 /***********************************************************************/
  13 //This is used by SharedJAbund and SharedSorAbund
  14 EstOutput UVEst::getUVest(vector<SharedRAbundVector*> shared) {
  15         try {
  16                 EstOutput results;
  17                 results.resize(2,0);
  18
  19                 int S12, Atotal, Btotal, f1A, f2A, f1B, f2B, sumSharedA, sumSharedB, sumSharedA1, sumSharedB1, tempA, tempB;
  20                 S12 = 0; Atotal = 0; Btotal = 0; f1A = 0; f2A = 0; f1B = 0; f2B = 0; sumSharedA = 0; sumSharedB = 0; sumSharedA1 = 0; sumSharedB1 = 0;
  21                 float Upart1, Upart2, Upart3, Vpart1, Vpart2, Vpart3, Uest, Vest;
  22                 Upart1 = 0.0; Upart2 = 0.0; Upart3 = 0.0; Vpart1 = 0.0; Vpart2 = 0.0; Vpart3 = 0.0;
  23
  24                 /*Xi, Yi = abundance of the ith shared OTU in A and B
  25                 ntotal, mtotal = total number of sequences sampled in A and B
  26                 I(•) = if the argument, •, is true then I(•) is 1; otherwise it is 0.
  27                 sumSharedA = the sum of all shared otus in A
  28                 sumSharedB = the sum of all shared otus in B
  29                 sumSharedA1 = the sum of all shared otus in A where B = 1
  30                 sumSharedB1 = the sum of all shared otus in B where A = 1 */
  31
  32                 for (int i = 0; i < shared[0]->getNumBins(); i++) {
  33                         //store in temps to avoid multiple repetitive function calls
  34                         tempA = shared[0]->getAbundance(i);
  35                         tempB = shared[1]->getAbundance(i);
  36
  37                         Atotal += tempA;
  38                         Btotal += tempB;
  39
  40                         if ((tempA != 0) && (tempB != 0)) {//they are shared
  41                                 sumSharedA += tempA;
  42                                 sumSharedB += tempB;
  43
  44                                 //does A have one or two
  45                                 if (tempA == 1)                 { f1A++; sumSharedB1 += tempB;}
  46                                 else if (tempA == 2)    { f2A++; }
  47
  48                                 //does B have one or two
  49                                 if (tempB == 1)                 { f1B++; sumSharedA1 += tempA;}
  50                                 else if (tempB == 2)    { f2B++; }
  51                         }
  52                 }
  53
  54                 Upart1 = sumSharedA / (float) Atotal;
  55                 Upart2 = ((Btotal - 1) * f1B) / (float) (Btotal * 2 * f2B);
  56                 Upart3 = sumSharedA1 / (float) Atotal;
  57
  58                 if (isnan(Upart1) || isinf(Upart1)) { Upart1 = 0; }
  59                 if (isnan(Upart2) || isinf(Upart2)) { Upart2 = 0; }
  60                 if (isnan(Upart3) || isinf(Upart3)) { Upart3 = 0; }
  61
  62                 Uest = Upart1 + (Upart2 * Upart3);
  63
  64                 Vpart1 = sumSharedB / (float) Btotal;
  65                 Vpart2 = ((Atotal - 1) * f1A) / (float) (Atotal * 2 * f2A);
  66                 Vpart3 = sumSharedB1 / (float) Btotal;
  67
  68                 if (isnan(Vpart1) || isinf(Vpart1)) { Vpart1 = 0; }
  69                 if (isnan(Vpart2) || isinf(Vpart2)) { Vpart2 = 0; }
  70                 if (isnan(Vpart3) || isinf(Vpart3)) { Vpart3 = 0; }
  71
  72                 Vest = Vpart1 + (Vpart2 * Vpart3);
  73
  74                 if (Uest > 1) { Uest = 1; }
  75                 if (Vest > 1) { Vest = 1; }
  76
  77                 results[0] = Uest;
  78                 results[1] = Vest;
  79
  80                 return results;
  81         }
  82         catch(exception& e) {
  83                 m->errorOut(e, "UVEst", "getUVest");
  84                 exit(1);
  85         }
  86 }
  87
  88 /***********************************************************************/