//**********************************************************************************************************************
vector<string> IndicatorCommand::setParameters(){
try {
- CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
- CommandParameter pdesign("design", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(pdesign);
- CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(pshared);
- CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none",false,false); parameters.push_back(prelabund);
- CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
- CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
- CommandParameter ptree("tree", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none",false,false); parameters.push_back(ptree);
- CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
- CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
- CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
+ CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
+ CommandParameter pdesign("design", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none","summary",false,false,true); parameters.push_back(pdesign);
+ CommandParameter pshared("shared", "InputTypes", "", "", "SharedRel", "SharedRel", "none","summary",false,false,true); parameters.push_back(pshared);
+ CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRel", "SharedRel", "none","summary",false,false); parameters.push_back(prelabund);
+ CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
+ CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
+ CommandParameter ptree("tree", "InputTypes", "", "", "TreeDesign", "TreeDesign", "none","tree-summary",false,false,true); parameters.push_back(ptree);
+ CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
+ CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
+ CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false); parameters.push_back(pprocessors);
vector<string> myArray;
for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
exit(1);
}
}
-
+//**********************************************************************************************************************
+string IndicatorCommand::getOutputPattern(string type) {
+ try {
+ string pattern = "";
+
+ if (type == "tree") { pattern = "[filename],indicator.tre"; }
+ else if (type == "summary") { pattern = "[filename],indicator.summary"; }
+ else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
+
+ return pattern;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "IndicatorCommand", "getOutputPattern");
+ exit(1);
+ }
+}
//**********************************************************************************************************************
IndicatorCommand::IndicatorCommand(){
try {
string groupfile = "";
m->setTreeFile(treefile);
Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
- treeMap = new TreeMap();
+ ct = new CountTable();
bool mismatch = false;
-
- for (int i = 0; i < m->Treenames.size(); i++) {
- //sanity check - is this a group that is not in the sharedfile?
+
+ set<string> nameMap;
+ map<string, string> groupMap;
+ set<string> gps;
+ for (int i = 0; i < m->Treenames.size(); i++) {
+ nameMap.insert(m->Treenames[i]);
+ //sanity check - is this a group that is not in the sharedfile?
+ if (i == 0) { gps.insert("Group1"); }
if (designfile == "") {
if (!(m->inUsersGroups(m->Treenames[i], m->getAllGroups()))) {
m->mothurOut("[ERROR]: " + m->Treenames[i] + " is not a group in your shared or relabund file."); m->mothurOutEndLine();
mismatch = true;
}
- treeMap->addSeq(m->Treenames[i], "Group1");
+ groupMap[m->Treenames[i]] = "Group1";
}else{
vector<string> myGroups; myGroups.push_back(m->Treenames[i]);
vector<string> myNames = designMap->getNamesSeqs(myGroups);
mismatch = true;
}
}
- treeMap->addSeq(m->Treenames[i], "Group1");
+ groupMap[m->Treenames[i]] = "Group1";
}
- }
+ }
+ ct->createTable(nameMap, groupMap, gps);
if ((designfile != "") && (m->Treenames.size() != Groups.size())) { cout << Groups.size() << '\t' << m->Treenames.size() << endl; m->mothurOut("[ERROR]: You design file does not match your tree, aborting."); m->mothurOutEndLine(); mismatch = true; }
if (designfile != "") { delete designMap; }
if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
- delete treeMap;
+ delete ct;
return 0;
}
read = new ReadNewickTree(treefile);
- int readOk = read->read(treeMap);
+ int readOk = read->read(ct);
- if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete treeMap; delete read; return 0; }
+ if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete ct; delete read; return 0; }
vector<Tree*> T = read->getTrees();
if (designfile != "") { delete designMap; }
if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
- for (int i = 0; i < T.size(); i++) { delete T[i]; } delete treeMap; return 0;
+ for (int i = 0; i < T.size(); i++) { delete T[i]; } delete ct; return 0;
}
- map<string, string> nameMap;
- T[0]->assembleTree(nameMap);
+ T[0]->assembleTree();
/***************************************************/
// create ouptut tree - respecting pickedGroups //
/***************************************************/
- Tree* outputTree = new Tree(m->getNumGroups(), treeMap);
+ Tree* outputTree = new Tree(m->getNumGroups(), ct);
outputTree->getSubTree(T[0], m->getGroups());
- outputTree->assembleTree(nameMap);
+ outputTree->assembleTree();
//no longer need original tree, we have output tree to use and label
for (int i = 0; i < T.size(); i++) { delete T[i]; }
if (designfile != "") { delete designMap; }
if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } }
else { for (int i = 0; i < lookupFloat.size(); i++) { delete lookupFloat[i]; } }
- delete outputTree; delete treeMap; return 0;
+ delete outputTree; delete ct; return 0;
}
/***************************************************/
// get indicator species values //
/***************************************************/
GetIndicatorSpecies(outputTree);
- delete outputTree; delete treeMap;
+ delete outputTree; delete ct;
}else { //run with design file only
//get indicator species
try {
string thisOutputDir = outputDir;
if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
- string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
+ map<string, string> variables;
+ variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName));
+ string outputFileName = getOutputFileName("summary", variables);
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);
- m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
+ m->mothurOutEndLine(); m->mothurOut("Species\tIndicator_Groups\tIndicatorValue\tpValue\n");
int numBins = 0;
if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
vector<float> indicatorValues; //size of numBins
vector<float> pValues;
+ vector<string> indicatorGroups;
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 != "") {
if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
- indicatorValues = getValues(groupings, randomGroupingsMap);
+ indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
}else {
if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
- indicatorValues = getValues(groupings, randomGroupingsMap);
+ indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
}
/******************************************************/
//output indicator values to table form //
/*****************************************************/
- out << "OTU\tIndicator_Value\tpValue" << endl;
+ out << "OTU\tIndicator_Groups\tIndicator_Value\tpValue" << endl;
for (int j = 0; j < indicatorValues.size(); j++) {
if (m->control_pressed) { out.close(); return 0; }
- out << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
+ out << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t';
if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; }
else { out << "<" << (1/(float)iters) << endl; }
if (pValues[j] <= 0.05) {
- cout << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
+ cout << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t';
string pValueString = "<" + toString((1/(float)iters));
if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
else { cout << "<" << (1/(float)iters); }
- m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
+ m->mothurOutJustToLog(m->currentBinLabels[j] + "\t" + indicatorGroups[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
m->mothurOutEndLine();
}
}
string thisOutputDir = outputDir;
if (outputDir == "") { thisOutputDir += m->hasPath(inputFileName); }
- string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName)) + "indicator.summary";
+ map<string, string> variables;
+ variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(inputFileName));
+ string outputFileName = getOutputFileName("summary",variables);
outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
ofstream out;
//print headings
out << "TreeNode\t";
- for (int i = 0; i < numBins; i++) { out << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; }
+ for (int i = 0; i < numBins; i++) { out << m->currentBinLabels[i] << "_IndGroups" << '\t' << m->currentBinLabels[i] << "_IndValue" << '\t' << "pValue" << '\t'; }
out << endl;
- m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicatorValue\tpValue\n");
+ m->mothurOutEndLine(); m->mothurOut("Node\tSpecies\tIndicator_Groups\tIndicatorValue\tpValue\n");
string treeOutputDir = outputDir;
if (outputDir == "") { treeOutputDir += m->hasPath(treefile); }
- string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
+ variables["[filename]"] = treeOutputDir + m->getRootName(m->getSimpleName(treefile));
+ string outputTreeFileName = getOutputFileName("tree", variables);
//create a map from tree node index to names of descendants, save time later to know which sharedRabund you need
vector<float> indicatorValues; //size of numBins
vector<float> pValues;
+ vector<string> indicatorGroups;
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 != "") {
if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
- indicatorValues = getValues(groupings, randomGroupingsMap);
+ indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
pValues = getPValues(groupings, randomGroupingsMap, lookup.size(), indicatorValues);
}else {
if (groupsAlreadyAdded.size() != lookupFloat.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
- indicatorValues = getValues(groupings, randomGroupingsMap);
+ indicatorValues = getValues(groupings, indicatorGroups, randomGroupingsMap);
pValues = getPValues(groupings, randomGroupingsMap, lookupFloat.size(), indicatorValues);
}
if (m->control_pressed) { out.close(); return 0; }
if (pValues[j] < (1/(float)iters)) {
- out << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
+ out << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t' << '<' << (1/(float)iters) << '\t';
}else {
- out << indicatorValues[j] << '\t' << pValues[j] << '\t';
+ out << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t' << pValues[j] << '\t';
}
if (pValues[j] <= 0.05) {
- cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorValues[j] << '\t';
+ cout << i+1 << '\t' << m->currentBinLabels[j] << '\t' << indicatorGroups[j] << '\t' << indicatorValues[j] << '\t';
string pValueString = "<" + toString((1/(float)iters));
if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];}
else { cout << "<" << (1/(float)iters); }
- m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
+ m->mothurOutJustToLog(toString(i) + "\t" + m->currentBinLabels[j] + "\t" + indicatorGroups[j] + "\t" + toString(indicatorValues[j]) + "\t" + pValueString);
m->mothurOutEndLine();
}
}
}
}
//**********************************************************************************************************************
-vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
+vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundFloatVector*> >& groupings, vector<string>& indicatorGroupings, map< vector<int>, vector<int> > groupingsMap){
try {
vector<float> values;
map< vector<int>, vector<int> >::iterator it;
-
+
+ indicatorGroupings.clear();
+
+ //create grouping strings
+ vector<string> groupingsGroups;
+ for (int j = 0; j < groupings.size(); j++) {
+ string tempGrouping = "";
+ for (int k = 0; k < groupings[j].size()-1; k++) {
+ tempGrouping += groupings[j][k]->getGroup() + "-";
+ }
+ tempGrouping += groupings[j][groupings[j].size()-1]->getGroup();
+ groupingsGroups.push_back(tempGrouping);
+ }
+
+
//for each otu
for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
}
float maxIndVal = 0.0;
+ string maxGrouping = "";
for (int j = 0; j < terms.size(); j++) {
float thisAij = (terms[j] / AijDenominator); //relative abundance
float thisValue = thisAij * Bij[j] * 100.0;
//save largest
- if (thisValue > maxIndVal) { maxIndVal = thisValue; }
+ if (thisValue > maxIndVal) { maxIndVal = thisValue; maxGrouping = groupingsGroups[j]; }
}
values.push_back(maxIndVal);
+ indicatorGroupings.push_back(maxGrouping);
}
return values;
}
//**********************************************************************************************************************
//same as above, just data type difference
-vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, map< vector<int>, vector<int> > groupingsMap){
+vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >& groupings, vector<string>& indicatorGroupings, map< vector<int>, vector<int> > groupingsMap){
try {
vector<float> values;
-
- /*for (int j = 0; j < groupings.size(); j++) {
- cout << "grouping " << j << endl;
- for (int k = 0; k < groupings[j].size(); k++) {
- cout << groupings[j][k]->getGroup() << endl;
- }
- }*/
map< vector<int>, vector<int> >::iterator it;
+
+ indicatorGroupings.clear();
+
+ //create grouping strings
+ vector<string> groupingsGroups;
+ for (int j = 0; j < groupings.size(); j++) {
+ string tempGrouping = "";
+ for (int k = 0; k < groupings[j].size()-1; k++) {
+ tempGrouping += groupings[j][k]->getGroup() + "-";
+ }
+ tempGrouping += groupings[j][groupings[j].size()-1]->getGroup();
+ groupingsGroups.push_back(tempGrouping);
+ }
+
//for each otu
for (int i = 0; i < groupings[0][0]->getNumBins(); i++) {
}
float maxIndVal = 0.0;
+ string maxGrouping = "";
for (int j = 0; j < terms.size(); j++) {
float thisAij = (terms[j] / AijDenominator); //relative abundance
float thisValue = thisAij * Bij[j] * 100.0;
//save largest
- if (thisValue > maxIndVal) { maxIndVal = thisValue; }
+ if (thisValue > maxIndVal) { maxIndVal = thisValue; maxGrouping = groupingsGroups[j]; }
}
values.push_back(maxIndVal);
+ indicatorGroupings.push_back(maxGrouping);
}
return values;
try {
vector<float> pvalues;
pvalues.resize(indicatorValues.size(), 0);
+ vector<string> notUsedGroupings; //we dont care about the grouping for the pvalues since they are randomized, but we need to pass the function something to make it work.
for(int i=0;i<numIters;i++){
if (m->control_pressed) { break; }
groupingsMap = randomizeGroupings(groupings, num);
- vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
+ vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap);
for (int j = 0; j < indicatorValues.size(); j++) {
if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
try {
vector<float> pvalues;
pvalues.resize(indicatorValues.size(), 0);
+ vector<string> notUsedGroupings; //we dont care about the grouping for the pvalues since they are randomized, but we need to pass the function something to make it work.
for(int i=0;i<numIters;i++){
if (m->control_pressed) { break; }
groupingsMap = randomizeGroupings(groupings, num);
- vector<float> randomIndicatorValues = getValues(groupings, groupingsMap);
+ vector<float> randomIndicatorValues = getValues(groupings, notUsedGroupings, groupingsMap);
for (int j = 0; j < indicatorValues.size(); j++) {
if (randomIndicatorValues[j] >= indicatorValues[j]) { pvalues[j]++; }
projects / mothur.git / blobdiff