coverage.cpp

   1 /*
   2  *  coverage.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 3/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 /* This class library coverage at the given distances of the Cramer-von Mises statistic.
  11         you may refer to the "Integration of Microbial Ecology and Statistics: A Test To Compare Gene Libraries"
  12         paper in Applied and Environmental Microbiology, Sept. 2004, p. 5485-5492 0099-2240/04/$8.00+0
  13         DOI: 10.1128/AEM.70.9.5485-5492.2004 Copyright 2004 American Society for Microbiology for more information. */
  14
  15
  16 #include "coverage.h"
  17
  18 //**********************************************************************************************************************
  19 Coverage::Coverage() {
  20                 globaldata = GlobalData::getInstance();
  21                 numUserGroups = globaldata->Groups.size();
  22                 numGroups = globaldata->gGroupmap->getNumGroups();
  23 }
  24
  25 //**********************************************************************************************************************
  26 void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >& data, vector<float> dist) {
  27         try {
  28                 vector<float> min;
  29                 vector<string> groups;
  30
  31                 //initialize data
  32                 data.resize(dist.size());
  33                 for (int l = 0; l < data.size(); l++) {
  34                         data[l].resize(numGroups);
  35                         for (int k = 0; k < data[l].size(); k++) {
  36                                 data[l][k].push_back(0.0);
  37                         }
  38                 }
  39
  40                 /**************************************/
  41                 //get the minimums for each comparision
  42                 /**************************************/
  43                 int count = 0;
  44
  45                 for (int i = 0; i < numGroups; i++) {
  46                         for (int j = 0; j < numGroups; j++) {
  47
  48                                 //is this "box" one the user wants analyzed?
  49                                 if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) {
  50
  51                                         min = matrix->getMins(count);  //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
  52
  53                                         //find the coverage at this distance
  54                                         sort(min.begin(), min.end());
  55
  56                                         //loop through each distance and fill data
  57                                         for (int k = 0; k < data.size(); k++) {
  58
  59                                                 int index = -1;
  60                                                 //find index in min where value is higher than d
  61                                                 for (int m = 0; m < min.size(); m++) {
  62                                                         if (min[m] > dist[k])   { index = m; break;     }
  63                                                 }
  64
  65                                                 // if you don't find one than all the mins are less than d
  66                                                 if (index == -1) { index = min.size();  }
  67
  68                                                 //save value in data
  69                                                 data[k][i][j] = 1.0 - ((min.size()-index)/(float)min.size());
  70
  71                                         }
  72                                 }
  73                                 count++;
  74                         }
  75                 }
  76
  77         }
  78         catch(exception& e) {
  79                 cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  80                 exit(1);
  81         }
  82         catch(...) {
  83                 cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  84                 exit(1);
  85         }
  86 }
  87
  88 //**********************************************************************************************************************
  89 void Coverage::getValues(FullMatrix* matrix, vector< vector< vector<float> > >& data, vector<float> dist, string mode) {
  90         try {
  91                 vector<float> min1;
  92                 vector<float> min2;
  93                 vector<string> groups;
  94
  95                 //initialize data
  96                 data.resize(dist.size());
  97                 for (int l = 0; l < data.size(); l++) {
  98                         data[l].resize(numGroups);
  99                         for (int k = 0; k < data[l].size(); k++) {
 100                                 data[l][k].push_back(0.0);
 101                         }
 102                 }
 103
 104                 int count = 0;
 105                 int count2 = 0;
 106
 107                 //for each box
 108                 for (int i = 0; i < numGroups; i++) {
 109                         for (int j = 0; j < numGroups; j++) {
 110
 111                                 if (i != j) {
 112                                         //is this "box" one the user wants analyzed?
 113                                         if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) {
 114
 115                                                 matrix->shuffle(globaldata->gGroupmap->namesOfGroups[i], globaldata->gGroupmap->namesOfGroups[j]);
 116
 117                                                 min1 = matrix->getMins(count);  //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
 118                                                 min2 = matrix->getMins(count2);  //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
 119
 120                                                 //find the coverage at this distance
 121                                                 sort(min1.begin(), min1.end());
 122
 123                                                 //find the coverage at this distance
 124                                                 sort(min2.begin(), min2.end());
 125
 126                                                 float distDiff = 0;
 127
 128                                                 //loop through each distance and fill data
 129                                                 for (int k = 0; k < data.size(); k++) {
 130                                                         //****** coverage of AA **********//
 131                                                         int index = -1;
 132                                                         //find index in min where value is higher than d
 133                                                         for (int m = 0; m < min1.size(); m++) {
 134                                                                         if (min1[m] > dist[k])  { index = m; break;     }
 135                                                         }
 136
 137                                                         // if you don't find one than all the mins are less than d
 138                                                         if (index == -1) { index = min1.size();  }
 139
 140                                                         //****** coverage of AB **********//
 141                                                         int index2 = -1;
 142                                                         //find index in min where value is higher than d
 143                                                         for (int m = 0; m < min2.size(); m++) {
 144                                                                         if (min2[m] > dist[k])  { index2 = m; break;    }
 145                                                         }
 146
 147                                                         // if you don't find one than all the mins are less than d
 148                                                         if (index2 == -1) { index2 = min2.size();  }
 149
 150                                                         //coverage of ii
 151                                                         float covII = 1.0 - ((min1.size()-index)/(float)min1.size());
 152
 153                                                         //coverage of ij
 154                                                         float covIJ = 1.0 - ((min2.size()-index2)/(float)min2.size());
 155
 156                                                         //save value in data (Caa - Cab)^2 * distDiff
 157                                                         data[k][i][j] = ((covII-covIJ) * (covII-covIJ)) * distDiff;
 158
 159                                                         //update distDiff
 160                                                         if (k < data.size() - 1) {
 161                                                                 distDiff = dist[k+1] - dist[k];
 162                                                         }
 163                                                 }
 164
 165                                                 //put matrix back to original
 166                                                 matrix->restore();
 167                                         }
 168                                 }
 169                                 count2++;
 170                         }
 171                         count += numGroups+1; //go from AA to BB to CC
 172                 }
 173
 174         }
 175         catch(exception& e) {
 176                 cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 177                 exit(1);
 178         }
 179         catch(...) {
 180                 cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 181                 exit(1);
 182         }
 183
 184 }
 185