2 * getrelabundcommand.cpp
5 * Created by westcott on 6/21/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "getrelabundcommand.h"
12 //**********************************************************************************************************************
13 vector<string> GetRelAbundCommand::setParameters(){
15 CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pshared);
16 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
17 CommandParameter pscale("scale", "Multiple", "totalgroup-totalotu-averagegroup-averageotu", "totalgroup", "", "", "",false,false); parameters.push_back(pscale);
18 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
19 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
20 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
22 vector<string> myArray;
23 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
27 m->errorOut(e, "GetRelAbundCommand", "setParameters");
31 //**********************************************************************************************************************
32 string GetRelAbundCommand::getHelpString(){
34 string helpString = "";
35 helpString += "The get.relabund command parameters are shared, groups, scale and label. shared is required, unless you have a valid current file.\n";
36 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";
37 helpString += "The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n";
38 helpString += "The scale parameter allows you to select what scale you would like to use. Choices are totalgroup, totalotu, averagegroup, averageotu, default is totalgroup.\n";
39 helpString += "The get.relabund command should be in the following format: get.relabund(groups=yourGroups, label=yourLabels).\n";
40 helpString += "Example get.relabund(groups=A-B-C, scale=averagegroup).\n";
41 helpString += "The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n";
42 helpString += "The get.relabund command outputs a .relabund file.\n";
43 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
47 m->errorOut(e, "GetRelAbundCommand", "getHelpString");
51 //**********************************************************************************************************************
52 string GetRelAbundCommand::getOutputFileNameTag(string type, string inputName=""){
54 string outputFileName = "";
55 map<string, vector<string> >::iterator it;
57 //is this a type this command creates
58 it = outputTypes.find(type);
59 if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
61 if (type == "relabund") { outputFileName = "relabund" ; }
62 else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
64 return outputFileName;
67 m->errorOut(e, "GetRelAbundCommand", "getOutputFileNameTag");
71 //**********************************************************************************************************************
72 GetRelAbundCommand::GetRelAbundCommand(){
74 abort = true; calledHelp = true;
76 vector<string> tempOutNames;
77 outputTypes["relabund"] = tempOutNames;
80 m->errorOut(e, "GetRelAbundCommand", "GetRelAbundCommand");
84 //**********************************************************************************************************************
85 GetRelAbundCommand::GetRelAbundCommand(string option) {
87 abort = false; calledHelp = false;
90 //allow user to run help
91 if(option == "help") { help(); abort = true; calledHelp = true; }
92 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
95 vector<string> myArray = setParameters();
97 OptionParser parser(option);
98 map<string,string> parameters = parser.getParameters();
99 map<string,string>::iterator it;
101 ValidParameters validParameter;
103 //check to make sure all parameters are valid for command
104 for (it = parameters.begin(); it != parameters.end(); it++) {
105 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
108 //initialize outputTypes
109 vector<string> tempOutNames;
110 outputTypes["relabund"] = tempOutNames;
112 //if the user changes the input directory command factory will send this info to us in the output parameter
113 string inputDir = validParameter.validFile(parameters, "inputdir", false);
114 if (inputDir == "not found"){ inputDir = ""; }
117 it = parameters.find("shared");
118 //user has given a template file
119 if(it != parameters.end()){
120 path = m->hasPath(it->second);
121 //if the user has not given a path then, add inputdir. else leave path alone.
122 if (path == "") { parameters["shared"] = inputDir + it->second; }
127 sharedfile = validParameter.validFile(parameters, "shared", true);
128 if (sharedfile == "not open") { sharedfile = ""; abort = true; }
129 else if (sharedfile == "not found") {
130 //if there is a current shared file, use it
131 sharedfile = m->getSharedFile();
132 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
133 else { m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
134 }else { m->setSharedFile(sharedfile); }
137 //if the user changes the output directory command factory will send this info to us in the output parameter
138 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(sharedfile); }
140 //check for optional parameter and set defaults
141 // ...at some point should added some additional type checking...
142 label = validParameter.validFile(parameters, "label", false);
143 if (label == "not found") { label = ""; }
145 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
146 else { allLines = 1; }
149 groups = validParameter.validFile(parameters, "groups", false);
150 if (groups == "not found") { groups = ""; pickedGroups = false; }
153 m->splitAtDash(groups, Groups);
154 m->setGroups(Groups);
157 scale = validParameter.validFile(parameters, "scale", false); if (scale == "not found") { scale = "totalgroup"; }
159 if ((scale != "totalgroup") && (scale != "totalotu") && (scale != "averagegroup") && (scale != "averageotu")) {
160 m->mothurOut(scale + " is not a valid scaling option for the get.relabund command. Choices are totalgroup, totalotu, averagegroup, averageotu."); m->mothurOutEndLine(); abort = true;
165 catch(exception& e) {
166 m->errorOut(e, "GetRelAbundCommand", "GetRelAbundCommand");
170 //**********************************************************************************************************************
172 int GetRelAbundCommand::execute(){
175 if (abort == true) { if (calledHelp) { return 0; } return 2; }
177 string outputFileName = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + getOutputFileNameTag("relabund");
179 m->openOutputFile(outputFileName, out);
180 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
182 input = new InputData(sharedfile, "sharedfile");
183 lookup = input->getSharedRAbundVectors();
184 string lastLabel = lookup[0]->getLabel();
186 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
187 set<string> processedLabels;
188 set<string> userLabels = labels;
190 //as long as you are not at the end of the file or done wih the lines you want
191 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
193 if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } m->clearGroups(); delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
195 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
197 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
198 if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
199 getRelAbundance(lookup, out);
201 processedLabels.insert(lookup[0]->getLabel());
202 userLabels.erase(lookup[0]->getLabel());
205 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
206 string saveLabel = lookup[0]->getLabel();
208 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
209 lookup = input->getSharedRAbundVectors(lastLabel);
210 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
211 if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
212 getRelAbundance(lookup, out);
214 processedLabels.insert(lookup[0]->getLabel());
215 userLabels.erase(lookup[0]->getLabel());
217 //restore real lastlabel to save below
218 lookup[0]->setLabel(saveLabel);
221 lastLabel = lookup[0]->getLabel();
222 //prevent memory leak
223 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; }
225 if (m->control_pressed) { outputTypes.clear(); m->clearGroups(); delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
227 //get next line to process
228 lookup = input->getSharedRAbundVectors();
231 if (m->control_pressed) { outputTypes.clear(); m->clearGroups(); delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
233 //output error messages about any remaining user labels
234 set<string>::iterator it;
235 bool needToRun = false;
236 for (it = userLabels.begin(); it != userLabels.end(); it++) {
237 m->mothurOut("Your file does not include the label " + *it);
238 if (processedLabels.count(lastLabel) != 1) {
239 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
242 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
246 //run last label if you need to
247 if (needToRun == true) {
248 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
249 lookup = input->getSharedRAbundVectors(lastLabel);
251 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
252 if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
253 getRelAbundance(lookup, out);
255 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
258 //reset groups parameter
263 if (m->control_pressed) { outputTypes.clear(); m->mothurRemove(outputFileName); return 0;}
265 m->mothurOutEndLine();
266 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
267 m->mothurOut(outputFileName); m->mothurOutEndLine(); outputNames.push_back(outputFileName); outputTypes["relabund"].push_back(outputFileName);
268 m->mothurOutEndLine();
270 //set relabund file as new current relabundfile
272 itTypes = outputTypes.find("relabund");
273 if (itTypes != outputTypes.end()) {
274 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRelAbundFile(current); }
279 catch(exception& e) {
280 m->errorOut(e, "GetRelAbundCommand", "execute");
284 //**********************************************************************************************************************
286 int GetRelAbundCommand::getRelAbundance(vector<SharedRAbundVector*>& thisLookUp, ofstream& out){
289 for (int i = 0; i < thisLookUp.size(); i++) {
290 out << thisLookUp[i]->getLabel() << '\t' << thisLookUp[i]->getGroup() << '\t' << thisLookUp[i]->getNumBins() << '\t';
292 for (int j = 0; j < thisLookUp[i]->getNumBins(); j++) {
294 if (m->control_pressed) { return 0; }
296 int abund = thisLookUp[i]->getAbundance(j);
298 float relabund = 0.0;
300 if (scale == "totalgroup") {
301 relabund = abund / (float) thisLookUp[i]->getNumSeqs();
302 }else if (scale == "totalotu") {
303 //calc the total in this otu
305 for (int l = 0; l < thisLookUp.size(); l++) { totalOtu += thisLookUp[l]->getAbundance(j); }
306 relabund = abund / (float) totalOtu;
307 }else if (scale == "averagegroup") {
308 relabund = abund / (float) (thisLookUp[i]->getNumSeqs() / (float) thisLookUp[i]->getNumBins());
309 }else if (scale == "averageotu") {
310 //calc the total in this otu
312 for (int l = 0; l < thisLookUp.size(); l++) { totalOtu += thisLookUp[l]->getAbundance(j); }
313 float averageOtu = totalOtu / (float) thisLookUp.size();
315 relabund = abund / (float) averageOtu;
316 }else{ m->mothurOut(scale + " is not a valid scaling option."); m->mothurOutEndLine(); m->control_pressed = true; return 0; }
318 out << relabund << '\t';
325 catch(exception& e) {
326 m->errorOut(e, "GetRelAbundCommand", "getRelAbundance");
330 //**********************************************************************************************************************