1 #ifndef INDICATORCOMMAND_H
2 #define INDICATORCOMMAND_H
8 * Created by westcott on 11/12/10.
9 * Copyright 2010 Schloss Lab. All rights reserved.
13 #include "command.hpp"
15 #include "counttable.h"
16 #include "sharedrabundvector.h"
17 #include "sharedrabundfloatvector.h"
18 #include "inputdata.h"
20 class IndicatorCommand : public Command {
22 IndicatorCommand(string);
26 vector<string> setParameters();
27 string getCommandName() { return "indicator"; }
28 string getCommandCategory() { return "Hypothesis Testing"; }
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"; }
36 void help() { m->mothurOut(getHelpString()); }
42 string treefile, sharedfile, relabundfile, groups, label, inputFileName, outputDir, designfile;
44 int iters, processors;
45 vector<string> outputNames, Groups;
46 vector<SharedRAbundVector*> lookup;
47 vector<SharedRAbundFloatVector*> lookupFloat;
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> >);
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);
61 vector<float> driver(vector< vector<SharedRAbundFloatVector*> >&, int, vector<float>, int);
62 vector<float> driver(vector< vector<SharedRAbundVector*> >&, int, vector<float>, int);
64 vector<float> getPValues(vector< vector<SharedRAbundFloatVector*> >&, int, vector<float>);
65 vector<float> getPValues(vector< vector<SharedRAbundVector*> >&, int, vector<float>);
70 /**************************************************************************************************/
72 struct indicatorData {
73 vector< vector<SharedRAbundFloatVector*> > groupings;
76 vector<float> indicatorValues;
77 vector<float> pvalues;
80 indicatorData(MothurOut* mout, int it, vector< vector<SharedRAbundFloatVector*> > ng, int n, vector<float> iv) {
88 /**************************************************************************************************/
89 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
91 static DWORD WINAPI MyIndicatorThreadFunction(LPVOID lpParam){
92 indicatorData* pDataArray;
93 pDataArray = (indicatorData*)lpParam;
97 pDataArray->pvalues.resize(pDataArray->indicatorValues.size(), 0);
99 for(int i=0;i<pDataArray->iters;i++){
100 if (pDataArray->m->control_pressed) { break; }
102 //groupingsMap = randomizeGroupings(groupings, num);
103 ///////////////////////////////////////////////////////////////////////
104 map< vector<int>, vector<int> > randomGroupings;
106 for (int j = 0; j < pDataArray->num; j++) {
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;
119 from.push_back(z); from.push_back(a);
120 to.push_back(x); to.push_back(b);
122 randomGroupings[from] = to;
124 ///////////////////////////////////////////////////////////////////////
126 //vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, randomGroupings);
127 ///////////////////////////////////////////////////////////////////////
128 vector<float> randomIndicatorValues;
129 map< vector<int>, vector<int> >::iterator it;
132 for (int i = 0; i < pDataArray->groupings[0][0]->getNumBins(); i++) {
134 if (pDataArray->m->control_pressed) { return 0; }
137 float AijDenominator = 0.0;
140 //get overall abundance of each grouping
141 for (int j = 0; j < pDataArray->groupings.size(); j++) {
143 float totalAbund = 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);
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++; }
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++; }
160 float Aij = (totalAbund / (float) pDataArray->groupings[j].size());
161 terms.push_back(Aij);
163 //percentage of sites represented
164 Bij.push_back(numNotZero / (float) pDataArray->groupings[j].size());
166 AijDenominator += Aij;
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;
175 if (thisValue > maxIndVal) { maxIndVal = thisValue; }
178 randomIndicatorValues.push_back(maxIndVal);
181 ///////////////////////////////////////////////////////////////////////
183 for (int j = 0; j < pDataArray->indicatorValues.size(); j++) {
184 if (randomIndicatorValues[j] >= pDataArray->indicatorValues[j]) { pDataArray->pvalues[j]++; }
191 catch(exception& e) {
192 pDataArray->m->errorOut(e, "IndicatorCommand", "MyIndicatorThreadFunction");