]> git.donarmstrong.com Git - mothur.git/blob - cooccurrencecommand.cpp
replaced p-value with zscore for lookup in ztable.
[mothur.git] / cooccurrencecommand.cpp
1 /*
2  *  cooccurrencecommand.cpp
3  *  Mothur
4  *
5  *  Created by kiverson on 1/2/12.
6  *  Copyright 2012 Schloss Lab. All rights reserved.
7  *
8  */
9
10 #include "cooccurrencecommand.h"
11
12 //**********************************************************************************************************************
13 vector<string> CooccurrenceCommand::setParameters() {   
14         try { 
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);
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, "CooccurrenceCommand", "setParameters");
30                 exit(1);
31         }
32 }
33 //**********************************************************************************************************************
34 string CooccurrenceCommand::getHelpString(){    
35         try {
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";
46                 return helpString;
47         }
48         catch(exception& e) {
49                 m->errorOut(e, "CooccurrenceCommand", "getHelpString");
50                 exit(1);
51         }
52 }
53 //**********************************************************************************************************************
54 CooccurrenceCommand::CooccurrenceCommand(){     
55         try {
56                 abort = true; calledHelp = true; 
57                 setParameters();
58         vector<string> tempOutNames;
59                 outputTypes["summary"] = tempOutNames;
60
61         }
62         catch(exception& e) {
63                 m->errorOut(e, "CooccurrenceCommand", "CooccurrenceCommand");
64                 exit(1);
65         }
66 }
67 //**********************************************************************************************************************
68 CooccurrenceCommand::CooccurrenceCommand(string option) {
69         try {
70                 abort = false; calledHelp = false;   
71                 allLines = 1;
72                                 
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;}
76                 
77                 else {
78                         vector<string> myArray = setParameters();
79                         
80                         OptionParser parser(option);
81                         map<string,string> parameters = parser.getParameters();
82                         map<string,string>::iterator it;
83                         
84                         ValidParameters validParameter;
85                         
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;  }
89                         }
90
91                         
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 = "";          }
95                         else {
96                                 string path;
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;           }
103                                 }
104                         }
105                 
106             vector<string> tempOutNames;
107             outputTypes["summary"] = tempOutNames;
108                 
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 = ""; }
113                         else { 
114                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
115                                 else { allLines = 1;  }
116                         }
117                         
118                         //get shared file
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); }
127                         
128                         
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);             }
131
132                         
133                         metric = validParameter.validFile(parameters, "metric", false);                         if (metric == "not found") { metric = "cscore"; }
134                         
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; 
137                         }
138                         
139                         matrix = validParameter.validFile(parameters, "matrixmodel", false);                            if (matrix == "not found") { matrix = "sim2"; }
140                         
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; 
143                         }
144             
145             groups = validParameter.validFile(parameters, "groups", false);                     
146                         if (groups == "not found") { groups = "";   }
147                         else { 
148                                 m->splitAtDash(groups, Groups); 
149                         }                       
150                         m->setGroups(Groups);
151             
152             string temp = validParameter.validFile(parameters, "iters", false);                 if (temp == "not found") { temp = "1000"; }
153                         m->mothurConvert(temp, runs); 
154
155                 }
156
157         }
158         catch(exception& e) {
159                 m->errorOut(e, "CooccurrenceCommand", "CooccurrenceCommand");
160                 exit(1);
161         }
162 }
163 //**********************************************************************************************************************
164
165 int CooccurrenceCommand::execute(){
166         try {
167         
168                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
169                 
170                 InputData* input = new InputData(sharedfile, "sharedfile");
171                 vector<SharedRAbundVector*> lookup = input->getSharedRAbundVectors();
172                 string lastLabel = lookup[0]->getLabel();
173                 
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;
177
178         ofstream out;
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\tzScore\n";
184
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))) {
187                         
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; }
189         
190                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){                  
191
192                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
193                                 
194                                 getCooccurrence(lookup, out);
195                                 
196                                 processedLabels.insert(lookup[0]->getLabel());
197                                 userLabels.erase(lookup[0]->getLabel());
198                         }
199                         
200                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
201                                 string saveLabel = lookup[0]->getLabel();
202                         
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);
207                                 
208                                 processedLabels.insert(lookup[0]->getLabel());
209                                 userLabels.erase(lookup[0]->getLabel());
210                                 
211                                 //restore real lastlabel to save below
212                                 lookup[0]->setLabel(saveLabel);
213                         }
214                         
215                         lastLabel = lookup[0]->getLabel();
216                         //prevent memory leak
217                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
218                         
219                         if (m->control_pressed) {  outputTypes.clear(); delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
220
221                         //get next line to process
222                         lookup = input->getSharedRAbundVectors();                               
223                 }
224                 
225                 if (m->control_pressed) { delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
226
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();
234                                 needToRun = true;
235                         }else {
236                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
237                         }
238                 }
239         
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);
244                         
245                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
246                         
247                         getCooccurrence(lookup, out);
248                         
249                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
250                 }
251         
252         out.close(); 
253         
254                 //reset groups parameter 
255                 delete input; 
256         m->clearGroups(); 
257
258         m->mothurOutEndLine();
259                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
260                 m->mothurOut(outputFileName); m->mothurOutEndLine();    
261                 m->mothurOutEndLine();
262         
263                 return 0;
264         }
265         catch(exception& e) {
266                 m->errorOut(e, "CooccurrenceCommand", "execute");
267                 exit(1);
268         }
269 }
270 //**********************************************************************************************************************
271
272 int CooccurrenceCommand::getCooccurrence(vector<SharedRAbundVector*>& thisLookUp, ofstream& out){
273     try {
274         int numOTUS = thisLookUp[0]->getNumBins();
275         vector< vector<int> > co_matrix; co_matrix.resize(thisLookUp[0]->getNumBins());
276         for (int i = 0; i < thisLookUp[0]->getNumBins(); i++) { co_matrix[i].resize((thisLookUp.size()), 0); }
277         vector<int> columntotal; columntotal.resize(thisLookUp.size(), 0);
278         vector<int> rowtotal; rowtotal.resize(numOTUS, 0);
279         
280         for (int i = 0; i < thisLookUp.size(); i++) { //nrows in the shared file
281             for (int j = 0; j < thisLookUp[i]->getNumBins(); j++) { //cols of original shared file
282                 if (m->control_pressed) { return 0; }
283                 int abund = thisLookUp[i]->getAbundance(j);
284                 
285                 if(abund > 0) {
286                     co_matrix[j][i] = 1;
287                     rowtotal[j]++;
288                     columntotal[i]++;
289                 }
290             }
291         }
292         
293         //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.
294         //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!
295         int nrows = numOTUS;//rows of inital matrix
296         int ncols = thisLookUp.size();//groups
297         double initscore = 0.0;
298         
299         vector<double> stats;
300         double probabilityMatrix[ncols * nrows];
301         vector<vector<int> > nullmatrix(nrows, vector<int>(ncols, 0));
302         
303         TrialSwap2 trial;
304         
305         int n = accumulate( columntotal.begin(), columntotal.end(), 0 );
306         
307         //============================================================
308         
309         //generate a probability matrix. Only do this once.
310         float start = 0.0;
311         
312         if (matrix == "sim1") {
313             for(int i=0;i<nrows;i++) {
314                 for(int j=0;j<ncols;j++) {
315                     probabilityMatrix[ncols * i + j] = start + 1/double(nrows*ncols);
316                     start = start + 1/double(nrows*ncols);
317                 }
318             }
319         }
320         //don't need a prob matrix because we just shuffle the rows, may use this in the future
321         else if (matrix == "sim2") { }
322 //            for(int i=0;i<nrows;i++) {
323 //                start = 0.0;
324 //                for(int j=0;j<ncols;j++) {
325 //                    probabilityMatrix[ncols * i + j] = start + 1/double(ncols);
326 //                    start = start + 1/double(ncols);
327 //                }
328 //            }
329 //        }
330         
331         else if (matrix == "sim3") {
332             for(int j=0;j<ncols;j++) {
333                 start = 0.0;
334                 for(int i=0;i<nrows;i++) {
335                     probabilityMatrix[ncols * i + j] = start + 1/double(nrows);
336                     start = start + 1/double(nrows);
337                 }
338             }
339         }
340         
341         else if (matrix == "sim4") {
342             for(int i=0;i<nrows;i++) {
343                 start = 0.0;
344                 for(int j=0;j<ncols;j++) {
345                     probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n);
346                     start = start + columntotal[j]/double(n);
347                 }
348             }
349         }
350         
351         else if (matrix == "sim5") {
352             for(int j=0;j<ncols;j++) {
353                 start = 0.0;
354                 for(int i=0;i<nrows;i++) {
355                     probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n);
356                     start = start + rowtotal[i]/double(n);
357                 }
358             }
359         }
360         
361         else if (matrix == "sim6") {
362             for(int i=0;i<nrows;i++) {
363                 for(int j=0;j<ncols;j++) {
364                     probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n*nrows);
365                     start = start + columntotal[j]/double(n*nrows);
366                 }
367             }
368         }
369         
370         
371         else if (matrix == "sim7") {
372             for(int i=0;i<nrows;i++) {
373                 for(int j=0;j<ncols;j++) {
374                     probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n*ncols);
375                     start = start + rowtotal[i]/double(n*ncols);
376                 }
377             }
378         }
379         
380         else if (matrix == "sim8") {
381             for(int i=0;i<nrows;i++) {
382                 for(int j=0;j<ncols;j++) {
383                     probabilityMatrix[ncols * i + j] = start + (rowtotal[i]*columntotal[j])/double(n*n);
384                     start = start + (rowtotal[i]*columntotal[j])/double(n*n);
385                 }
386             }
387         }
388         else if (matrix == "sim9" || matrix == "sim2") { }
389         else {
390             m->mothurOut("[ERROR]: No model selected! \n");
391             m->control_pressed = true;
392         }
393         
394         
395         //co_matrix is the transposed shared file, initmatrix is the original shared file
396         if (metric == "cscore") { initscore = trial.calc_c_score(co_matrix, rowtotal, ncols, nrows); }
397         else if (metric == "checker") { initscore = trial.calc_checker(co_matrix, rowtotal, ncols, nrows); }
398         else if (metric == "vratio") { initscore = trial.calc_vratio(nrows, ncols, rowtotal, columntotal); }
399         else if (metric == "combo") { initscore = trial.calc_combo(nrows, ncols, co_matrix); }
400         else { m->mothurOut("[ERROR]: No metric selected!\n"); m->control_pressed = true; return 1; }
401         
402         m->mothurOut("Initial c score: " + toString(initscore)); m->mothurOutEndLine();
403         
404         double previous;
405         double current;
406         double randnum;
407         int count;
408
409         //burn-in for sim9    
410         if(matrix == "sim9") {
411             for(int i=0;i<10000;i++) trial.swap_checkerboards (co_matrix, ncols, nrows);
412         }
413
414         //populate null matrix from probability matrix, do this a lot.
415         for(int k=0;k<runs;k++){
416             nullmatrix.clear();
417             //zero-fill the null matrix
418             nullmatrix.assign(nrows, vector<int>(ncols, 0));
419             
420             if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
421                 count = 0;
422                 while(count < n) {
423                     if (m->control_pressed) { return 0; }
424                 nextnum2:
425                     previous = 0.0;
426                     randnum = rand() / double(RAND_MAX);
427                     for(int i=0;i<nrows;i++) {
428                         for(int j=0;j<ncols;j++) {
429                             current = probabilityMatrix[ncols * i + j];
430                             if(randnum <= current && randnum > previous) {
431                                 nullmatrix[i][j] = 1;
432                                 count++;
433                                 if (count > n) break;
434                                 else
435                                     goto nextnum2;
436                             }
437                             previous = current;
438                         }
439                     }
440                 }
441             }
442             
443             else if (matrix == "sim2") {
444                 for(int i=0;i<nrows;i++) {
445                     random_shuffle( co_matrix[i].begin(), co_matrix[i].end() ); 
446                 }
447                 //do this for the scoring since those all have nullmatrix as a parameter
448                 //nullmatrix gets cleared at the begining of each run
449                 nullmatrix = co_matrix;
450             }
451             
452             else if(matrix == "sim4") {
453                 for(int i=0;i<nrows;i++) {
454                     count = 0;
455                     while(count < rowtotal[i]) {
456                         previous = 0.0;
457                         if (m->control_pressed) { return 0; }
458                         randnum = rand() / double(RAND_MAX);
459                         for(int j=0;j<ncols;j++) {
460                             current = probabilityMatrix[ncols * i + j];
461                             if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
462                                 nullmatrix[i][j] = 1;
463                                 count++;
464                                 previous = 0.0;
465                                 break;
466                             }
467                             previous = current;
468                         }
469                     }
470                 }
471             }
472             
473             else if(matrix == "sim3" || matrix == "sim5") {
474                 //columns
475                 for(int j=0;j<ncols;j++) {
476                     count = 0;
477                     while(count < columntotal[j]) {
478                         if (m->control_pressed) { return 0; }
479                         randnum = rand() / double(RAND_MAX);
480                         for(int i=0;i<nrows;i++) {
481                             current = probabilityMatrix[ncols * i + j];
482                             if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
483                                 nullmatrix[i][j] = 1;
484                                 count++;
485                                 previous = 0.0;
486                                 break;
487                             }
488                             previous = current;
489                         }
490                     }
491                 }
492             }
493             
494             //swap_checkerboards takes the original matrix and swaps checkerboards
495             else if(matrix == "sim9") {
496                 trial.swap_checkerboards (co_matrix, ncols, nrows);
497                 nullmatrix = co_matrix;
498             }
499             else {
500                 m->mothurOut("[ERROR]: No null model selected!\n\n"); m->control_pressed = true;
501                 return 1;
502             }
503             
504             //run metric on null matrix and add score to the stats vector
505             if (metric == "cscore"){
506                 stats.push_back(trial.calc_c_score(nullmatrix, rowtotal, ncols, nrows));
507             }
508             else if (metric == "checker") {
509                 stats.push_back(trial.calc_checker(nullmatrix, rowtotal, ncols, nrows));
510             }
511             else if (metric == "vratio") {
512                 stats.push_back(trial.calc_vratio(nrows, ncols, rowtotal, columntotal));
513             }
514             else if (metric == "combo") {
515                 stats.push_back(trial.calc_combo(nrows, ncols, nullmatrix));
516             }
517             else {
518                 m->mothurOut("[ERROR]: No metric selected!\n\n"); m->control_pressed = true;
519                 return 1;
520             }
521             
522         }
523         
524         
525         
526         double total = 0.0;
527         for (int i=0; i<stats.size();i++) { total+=stats[i]; }
528         
529         double nullMean = double (total/(double)stats.size());
530         
531         m->mothurOutEndLine(); m->mothurOut("average metric score: " + toString(nullMean)); m->mothurOutEndLine();
532         
533         //calc_p_value is not a statistical p-value, it's just the average that are either > or < the initscore.
534         //All it does is show what is expected in a competitively structured community
535         //zscore is output so p-value can be looked up in a ztable
536         double pvalue = 0.0;
537         if (metric == "cscore" || metric == "checker") { pvalue = trial.calc_pvalue_greaterthan (stats, initscore); }
538         else{ pvalue = trial.calc_pvalue_lessthan (stats, initscore); }
539
540         double sd = trial.getSD(runs, stats, nullMean);
541
542         double zscore = trial.get_zscore(sd, nullMean, initscore);
543         
544         m->mothurOut("zscore: " + toString(zscore)); m->mothurOutEndLine();
545         m->mothurOut("standard deviation: " + toString(sd)); m->mothurOutEndLine();
546         out << metric << '\t' << thisLookUp[0]->getLabel() << '\t' << nullMean << '\t' << zscore << endl;
547         
548         return 0;
549     }
550     catch(exception& e) {
551         m->errorOut(e, "CooccurrenceCommand", "Cooccurrence");
552         exit(1);
553     }
554 }
555 //**********************************************************************************************************************
556
557