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\tpValue\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> > 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);
 281
 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);
 286
 287                                 if(abund > 0) {
 288                                     initmatrix[i][j] = 1;
 289                     co_matrix[j][i] = 1;
 290                     rowtotal[j]++;
 291                     columntotal[i]++;
 292                                 }
 293                         }
 294         }
 295
 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;
 301
 302         vector<int> columntotal; columntotal.resize(ncols, 0);
 303         vector<int> rowtotal; rowtotal.resize(nrows, 0);
 304         vector<double> stats;
 305         double probabilityMatrix[ncols * nrows];
 306         vector<vector<int> > nullmatrix(nrows, vector<int>(ncols, 0));
 307
 308         TrialSwap2 trial;
 309
 310         n = accumulate( columntotal.begin(), columntotal.end(), 0 );
 311
 312         //============================================================
 313
 314         //generate a probability matrix. Only do this once.
 315         float start = 0.0;
 316
 317         if (matrix == "sim1") {
 318             for(int i=0;i<nrows;i++) {
 319                 for(int j=0;j<ncols;j++) {
 320                     probabilityMatrix[ncols * i + j] = start + 1/double(nrows*ncols);
 321                     start = start + 1/double(nrows*ncols);
 322                 }
 323             }
 324         }
 325         else if (matrix == "sim2") {
 326             for(int i=0;i<nrows;i++) {
 327                 start = 0.0;
 328                 for(int j=0;j<ncols;j++) {
 329                     probabilityMatrix[ncols * i + j] = start + 1/double(ncols);
 330                     start = start + 1/double(ncols);
 331                 }
 332             }
 333         }
 334
 335         else if (matrix == "sim3") {
 336             for(int j=0;j<ncols;j++) {
 337                 start = 0.0;
 338                 for(int i=0;i<nrows;i++) {
 339                     probabilityMatrix[ncols * i + j] = start + 1/double(nrows);
 340                     start =  start + 1/double(nrows);
 341                 }
 342             }
 343         }
 344
 345         else if (matrix == "sim4") {
 346             for(int i=0;i<nrows;i++) {
 347                 start = 0.0;
 348                 for(int j=0;j<ncols;j++) {
 349                     probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n);
 350                     start =  start + columntotal[j]/double(n);
 351                 }
 352             }
 353         }
 354
 355         else if (matrix == "sim5") {
 356             for(int j=0;j<ncols;j++) {
 357                 start = 0.0;
 358                 for(int i=0;i<nrows;i++) {
 359                     probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n);
 360                     start =  start + rowtotal[i]/double(n);
 361                 }
 362             }
 363         }
 364
 365         else if (matrix == "sim6") {
 366             for(int i=0;i<nrows;i++) {
 367                 for(int j=0;j<ncols;j++) {
 368                     probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n*nrows);
 369                     start =  start + columntotal[j]/double(n*nrows);
 370                 }
 371             }
 372         }
 373
 374
 375         else if (matrix == "sim7") {
 376             for(int i=0;i<nrows;i++) {
 377                 for(int j=0;j<ncols;j++) {
 378                     probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n*ncols);
 379                     start =  start + rowtotal[i]/double(n*ncols);
 380                 }
 381             }
 382         }
 383
 384         else if (matrix == "sim8") {
 385             for(int i=0;i<nrows;i++) {
 386                 for(int j=0;j<ncols;j++) {
 387                     probabilityMatrix[ncols * i + j] = start + (rowtotal[i]*columntotal[j])/double(n*n);
 388                     start =  start + (rowtotal[i]*columntotal[j])/double(n*n);
 389                 }
 390             }
 391         }
 392
 393         else {
 394             if(sim != 9) {
 395                 m->mothurOut("[ERROR]: No model selected! \n");
 396                 m->control_pressed = true;
 397             }
 398         }
 399
 400
 401
 402         if (metric == "cscore")         { initscore = trial.calc_c_score(initmatrix, rowtotal, ncols, nrows);    }
 403         else if (metric == "checker")   { initscore = trial.calc_checker(initmatrix, rowtotal, ncols, nrows);    }
 404         else if (metric == "vratio")    { initscore = trial.calc_vratio(nrows, ncols, rowtotal, columntotal);   }
 405         else if (metric == "combo")     { initscore = trial.calc_combo(nrows, ncols, initmatrix);                }
 406         else                            {  m->mothurOut("[ERROR]: No metric selected!\n");  m->control_pressed = true; return 1;            }
 407
 408         m->mothurOut("Initial c score: " + toString(initscore)); m->mothurOutEndLine();
 409
 410         double previous;
 411         double current;
 412         double randnum;
 413         int count;
 414
 415         //burn-in
 416         for(int i=0;i<10000;i++){
 417             nullmatrix.clear();
 418             //zero-fill the null matrix
 419             nullmatrix.assign(nrows, vector<int>(ncols, 0));
 420
 421             if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
 422                 count = 0;
 423                 while(count < n) {
 424                     nextnum:
 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 nextnum;
 436                             }
 437                             previous = current;
 438                         }
 439                     }
 440                 }
 441             }
 442
 443             else if(matrix == "sim2" || matrix == "sim4") {
 444                     for(int i=0;i<nrows;i++) {
 445                         previous = 0.0;
 446                         count = 0;
 447                         while(count < rowtotal[i]) {
 448                             randnum = rand() / double(RAND_MAX);
 449                             for(int j=0;j<ncols;j++) {
 450                                 current = probabilityMatrix[ncols * i + j];
 451                                 if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
 452                                     nullmatrix[i][j] = 1;
 453                                     count++;
 454                                     previous = 0.0;
 455                                     break;
 456                                 }
 457                                 previous = current;
 458                             }
 459                         }
 460                     }
 461             }
 462
 463             else if(matrix == "sim3" || matrix == "sim5") {
 464                 //columns
 465                 for(int j=0;j<ncols;j++) {
 466                     count = 0;
 467                     while(count < columntotal[j]) {
 468                         randnum = rand() / double(RAND_MAX);
 469                         for(int i=0;i<nrows;i++) {
 470                             current = probabilityMatrix[ncols * i + j];
 471                             if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
 472                                 nullmatrix[i][j] = 1;
 473                                 count++;
 474                                 previous = 0.0;
 475                                 break;
 476                             }
 477                             previous = current;
 478                         }
 479                     }
 480                 }
 481             }
 482
 483         }
 484
 485         //populate null matrix from probability matrix, do this a lot.
 486         for(int i=0;i<runs;i++){
 487             nullmatrix.clear();
 488             //zero-fill the null matrix
 489             nullmatrix.assign(nrows, vector<int>(ncols, 0));
 490
 491             if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
 492                 count = 0;
 493                 while(count < n) {
 494                     nextnum:
 495                     previous = 0.0;
 496                     randnum = rand() / double(RAND_MAX);
 497                     for(int i=0;i<nrows;i++) {
 498                         for(int j=0;j<ncols;j++) {
 499                             current = probabilityMatrix[ncols * i + j];
 500                             if(randnum <= current && randnum > previous) {
 501                                 nullmatrix[i][j] = 1;
 502                                 count++;
 503                                 if (count > n) break;
 504                                 else
 505                                 goto nextnum;
 506                             }
 507                             previous = current;
 508                         }
 509                     }
 510                 }
 511             }
 512
 513             else if(matrix == "sim2" || matrix == "sim4") {
 514                     for(int i=0;i<nrows;i++) {
 515                         previous = 0.0;
 516                         count = 0;
 517                         while(count < rowtotal[i]) {
 518                             randnum = rand() / double(RAND_MAX);
 519                             for(int j=0;j<ncols;j++) {
 520                                 current = probabilityMatrix[ncols * i + j];
 521                                 if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
 522                                     nullmatrix[i][j] = 1;
 523                                     count++;
 524                                     previous = 0.0;
 525                                     break;
 526                                 }
 527                                 previous = current;
 528                             }
 529                         }
 530                     }
 531             }
 532
 533             else if(matrix == "sim3" || matrix == "sim5") {
 534                 //columns
 535                 for(int j=0;j<ncols;j++) {
 536                     count = 0;
 537                     while(count < columntotal[j]) {
 538                         randnum = rand() / double(RAND_MAX);
 539                         for(int i=0;i<nrows;i++) {
 540                             current = probabilityMatrix[ncols * i + j];
 541                             if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
 542                                 nullmatrix[i][j] = 1;
 543                                 count++;
 544                                 previous = 0.0;
 545                                 break;
 546                             }
 547                             previous = current;
 548                         }
 549                     }
 550                 }
 551             }
 552
 553             //swap_checkerboards takes the original matrix and swaps checkerboards
 554             else if(matrix == "sim9") {
 555                 trial.swap_checkerboards (initmatrix, rowtotal, columntotal, ncols, nrows);
 556             }
 557             else {
 558                 cout << "[ERROR]: No null model selected!\n" << endl;
 559                 return 1;
 560             }
 561
 562             //run metric on null matrix and add score to the stats vector
 563             if (metric == "cscore"){
 564                 stats.push_back(trial.calc_c_score(nullmatrix, rowtotal, ncols, nrows));
 565             }
 566             else if (metric == "checker") {
 567                 stats.push_back(trial.calc_checker(nullmatrix, rowtotal, ncols, nrows));
 568             }
 569             else if (metric == "vratio") {
 570                 stats.push_back(trial.calc_vratio(nrows, ncols, rowtotal, columntotal));
 571             }
 572             else if (metric == "combo") {
 573                 stats.push_back(trial.calc_combo(nrows, ncols, nullmatrix));
 574             }
 575             else {
 576                 cout << "[ERROR]: No metric selected!\n" << endl;
 577                 return 1;
 578             }
 579
 580         }
 581
 582
 583
 584         double total = 0.0;
 585         for (int i=0; i<stats.size();i++)   {   total+=stats[i];   }
 586
 587         double nullMean = double (total/(double)stats.size());
 588
 589         m->mothurOutEndLine(); m->mothurOut("average metric score: " + toString(nullMean)); m->mothurOutEndLine();
 590
 591         double pvalue = 0.0;
 592         if (metric == "cscore" || metric == "checker") {    pvalue = trial.calc_pvalue_greaterthan (stats, initscore);   }
 593         else{   pvalue = trial.calc_pvalue_lessthan (stats, initscore); }
 594
 595         m->mothurOut("pvalue: " + toString(pvalue)); m->mothurOutEndLine();
 596         out << metric << '\t' << thisLookUp[0]->getLabel() << '\t' << nullMean << '\t' << pvalue << endl;
 597
 598         return 0;
 599         }
 600         catch(exception& e) {
 601                 m->errorOut(e, "CooccurrenceCommand", "Cooccurrence");
 602                 exit(1);
 603         }
 604 }
 605 //**********************************************************************************************************************
 606
 607