5 * Created by westcott on 9/16/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "metastatscommand.h"
11 #include "sharedutilities.h"
12 #include "sharedrabundfloatvector.h"
15 //**********************************************************************************************************************
16 vector<string> MetaStatsCommand::setParameters(){
18 CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","metastats",false,true,true); parameters.push_back(pshared);
19 CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(pdesign);
20 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
21 CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
22 CommandParameter pthreshold("threshold", "Number", "", "0.05", "", "", "","",false,false); parameters.push_back(pthreshold);
23 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
24 CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
25 CommandParameter psets("sets", "String", "", "", "", "", "","",false,false); parameters.push_back(psets);
26 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
27 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
29 vector<string> myArray;
30 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
34 m->errorOut(e, "MetaStatsCommand", "setParameters");
38 //**********************************************************************************************************************
39 string MetaStatsCommand::getHelpString(){
41 string helpString = "";
42 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";
43 helpString += "The metastats command outputs a .metastats file. \n";
44 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";
45 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";
46 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";
47 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";
48 helpString += "The iters parameter allows you to set number of bootstrap permutations for estimating null distribution of t statistic. The default is 1000. \n";
49 helpString += "The threshold parameter allows you to set the significance level to reject null hypotheses (default 0.05).\n";
50 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";
51 helpString += "The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n";
52 helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n";
53 helpString += "The metastats command should be in the following format: metastats(design=yourDesignFile).\n";
54 helpString += "Example metastats(design=temp.design, groups=A-B-C).\n";
55 helpString += "The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n";
56 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
60 m->errorOut(e, "MetaStatsCommand", "getHelpString");
64 //**********************************************************************************************************************
65 string MetaStatsCommand::getOutputPattern(string type) {
69 if (type == "metastats") { pattern = "[filename],[distance],[group],metastats"; }
70 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
75 m->errorOut(e, "MetaStatsCommand", "getOutputPattern");
79 //**********************************************************************************************************************
80 MetaStatsCommand::MetaStatsCommand(){
82 abort = true; calledHelp = true;
84 vector<string> tempOutNames;
85 outputTypes["metastats"] = tempOutNames;
88 m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
92 //**********************************************************************************************************************
94 MetaStatsCommand::MetaStatsCommand(string option) {
96 abort = false; calledHelp = false;
100 //allow user to run help
101 if(option == "help") { help(); abort = true; calledHelp = true; }
102 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
105 vector<string> myArray = setParameters();
107 OptionParser parser(option);
108 map<string,string> parameters = parser.getParameters();
110 ValidParameters validParameter;
112 //check to make sure all parameters are valid for command
113 map<string,string>::iterator it;
114 for (it = parameters.begin(); it != parameters.end(); it++) {
115 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
118 //initialize outputTypes
119 vector<string> tempOutNames;
120 outputTypes["metastats"] = tempOutNames;
123 //if the user changes the input directory command factory will send this info to us in the output parameter
124 string inputDir = validParameter.validFile(parameters, "inputdir", false);
125 if (inputDir == "not found"){ inputDir = ""; }
128 it = parameters.find("design");
129 //user has given a template file
130 if(it != parameters.end()){
131 path = m->hasPath(it->second);
132 //if the user has not given a path then, add inputdir. else leave path alone.
133 if (path == "") { parameters["design"] = inputDir + it->second; }
136 it = parameters.find("shared");
137 //user has given a template file
138 if(it != parameters.end()){
139 path = m->hasPath(it->second);
140 //if the user has not given a path then, add inputdir. else leave path alone.
141 if (path == "") { parameters["shared"] = inputDir + it->second; }
146 //check for required parameters
147 sharedfile = validParameter.validFile(parameters, "shared", true);
148 if (sharedfile == "not open") { abort = true; }
149 else if (sharedfile == "not found") { //if there is a current shared file, use it
150 sharedfile = m->getSharedFile();
151 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
152 else { m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
153 }else { m->setSharedFile(sharedfile); }
155 //check for required parameters
156 designfile = validParameter.validFile(parameters, "design", true);
157 if (designfile == "not open") { abort = true; }
158 else if (designfile == "not found") {
159 //if there is a current design file, use it
160 designfile = m->getDesignFile();
161 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
162 else { m->mothurOut("You have no current designfile and the design parameter is required."); m->mothurOutEndLine(); abort = true; }
163 }else { m->setDesignFile(designfile); }
165 //if the user changes the output directory command factory will send this info to us in the output parameter
166 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
168 outputDir += m->hasPath(sharedfile); //if user entered a file with a path then preserve it
171 //check for optional parameter and set defaults
172 // ...at some point should added some additional type checking...
173 label = validParameter.validFile(parameters, "label", false);
174 if (label == "not found") { label = ""; }
176 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
177 else { allLines = 1; }
180 groups = validParameter.validFile(parameters, "groups", false);
181 if (groups == "not found") { groups = ""; pickedGroups = false; }
184 m->splitAtDash(groups, Groups);
185 m->setGroups(Groups);
188 sets = validParameter.validFile(parameters, "sets", false);
189 if (sets == "not found") { sets = ""; }
191 m->splitAtDash(sets, Sets);
195 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
196 m->mothurConvert(temp, iters);
198 temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "0.05"; }
199 m->mothurConvert(temp, threshold);
201 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
202 m->setProcessors(temp);
203 m->mothurConvert(temp, processors);
207 catch(exception& e) {
208 m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
212 //**********************************************************************************************************************
214 int MetaStatsCommand::execute(){
217 if (abort == true) { if (calledHelp) { return 0; } return 2; }
219 //just used to convert files to test metastats online
220 /****************************************************/
221 bool convertInputToShared = false;
222 convertSharedToInput = false;
223 if (convertInputToShared) { convertToShared(sharedfile); return 0; }
224 /****************************************************/
226 designMap = new GroupMap(designfile);
227 designMap->readDesignMap();
229 input = new InputData(sharedfile, "sharedfile");
230 lookup = input->getSharedRAbundVectors();
231 string lastLabel = lookup[0]->getLabel();
233 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
234 set<string> processedLabels;
235 set<string> userLabels = labels;
237 //setup the pairwise comparions of sets for metastats
238 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
239 //make sure sets are all in designMap
240 SharedUtil* util = new SharedUtil();
241 vector<string> dGroups = designMap->getNamesOfGroups();
242 util->setGroups(Sets, dGroups);
245 int numGroups = Sets.size();
246 for (int a=0; a<numGroups; a++) {
247 for (int l = 0; l < a; l++) {
248 vector<string> groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]);
249 namesOfGroupCombos.push_back(groups);
255 if (numGroups == 2) { processors = 1; }
256 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; }
259 int numPairs = namesOfGroupCombos.size();
260 int numPairsPerProcessor = numPairs / processors;
262 for (int i = 0; i < processors; i++) {
263 int startPos = i * numPairsPerProcessor;
264 if(i == processors - 1){
265 numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
267 lines.push_back(linePair(startPos, numPairsPerProcessor));
271 //as long as you are not at the end of the file or done wih the lines you want
272 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
274 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; }
276 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
278 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
281 processedLabels.insert(lookup[0]->getLabel());
282 userLabels.erase(lookup[0]->getLabel());
285 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
286 string saveLabel = lookup[0]->getLabel();
288 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
289 lookup = input->getSharedRAbundVectors(lastLabel);
290 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
294 processedLabels.insert(lookup[0]->getLabel());
295 userLabels.erase(lookup[0]->getLabel());
297 //restore real lastlabel to save below
298 lookup[0]->setLabel(saveLabel);
301 lastLabel = lookup[0]->getLabel();
302 //prevent memory leak
303 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; }
305 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; }
307 //get next line to process
308 lookup = input->getSharedRAbundVectors();
311 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; }
313 //output error messages about any remaining user labels
314 set<string>::iterator it;
315 bool needToRun = false;
316 for (it = userLabels.begin(); it != userLabels.end(); it++) {
317 m->mothurOut("Your file does not include the label " + *it);
318 if (processedLabels.count(lastLabel) != 1) {
319 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
322 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
326 //run last label if you need to
327 if (needToRun == true) {
328 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
329 lookup = input->getSharedRAbundVectors(lastLabel);
331 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
335 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
338 //reset groups parameter
343 if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0;}
345 m->mothurOutEndLine();
346 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
347 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
348 m->mothurOutEndLine();
352 catch(exception& e) {
353 m->errorOut(e, "MetaStatsCommand", "execute");
357 //**********************************************************************************************************************
359 int MetaStatsCommand::process(vector<SharedRAbundVector*>& thisLookUp){
364 driver(0, namesOfGroupCombos.size(), thisLookUp);
367 vector<int> processIDS;
368 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
369 //loop through and create all the processes you want
370 while (process != processors) {
374 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
377 driver(lines[process].start, lines[process].num, thisLookUp);
380 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
381 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
387 driver(lines[0].start, lines[0].num, thisLookUp);
389 //force parent to wait until all the processes are done
390 for (int i=0;i<(processors-1);i++) {
391 int temp = processIDS[i];
396 //////////////////////////////////////////////////////////////////////////////////////////////////////
397 //Windows version shared memory, so be careful when passing variables through the summarySharedData struct.
398 //Above fork() will clone, so memory is separate, but that's not the case with windows,
399 //Taking advantage of shared memory to pass results vectors.
400 //////////////////////////////////////////////////////////////////////////////////////////////////////
402 vector<metastatsData*> pDataArray;
403 DWORD dwThreadIdArray[processors-1];
404 HANDLE hThreadArray[processors-1];
406 //Create processor worker threads.
407 for( int i=1; i<processors; i++ ){
409 //make copy of lookup so we don't get access violations
410 vector<SharedRAbundVector*> newLookup;
411 vector<string> designMapGroups;
412 for (int k = 0; k < thisLookUp.size(); k++) {
413 SharedRAbundVector* temp = new SharedRAbundVector();
414 temp->setLabel(thisLookUp[k]->getLabel());
415 temp->setGroup(thisLookUp[k]->getGroup());
416 newLookup.push_back(temp);
417 designMapGroups.push_back(designMap->getGroup(thisLookUp[k]->getGroup()));
421 for (int k = 0; k < thisLookUp[0]->getNumBins(); k++) {
422 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
423 for (int j = 0; j < thisLookUp.size(); j++) { newLookup[j]->push_back(thisLookUp[j]->getAbundance(k), thisLookUp[j]->getGroup()); }
426 // Allocate memory for thread data.
427 metastatsData* tempSum = new metastatsData(sharedfile, outputDir, m, lines[i].start, lines[i].num, namesOfGroupCombos, newLookup, designMapGroups, iters, threshold);
428 pDataArray.push_back(tempSum);
429 processIDS.push_back(i);
431 hThreadArray[i-1] = CreateThread(NULL, 0, MyMetastatsThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
435 driver(lines[0].start, lines[0].num, thisLookUp);
437 //Wait until all threads have terminated.
438 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
440 //Close all thread handles and free memory allocations.
441 for(int i=0; i < pDataArray.size(); i++){
442 if (pDataArray[i]->count != (pDataArray[i]->num)) {
443 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;
445 for (int j = 0; j < pDataArray[i]->thisLookUp.size(); j++) { delete pDataArray[i]->thisLookUp[j]; }
446 for (int j = 0; j < pDataArray[i]->outputNames.size(); j++) {
447 outputNames.push_back(pDataArray[i]->outputNames[j]);
448 outputTypes["metastats"].push_back(pDataArray[i]->outputNames[j]);
451 CloseHandle(hThreadArray[i]);
452 delete pDataArray[i];
461 catch(exception& e) {
462 m->errorOut(e, "MetaStatsCommand", "process");
466 //**********************************************************************************************************************
467 int MetaStatsCommand::driver(int start, int num, vector<SharedRAbundVector*>& thisLookUp) {
471 for (int c = start; c < (start+num); c++) {
474 string setA = namesOfGroupCombos[c][0];
475 string setB = namesOfGroupCombos[c][1];
478 map<string, string> variables;
479 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
480 variables["[distance]"] = thisLookUp[0]->getLabel();
481 variables["[group]"] = setA + "-" + setB;
482 string outputFileName = getOutputFileName("metastats",variables);
483 outputNames.push_back(outputFileName); outputTypes["metastats"].push_back(outputFileName);
484 //int nameLength = outputFileName.length();
485 //char * output = new char[nameLength];
486 //strcpy(output, outputFileName.c_str());
488 //build matrix from shared rabunds
490 //data = new double*[thisLookUp[0]->getNumBins()];
492 vector< vector<double> > data2; data2.resize(thisLookUp[0]->getNumBins());
494 vector<SharedRAbundVector*> subset;
497 for (int i = 0; i < thisLookUp.size(); i++) {
498 string thisGroup = thisLookUp[i]->getGroup();
500 //is this group for a set we want to compare??
501 //sorting the sets by putting setB at the back and setA in the front
502 if ((designMap->getGroup(thisGroup) == setB)) {
503 subset.push_back(thisLookUp[i]);
505 }else if ((designMap->getGroup(thisGroup) == setA)) {
506 subset.insert(subset.begin()+setACount, thisLookUp[i]);
511 if ((setACount == 0) || (setBCount == 0)) {
512 m->mothurOut("Missing shared info for " + setA + " or " + setB + ". Skipping comparison."); m->mothurOutEndLine();
513 outputNames.pop_back();
517 for (int j = 0; j < thisLookUp[0]->getNumBins(); j++) {
518 //data[j] = new double[subset.size()];
519 data2[j].resize(subset.size(), 0.0);
521 for (int i = 0; i < subset.size(); i++) {
522 data2[j][i] = (subset[i]->getAbundance(j));
526 m->mothurOut("Comparing " + setA + " and " + setB + "..."); m->mothurOutEndLine();
527 //metastat_main(output, thisLookUp[0]->getNumBins(), subset.size(), threshold, iters, data, setACount);
528 if (convertSharedToInput) { convertToInput(subset, outputFileName); }
530 m->mothurOutEndLine();
531 MothurMetastats mothurMeta(threshold, iters);
532 mothurMeta.runMetastats(outputFileName , data2, setACount);
533 m->mothurOutEndLine();
534 m->mothurOutEndLine();
539 //for(int i = 0; i < thisLookUp[0]->getNumBins(); i++) { delete[] data[i]; }
546 catch(exception& e) {
547 m->errorOut(e, "MetaStatsCommand", "driver");
551 //**********************************************************************************************************************
552 /*Metastats files look like:
553 13_0 14_0 13_52 14_52 70S 71S 72S M1 M2 M3 C11 C12 C21 C15 C16 C19 C3 C4 C9
554 Alphaproteobacteria 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0
555 Mollicutes 0 0 2 0 0 59 5 11 4 1 0 2 8 1 0 1 0 3 0
556 Verrucomicrobiae 0 0 0 0 0 1 6 0 0 0 0 0 0 0 0 0 0 0 0
557 Deltaproteobacteria 0 0 0 0 0 6 1 0 1 0 1 1 7 0 0 0 0 0 0
558 Cyanobacteria 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
559 Epsilonproteobacteria 0 0 0 0 0 0 0 0 6 0 0 3 1 0 0 0 0 0 0
560 Clostridia 75 65 207 226 801 280 267 210 162 197 81 120 106 148 120 94 84 98 121
561 Bacilli 3 2 16 8 21 52 31 70 46 65 4 28 5 23 62 26 20 30 25
562 Bacteroidetes (class) 21 25 22 64 226 193 296 172 98 55 19 149 201 85 50 76 113 92 82
563 Gammaproteobacteria 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 0
564 TM7_genera_incertae_sedis 0 0 0 0 0 0 0 0 1 0 1 2 0 2 0 0 0 0 0
565 Actinobacteria (class) 1 1 1 2 0 0 0 9 3 7 1 1 1 3 1 2 1 2 3
566 Betaproteobacteria 0 0 3 3 0 0 9 1 1 0 1 2 3 1 1 0 0 0 0
568 //this function is just used to convert files to test the differences between the metastats version and mothurs version
569 int MetaStatsCommand::convertToShared(string filename) {
572 m->openInputFile(filename, in);
574 string header = m->getline(in); m->gobble(in);
576 vector<string> groups = m->splitWhiteSpace(header);
577 vector<SharedRAbundFloatVector*> newLookup;
578 cout << groups.size() << endl;
579 for (int i = 0; i < groups.size(); i++) {
580 cout << "creating group " << groups[i] << endl;
581 SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
582 temp->setLabel("0.03");
583 temp->setGroup(groups[i]);
584 newLookup.push_back(temp);
589 if (m->control_pressed) { break; }
592 in >> otuname; m->gobble(in);
594 cout << otuname << endl;
595 for (int i = 0; i < groups.size(); i++) {
597 in >> temp; m->gobble(in);
598 newLookup[i]->push_back(temp, groups[i]);
605 m->openOutputFile(filename+".shared", out);
607 out << "label\tgroup\tnumOTUs\t";
609 string snumBins = toString(otuCount);
610 for (int i = 0; i < otuCount; i++) {
611 string binLabel = "Otu";
612 string sbinNumber = toString(i+1);
613 if (sbinNumber.length() < snumBins.length()) {
614 int diff = snumBins.length() - sbinNumber.length();
615 for (int h = 0; h < diff; h++) { binLabel += "0"; }
617 binLabel += sbinNumber;
618 out << binLabel << '\t';
622 for (int i = 0; i < groups.size(); i++) {
623 out << "0.03" << '\t' << groups[i] << '\t';
624 newLookup[i]->print(out);
628 cout << filename+".shared" << endl;
632 catch(exception& e) {
633 m->errorOut(e, "MetaStatsCommand", "convertToShared");
637 //**********************************************************************************************************************
639 int MetaStatsCommand::convertToInput(vector<SharedRAbundVector*>& subset, string thisfilename) {
642 m->openOutputFile(thisfilename+".matrix", out);
645 for (int i = 0; i < subset.size()-1; i++) {
646 out << subset[i]->getGroup() << '\t';
648 out << subset[subset.size()-1]->getGroup() << endl;
650 for (int i = 0; i < subset[0]->getNumBins(); i++) {
651 out << m->currentSharedBinLabels[i] << '\t';
652 for (int j = 0; j < subset.size()-1; j++) {
653 out << subset[j]->getAbundance(i) << '\t';
655 out << subset[subset.size()-1]->getAbundance(i) << endl;
659 cout << thisfilename+".matrix" << endl;
663 catch(exception& e) {
664 m->errorOut(e, "MetaStatsCommand", "convertToInput");
669 //**********************************************************************************************************************