]> git.donarmstrong.com Git - mothur.git/blob - otuassociationcommand.cpp
re-wrote co-occurrence command
[mothur.git] / otuassociationcommand.cpp
1 /*
2  *  otuassociationcommand.cpp
3  *  Mothur
4  *
5  *  Created by westcott on 1/19/12.
6  *  Copyright 2012 Schloss Lab. All rights reserved.
7  *
8  */
9
10 #include "otuassociationcommand.h"
11 #include "linearalgebra.h"
12
13 //**********************************************************************************************************************
14 vector<string> OTUAssociationCommand::setParameters(){  
15         try {
16                 CommandParameter pshared("shared", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(pshared);
17                 CommandParameter prelabund("relabund", "InputTypes", "", "", "SharedRelMeta", "SharedRelMeta", "none",false,false); parameters.push_back(prelabund);
18                 CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
19                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
20                 CommandParameter pmethod("method", "Multiple", "pearson-spearman-kendall", "pearson", "", "", "",false,false); parameters.push_back(pmethod);
21                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
22                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
23                 
24                 vector<string> myArray;
25                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
26                 return myArray;
27         }
28         catch(exception& e) {
29                 m->errorOut(e, "OTUAssociationCommand", "setParameters");
30                 exit(1);
31         }
32 }
33 //**********************************************************************************************************************
34 string OTUAssociationCommand::getHelpString(){  
35         try {
36                 string helpString = "";
37                 helpString += "The otu.association command reads a shared or relabund file and calculates the correlation coefficients between otus.\n";
38                 helpString += "The otu.association command parameters are shared, relabund, groups, method and label.  The shared or relabund parameter is required.\n";
39                 helpString += "The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes.\n";
40                 helpString += "The label parameter allows you to select what distances level you would like used, and are also separated by dashes.\n";
41                 helpString += "The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n";
42                 helpString += "The otu.association command should be in the following format: otu.association(shared=yourSharedFile, method=yourMethod).\n";
43                 helpString += "Example otu.association(shared=genus.pool.shared, method=kendall).\n";
44                 helpString += "The otu.association command outputs a .otu.corr file.\n";
45                 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
46                 return helpString;
47         }
48         catch(exception& e) {
49                 m->errorOut(e, "OTUAssociationCommand", "getHelpString");
50                 exit(1);
51         }
52 }
53 //**********************************************************************************************************************
54 OTUAssociationCommand::OTUAssociationCommand(){ 
55         try {
56                 abort = true; calledHelp = true; 
57                 setParameters();
58                 vector<string> tempOutNames;
59                 outputTypes["otu.corr"] = tempOutNames;
60         }
61         catch(exception& e) {
62                 m->errorOut(e, "OTUAssociationCommand", "OTUAssociationCommand");
63                 exit(1);
64         }
65 }
66 //**********************************************************************************************************************
67 OTUAssociationCommand::OTUAssociationCommand(string option)  {
68         try {
69                 abort = false; calledHelp = false;   
70                 allLines = 1;
71                 
72                 //allow user to run help
73                 if(option == "help") { help(); abort = true; calledHelp = true; }
74                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
75                 
76                 else {
77                         vector<string> myArray = setParameters();
78                         
79                         OptionParser parser(option);
80                         map<string, string> parameters = parser.getParameters();
81                         
82                         ValidParameters validParameter;
83                         map<string, string>::iterator it;
84                         
85                         //check to make sure all parameters are valid for command
86                         for (it = parameters.begin(); it != parameters.end(); it++) { 
87                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
88                         }
89                         
90                         vector<string> tempOutNames;
91                         outputTypes["otu.corr"] = tempOutNames;
92                         
93                         //if the user changes the input directory command factory will send this info to us in the output parameter 
94                         string inputDir = validParameter.validFile(parameters, "inputdir", false);              
95                         if (inputDir == "not found"){   inputDir = "";          }
96                         else {
97                                 string path;
98                                 it = parameters.find("shared");
99                                 //user has given a template file
100                                 if(it != parameters.end()){ 
101                                         path = m->hasPath(it->second);
102                                         //if the user has not given a path then, add inputdir. else leave path alone.
103                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
104                                 }
105                                 
106                                 it = parameters.find("relabund");
107                                 //user has given a template file
108                                 if(it != parameters.end()){ 
109                                         path = m->hasPath(it->second);
110                                         //if the user has not given a path then, add inputdir. else leave path alone.
111                                         if (path == "") {       parameters["relabund"] = inputDir + it->second;         }
112                                 }
113                         }
114                         
115                         
116                         //check for required parameters                 
117                         sharedfile = validParameter.validFile(parameters, "shared", true);
118                         if (sharedfile == "not open") { abort = true; }
119                         else if (sharedfile == "not found") { sharedfile = ""; }
120                         else { inputFileName = sharedfile; m->setSharedFile(sharedfile); }
121                         
122                         relabundfile = validParameter.validFile(parameters, "relabund", true);
123                         if (relabundfile == "not open") { abort = true; }
124                         else if (relabundfile == "not found") { relabundfile = ""; }
125                         else { inputFileName = relabundfile; m->setRelAbundFile(relabundfile); }
126                         
127                         groups = validParameter.validFile(parameters, "groups", false);                 
128                         if (groups == "not found") { groups = "";  pickedGroups = false;  }
129                         else { 
130                                 pickedGroups = true;
131                                 m->splitAtDash(groups, Groups); 
132                         }                       
133                         m->setGroups(Groups);
134                         
135                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(inputFileName);  }
136                         
137                         label = validParameter.validFile(parameters, "label", false);                   
138                         if (label == "not found") { label = ""; }
139                         else { 
140                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
141                                 else { allLines = 1;  }
142                         }
143                         
144                         if ((relabundfile == "") && (sharedfile == "")) { 
145                                 //is there are current file available for any of these?
146                                 //give priority to shared, then relabund
147                                 //if there is a current shared file, use it
148                                 sharedfile = m->getSharedFile(); 
149                                 if (sharedfile != "") { inputFileName = sharedfile; m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
150                                 else { 
151                                         relabundfile = m->getRelAbundFile(); 
152                                         if (relabundfile != "") { inputFileName = relabundfile;  m->mothurOut("Using " + relabundfile + " as input file for the relabund parameter."); m->mothurOutEndLine(); }
153                                         else { 
154                                                 m->mothurOut("You must provide either a shared or relabund file."); m->mothurOutEndLine(); abort = true; 
155                                         }
156                                 }
157                         }       
158                         
159                         
160                         if ((relabundfile != "") && (sharedfile != "")) { m->mothurOut("You may only use one of the following : shared or relabund file."); m->mothurOutEndLine(); abort = true;  }
161                         
162                         method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "pearson";             }
163                         
164                         if ((method != "pearson") && (method != "spearman") && (method != "kendall")) { m->mothurOut(method + " is not a valid method. Valid methods are pearson, spearman, and kendall."); m->mothurOutEndLine(); abort = true; }
165                         
166                 }
167         }
168         catch(exception& e) {
169                 m->errorOut(e, "OTUAssociationCommand", "OTUAssociationCommand");               
170                 exit(1);
171         }
172 }
173 //**********************************************************************************************************************
174
175 int OTUAssociationCommand::execute(){
176         try {
177                 
178                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
179                 
180                 //function are identical just different datatypes
181                 if (sharedfile != "")                   {  processShared();             } 
182                 else if (relabundfile != "")    {  processRelabund();   }
183                                 
184                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); } return 0; }
185                 
186                 m->mothurOutEndLine();
187                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
188                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
189                 m->mothurOutEndLine();
190                 
191                 return 0;
192         }
193         catch(exception& e) {
194                 m->errorOut(e, "OTUAssociationCommand", "execute");     
195                 exit(1);
196         }
197 }
198 //**********************************************************************************************************************
199 int OTUAssociationCommand::processShared(){
200         try {
201                 InputData* input = new InputData(sharedfile, "sharedfile");
202                 vector<SharedRAbundVector*> lookup = input->getSharedRAbundVectors();
203                 string lastLabel = lookup[0]->getLabel();
204                 
205                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
206                 set<string> processedLabels;
207                 set<string> userLabels = labels;
208                 
209                 //as long as you are not at the end of the file or done wih the lines you want
210                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
211                         
212                         if (m->control_pressed) {  delete input; return 0;  }
213                         
214                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){  
215                                 processedLabels.insert(lookup[0]->getLabel());
216                                 userLabels.erase(lookup[0]->getLabel());
217                                 
218                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
219                                 process(lookup);
220                         }
221                         
222                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
223                                 string saveLabel = lookup[0]->getLabel();
224                                 
225                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
226                                 lookup = input->getSharedRAbundVectors(lastLabel);
227                                 
228                                 processedLabels.insert(lookup[0]->getLabel());
229                                 userLabels.erase(lookup[0]->getLabel());
230                                 
231                                 //restore real lastlabel to save below
232                                 lookup[0]->setLabel(saveLabel);
233                                 
234                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
235                                 process(lookup);
236                         }
237                         
238                         lastLabel = lookup[0]->getLabel();                      
239                         
240                         //get next line to process
241                         //prevent memory leak
242                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
243                         lookup = input->getSharedRAbundVectors();
244                 }
245                 
246                 
247                 if (m->control_pressed) { delete input; return 0;  }
248                 
249                 //output error messages about any remaining user labels
250                 set<string>::iterator it;
251                 bool needToRun = false;
252                 for (it = userLabels.begin(); it != userLabels.end(); it++) {  
253                         m->mothurOut("Your file does not include the label " + *it); 
254                         if (processedLabels.count(lastLabel) != 1) {
255                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
256                                 needToRun = true;
257                         }else {
258                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
259                         }
260                 }
261                 
262                 //run last label if you need to
263                 if (needToRun == true)  {
264                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } } 
265                         lookup = input->getSharedRAbundVectors(lastLabel);
266                         
267                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
268                         process(lookup);
269                 }       
270                 
271                 delete input;
272                 
273                 return 0;
274         }
275         catch(exception& e) {
276                 m->errorOut(e, "OTUAssociationCommand", "processShared");       
277                 exit(1);
278         }
279 }
280 //**********************************************************************************************************************
281 int OTUAssociationCommand::process(vector<SharedRAbundVector*>& lookup){
282         try {
283                 
284                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + "." + method + ".otu.corr";
285                 outputNames.push_back(outputFileName); outputTypes["shared"].push_back(outputFileName);
286                 
287                 ofstream out;
288                 m->openOutputFile(outputFileName, out);
289                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
290                 
291                 //column headings
292                 out << "OTUA\tOTUB\t" << method << "Coef\tSignificance\n";
293                 
294                 vector< vector<double> > xy; xy.resize(lookup[0]->getNumBins());
295                 for (int i = 0; i < lookup[0]->getNumBins(); i++) { for (int j = 0; j < lookup.size(); j++) { xy[i].push_back(lookup[j]->getAbundance(i)); } }
296                 
297                 LinearAlgebra linear;
298                 for (int i = 0; i < xy.size(); i++) {
299                         
300                         for (int k = 0; k < i; k++) {
301                                 
302                                 if (m->control_pressed) { out.close(); return 0; }
303
304                                 double coef = 0.0;
305                                 double sig = 0.0;
306                                 if (method == "spearman")               {   coef = linear.calcSpearman(xy[i], xy[k], sig);      }
307                                 else if (method == "pearson")   {       coef = linear.calcPearson(xy[i], xy[k], sig);   }
308                                 else if (method == "kendall")   {       coef = linear.calcKendall(xy[i], xy[k], sig);   }                   
309                                 else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
310                         
311                 out << m->binLabelsInFile[i] << '\t' << m->binLabelsInFile[k] << '\t' << coef << '\t' << sig << endl;
312                         }
313                 }
314                 
315                 out.close();
316                 
317                
318                 return 0;
319                 
320         }
321         catch(exception& e) {
322                 m->errorOut(e, "OTUAssociationCommand", "process");     
323                 exit(1);
324         }
325 }
326 //**********************************************************************************************************************
327 int OTUAssociationCommand::processRelabund(){
328         try {
329                 InputData* input = new InputData(relabundfile, "relabund");
330                 vector<SharedRAbundFloatVector*> lookup = input->getSharedRAbundFloatVectors();
331                 string lastLabel = lookup[0]->getLabel();
332                 
333                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
334                 set<string> processedLabels;
335                 set<string> userLabels = labels;
336                 
337                 //as long as you are not at the end of the file or done wih the lines you want
338                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
339                         
340                         if (m->control_pressed) {  delete input; return 0;  }
341                         
342                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){  
343                                 processedLabels.insert(lookup[0]->getLabel());
344                                 userLabels.erase(lookup[0]->getLabel());
345                                 
346                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
347                                 process(lookup);
348                         }
349                         
350                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
351                                 string saveLabel = lookup[0]->getLabel();
352                                 
353                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
354                                 lookup = input->getSharedRAbundFloatVectors(lastLabel);
355                                 
356                                 processedLabels.insert(lookup[0]->getLabel());
357                                 userLabels.erase(lookup[0]->getLabel());
358                                 
359                                 //restore real lastlabel to save below
360                                 lookup[0]->setLabel(saveLabel);
361                                 
362                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
363                                 process(lookup);
364                         }
365                         
366                         lastLabel = lookup[0]->getLabel();                      
367                         
368                         //get next line to process
369                         //prevent memory leak
370                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
371                         lookup = input->getSharedRAbundFloatVectors();
372                 }
373                 
374                 
375                 if (m->control_pressed) { delete input; return 0;  }
376                 
377                 //output error messages about any remaining user labels
378                 set<string>::iterator it;
379                 bool needToRun = false;
380                 for (it = userLabels.begin(); it != userLabels.end(); it++) {  
381                         m->mothurOut("Your file does not include the label " + *it); 
382                         if (processedLabels.count(lastLabel) != 1) {
383                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
384                                 needToRun = true;
385                         }else {
386                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
387                         }
388                 }
389                 
390                 //run last label if you need to
391                 if (needToRun == true)  {
392                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       } } 
393                         lookup = input->getSharedRAbundFloatVectors(lastLabel);
394                         
395                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
396                         process(lookup);
397                 }       
398                 
399                 delete input;
400                 
401                 return 0;
402         }
403         catch(exception& e) {
404                 m->errorOut(e, "OTUAssociationCommand", "processRelabund");     
405                 exit(1);
406         }
407 }
408 //**********************************************************************************************************************
409 int OTUAssociationCommand::process(vector<SharedRAbundFloatVector*>& lookup){
410         try {
411                 
412                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(inputFileName)) + lookup[0]->getLabel() + "." + method + ".otu.corr";
413                 outputNames.push_back(outputFileName); outputTypes["shared"].push_back(outputFileName);
414                 
415                 ofstream out;
416                 m->openOutputFile(outputFileName, out);
417                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
418                 
419                 //column headings
420                 out << "OTUA\tOTUB\t" << method << "Coef\tSignificance\n";
421                 
422                 vector< vector<double> > xy; xy.resize(lookup[0]->getNumBins());
423                 for (int i = 0; i < lookup[0]->getNumBins(); i++) { for (int j = 0; j < lookup.size(); j++) { xy[i].push_back(lookup[j]->getAbundance(i)); } }
424                 
425                 LinearAlgebra linear;
426                 for (int i = 0; i < xy.size(); i++) {
427                         
428                         for (int k = 0; k < i; k++) {
429                                 
430                                 if (m->control_pressed) { out.close(); return 0; }
431                                 
432                                 double coef = 0.0;
433                                 double sig = 0.0;
434                                 if (method == "spearman")               {   coef = linear.calcSpearman(xy[i], xy[k], sig);      }
435                                 else if (method == "pearson")   {       coef = linear.calcPearson(xy[i], xy[k], sig);   }
436                                 else if (method == "kendall")   {       coef = linear.calcKendall(xy[i], xy[k], sig);   }                   
437                                 else { m->mothurOut("[ERROR]: invalid method, choices are spearman, pearson or kendall."); m->mothurOutEndLine(); m->control_pressed = true; }
438                                 
439                 out << m->binLabelsInFile[i] << '\t' << m->binLabelsInFile[k] << '\t' << coef << '\t' << sig << endl; 
440                         }
441                 }
442                 
443                 out.close();
444                 
445                 return 0;
446                 
447         }
448         catch(exception& e) {
449                 m->errorOut(e, "OTUAssociationCommand", "process");     
450                 exit(1);
451         }
452 }
453 /*****************************************************************/
454
455
456
457
458
459
460
461