2 * cooccurrencecommand.cpp
5 * Created by kiverson on 1/2/12.
6 * Copyright 2012 Schloss Lab. All rights reserved.
10 #include "cooccurrencecommand.h"
12 //**********************************************************************************************************************
13 vector<string> CooccurrenceCommand::setParameters() {
15 CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pshared);
16 CommandParameter pmetric("metric", "Multiple", "cscore-checker-combo-vratio", "cscore", "", "", "",false,false); parameters.push_back(pmetric);
17 CommandParameter pmatrix("matrixmodel", "Multiple", "sim1-sim2-sim3-sim4-sim5-sim6-sim7-sim8-sim9", "sim2", "", "", "",false,false); parameters.push_back(pmatrix);
18 CommandParameter pruns("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(pruns);
19 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
20 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
21 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
22 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
24 vector<string> myArray;
25 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
29 m->errorOut(e, "CooccurrenceCommand", "setParameters");
33 //**********************************************************************************************************************
34 string CooccurrenceCommand::getHelpString(){
36 string helpString = "The cooccurrence command calculates four metrics and tests their significance to assess whether presence-absence patterns are different than what one would expect by chance.";
37 helpString += "The cooccurrence command parameters are shared, metric, matrixmodel, iters, label and groups.";
38 helpString += "The matrixmodel parameter options are sim1, sim2, sim3, sim4, sim5, sim6, sim7, sim8 and sim9. Default=sim2";
39 helpString += "The metric parameter options are cscore, checker, combo and vratio. Default=cscore";
40 helpString += "The label parameter is used to analyze specific labels in your input.\n";
41 helpString += "The groups parameter allows you to specify which of the groups you would like analyzed.\n";
42 helpString += "The cooccurrence command should be in the following format: \n";
43 helpString += "cooccurrence(shared=yourSharedFile) \n";
44 helpString += "Example cooccurrence(shared=final.an.shared).\n";
45 helpString += "Note: No spaces between parameter labels (i.e. shared), '=' and parameters (i.e.yourShared).\n";
49 m->errorOut(e, "CooccurrenceCommand", "getHelpString");
53 //**********************************************************************************************************************
54 CooccurrenceCommand::CooccurrenceCommand(){
56 abort = true; calledHelp = true;
58 vector<string> tempOutNames;
59 outputTypes["summary"] = tempOutNames;
63 m->errorOut(e, "CooccurrenceCommand", "CooccurrenceCommand");
67 //**********************************************************************************************************************
68 CooccurrenceCommand::CooccurrenceCommand(string option) {
70 abort = false; calledHelp = false;
73 //allow user to run help
74 if(option == "help") { help(); abort = true; calledHelp = true; }
75 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
78 vector<string> myArray = setParameters();
80 OptionParser parser(option);
81 map<string,string> parameters = parser.getParameters();
82 map<string,string>::iterator it;
84 ValidParameters validParameter;
86 //check to make sure all parameters are valid for command
87 for (it = parameters.begin(); it != parameters.end(); it++) {
88 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
92 //if the user changes the input directory command factory will send this info to us in the output parameter
93 string inputDir = validParameter.validFile(parameters, "inputdir", false);
94 if (inputDir == "not found"){ inputDir = ""; }
97 it = parameters.find("shared");
98 //user has given a template file
99 if(it != parameters.end()){
100 path = m->hasPath(it->second);
101 //if the user has not given a path then, add inputdir. else leave path alone.
102 if (path == "") { parameters["shared"] = inputDir + it->second; }
106 vector<string> tempOutNames;
107 outputTypes["summary"] = tempOutNames;
109 //check for optional parameter and set defaults
110 // ...at some point should added some additional type checking...
111 label = validParameter.validFile(parameters, "label", false);
112 if (label == "not found") { label = ""; }
114 if(label != "all") { m->splitAtDash(label, labels); allLines = 0; }
115 else { allLines = 1; }
119 sharedfile = validParameter.validFile(parameters, "shared", true);
120 if (sharedfile == "not open") { sharedfile = ""; abort = true; }
121 else if (sharedfile == "not found") {
122 //if there is a current shared file, use it
123 sharedfile = m->getSharedFile();
124 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
125 else { m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
126 }else { m->setSharedFile(sharedfile); }
129 //if the user changes the output directory command factory will send this info to us in the output parameter
130 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(sharedfile); }
133 metric = validParameter.validFile(parameters, "metric", false); if (metric == "not found") { metric = "cscore"; }
135 if ((metric != "cscore") && (metric != "checker") && (metric != "combo") && (metric != "vratio")) {
136 m->mothurOut("[ERROR]: " + metric + " is not a valid metric option for the cooccurrence command. Choices are cscore, checker, combo, vratio."); m->mothurOutEndLine(); abort = true;
139 matrix = validParameter.validFile(parameters, "matrixmodel", false); if (matrix == "not found") { matrix = "sim2"; }
141 if ((matrix != "sim1") && (matrix != "sim2") && (matrix != "sim3") && (matrix != "sim4") && (matrix != "sim5" ) && (matrix != "sim6" ) && (matrix != "sim7" ) && (matrix != "sim8" ) && (matrix != "sim9" )) {
142 m->mothurOut("[ERROR]: " + matrix + " is not a valid matrix option for the cooccurrence command. Choices are sim1, sim2, sim3, sim4, sim5, sim6, sim7, sim8, sim9."); m->mothurOutEndLine(); abort = true;
145 groups = validParameter.validFile(parameters, "groups", false);
146 if (groups == "not found") { groups = ""; }
148 m->splitAtDash(groups, Groups);
150 m->setGroups(Groups);
152 string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
153 m->mothurConvert(temp, runs);
158 catch(exception& e) {
159 m->errorOut(e, "CooccurrenceCommand", "CooccurrenceCommand");
163 //**********************************************************************************************************************
165 int CooccurrenceCommand::execute(){
168 if (abort == true) { if (calledHelp) { return 0; } return 2; }
170 InputData* input = new InputData(sharedfile, "sharedfile");
171 vector<SharedRAbundVector*> lookup = input->getSharedRAbundVectors();
172 string lastLabel = lookup[0]->getLabel();
174 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
175 set<string> processedLabels;
176 set<string> userLabels = labels;
179 string outputFileName = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + "cooccurence.summary";
180 m->openOutputFile(outputFileName, out);
181 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
182 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
183 out << "metric\tlabel\tScore\tpValue\n";
185 //as long as you are not at the end of the file or done wih the lines you want
186 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
188 if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
190 if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
192 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
194 getCooccurrence(lookup, out);
196 processedLabels.insert(lookup[0]->getLabel());
197 userLabels.erase(lookup[0]->getLabel());
200 if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
201 string saveLabel = lookup[0]->getLabel();
203 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
204 lookup = input->getSharedRAbundVectors(lastLabel);
205 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
206 getCooccurrence(lookup, out);
208 processedLabels.insert(lookup[0]->getLabel());
209 userLabels.erase(lookup[0]->getLabel());
211 //restore real lastlabel to save below
212 lookup[0]->setLabel(saveLabel);
215 lastLabel = lookup[0]->getLabel();
216 //prevent memory leak
217 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; }
219 if (m->control_pressed) { outputTypes.clear(); delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
221 //get next line to process
222 lookup = input->getSharedRAbundVectors();
225 if (m->control_pressed) { delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
227 //output error messages about any remaining user labels
228 set<string>::iterator it;
229 bool needToRun = false;
230 for (it = userLabels.begin(); it != userLabels.end(); it++) {
231 m->mothurOut("Your file does not include the label " + *it);
232 if (processedLabels.count(lastLabel) != 1) {
233 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
236 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
240 //run last label if you need to
241 if (needToRun == true) {
242 for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
243 lookup = input->getSharedRAbundVectors(lastLabel);
245 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
247 getCooccurrence(lookup, out);
249 for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; }
254 //reset groups parameter
258 m->mothurOutEndLine();
259 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
260 m->mothurOut(outputFileName); m->mothurOutEndLine();
261 m->mothurOutEndLine();
265 catch(exception& e) {
266 m->errorOut(e, "CooccurrenceCommand", "execute");
270 //**********************************************************************************************************************
272 int CooccurrenceCommand::getCooccurrence(vector<SharedRAbundVector*>& thisLookUp, ofstream& out){
274 int numOTUS = thisLookUp[0]->getNumBins();
275 vector< vector<int> > initmatrix; initmatrix.resize(thisLookUp.size());
276 vector< vector<int> > co_matrix; co_matrix.resize(thisLookUp[0]->getNumBins());
277 for (int i = 0; i < thisLookUp[0]->getNumBins(); i++) { co_matrix[i].resize((thisLookUp.size()), 0); }
278 for (int i = 0; i < thisLookUp.size(); i++) { initmatrix[i].resize((thisLookUp[i]->getNumBins()), 0); }
279 vector<int> columntotal; columntotal.resize(thisLookUp.size(), 0);
280 vector<int> rowtotal; rowtotal.resize(numOTUS, 0);
282 for (int i = 0; i < thisLookUp.size(); i++) {
283 for (int j = 0; j < thisLookUp[i]->getNumBins(); j++) {
284 if (m->control_pressed) { return 0; }
285 int abund = thisLookUp[i]->getAbundance(j);
288 initmatrix[i][j] = 1;
296 //nrows is ncols of inital matrix. All the functions need this value. They assume the transposition has already taken place and nrows and ncols refer to that matrix.
297 //comatrix and initmatrix are still vectors of vectors of ints as in the original script. The abundancevector is only what was read in ie not a co-occurrence matrix!
298 int nrows = numOTUS;//rows of inital matrix
299 int ncols = thisLookUp.size();//groups
300 double initscore = 0.0;
302 vector<int> initcolumntotal; initcolumntotal.resize(ncols, 0);
303 vector<int> initrowtotal; initrowtotal.resize(nrows, 0);
304 vector<double> stats;
308 initcolumntotal = rowtotal;
309 initrowtotal = columntotal;
310 trial.update_row_col_totals(co_matrix, rowtotal, columntotal);
312 if (metric == "cscore") { initscore = trial.calc_c_score(co_matrix, rowtotal); }
313 else if (metric == "checker") { initscore = trial.calc_checker(co_matrix, rowtotal); }
314 else if (metric == "vratio") { initscore = trial.calc_vratio(rowtotal, columntotal); }
315 else if (metric == "combo") { initscore = trial.calc_combo(co_matrix); }
316 else { m->mothurOut("[ERROR]: No metric selected!\n"); m->control_pressed = true; return 1; }
318 m->mothurOut("Initial c score: " + toString(initscore)); m->mothurOutEndLine();
321 for(int i=0;i<10000;i++) {
322 if (m->control_pressed) { return 0; }
323 if (matrix == "sim1") {
324 trial.sim1(co_matrix);
325 }else if (matrix == "sim2") {
326 trial.sim2(co_matrix);
327 }else if (matrix == "sim3") {
328 trial.sim3(initmatrix);
329 co_matrix = initmatrix;
330 }else if (matrix == "sim4") {
331 trial.sim4(columntotal, rowtotal, co_matrix);
332 }else if (matrix == "sim5") {
333 trial.sim5(initcolumntotal, initrowtotal, initmatrix);
334 trial.transpose_matrix(initmatrix,co_matrix);
335 }else if (matrix == "sim6") {
336 trial.sim6(columntotal, co_matrix);
337 }else if (matrix == "sim7") {
338 trial.sim7(initcolumntotal, initmatrix);
339 co_matrix = initmatrix;
340 }else if (matrix == "sim8") {
341 trial.sim8(columntotal, rowtotal, co_matrix);
342 }else if (matrix == "sim9") {
343 trial.swap_checkerboards (co_matrix);
345 m->mothurOut("[ERROR]: No model selected! \n");
346 m->control_pressed = true;
351 for(int i=0;i<runs;i++) {
352 if (m->control_pressed) { return 0; }
353 //calc metric of nullmatrix
354 if (matrix == "sim1") {
355 trial.sim1(co_matrix);
356 }else if (matrix == "sim2") {
357 trial.sim2(co_matrix);
358 }else if (matrix == "sim3") {
359 trial.sim3(initmatrix);
360 co_matrix = initmatrix;
361 }else if (matrix == "sim4") {
362 trial.sim4(columntotal, rowtotal, co_matrix);
363 }else if (matrix == "sim5") {
364 trial.sim5(initcolumntotal, initrowtotal, initmatrix);
365 trial.transpose_matrix(initmatrix,co_matrix);
366 }else if (matrix == "sim6") {
367 trial.sim6(columntotal, co_matrix);
368 }else if (matrix == "sim7") {
369 trial.sim7(initcolumntotal, initmatrix);
370 co_matrix = initmatrix;
371 }else if (matrix == "sim8") {
372 trial.sim8(columntotal, rowtotal, co_matrix);
373 }else if (matrix == "sim9") {
374 trial.swap_checkerboards (co_matrix);
376 m->mothurOut("[ERROR]: No model selected! \n");
377 m->control_pressed = true;
381 trial.update_row_col_totals(co_matrix, rowtotal, columntotal);
383 if (metric == "cscore") {
384 stats.push_back(trial.calc_c_score(co_matrix, rowtotal));
385 }else if (metric == "checker") {
386 stats.push_back(trial.calc_checker(co_matrix, rowtotal));
387 }else if (metric == "vratio") {
388 stats.push_back(trial.calc_vratio(rowtotal, columntotal));
389 }else if (metric == "combo") {
390 stats.push_back(trial.calc_combo(co_matrix));
392 m->mothurOut("[ERROR]: No metric selected!\n");
393 m->control_pressed = true;
400 for (int i=0; i<stats.size();i++) { total+=stats[i]; }
402 double nullMean = double (total/(double)stats.size());
404 m->mothurOutEndLine(); m->mothurOut("average metric score: " + toString(nullMean)); m->mothurOutEndLine();
407 if (metric == "cscore" || metric == "checker") { pvalue = trial.calc_pvalue_greaterthan (stats, initscore); }
408 else{ pvalue = trial.calc_pvalue_lessthan (stats, initscore); }
410 m->mothurOut("pvalue: " + toString(pvalue)); m->mothurOutEndLine();
411 out << metric << '\t' << thisLookUp[0]->getLabel() << '\t' << nullMean << '\t' << pvalue << endl;
415 catch(exception& e) {
416 m->errorOut(e, "CooccurrenceCommand", "Cooccurrence");
420 //**********************************************************************************************************************