+//divide shared or relabund file by groupings in the design file
+//report all otu values to file
+int IndicatorCommand::GetIndicatorSpecies(){
+ try {
+ string thisOutputDir = outputDir;
+ if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
+ string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
+ outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
+
+ ofstream out;
+ m->openOutputFile(outputFileName, out);
+ out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
+
+ int numBins = 0;
+ if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
+ else { numBins = lookupFloat[0]->getNumBins(); }
+
+ if (m->control_pressed) { out.close(); return 0; }
+
+ /*****************************************************/
+ //create vectors containing rabund info //
+ /*****************************************************/
+
+ vector<float> indicatorValues; //size of numBins
+ vector<float> pValues;
+ map< vector<int>, vector<int> > randomGroupingsMap; //maps location in groupings to location in groupings, ie, [0][0] -> [1][2]. This is so we don't have to actually move the sharedRabundVectors.
+
+ if (sharedfile != "") {
+ vector< vector<SharedRAbundVector*> > groupings;
+ set<string> groupsAlreadyAdded;
+ vector<SharedRAbundVector*> subset;
+
+ //for each grouping
+ for (int i = 0; i < designMap->namesOfGroups.size(); i++) {
+
+ for (int k = 0; k < lookup.size(); k++) {
+ //are you from this grouping?
+ if (designMap->getGroup(lookup[k]->getGroup()) == designMap->namesOfGroups[i]) {
+ subset.push_back(lookup[k]);
+ groupsAlreadyAdded.insert(lookup[k]->getGroup());
+ }
+ }
+ if (subset.size() != 0) { groupings.push_back(subset); }
+ subset.clear();
+ }
+
+ if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
+
+ indicatorValues = getValues(groupings, randomGroupingsMap);
+
+ pValues.resize(indicatorValues.size(), 0);
+ for(int i=0;i<iters;i++){
+ if (m->control_pressed) { break; }
+ randomGroupingsMap = randomizeGroupings(groupings, lookup.size());
+ vector<float> randomIndicatorValues = getValues(groupings, randomGroupingsMap);
+
+ for (int j = 0; j < indicatorValues.size(); j++) {
+ if (randomIndicatorValues[j] >= indicatorValues[j]) { pValues[j]++; }
+ }
+ }
+
+ for (int i = 0; i < pValues.size(); i++) { pValues[i] /= (double)iters; }
+
+ }else {
+ vector< vector<SharedRAbundFloatVector*> > groupings;
+ set<string> groupsAlreadyAdded;
+ vector<SharedRAbundFloatVector*> subset;
+
+ //for each grouping
+ for (int i = 0; i < designMap->namesOfGroups.size(); i++) {
+ for (int k = 0; k < lookupFloat.size(); k++) {
+ //are you from this grouping?
+ if (designMap->getGroup(lookupFloat[k]->getGroup()) == designMap->namesOfGroups[i]) {
+ subset.push_back(lookupFloat[k]);
+ groupsAlreadyAdded.insert(lookupFloat[k]->getGroup());
+ }
+ }
+ if (subset.size() != 0) { groupings.push_back(subset); }
+ subset.clear();
+ }
+
+ if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
+
+ indicatorValues = getValues(groupings, randomGroupingsMap);
+
+ pValues.resize(indicatorValues.size(), 0);
+ for(int i=0;i<iters;i++){
+ if (m->control_pressed) { break; }
+ randomGroupingsMap = randomizeGroupings(groupings, lookupFloat.size());
+ vector<float> randomIndicatorValues = getValues(groupings, randomGroupingsMap);
+
+ for (int j = 0; j < indicatorValues.size(); j++) {
+ if (randomIndicatorValues[j] >= indicatorValues[j]) { pValues[j]++; }
+ }
+ }
+
+ for (int i = 0; i < pValues.size(); i++) { pValues[i] /= (double)iters; }
+ }
+
+ if (m->control_pressed) { out.close(); return 0; }
+
+
+ /******************************************************/
+ //output indicator values to table form //
+ /*****************************************************/
+ int indicatorOTU;
+ double maxValue, pvalue; maxValue=0.0; pvalue=0.0;
+ out << "OTU\tIndicator_Value\tpValue" << endl;
+ for (int j = 0; j < indicatorValues.size(); j++) {
+
+ if (m->control_pressed) { out.close(); return 0; }
+
+ out << (j+1) << '\t' << indicatorValues[j] << '\t';
+
+ if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
+ else { out << "<" << (1/(float)iters) << endl; }
+
+ if (maxValue < indicatorValues[j]) {
+ maxValue = indicatorValues[j];
+ indicatorOTU = j;
+ }
+ }
+
+ m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
+ cout << "OTU" << indicatorOTU+1 << '\t' << maxValue << '\t';
+ string pValueString = "<" + toString((1/(float)iters));
+ if (pValues[indicatorOTU] > (1/(float)iters)) { pValueString = toString(pValues[indicatorOTU]); cout << pValues[indicatorOTU];}
+ else { cout << "<" << (1/(float)iters); }
+ m->mothurOutJustToLog("OTU" + toString(indicatorOTU+1) + "\t" + toString(maxValue) + "\t" + pValueString);
+ m->mothurOutEndLine();
+
+ out.close();
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "IndicatorCommand", "GetIndicatorSpecies");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************