5 * Created by westcott on 9/16/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "metastatscommand.h"
11 #include "metastats.h"
13 //**********************************************************************************************************************
15 MetaStatsCommand::MetaStatsCommand(string option) {
17 globaldata = GlobalData::getInstance();
22 //allow user to run help
23 if(option == "help") { help(); abort = true; }
26 //valid paramters for this command
27 string AlignArray[] = {"groups","label","outputdir","iters","threshold","g","design","inputdir"};
28 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
30 OptionParser parser(option);
31 map<string,string> parameters = parser.getParameters();
33 ValidParameters validParameter;
35 //check to make sure all parameters are valid for command
36 map<string,string>::iterator it;
37 for (it = parameters.begin(); it != parameters.end(); it++) {
38 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
41 //if the user changes the output directory command factory will send this info to us in the output parameter
42 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
44 outputDir += m->hasPath(globaldata->inputFileName); //if user entered a file with a path then preserve it
47 //if the user changes the input directory command factory will send this info to us in the output parameter
48 string inputDir = validParameter.validFile(parameters, "inputdir", false);
49 if (inputDir == "not found"){ inputDir = ""; }
52 it = parameters.find("design");
53 //user has given a template file
54 if(it != parameters.end()){
55 path = m->hasPath(it->second);
56 //if the user has not given a path then, add inputdir. else leave path alone.
57 if (path == "") { parameters["design"] = inputDir + it->second; }
61 //check for required parameters
62 designfile = validParameter.validFile(parameters, "design", true);
63 if (designfile == "not open") { abort = true; }
64 else if (designfile == "not found") { designfile = ""; m->mothurOut("You must provide an design file."); m->mothurOutEndLine(); abort = true; }
66 //make sure the user has already run the read.otu command
67 if ((globaldata->getSharedFile() == "")) {
68 m->mothurOut("You must read a list and a group, or a shared file before you can use the normalize.shared command."); m->mothurOutEndLine(); abort = true;
71 //check for optional parameter and set defaults
72 // ...at some point should added some additional type checking...
73 label = validParameter.validFile(parameters, "label", false);
74 if (label == "not found") { label = ""; }
76 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
77 else { allLines = 1; }
80 //if the user has not specified any labels use the ones from read.otu
82 allLines = globaldata->allLines;
83 labels = globaldata->labels;
86 groups = validParameter.validFile(parameters, "groups", false);
87 if (groups == "not found") { groups = ""; pickedGroups = false; }
90 m->splitAtDash(groups, Groups);
91 globaldata->Groups = Groups;
94 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
97 temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "0.05"; }
98 convert(temp, threshold);
100 //temp = validParameter.validFile(parameters, "g", false); if (temp == "not found") { temp = "0"; }
105 catch(exception& e) {
106 m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
111 //**********************************************************************************************************************
113 void MetaStatsCommand::help(){
115 m->mothurOut("Paper reference goes here....\n");
116 m->mothurOut("The metastats command can only be executed after a successful read.otu command of a list and group or shared file.\n");
117 m->mothurOut("The metastats command outputs a .... file. \n");
118 m->mothurOut("The metastats command parameters are iters, threshold, groups, design and label. No parameters are required.\n");
119 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");
120 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");
121 m->mothurOut("The iters parameter ...\n");
122 m->mothurOut("The threshold parameter ...\n");
123 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");
124 m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n");
125 m->mothurOut("The metastats command should be in the following format: metastats(groups=yourGroups, label=yourLabels).\n");
126 m->mothurOut("Example metastats(groups=A-B-C).\n");
127 m->mothurOut("The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n");
128 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
131 catch(exception& e) {
132 m->errorOut(e, "MetaStatsCommand", "help");
137 //**********************************************************************************************************************
139 MetaStatsCommand::~MetaStatsCommand(){}
141 //**********************************************************************************************************************
143 int MetaStatsCommand::execute(){
146 if (abort == true) { return 0; }
148 designMap = new GroupMap(designfile);
149 designMap->readDesignMap();
151 read = new ReadOTUFile(globaldata->inputFileName);
152 read->read(&*globaldata);
153 input = globaldata->ginput;
154 lookup = input->getSharedRAbundVectors();
155 string lastLabel = lookup[0]->getLabel();
157 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
158 set<string> processedLabels;
159 set<string> userLabels = labels;
161 //as long as you are not at the end of the file or done wih the lines you want
162 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
164 if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } globaldata->Groups.clear(); delete read; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
166 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
168 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
171 processedLabels.insert(lookup[0]->getLabel());
172 userLabels.erase(lookup[0]->getLabel());
175 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
176 string saveLabel = lookup[0]->getLabel();
178 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
179 lookup = input->getSharedRAbundVectors(lastLabel);
180 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
184 processedLabels.insert(lookup[0]->getLabel());
185 userLabels.erase(lookup[0]->getLabel());
187 //restore real lastlabel to save below
188 lookup[0]->setLabel(saveLabel);
191 lastLabel = lookup[0]->getLabel();
192 //prevent memory leak
193 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; }
195 if (m->control_pressed) { globaldata->Groups.clear(); delete read; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
197 //get next line to process
198 lookup = input->getSharedRAbundVectors();
201 if (m->control_pressed) { globaldata->Groups.clear(); delete read; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
203 //output error messages about any remaining user labels
204 set<string>::iterator it;
205 bool needToRun = false;
206 for (it = userLabels.begin(); it != userLabels.end(); it++) {
207 m->mothurOut("Your file does not include the label " + *it);
208 if (processedLabels.count(lastLabel) != 1) {
209 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
212 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
216 //run last label if you need to
217 if (needToRun == true) {
218 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
219 lookup = input->getSharedRAbundVectors(lastLabel);
221 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
225 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
228 //reset groups parameter
229 globaldata->Groups.clear();
230 delete input; globaldata->ginput = NULL;
233 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0;}
235 m->mothurOutEndLine();
236 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
237 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
238 m->mothurOutEndLine();
242 catch(exception& e) {
243 m->errorOut(e, "MetaStatsCommand", "execute");
247 //**********************************************************************************************************************
249 int MetaStatsCommand::process(vector<SharedRAbundVector*>& thisLookUp){
251 if (pickedGroups) { eliminateZeroOTUS(thisLookUp); }
253 string outputFileName = outputDir + m->getRootName(m->getSimpleName(globaldata->inputFileName)) + thisLookUp[0]->getLabel() + ".metastats";
254 outputNames.push_back(outputFileName);
255 int nameLength = outputFileName.length();
256 char * output = new char[nameLength];
257 strcpy(output, outputFileName.c_str());
259 //build matrix from shared rabunds
261 data = new double*[thisLookUp[0]->getNumBins()];
263 for (int j = 0; j < thisLookUp[0]->getNumBins(); j++) {
264 data[j] = new double[thisLookUp.size()];
265 for (int i = 0; i < thisLookUp.size(); i++) {
266 data[j][i] = (thisLookUp[i]->getAbundance(j));
270 if (g == 0) { g = thisLookUp.size() / 2; }
272 metastat_main(output, thisLookUp[0]->getNumBins(), thisLookUp.size(), threshold, iters, data, g);
276 catch(exception& e) {
277 m->errorOut(e, "MetaStatsCommand", "normalize");
281 //**********************************************************************************************************************
282 int MetaStatsCommand::eliminateZeroOTUS(vector<SharedRAbundVector*>& thislookup) {
285 vector<SharedRAbundVector*> newLookup;
286 for (int i = 0; i < thislookup.size(); i++) {
287 SharedRAbundVector* temp = new SharedRAbundVector();
288 temp->setLabel(thislookup[i]->getLabel());
289 temp->setGroup(thislookup[i]->getGroup());
290 newLookup.push_back(temp);
294 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
295 if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) { delete newLookup[j]; } return 0; }
297 //look at each sharedRabund and make sure they are not all zero
299 for (int j = 0; j < thislookup.size(); j++) {
300 if (thislookup[j]->getAbundance(i) != 0) { allZero = false; break; }
303 //if they are not all zero add this bin
305 for (int j = 0; j < thislookup.size(); j++) {
306 newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
311 for (int j = 0; j < thislookup.size(); j++) { delete thislookup[j]; }
313 thislookup = newLookup;
318 catch(exception& e) {
319 m->errorOut(e, "MetaStatsCommand", "eliminateZeroOTUS");
323 //**********************************************************************************************************************