indicatorcommand.h

   1 #ifndef INDICATORCOMMAND_H
   2 #define INDICATORCOMMAND_H
   3
   4 /*
   5  *  indicatorcommand.h
   6  *  Mothur
   7  *
   8  *  Created by westcott on 11/12/10.
   9  *  Copyright 2010 Schloss Lab. All rights reserved.
  10  *
  11  */
  12
  13 #include "command.hpp"
  14 #include "readtree.h"
  15 #include "counttable.h"
  16 #include "sharedrabundvector.h"
  17 #include "sharedrabundfloatvector.h"
  18 #include "inputdata.h"
  19
  20 class IndicatorCommand : public Command {
  21 public:
  22         IndicatorCommand(string);
  23         IndicatorCommand();
  24         ~IndicatorCommand(){}
  25
  26         vector<string> setParameters();
  27         string getCommandName()                 { return "indicator";                           }
  28         string getCommandCategory()             { return "Hypothesis Testing";          }
  29
  30         string getHelpString();
  31     string getOutputPattern(string);
  32         string getCitation() { return "Dufrene M, Legendre P (1997). Species assemblages and indicator species: The need for a flexible asymmetrical approach. Ecol Monogr 67: 345-66.\n McCune B, Grace JB, Urban DL (2002). Analysis of ecological communities. MjM Software Design: Gleneden Beach, OR. \nLegendre P, Legendre L (1998). Numerical Ecology. Elsevier: New York. \nhttp://www.mothur.org/wiki/Indicator"; }
  33         string getDescription()         { return "calculate the indicator value for each OTU"; }
  34
  35         int execute();
  36         void help() { m->mothurOut(getHelpString()); }
  37
  38 private:
  39         ReadTree* read;
  40         CountTable* ct;
  41         GroupMap* designMap;
  42         string treefile, sharedfile, relabundfile, groups, label, inputFileName, outputDir, designfile;
  43         bool abort;
  44         int iters, processors;
  45         vector<string> outputNames, Groups;
  46         vector<SharedRAbundVector*> lookup;
  47         vector<SharedRAbundFloatVector*> lookupFloat;
  48
  49         int getShared();
  50         int getSharedFloat();
  51         int GetIndicatorSpecies(Tree*&);
  52         int GetIndicatorSpecies();
  53         set<string> getDescendantList(Tree*&, int, map<int, set<string> >, map<int, set<int> >&);
  54         vector<float> getValues(vector< vector<SharedRAbundVector*> >&, vector<string>&, map< vector<int>, vector<int> >);
  55         vector<float> getValues(vector< vector<SharedRAbundFloatVector*> >&, vector<string>&, map< vector<int>, vector<int> >);
  56
  57         map<int, float> getDistToRoot(Tree*&);
  58         map< vector<int>, vector<int> > randomizeGroupings(vector< vector<SharedRAbundVector*> >&, int);
  59         map< vector<int>, vector<int> > randomizeGroupings(vector< vector<SharedRAbundFloatVector*> >&, int);
  60
  61         vector<float> driver(vector< vector<SharedRAbundFloatVector*> >&, int, vector<float>, int);
  62         vector<float> driver(vector< vector<SharedRAbundVector*> >&, int, vector<float>, int);
  63
  64         vector<float> getPValues(vector< vector<SharedRAbundFloatVector*> >&, int, vector<float>);
  65         vector<float> getPValues(vector< vector<SharedRAbundVector*> >&, int, vector<float>);
  66
  67
  68 };
  69
  70 /**************************************************************************************************/
  71
  72 struct indicatorData {
  73     vector< vector<SharedRAbundFloatVector*> > groupings;
  74         MothurOut* m;
  75     int iters, num;
  76     vector<float> indicatorValues;
  77     vector<float> pvalues;
  78
  79         indicatorData(){}
  80         indicatorData(MothurOut* mout, int it, vector< vector<SharedRAbundFloatVector*> > ng, int n, vector<float> iv) {
  81                 m = mout;
  82         iters = it;
  83         groupings = ng;
  84         indicatorValues = iv;
  85         num = n;
  86     }
  87 };
  88 /**************************************************************************************************/
  89 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
  90 #else
  91 static DWORD WINAPI MyIndicatorThreadFunction(LPVOID lpParam){
  92         indicatorData* pDataArray;
  93         pDataArray = (indicatorData*)lpParam;
  94
  95         try {
  96
  97                 pDataArray->pvalues.resize(pDataArray->indicatorValues.size(), 0);
  98
  99                 for(int i=0;i<pDataArray->iters;i++){
 100                         if (pDataArray->m->control_pressed) { break; }
 101
 102                         //groupingsMap = randomizeGroupings(groupings, num);
 103             ///////////////////////////////////////////////////////////////////////
 104             map< vector<int>, vector<int> > randomGroupings;
 105
 106             for (int j = 0; j < pDataArray->num; j++) {
 107
 108                 //get random groups to swap to switch with
 109                 //generate random int between 0 and groupings.size()-1
 110                 int z = pDataArray->m->getRandomIndex(pDataArray->groupings.size()-1);
 111                 int x = pDataArray->m->getRandomIndex(pDataArray->groupings.size()-1);
 112                 int a = pDataArray->m->getRandomIndex(pDataArray->groupings[z].size()-1);
 113                 int b = pDataArray->m->getRandomIndex(pDataArray->groupings[x].size()-1);
 114                 //cout << i << '\t' << z << '\t' << x << '\t' << a << '\t' << b << endl;
 115
 116                 vector<int> from;
 117                 vector<int> to;
 118
 119                 from.push_back(z); from.push_back(a);
 120                 to.push_back(x); to.push_back(b);
 121
 122                 randomGroupings[from] = to;
 123             }
 124             ///////////////////////////////////////////////////////////////////////
 125
 126                         //vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, randomGroupings);
 127             ///////////////////////////////////////////////////////////////////////
 128             vector<float> randomIndicatorValues;
 129             map< vector<int>, vector<int> >::iterator it;
 130
 131             //for each otu
 132             for (int i = 0; i < pDataArray->groupings[0][0]->getNumBins(); i++) {
 133
 134                 if (pDataArray->m->control_pressed) { return 0; }
 135
 136                 vector<float> terms;
 137                 float AijDenominator = 0.0;
 138                 vector<float> Bij;
 139
 140                 //get overall abundance of each grouping
 141                 for (int j = 0; j < pDataArray->groupings.size(); j++) {
 142
 143                     float totalAbund = 0;
 144                     int numNotZero = 0;
 145                     for (int k = 0; k < pDataArray->groupings[j].size(); k++) {
 146                         vector<int> temp; temp.push_back(j); temp.push_back(k);
 147                         it = randomGroupings.find(temp);
 148
 149                         if (it == randomGroupings.end()) { //this one didnt get moved
 150                             totalAbund += pDataArray->groupings[j][k]->getAbundance(i);
 151                             if (pDataArray->groupings[j][k]->getAbundance(i) != 0.0) { numNotZero++; }
 152                         }else {
 153                             totalAbund += pDataArray->groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i);
 154                             if (pDataArray->groupings[(it->second)[0]][(it->second)[1]]->getAbundance(i) != 0.0) { numNotZero++; }
 155                         }
 156
 157                     }
 158
 159                     //mean abundance
 160                     float Aij = (totalAbund / (float) pDataArray->groupings[j].size());
 161                     terms.push_back(Aij);
 162
 163                     //percentage of sites represented
 164                     Bij.push_back(numNotZero / (float) pDataArray->groupings[j].size());
 165
 166                     AijDenominator += Aij;
 167                 }
 168
 169                 float maxIndVal = 0.0;
 170                 for (int j = 0; j < terms.size(); j++) {
 171                     float thisAij = (terms[j] / AijDenominator); //relative abundance
 172                     float thisValue = thisAij * Bij[j] * 100.0;
 173
 174                     //save largest
 175                     if (thisValue > maxIndVal) { maxIndVal = thisValue; }
 176                 }
 177
 178                 randomIndicatorValues.push_back(maxIndVal);
 179             }
 180
 181             ///////////////////////////////////////////////////////////////////////
 182
 183                         for (int j = 0; j < pDataArray->indicatorValues.size(); j++) {
 184                                 if (randomIndicatorValues[j] >= pDataArray->indicatorValues[j]) { pDataArray->pvalues[j]++; }
 185                         }
 186                 }
 187
 188         return 0;
 189
 190         }
 191         catch(exception& e) {
 192                 pDataArray->m->errorOut(e, "IndicatorCommand", "MyIndicatorThreadFunction");
 193                 exit(1);
 194         }
 195 }
 196 #endif
 197
 198
 199
 200 #endif
 201