1 /*
2 * metastatscommand.cpp
3 * Mothur
4 *
5 * Created by westcott on 9/16/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
7 *
8 */
9
10 #include "metastatscommand.h"
11
12
13 //CommandParameter(string n, string t, string o, string d, string only, string atLeast, string linked, string opt, bool m, bool r, bool i) : name(n), type(t), options(o), optionsDefault(d), chooseOnlyOneGroup(only), chooseAtLeastOneGroup(atLeast), linkedGroup(linked), outputTypes(opt),multipleSelectionAllowed(m), required(r), important(i) {}
14
15 //**********************************************************************************************************************
setParameters()16 vector<string> MetaStatsCommand::setParameters(){
17 try {
18 CommandParameter pshared("shared", "InputTypes", "", "", "shared-clr", "none", "none","metastats",false,false,true); parameters.push_back(pshared);
19 CommandParameter pclr("clr", "InputTypes", "", "", "shared-clr", "none", "none","metastats",false,false,true); parameters.push_back(pclr);
20 CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(pdesign);
21 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
22 CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
23 CommandParameter pthreshold("threshold", "Number", "", "0.05", "", "", "","",false,false); parameters.push_back(pthreshold);
24 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
25 CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
26 CommandParameter psets("sets", "String", "", "", "", "", "","",false,false); parameters.push_back(psets);
27 CommandParameter pseed("seed", "Number", "", "0", "", "", "","",false,false); parameters.push_back(pseed);
28 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
29 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
30
31 abort = false; calledHelp = false; allLines = true;
32
33 vector<string> tempOutNames;
34 outputTypes["metastats"] = tempOutNames;
35
36 vector<string> myArray;
37 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
38 return myArray;
39 }
40 catch(exception& e) {
41 m->errorOut(e, "MetaStatsCommand", "setParameters");
42 exit(1);
43 }
44 }
45 //**********************************************************************************************************************
getHelpString()46 string MetaStatsCommand::getHelpString(){
47 try {
48 string helpString = "";
49 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";
50 helpString += "The metastats command outputs a .metastats file. \n";
51 helpString += "The metastats command parameters are shared, clr, iters, threshold, groups, label, design, sets and processors. The shared or clr and design parameters are required, unless you have valid current files.\n";
52 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";
53 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";
54 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";
55 helpString += "The iters parameter allows you to set number of bootstrap permutations for estimating null distribution of t statistic. The default is 1000. \n";
56 helpString += "The threshold parameter allows you to set the significance level to reject null hypotheses (default 0.05).\n";
57 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";
58 helpString += "The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n";
59 helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n";
60 helpString += "The metastats command should be in the following format: metastats(design=yourDesignFile).\n";
61 helpString += "Example metastats(design=temp.design, groups=A-B-C).\n";
62 helpString += "The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n";
63
64 return helpString;
65 }
66 catch(exception& e) {
67 m->errorOut(e, "MetaStatsCommand", "getHelpString");
68 exit(1);
69 }
70 }
71 //**********************************************************************************************************************
getOutputPattern(string type)72 string MetaStatsCommand::getOutputPattern(string type) {
73 try {
74 string pattern = "";
75
76 if (type == "metastats") { pattern = "[filename],[distance],[group],metastats"; }
77 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->setControl_pressed(true); }
78
79 return pattern;
80 }
81 catch(exception& e) {
82 m->errorOut(e, "MetaStatsCommand", "getOutputPattern");
83 exit(1);
84 }
85 }
86 //**********************************************************************************************************************
MetaStatsCommand(string option)87 MetaStatsCommand::MetaStatsCommand(string option) : Command() {
88 try {
89 if(option == "help") { help(); abort = true; calledHelp = true; }
90 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
91 else if(option == "category") { abort = true; calledHelp = true; }
92
93 else {
94 OptionParser parser(option, setParameters());
95 map<string,string> parameters = parser.getParameters();
96
97 ValidParameters validParameter;
98 sharedfile = validParameter.validFile(parameters, "shared");
99 if (sharedfile == "not open") { abort = true; }
100 else if (sharedfile == "not found") { sharedfile = ""; }
101 else { current->setSharedFile(sharedfile); inputfile = sharedfile; format = "sharedfile"; }
102
103 clrfile = validParameter.validFile(parameters, "clr");
104 if (clrfile == "not open") { abort = true; }
105 else if (clrfile == "not found") { clrfile = ""; }
106 else {
107 current->setCLRFile(clrfile); inputfile = clrfile; format = "clrfile";
108 m->mothurOut("[NOTE]: When using a clr file mothur will run the fisher exact test with the floor of the values generated.\n");
109 }
110
111 if ((sharedfile == "") && (clrfile == "")) {
112 //is there are current file available for any of these?
113 //give priority to shared, then list, then rabund, then sabund
114 //if there is a current shared file, use it
115 sharedfile = current->getSharedFile();
116 if (sharedfile != "") { inputfile = sharedfile; format = "sharedfile"; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter.\n"); }
117 else {
118 clrfile = current->getCLRFile();
119 if (clrfile != "") { inputfile = clrfile; format = "clrfile"; m->mothurOut("Using " + clrfile + " as input file for the clr parameter.\n"); m->mothurOut("[NOTE]: When using a clr file mothur will run the fisher exact test with the floor of the values generated.\n"); }
120 else { m->mothurOut("No valid current files. You must provide a clrfile or shared file.\n"); abort = true; }
121 }
122 }
123
124 //check for required parameters
125 designfile = validParameter.validFile(parameters, "design");
126 if (designfile == "not open") { abort = true; }
127 else if (designfile == "not found") {
128 //if there is a current design file, use it
129 designfile = current->getDesignFile();
130 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter.\n"); }
131 else { m->mothurOut("You have no current designfile and the design parameter is required.\n"); abort = true; }
132 }else { current->setDesignFile(designfile); }
133
134
135 if (outputdir == ""){
136 outputdir += util.hasPath(inputfile);
137 }
138
139 //check for optional parameter and set defaults
140 // ...at some point should added some additional type checking...
141 label = validParameter.valid(parameters, "label");
142 if (label == "not found") { label = ""; }
143 else {
144 if(label != "all") { util.splitAtDash(label, labels); allLines = false; }
145 else { allLines = true; }
146 }
147
148 groups = validParameter.valid(parameters, "groups");
149 if (groups == "not found") { groups = ""; pickedGroups = false; }
150 else {
151 pickedGroups = true;
152 util.splitAtDash(groups, Groups);
153 if (Groups.size() != 0) { if (Groups[0]== "all") { Groups.clear(); } }
154 }
155
156 sets = validParameter.valid(parameters, "sets");
157 if (sets == "not found") { sets = ""; }
158 else {
159 util.splitAtDash(sets, Sets);
160 if (Sets.size() != 0) { if (Sets[0] != "all") { Groups.clear(); } }
161 }
162
163 string temp = validParameter.valid(parameters, "iters"); if (temp == "not found") { temp = "1000"; }
164 util.mothurConvert(temp, iters);
165
166 temp = validParameter.valid(parameters, "threshold"); if (temp == "not found") { temp = "0.05"; }
167 util.mothurConvert(temp, threshold);
168
169 temp = validParameter.valid(parameters, "processors"); if (temp == "not found"){ temp = current->getProcessors(); }
170 processors = current->setProcessors(temp);
171 }
172
173 }
174 catch(exception& e) {
175 m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
176 exit(1);
177 }
178 }
179 //**********************************************************************************************************************
180
execute()181 int MetaStatsCommand::execute(){
182 try {
183
184 if (abort) { if (calledHelp) { return 0; } return 2; }
185
186 DesignMap* designMap = new DesignMap(designfile); if (m->getControl_pressed()) { delete designMap; return 0; }
187
188 InputData input(inputfile, format, Groups);
189 set<string> processedLabels;
190 set<string> userLabels = labels;
191 string lastLabel = "";
192
193 SharedRAbundVectors* lookup = NULL; SharedCLRVectors* clr = NULL;
194
195 if (format == "sharedfile") {
196 lookup = util.getNextShared(input, allLines, userLabels, processedLabels, lastLabel);
197 Groups = lookup->getNamesGroups();
198 }else {
199 clr = util.getNextCLR(input, allLines, userLabels, processedLabels, lastLabel);
200 Groups = clr->getNamesGroups();
201 }
202
203 if (Sets.size() == 0) { Sets = designMap->getCategory(); }
204 int numGroups = (int)Sets.size();
205 for (int a=0; a<numGroups; a++) {
206 for (int l = 0; l < a; l++) {
207 vector<string> groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]);
208 namesOfGroupCombos.push_back(groups);
209 }
210 }
211
212 if (numGroups == 2) { processors = 1; }
213 else if (numGroups < 2) { m->mothurOut("[ERROR]: Not enough sets, I need at least 2 valid sets. Unable to complete command.\n"); m->setControl_pressed(true); }
214
215 while ((lookup != NULL) || (clr != NULL)){
216
217 if (m->getControl_pressed()) { if (lookup != NULL) { delete lookup; } if (clr != NULL) { delete clr; }break; }
218
219 if (format == "sharedfile") {
220 process(lookup, designMap); delete lookup;
221 lookup = util.getNextShared(input, allLines, userLabels, processedLabels, lastLabel);
222 }
223 else {
224 process(clr, designMap); delete clr;
225 clr = util.getNextCLR(input, allLines, userLabels, processedLabels, lastLabel);
226 }
227 }
228
229 delete designMap;
230 if (m->getControl_pressed()) { outputTypes.clear(); if (lookup != NULL) { delete lookup; } if (clr != NULL) { delete clr; } for (int i = 0; i < outputNames.size(); i++) { util.mothurRemove(outputNames[i]); } return 0; }
231
232 m->mothurOut("\nOutput File Names: \n");
233 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i] +"\n"); } m->mothurOutEndLine();
234
235 return 0;
236 }
237 catch(exception& e) {
238 m->errorOut(e, "MetaStatsCommand", "execute");
239 exit(1);
240 }
241 }
242 /**************************************************************************************************/
243 struct metastatsData {
244 SharedRAbundVectors* thisLookUp;
245 SharedCLRVectors* thisCLR;
246 vector< vector<string> > namesOfGroupCombos;
247 vector<string> designMapGroups, outputNames;
248 int start, num, iters, count;
249 float threshold;
250 Utils util;
251 MothurOut* m;
252
metastatsDatametastatsData253 metastatsData(){}
metastatsDatametastatsData254 metastatsData(int st, int en, vector<string> on, vector< vector<string> > ns, SharedRAbundVectors*& lu, vector<string> dg, int i, float thr) {
255 m = MothurOut::getInstance();
256 outputNames = on;
257 start = st;
258 num = en;
259 namesOfGroupCombos = ns;
260 thisLookUp = lu;
261 designMapGroups = dg;
262 iters = i;
263 threshold = thr;
264 count=0;
265 thisCLR = NULL;
266 }
267
metastatsDatametastatsData268 metastatsData(int st, int en, vector<string> on, vector< vector<string> > ns, SharedCLRVectors*& lu, vector<string> dg, int i, float thr) {
269 m = MothurOut::getInstance();
270 outputNames = on;
271 start = st;
272 num = en;
273 namesOfGroupCombos = ns;
274 thisCLR = lu;
275 designMapGroups = dg;
276 iters = i;
277 threshold = thr;
278 count=0;
279 thisLookUp = NULL;
280 }
281 };
282 //**********************************************************************************************************************
driverShared(metastatsData * params)283 int driverShared(metastatsData* params) {
284 try {
285
286 vector<string> thisLookupNames = params->thisLookUp->getNamesGroups();
287 vector<SharedRAbundVector*> thisLookupRabunds = params->thisLookUp->getSharedRAbundVectors();
288
289 //for each combo
290 for (int c = params->start; c < (params->start+params->num); c++) {
291
292 //get set names
293 string setA = params->namesOfGroupCombos[c][0];
294 string setB = params->namesOfGroupCombos[c][1];
295 string outputFileName = params->outputNames[c];
296
297 vector< vector<double> > data2; data2.resize(params->thisLookUp->getNumBins());
298
299 vector<SharedRAbundVector*> subset;
300 vector<string> subsetGroups;
301 int setACount = 0; int setBCount = 0;
302 for (int i = 0; i < params->thisLookUp->size(); i++) {
303 string thisGroup = thisLookupNames[i];
304 if (params->designMapGroups[i] == setB) {
305 subset.push_back(thisLookupRabunds[i]);
306 subsetGroups.push_back(thisGroup);
307 setBCount++;
308 }else if (params->designMapGroups[i] == setA) {
309 subset.insert(subset.begin()+setACount, thisLookupRabunds[i]);
310 subsetGroups.insert(subsetGroups.begin()+setACount, thisGroup);
311 setACount++;
312 }
313 }
314
315 if ((setACount == 0) || (setBCount == 0)) { params->m->mothurOut("Missing shared info for " + setA + " or " + setB + ". Skipping comparison.\n"); }
316 else {
317 for (int j = 0; j < params->thisLookUp->getNumBins(); j++) {
318 data2[j].resize(subset.size(), 0.0);
319 for (int i = 0; i < subset.size(); i++) { data2[j][i] = (subset[i]->get(j)); }
320 }
321
322 params->m->mothurOut("\nComparing " + setA + " and " + setB + "...\n");
323 MothurMetastats mothurMeta(params->threshold, params->iters);
324 mothurMeta.runMetastats(outputFileName , data2, setACount, params->thisLookUp->getOTUNames(), true);
325 params->m->mothurOutEndLine();
326 }
327 }
328
329 for(int i = 0; i < thisLookupRabunds.size(); i++) { delete thisLookupRabunds[i]; }
330
331 return 0;
332
333 }
334 catch(exception& e) {
335 params->m->errorOut(e, "MetaStatsCommand", "driver");
336 exit(1);
337 }
338 }
339 //**********************************************************************************************************************
process(SharedRAbundVectors * & thisLookUp,DesignMap * & designMap)340 int MetaStatsCommand::process(SharedRAbundVectors*& thisLookUp, DesignMap*& designMap){
341 try {
342 vector<linePair> lines;
343 vector<std::thread*> workerThreads;
344 vector<metastatsData*> data;
345
346 int remainingPairs = namesOfGroupCombos.size();
347 int startIndex = 0;
348 vector<string> thisLabelsOutputFiles;
349 for (int remainingProcessors = processors; remainingProcessors > 0; remainingProcessors--) {
350 int numPairs = remainingPairs; //case for last processor
351 if (remainingProcessors != 1) { numPairs = ceil(remainingPairs / remainingProcessors); }
352 lines.push_back(linePair(startIndex, numPairs)); //startIndex, numPairs
353
354 for (int i = startIndex; i < startIndex+numPairs; i++) {
355 //get set names
356 string setA = namesOfGroupCombos[i][0];
357 string setB = namesOfGroupCombos[i][1];
358
359 //get filename
360 map<string, string> variables;
361 variables["[filename]"] = outputdir + util.getRootName(util.getSimpleName(inputfile));
362 variables["[distance]"] = thisLookUp->getLabel();
363 variables["[group]"] = setA + "-" + setB;
364 string outputFileName = getOutputFileName("metastats",variables);
365 outputNames.push_back(outputFileName); outputTypes["metastats"].push_back(outputFileName);
366 thisLabelsOutputFiles.push_back(outputFileName);
367 }
368
369 startIndex = startIndex + numPairs;
370 remainingPairs = remainingPairs - numPairs;
371 }
372
373 vector<string> designMapGroups = thisLookUp->getNamesGroups();
374 for (int j = 0; j < designMapGroups.size(); j++) { designMapGroups[j] = designMap->get(designMapGroups[j]); }
375
376 //Lauch worker threads
377 for (int i = 0; i < processors-1; i++) {
378 //make copy of lookup so we don't get access violations
379 SharedRAbundVectors* newLookup = new SharedRAbundVectors(*thisLookUp);
380
381 metastatsData* dataBundle = new metastatsData(lines[i+1].start, lines[i+1].end, thisLabelsOutputFiles, namesOfGroupCombos, newLookup, designMapGroups, iters, threshold);
382 data.push_back(dataBundle);
383
384 std::thread* thisThread = new std::thread(driverShared, dataBundle);
385 workerThreads.push_back(thisThread);
386 }
387
388 metastatsData* dataBundle = new metastatsData(lines[0].start, lines[0].end, thisLabelsOutputFiles, namesOfGroupCombos, thisLookUp, designMapGroups, iters, threshold);
389 driverShared(dataBundle);
390
391 for (int i = 0; i < processors-1; i++) {
392 workerThreads[i]->join();
393
394 delete data[i]->thisLookUp;
395 delete data[i];
396 delete workerThreads[i];
397 }
398 delete dataBundle;
399 return 0;
400 }
401 catch(exception& e) {
402 m->errorOut(e, "MetaStatsCommand", "process");
403 exit(1);
404 }
405 }
406 //**********************************************************************************************************************
driverCLR(metastatsData * params)407 int driverCLR(metastatsData* params) {
408 try {
409
410 vector<string> thisLookupNames = params->thisCLR->getNamesGroups();
411 vector<SharedCLRVector*> thisCLRRabunds = params->thisCLR->getSharedCLRVectors();
412
413 //for each combo
414 for (int c = params->start; c < (params->start+params->num); c++) {
415
416 //get set names
417 string setA = params->namesOfGroupCombos[c][0];
418 string setB = params->namesOfGroupCombos[c][1];
419 string outputFileName = params->outputNames[c];
420
421 vector< vector<double> > data2; data2.resize(params->thisCLR->getNumBins());
422
423 vector<SharedCLRVector*> subset;
424 int setACount = 0; int setBCount = 0;
425 for (int i = 0; i < params->thisCLR->size(); i++) {
426 string thisGroup = thisLookupNames[i];
427 if (params->designMapGroups[i] == setB) {
428 subset.push_back(thisCLRRabunds[i]);
429 setBCount++;
430 }else if (params->designMapGroups[i] == setA) {
431 subset.insert(subset.begin()+setACount, thisCLRRabunds[i]);
432 setACount++;
433 }
434 }
435
436 if ((setACount == 0) || (setBCount == 0)) { params->m->mothurOut("Missing shared info for " + setA + " or " + setB + ". Skipping comparison.\n"); }
437 else {
438 for (int j = 0; j < params->thisCLR->getNumBins(); j++) {
439 data2[j].resize(subset.size(), 0.0);
440 for (int i = 0; i < subset.size(); i++) { data2[j][i] = (subset[i]->get(j)); }
441 }
442
443 params->m->mothurOut("\nComparing " + setA + " and " + setB + "...\n");
444 MothurMetastats mothurMeta(params->threshold, params->iters);
445 mothurMeta.runMetastats(outputFileName, data2, setACount, params->thisCLR->getOTUNames(), false);
446 params->m->mothurOutEndLine();
447 }
448 }
449
450 for(int i = 0; i < thisCLRRabunds.size(); i++) { delete thisCLRRabunds[i]; }
451
452 return 0;
453
454 }
455 catch(exception& e) {
456 params->m->errorOut(e, "MetaStatsCommand", "driver");
457 exit(1);
458 }
459 }
460 //**********************************************************************************************************************
process(SharedCLRVectors * & thisCLR,DesignMap * & designMap)461 int MetaStatsCommand::process(SharedCLRVectors*& thisCLR, DesignMap*& designMap){
462 try {
463 vector<linePair> lines;
464 vector<std::thread*> workerThreads;
465 vector<metastatsData*> data;
466
467 int remainingPairs = namesOfGroupCombos.size();
468 int startIndex = 0;
469 vector<string> thisLabelsOutputFiles;
470 for (int remainingProcessors = processors; remainingProcessors > 0; remainingProcessors--) {
471 int numPairs = remainingPairs; //case for last processor
472 if (remainingProcessors != 1) { numPairs = ceil(remainingPairs / remainingProcessors); }
473 lines.push_back(linePair(startIndex, numPairs)); //startIndex, numPairs
474
475 for (int i = startIndex; i < startIndex+numPairs; i++) {
476 //get set names
477 string setA = namesOfGroupCombos[i][0];
478 string setB = namesOfGroupCombos[i][1];
479
480 //get filename
481 map<string, string> variables;
482 variables["[filename]"] = outputdir + util.getRootName(util.getSimpleName(inputfile));
483 variables["[distance]"] = thisCLR->getLabel();
484 variables["[group]"] = setA + "-" + setB;
485 string outputFileName = getOutputFileName("metastats",variables);
486 outputNames.push_back(outputFileName); outputTypes["metastats"].push_back(outputFileName);
487 thisLabelsOutputFiles.push_back(outputFileName);
488 }
489
490 startIndex = startIndex + numPairs;
491 remainingPairs = remainingPairs - numPairs;
492 }
493
494 vector<string> designMapGroups = thisCLR->getNamesGroups();
495 for (int j = 0; j < designMapGroups.size(); j++) { designMapGroups[j] = designMap->get(designMapGroups[j]); }
496
497 //Lauch worker threads
498 for (int i = 0; i < processors-1; i++) {
499 //make copy of lookup so we don't get access violations
500 SharedCLRVectors* newLookup = new SharedCLRVectors(*thisCLR);
501
502 metastatsData* dataBundle = new metastatsData(lines[i+1].start, lines[i+1].end, thisLabelsOutputFiles, namesOfGroupCombos, newLookup, designMapGroups, iters, threshold);
503 data.push_back(dataBundle);
504
505 std::thread* thisThread = new std::thread(driverCLR, dataBundle);
506 workerThreads.push_back(thisThread);
507 }
508
509 metastatsData* dataBundle = new metastatsData(lines[0].start, lines[0].end, thisLabelsOutputFiles, namesOfGroupCombos, thisCLR, designMapGroups, iters, threshold);
510 driverCLR(dataBundle);
511
512 for (int i = 0; i < processors-1; i++) {
513 workerThreads[i]->join();
514
515 delete data[i]->thisLookUp;
516 delete data[i];
517 delete workerThreads[i];
518 }
519 delete dataBundle;
520 return 0;
521 }
522 catch(exception& e) {
523 m->errorOut(e, "MetaStatsCommand", "process");
524 exit(1);
525 }
526 }
527 //**********************************************************************************************************************/
528
529
530