*/
#include "metastatscommand.h"
-#include "metastats.h"
#include "sharedutilities.h"
+
//**********************************************************************************************************************
-vector<string> MetaStatsCommand::getValidParameters(){
+vector<string> MetaStatsCommand::setParameters(){
try {
- string Array[] = {"groups","label","outputdir","iters","threshold","g","design","sets","processors","inputdir"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","metastats",false,true,true); parameters.push_back(pshared);
+ CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(pdesign);
+ CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
+ CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
+ CommandParameter pthreshold("threshold", "Number", "", "0.05", "", "", "","",false,false); parameters.push_back(pthreshold);
+ CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
+ CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
+ CommandParameter psets("sets", "String", "", "", "", "", "","",false,false); parameters.push_back(psets);
+ CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
+ CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
+
+ vector<string> myArray;
+ for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
return myArray;
}
catch(exception& e) {
- m->errorOut(e, "MetaStatsCommand", "getValidParameters");
+ m->errorOut(e, "MetaStatsCommand", "setParameters");
exit(1);
}
}
//**********************************************************************************************************************
-MetaStatsCommand::MetaStatsCommand(){
+string MetaStatsCommand::getHelpString(){
try {
- abort = true;
- //initialize outputTypes
- vector<string> tempOutNames;
- outputTypes["metastats"] = tempOutNames;
+ string helpString = "";
+ helpString += "This command is based on the Metastats program, White, J.R., Nagarajan, N. & Pop, M. Statistical methods for detecting differentially abundant features in clinical metagenomic samples. PLoS Comput Biol 5, e1000352 (2009).\n";
+ helpString += "The metastats command outputs a .metastats file. \n";
+ helpString += "The metastats command parameters are shared, iters, threshold, groups, label, design, sets and processors. The shared and design parameters are required, unless you have valid current files.\n";
+ helpString += "The design parameter allows you to assign your groups to sets when you are running metastat. mothur will run all pairwise comparisons of the sets. It is required. \n";
+ helpString += "The design file looks like the group file. It is a 2 column tab delimited file, where the first column is the group name and the second column is the set the group belongs to.\n";
+ helpString += "The sets parameter allows you to specify which of the sets in your designfile you would like to analyze. The set names are separated by dashes. THe default is all sets in the designfile.\n";
+ helpString += "The iters parameter allows you to set number of bootstrap permutations for estimating null distribution of t statistic. The default is 1000. \n";
+ helpString += "The threshold parameter allows you to set the significance level to reject null hypotheses (default 0.05).\n";
+ helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like included. The group names are separated by dashes.\n";
+ helpString += "The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n";
+ helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n";
+ helpString += "The metastats command should be in the following format: metastats(design=yourDesignFile).\n";
+ helpString += "Example metastats(design=temp.design, groups=A-B-C).\n";
+ helpString += "The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n";
+ helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
+ return helpString;
}
catch(exception& e) {
- m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
+ m->errorOut(e, "MetaStatsCommand", "getHelpString");
exit(1);
}
}
//**********************************************************************************************************************
-vector<string> MetaStatsCommand::getRequiredParameters(){
- try {
- string Array[] = {"design"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
- return myArray;
- }
- catch(exception& e) {
- m->errorOut(e, "MetaStatsCommand", "getRequiredParameters");
- exit(1);
- }
+string MetaStatsCommand::getOutputPattern(string type) {
+ try {
+ string pattern = "";
+
+ if (type == "metastats") { pattern = "[filename],[distance],[group],metastats"; }
+ else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
+
+ return pattern;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "MetaStatsCommand", "getOutputPattern");
+ exit(1);
+ }
}
//**********************************************************************************************************************
-vector<string> MetaStatsCommand::getRequiredFiles(){
+MetaStatsCommand::MetaStatsCommand(){
try {
- string Array[] = {"shared"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
- return myArray;
+ abort = true; calledHelp = true;
+ setParameters();
+ vector<string> tempOutNames;
+ outputTypes["metastats"] = tempOutNames;
}
catch(exception& e) {
- m->errorOut(e, "MetaStatsCommand", "getRequiredFiles");
+ m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
exit(1);
}
}
MetaStatsCommand::MetaStatsCommand(string option) {
try {
- globaldata = GlobalData::getInstance();
- abort = false;
+ abort = false; calledHelp = false;
allLines = 1;
- labels.clear();
+
//allow user to run help
- if(option == "help") { help(); abort = true; }
+ if(option == "help") { help(); abort = true; calledHelp = true; }
+ else if(option == "citation") { citation(); abort = true; calledHelp = true;}
else {
- //valid paramters for this command
- string AlignArray[] = {"groups","label","outputdir","iters","threshold","g","design","sets","processors","inputdir"};
- vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
+ vector<string> myArray = setParameters();
OptionParser parser(option);
map<string,string> parameters = parser.getParameters();
vector<string> tempOutNames;
outputTypes["metastats"] = tempOutNames;
- //if the user changes the output directory command factory will send this info to us in the output parameter
- outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
- outputDir = "";
- outputDir += m->hasPath(globaldata->inputFileName); //if user entered a file with a path then preserve it
- }
//if the user changes the input directory command factory will send this info to us in the output parameter
string inputDir = validParameter.validFile(parameters, "inputdir", false);
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["design"] = inputDir + it->second; }
}
+
+ it = parameters.find("shared");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["shared"] = inputDir + it->second; }
+ }
+
}
+ //check for required parameters
+ sharedfile = validParameter.validFile(parameters, "shared", true);
+ if (sharedfile == "not open") { abort = true; }
+ else if (sharedfile == "not found") { //if there is a current shared file, use it
+ sharedfile = m->getSharedFile();
+ if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
+ else { m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
+ }else { m->setSharedFile(sharedfile); }
+
//check for required parameters
designfile = validParameter.validFile(parameters, "design", true);
if (designfile == "not open") { abort = true; }
- else if (designfile == "not found") { designfile = ""; m->mothurOut("You must provide an design file."); m->mothurOutEndLine(); abort = true; }
+ else if (designfile == "not found") {
+ //if there is a current design file, use it
+ designfile = m->getDesignFile();
+ if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
+ else { m->mothurOut("You have no current designfile and the design parameter is required."); m->mothurOutEndLine(); abort = true; }
+ }else { m->setDesignFile(designfile); }
- //make sure the user has already run the read.otu command
- if ((globaldata->getSharedFile() == "")) {
- m->mothurOut("You must read a list and a group, or a shared file before you can use the metastats command."); m->mothurOutEndLine(); abort = true;
+ //if the user changes the output directory command factory will send this info to us in the output parameter
+ outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
+ outputDir = "";
+ outputDir += m->hasPath(sharedfile); //if user entered a file with a path then preserve it
}
//check for optional parameter and set defaults
else { allLines = 1; }
}
- //if the user has not specified any labels use the ones from read.otu
- if (label == "") {
- allLines = globaldata->allLines;
- labels = globaldata->labels;
- }
-
groups = validParameter.validFile(parameters, "groups", false);
if (groups == "not found") { groups = ""; pickedGroups = false; }
else {
pickedGroups = true;
m->splitAtDash(groups, Groups);
- globaldata->Groups = Groups;
+ m->setGroups(Groups);
}
sets = validParameter.validFile(parameters, "sets", false);
string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
- convert(temp, iters);
+ m->mothurConvert(temp, iters);
temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "0.05"; }
- convert(temp, threshold);
+ m->mothurConvert(temp, threshold);
- temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
- convert(temp, processors);
+ temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
+ m->setProcessors(temp);
+ m->mothurConvert(temp, processors);
}
}
exit(1);
}
}
-
-//**********************************************************************************************************************
-
-void MetaStatsCommand::help(){
- try {
- m->mothurOut("This command is based on the Metastats program, White, J.R., Nagarajan, N. & Pop, M. Statistical methods for detecting differentially abundant features in clinical metagenomic samples. PLoS Comput Biol 5, e1000352 (2009).\n");
- m->mothurOut("The metastats command can only be executed after a successful read.otu command of a list and group or shared file.\n");
- m->mothurOut("The metastats command outputs a .metastats file. \n");
- m->mothurOut("The metastats command parameters are iters, threshold, groups, label, design, sets and processors. The design parameter is required.\n");
- m->mothurOut("The design parameter allows you to assign your groups to sets when you are running metastat. mothur will run all pairwise comparisons of the sets. It is required. \n");
- m->mothurOut("The design file looks like the group file. It is a 2 column tab delimited file, where the first column is the group name and the second column is the set the group belongs to.\n");
- m->mothurOut("The sets parameter allows you to specify which of the sets in your designfile you would like to analyze. The set names are separated by dashes. THe default is all sets in the designfile.\n");
- m->mothurOut("The iters parameter allows you to set number of bootstrap permutations for estimating null distribution of t statistic. The default is 1000. \n");
- m->mothurOut("The threshold parameter allows you to set the significance level to reject null hypotheses (default 0.05).\n");
- m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included. The group names are separated by dashes.\n");
- m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n");
- m->mothurOut("The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n");
- m->mothurOut("The metastats command should be in the following format: metastats(design=yourDesignFile).\n");
- m->mothurOut("Example metastats(design=temp.design, groups=A-B-C).\n");
- m->mothurOut("The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n");
- m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
-
- }
- catch(exception& e) {
- m->errorOut(e, "MetaStatsCommand", "help");
- exit(1);
- }
-}
-
-//**********************************************************************************************************************
-
-MetaStatsCommand::~MetaStatsCommand(){}
-
//**********************************************************************************************************************
int MetaStatsCommand::execute(){
try {
- if (abort == true) { return 0; }
+ if (abort == true) { if (calledHelp) { return 0; } return 2; }
+
+ //just used to convert files to test metastats online
+ /****************************************************/
+ bool convertInputToShared = false;
+ convertSharedToInput = false;
+ if (convertInputToShared) { convertToShared(sharedfile); return 0; }
+ /****************************************************/
+
designMap = new GroupMap(designfile);
designMap->readDesignMap();
-
- read = new ReadOTUFile(globaldata->inputFileName);
- read->read(&*globaldata);
- input = globaldata->ginput;
+
+ input = new InputData(sharedfile, "sharedfile");
lookup = input->getSharedRAbundVectors();
string lastLabel = lookup[0]->getLabel();
//setup the pairwise comparions of sets for metastats
//calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
//make sure sets are all in designMap
- SharedUtil* util = new SharedUtil();
- util->setGroups(Sets, designMap->namesOfGroups);
+ SharedUtil* util = new SharedUtil();
+ vector<string> dGroups = designMap->getNamesOfGroups();
+ util->setGroups(Sets, dGroups);
delete util;
int numGroups = Sets.size();
//only 1 combo
if (numGroups == 2) { processors = 1; }
else if (numGroups < 2) { m->mothurOut("Not enough sets, I need at least 2 valid sets. Unable to complete command."); m->mothurOutEndLine(); m->control_pressed = true; }
-
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- if(processors != 1){
- int numPairs = namesOfGroupCombos.size();
- int numPairsPerProcessor = numPairs / processors;
+
+ if(processors != 1){
+ int numPairs = namesOfGroupCombos.size();
+ int numPairsPerProcessor = numPairs / processors;
- for (int i = 0; i < processors; i++) {
- int startPos = i * numPairsPerProcessor;
- if(i == processors - 1){
- numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
- }
- lines.push_back(linePair(startPos, numPairsPerProcessor));
- }
- }
- #endif
+ for (int i = 0; i < processors; i++) {
+ int startPos = i * numPairsPerProcessor;
+ if(i == processors - 1){
+ numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
+ }
+ lines.push_back(linePair(startPos, numPairsPerProcessor));
+ }
+ }
//as long as you are not at the end of the file or done wih the lines you want
while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
- if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } globaldata->Groups.clear(); delete read; delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } m->clearGroups(); delete input; delete designMap; for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
//prevent memory leak
for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; }
- if (m->control_pressed) { outputTypes.clear(); globaldata->Groups.clear(); delete read; delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { outputTypes.clear(); m->clearGroups(); delete input; delete designMap; for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
//get next line to process
lookup = input->getSharedRAbundVectors();
}
- if (m->control_pressed) { outputTypes.clear(); globaldata->Groups.clear(); delete read; delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { outputTypes.clear(); m->clearGroups(); delete input; delete designMap; for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
//output error messages about any remaining user labels
set<string>::iterator it;
}
//reset groups parameter
- globaldata->Groups.clear();
- delete input; globaldata->ginput = NULL;
- delete read;
+ m->clearGroups();
+ delete input;
delete designMap;
- if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0;}
+ if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0;}
m->mothurOutEndLine();
m->mothurOut("Output File Names: "); m->mothurOutEndLine();
int MetaStatsCommand::process(vector<SharedRAbundVector*>& thisLookUp){
try {
- if (pickedGroups) { eliminateZeroOTUS(thisLookUp); }
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+
if(processors == 1){
driver(0, namesOfGroupCombos.size(), thisLookUp);
}else{
int process = 1;
vector<int> processIDS;
-
+ #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
int temp = processIDS[i];
wait(&temp);
}
- }
- #else
- driver(0, namesOfGroupCombos.size(), thisLookUp);
- #endif
+ #else
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+ //Windows version shared memory, so be careful when passing variables through the summarySharedData struct.
+ //Above fork() will clone, so memory is separate, but that's not the case with windows,
+ //Taking advantage of shared memory to pass results vectors.
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ vector<metastatsData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+
+ //make copy of lookup so we don't get access violations
+ vector<SharedRAbundVector*> newLookup;
+ vector<string> designMapGroups;
+ for (int k = 0; k < thisLookUp.size(); k++) {
+ SharedRAbundVector* temp = new SharedRAbundVector();
+ temp->setLabel(thisLookUp[k]->getLabel());
+ temp->setGroup(thisLookUp[k]->getGroup());
+ newLookup.push_back(temp);
+ designMapGroups.push_back(designMap->getGroup(thisLookUp[k]->getGroup()));
+ }
+
+ //for each bin
+ for (int k = 0; k < thisLookUp[0]->getNumBins(); k++) {
+ if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
+ for (int j = 0; j < thisLookUp.size(); j++) { newLookup[j]->push_back(thisLookUp[j]->getAbundance(k), thisLookUp[j]->getGroup()); }
+ }
+
+ // Allocate memory for thread data.
+ metastatsData* tempSum = new metastatsData(sharedfile, outputDir, m, lines[i].start, lines[i].num, namesOfGroupCombos, newLookup, designMapGroups, iters, threshold);
+ pDataArray.push_back(tempSum);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyMetastatsThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ //do my part
+ driver(lines[0].start, lines[0].num, thisLookUp);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ if (pDataArray[i]->count != (pDataArray[i]->num)) {
+ m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->num) + " groups assigned to it, quitting. \n"); m->control_pressed = true;
+ }
+ for (int j = 0; j < pDataArray[i]->thisLookUp.size(); j++) { delete pDataArray[i]->thisLookUp[j]; }
+ for (int j = 0; j < pDataArray[i]->outputNames.size(); j++) {
+ outputNames.push_back(pDataArray[i]->outputNames[j]);
+ outputTypes["metastats"].push_back(pDataArray[i]->outputNames[j]);
+ }
+
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+ #endif
+ }
+
return 0;
}
for (int c = start; c < (start+num); c++) {
//get set names
- string setA = namesOfGroupCombos[c][0];
+ string setA = namesOfGroupCombos[c][0];
string setB = namesOfGroupCombos[c][1];
-
+
//get filename
- string outputFileName = outputDir + m->getRootName(m->getSimpleName(globaldata->inputFileName)) + thisLookUp[0]->getLabel() + "." + setA + "-" + setB + ".metastats";
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
+ variables["[distance]"] = thisLookUp[0]->getLabel();
+ variables["[group]"] = setA + "-" + setB;
+ string outputFileName = getOutputFileName("metastats",variables);
outputNames.push_back(outputFileName); outputTypes["metastats"].push_back(outputFileName);
- int nameLength = outputFileName.length();
- char * output = new char[nameLength];
- strcpy(output, outputFileName.c_str());
+ //int nameLength = outputFileName.length();
+ //char * output = new char[nameLength];
+ //strcpy(output, outputFileName.c_str());
//build matrix from shared rabunds
- double** data;
- data = new double*[thisLookUp[0]->getNumBins()];
+ //double** data;
+ //data = new double*[thisLookUp[0]->getNumBins()];
+
+ vector< vector<double> > data2; data2.resize(thisLookUp[0]->getNumBins());
vector<SharedRAbundVector*> subset;
int setACount = 0;
m->mothurOut("Missing shared info for " + setA + " or " + setB + ". Skipping comparison."); m->mothurOutEndLine();
outputNames.pop_back();
}else {
+
//fill data
for (int j = 0; j < thisLookUp[0]->getNumBins(); j++) {
- data[j] = new double[subset.size()];
+ //data[j] = new double[subset.size()];
+ data2[j].resize(subset.size(), 0.0);
+
for (int i = 0; i < subset.size(); i++) {
- data[j][i] = (subset[i]->getAbundance(j));
+ data2[j][i] = (subset[i]->getAbundance(j));
}
}
m->mothurOut("Comparing " + setA + " and " + setB + "..."); m->mothurOutEndLine();
- metastat_main(output, thisLookUp[0]->getNumBins(), subset.size(), threshold, iters, data, setACount);
+ //metastat_main(output, thisLookUp[0]->getNumBins(), subset.size(), threshold, iters, data, setACount);
+ if (convertSharedToInput) { convertToInput(subset, outputFileName); }
+ m->mothurOutEndLine();
+ MothurMetastats mothurMeta(threshold, iters);
+ mothurMeta.runMetastats(outputFileName , data2, setACount);
+ m->mothurOutEndLine();
m->mothurOutEndLine();
}
//free memory
- delete output;
- for(int i = 0; i < thisLookUp[0]->getNumBins(); i++) { delete[] data[i]; }
- delete[] data;
+ //delete output;
+ //for(int i = 0; i < thisLookUp[0]->getNumBins(); i++) { delete[] data[i]; }
+ //delete[] data;
}
return 0;
}
}
//**********************************************************************************************************************
-int MetaStatsCommand::eliminateZeroOTUS(vector<SharedRAbundVector*>& thislookup) {
+/*Metastats files look like:
+ 13_0 14_0 13_52 14_52 70S 71S 72S M1 M2 M3 C11 C12 C21 C15 C16 C19 C3 C4 C9
+ Alphaproteobacteria 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0
+ Mollicutes 0 0 2 0 0 59 5 11 4 1 0 2 8 1 0 1 0 3 0
+ Verrucomicrobiae 0 0 0 0 0 1 6 0 0 0 0 0 0 0 0 0 0 0 0
+ Deltaproteobacteria 0 0 0 0 0 6 1 0 1 0 1 1 7 0 0 0 0 0 0
+ Cyanobacteria 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
+ Epsilonproteobacteria 0 0 0 0 0 0 0 0 6 0 0 3 1 0 0 0 0 0 0
+ Clostridia 75 65 207 226 801 280 267 210 162 197 81 120 106 148 120 94 84 98 121
+ Bacilli 3 2 16 8 21 52 31 70 46 65 4 28 5 23 62 26 20 30 25
+ Bacteroidetes (class) 21 25 22 64 226 193 296 172 98 55 19 149 201 85 50 76 113 92 82
+ Gammaproteobacteria 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 0
+ TM7_genera_incertae_sedis 0 0 0 0 0 0 0 0 1 0 1 2 0 2 0 0 0 0 0
+ Actinobacteria (class) 1 1 1 2 0 0 0 9 3 7 1 1 1 3 1 2 1 2 3
+ Betaproteobacteria 0 0 3 3 0 0 9 1 1 0 1 2 3 1 1 0 0 0 0
+*/
+//this function is just used to convert files to test the differences between the metastats version and mothurs version
+int MetaStatsCommand::convertToShared(string filename) {
try {
-
- vector<SharedRAbundVector*> newLookup;
- for (int i = 0; i < thislookup.size(); i++) {
- SharedRAbundVector* temp = new SharedRAbundVector();
- temp->setLabel(thislookup[i]->getLabel());
- temp->setGroup(thislookup[i]->getGroup());
- newLookup.push_back(temp);
- }
-
- //for each bin
- for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
- if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
-
- //look at each sharedRabund and make sure they are not all zero
- bool allZero = true;
- for (int j = 0; j < thislookup.size(); j++) {
- if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
- }
-
- //if they are not all zero add this bin
- if (!allZero) {
- for (int j = 0; j < thislookup.size(); j++) {
- newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
- }
- }
- }
-
- for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
-
- thislookup = newLookup;
-
- return 0;
-
+ ifstream in;
+ m->openInputFile(filename, in);
+
+ string header = m->getline(in); m->gobble(in);
+
+ vector<string> groups = m->splitWhiteSpace(header);
+ vector<SharedRAbundVector*> newLookup;
+ for (int i = 0; i < groups.size(); i++) {
+ SharedRAbundVector* temp = new SharedRAbundVector();
+ temp->setLabel("0.03");
+ temp->setGroup(groups[i]);
+ newLookup.push_back(temp);
+ }
+
+ int otuCount = 0;
+ while (!in.eof()) {
+ if (m->control_pressed) { break; }
+
+ string otuname;
+ in >> otuname; m->gobble(in);
+ otuCount++;
+
+ for (int i = 0; i < groups.size(); i++) {
+ int temp;
+ in >> temp; m->gobble(in);
+ newLookup[i]->push_back(temp, groups[i]);
+ }
+ m->gobble(in);
+ }
+ in.close();
+
+ ofstream out;
+ m->openOutputFile(filename+".shared", out);
+
+ out << "label\tgroup\tnumOTUs\t";
+
+ string snumBins = toString(otuCount);
+ for (int i = 0; i < otuCount; i++) {
+ string binLabel = "Otu";
+ string sbinNumber = toString(i+1);
+ if (sbinNumber.length() < snumBins.length()) {
+ int diff = snumBins.length() - sbinNumber.length();
+ for (int h = 0; h < diff; h++) { binLabel += "0"; }
+ }
+ binLabel += sbinNumber;
+ out << binLabel << '\t';
+ }
+ out << endl;
+
+ for (int i = 0; i < groups.size(); i++) {
+ out << "0.03" << '\t' << groups[i] << '\t';
+ newLookup[i]->print(out);
+ }
+ out.close();
+
+ cout << filename+".shared" << endl;
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "MetaStatsCommand", "convertToShared");
+ exit(1);
}
+}
+//**********************************************************************************************************************
+
+int MetaStatsCommand::convertToInput(vector<SharedRAbundVector*>& subset, string thisfilename) {
+ try {
+ ofstream out;
+ m->openOutputFile(thisfilename+".matrix", out);
+
+ out << "\t";
+ for (int i = 0; i < subset.size()-1; i++) {
+ out << subset[i]->getGroup() << '\t';
+ }
+ out << subset[subset.size()-1]->getGroup() << endl;
+
+ for (int i = 0; i < subset[0]->getNumBins(); i++) {
+ out << m->currentSharedBinLabels[i] << '\t';
+ for (int j = 0; j < subset.size()-1; j++) {
+ out << subset[j]->getAbundance(i) << '\t';
+ }
+ out << subset[subset.size()-1]->getAbundance(i) << endl;
+ }
+ out.close();
+
+ return 0;
+ }
catch(exception& e) {
- m->errorOut(e, "MetaStatsCommand", "eliminateZeroOTUS");
+ m->errorOut(e, "MetaStatsCommand", "convertToInput");
exit(1);
}
}
+
//**********************************************************************************************************************
projects / mothur.git / blobdiff