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 string CooccurrenceCommand::getOutputFileNameTag(string type, string inputName=""){
  55         try {
  56         string outputFileName = "";
  57                 map<string, vector<string> >::iterator it;
  58
  59         //is this a type this command creates
  60         it = outputTypes.find(type);
  61         if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
  62         else {
  63             if (type == "summary") {  outputFileName =  "cooccurence.summary"; }
  64             else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
  65         }
  66         return outputFileName;
  67         }
  68         catch(exception& e) {
  69                 m->errorOut(e, "CooccurrenceCommand", "getOutputFileNameTag");
  70                 exit(1);
  71         }
  72 }
  73 //**********************************************************************************************************************
  74 CooccurrenceCommand::CooccurrenceCommand(){
  75         try {
  76                 abort = true; calledHelp = true;
  77                 setParameters();
  78         vector<string> tempOutNames;
  79                 outputTypes["summary"] = tempOutNames;
  80
  81         }
  82         catch(exception& e) {
  83                 m->errorOut(e, "CooccurrenceCommand", "CooccurrenceCommand");
  84                 exit(1);
  85         }
  86 }
  87 //**********************************************************************************************************************
  88 CooccurrenceCommand::CooccurrenceCommand(string option) {
  89         try {
  90                 abort = false; calledHelp = false;
  91                 allLines = 1;
  92
  93                 //allow user to run help
  94                 if(option == "help") { help(); abort = true; calledHelp = true; }
  95                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  96
  97                 else {
  98                         vector<string> myArray = setParameters();
  99
 100                         OptionParser parser(option);
 101                         map<string,string> parameters = parser.getParameters();
 102                         map<string,string>::iterator it;
 103
 104                         ValidParameters validParameter;
 105
 106                         //check to make sure all parameters are valid for command
 107                         for (it = parameters.begin(); it != parameters.end(); it++) {
 108                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 109                         }
 110
 111
 112                         //if the user changes the input directory command factory will send this info to us in the output parameter
 113                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 114                         if (inputDir == "not found"){   inputDir = "";          }
 115                         else {
 116                                 string path;
 117                                 it = parameters.find("shared");
 118                                 //user has given a template file
 119                                 if(it != parameters.end()){
 120                                         path = m->hasPath(it->second);
 121                                         //if the user has not given a path then, add inputdir. else leave path alone.
 122                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 123                                 }
 124                         }
 125
 126             vector<string> tempOutNames;
 127             outputTypes["summary"] = tempOutNames;
 128
 129                 //check for optional parameter and set defaults
 130                         // ...at some point should added some additional type checking...
 131                         label = validParameter.validFile(parameters, "label", false);
 132                         if (label == "not found") { label = ""; }
 133                         else {
 134                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 135                                 else { allLines = 1;  }
 136                         }
 137
 138                         //get shared file
 139                         sharedfile = validParameter.validFile(parameters, "shared", true);
 140                         if (sharedfile == "not open") { sharedfile = ""; abort = true; }
 141                         else if (sharedfile == "not found") {
 142                                 //if there is a current shared file, use it
 143                                 sharedfile = m->getSharedFile();
 144                                 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 145                                 else {  m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
 146                         }else { m->setSharedFile(sharedfile); }
 147
 148
 149                         //if the user changes the output directory command factory will send this info to us in the output parameter
 150                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(sharedfile);             }
 151
 152
 153                         metric = validParameter.validFile(parameters, "metric", false);                         if (metric == "not found") { metric = "cscore"; }
 154
 155                         if ((metric != "cscore") && (metric != "checker") && (metric != "combo") && (metric != "vratio")) {
 156                                 m->mothurOut("[ERROR]: " + metric + " is not a valid metric option for the cooccurrence command. Choices are cscore, checker, combo, vratio."); m->mothurOutEndLine(); abort = true;
 157                         }
 158
 159                         matrix = validParameter.validFile(parameters, "matrixmodel", false);                            if (matrix == "not found") { matrix = "sim2"; }
 160
 161                         if ((matrix != "sim1") && (matrix != "sim2") && (matrix != "sim3") && (matrix != "sim4") && (matrix != "sim5" ) && (matrix != "sim6" ) && (matrix != "sim7" ) && (matrix != "sim8" ) && (matrix != "sim9" )) {
 162                                 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;
 163                         }
 164
 165             groups = validParameter.validFile(parameters, "groups", false);
 166                         if (groups == "not found") { groups = "";   }
 167                         else {
 168                                 m->splitAtDash(groups, Groups);
 169                         }
 170                         m->setGroups(Groups);
 171
 172             string temp = validParameter.validFile(parameters, "iters", false);                 if (temp == "not found") { temp = "1000"; }
 173                         m->mothurConvert(temp, runs);
 174
 175                 }
 176
 177         }
 178         catch(exception& e) {
 179                 m->errorOut(e, "CooccurrenceCommand", "CooccurrenceCommand");
 180                 exit(1);
 181         }
 182 }
 183 //**********************************************************************************************************************
 184
 185 int CooccurrenceCommand::execute(){
 186         try {
 187
 188                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 189
 190                 InputData* input = new InputData(sharedfile, "sharedfile");
 191                 vector<SharedRAbundVector*> lookup = input->getSharedRAbundVectors();
 192                 string lastLabel = lookup[0]->getLabel();
 193
 194                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 195                 set<string> processedLabels;
 196                 set<string> userLabels = labels;
 197
 198         ofstream out;
 199                 string outputFileName = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + getOutputFileNameTag("summary");
 200         m->openOutputFile(outputFileName, out);
 201         outputNames.push_back(outputFileName);  outputTypes["summary"].push_back(outputFileName);
 202         out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 203         out << "metric\tlabel\tScore\tzScore\tstandardDeviation\n";
 204
 205                 //as long as you are not at the end of the file or done wih the lines you want
 206                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 207
 208                         if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
 209
 210                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 211
 212                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 213
 214                                 getCooccurrence(lookup, out);
 215
 216                                 processedLabels.insert(lookup[0]->getLabel());
 217                                 userLabels.erase(lookup[0]->getLabel());
 218                         }
 219
 220                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 221                                 string saveLabel = lookup[0]->getLabel();
 222
 223                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 224                                 lookup = input->getSharedRAbundVectors(lastLabel);
 225                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 226                                 getCooccurrence(lookup, out);
 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
 235                         lastLabel = lookup[0]->getLabel();
 236                         //prevent memory leak
 237                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 238
 239                         if (m->control_pressed) {  outputTypes.clear(); delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
 240
 241                         //get next line to process
 242                         lookup = input->getSharedRAbundVectors();
 243                 }
 244
 245                 if (m->control_pressed) { delete input; out.close(); m->mothurRemove(outputFileName); return 0; }
 246
 247                 //output error messages about any remaining user labels
 248                 set<string>::iterator it;
 249                 bool needToRun = false;
 250                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 251                         m->mothurOut("Your file does not include the label " + *it);
 252                         if (processedLabels.count(lastLabel) != 1) {
 253                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 254                                 needToRun = true;
 255                         }else {
 256                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 257                         }
 258                 }
 259
 260                 //run last label if you need to
 261                 if (needToRun == true)  {
 262                         for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 263                         lookup = input->getSharedRAbundVectors(lastLabel);
 264
 265                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 266
 267                         getCooccurrence(lookup, out);
 268
 269                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 270                 }
 271
 272         out.close();
 273
 274                 //reset groups parameter
 275                 delete input;
 276         m->clearGroups();
 277
 278         m->mothurOutEndLine();
 279                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 280                 m->mothurOut(outputFileName); m->mothurOutEndLine();
 281                 m->mothurOutEndLine();
 282
 283                 return 0;
 284         }
 285         catch(exception& e) {
 286                 m->errorOut(e, "CooccurrenceCommand", "execute");
 287                 exit(1);
 288         }
 289 }
 290 //**********************************************************************************************************************
 291
 292 int CooccurrenceCommand::getCooccurrence(vector<SharedRAbundVector*>& thisLookUp, ofstream& out){
 293     try {
 294         int numOTUS = thisLookUp[0]->getNumBins();
 295
 296         if(numOTUS < 2) {
 297             m->mothurOut("Not enough OTUs for co-occurrence analysis, skipping"); m->mothurOutEndLine();
 298             return 0;
 299         }
 300
 301         vector< vector<int> > co_matrix; co_matrix.resize(thisLookUp[0]->getNumBins());
 302         for (int i = 0; i < thisLookUp[0]->getNumBins(); i++) { co_matrix[i].resize((thisLookUp.size()), 0); }
 303         vector<int> columntotal; columntotal.resize(thisLookUp.size(), 0);
 304         vector<int> rowtotal; rowtotal.resize(numOTUS, 0);
 305
 306         for (int i = 0; i < thisLookUp.size(); i++) { //nrows in the shared file
 307             for (int j = 0; j < thisLookUp[i]->getNumBins(); j++) { //cols of original shared file
 308                 if (m->control_pressed) { return 0; }
 309                 int abund = thisLookUp[i]->getAbundance(j);
 310
 311                 if(abund > 0) {
 312                     co_matrix[j][i] = 1;
 313                     rowtotal[j]++;
 314                     columntotal[i]++;
 315                 }
 316             }
 317         }
 318
 319         //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.
 320         //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!
 321         int nrows = numOTUS;//rows of inital matrix
 322         int ncols = thisLookUp.size();//groups
 323         double initscore = 0.0;
 324
 325         vector<double> stats;
 326         vector<double> probabilityMatrix; probabilityMatrix.resize(ncols * nrows, 0);
 327         vector<vector<int> > nullmatrix(nrows, vector<int>(ncols, 0));
 328
 329         TrialSwap2 trial;
 330
 331         int n = accumulate( columntotal.begin(), columntotal.end(), 0 );
 332
 333         //============================================================
 334
 335         //generate a probability matrix. Only do this once.
 336         float start = 0.0;
 337
 338         if (matrix == "sim1") {
 339             for(int i=0;i<nrows;i++) {
 340                 for(int j=0;j<ncols;j++) {
 341                     probabilityMatrix[ncols * i + j] = start + 1/double(nrows*ncols);
 342                     start = start + 1/double(nrows*ncols);
 343                 }
 344             }
 345         }
 346         //don't need a prob matrix because we just shuffle the rows, may use this in the future
 347         else if (matrix == "sim2") { }
 348 //            for(int i=0;i<nrows;i++) {
 349 //                start = 0.0;
 350 //                for(int j=0;j<ncols;j++) {
 351 //                    probabilityMatrix[ncols * i + j] = start + 1/double(ncols);
 352 //                    start = start + 1/double(ncols);
 353 //                }
 354 //            }
 355 //        }
 356
 357         else if (matrix == "sim3") {
 358             for(int j=0;j<ncols;j++) {
 359                 start = 0.0;
 360                 for(int i=0;i<nrows;i++) {
 361                     probabilityMatrix[ncols * i + j] = start + 1/double(nrows);
 362                     start = start + 1/double(nrows);
 363                 }
 364             }
 365         }
 366
 367         else if (matrix == "sim4") {
 368             for(int i=0;i<nrows;i++) {
 369                 start = 0.0;
 370                 for(int j=0;j<ncols;j++) {
 371                     probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n);
 372                     start = start + columntotal[j]/double(n);
 373                 }
 374             }
 375         }
 376
 377         else if (matrix == "sim5") {
 378             for(int j=0;j<ncols;j++) {
 379                 start = 0.0;
 380                 for(int i=0;i<nrows;i++) {
 381                     probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n);
 382                     start = start + rowtotal[i]/double(n);
 383                 }
 384             }
 385         }
 386
 387         else if (matrix == "sim6") {
 388             for(int i=0;i<nrows;i++) {
 389                 for(int j=0;j<ncols;j++) {
 390                     probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n*nrows);
 391                     start = start + columntotal[j]/double(n*nrows);
 392                 }
 393             }
 394         }
 395
 396
 397         else if (matrix == "sim7") {
 398             for(int i=0;i<nrows;i++) {
 399                 for(int j=0;j<ncols;j++) {
 400                     probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n*ncols);
 401                     start = start + rowtotal[i]/double(n*ncols);
 402                 }
 403             }
 404         }
 405
 406         else if (matrix == "sim8") {
 407             for(int i=0;i<nrows;i++) {
 408                 for(int j=0;j<ncols;j++) {
 409                     probabilityMatrix[ncols * i + j] = start + (rowtotal[i]*columntotal[j])/double(n*n);
 410                     start = start + (rowtotal[i]*columntotal[j])/double(n*n);
 411                 }
 412             }
 413         }
 414         else if (matrix == "sim9" || matrix == "sim2") { }
 415         else {
 416             m->mothurOut("[ERROR]: No model selected! \n");
 417             m->control_pressed = true;
 418         }
 419
 420
 421         //co_matrix is the transposed shared file, initmatrix is the original shared file
 422         if (metric == "cscore") { initscore = trial.calc_c_score(co_matrix, rowtotal, ncols, nrows); }
 423         else if (metric == "checker") { initscore = trial.calc_checker(co_matrix, rowtotal, ncols, nrows); }
 424         else if (metric == "vratio") { initscore = trial.calc_vratio(nrows, ncols, rowtotal, columntotal); }
 425         else if (metric == "combo") { initscore = trial.calc_combo(nrows, ncols, co_matrix); }
 426         else { m->mothurOut("[ERROR]: No metric selected!\n"); m->control_pressed = true; return 1; }
 427
 428         m->mothurOut("Initial c score: " + toString(initscore)); m->mothurOutEndLine();
 429
 430         double previous;
 431         double current;
 432         double randnum;
 433         int count;
 434
 435         //burn-in for sim9
 436         if(matrix == "sim9") {
 437             for(int i=0;i<10000;i++) trial.swap_checkerboards (co_matrix, ncols, nrows);
 438         }
 439
 440         //populate null matrix from probability matrix, do this a lot.
 441         for(int k=0;k<runs;k++){
 442             nullmatrix.clear();
 443             //zero-fill the null matrix
 444             nullmatrix.assign(nrows, vector<int>(ncols, 0));
 445
 446             if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
 447                 count = 0;
 448                 while(count < n) {
 449                     if (m->control_pressed) { return 0; }
 450                 nextnum2:
 451                     previous = 0.0;
 452                     randnum = rand() / double(RAND_MAX);
 453                     for(int i=0;i<nrows;i++) {
 454                         for(int j=0;j<ncols;j++) {
 455                             current = probabilityMatrix[ncols * i + j];
 456                             if(randnum <= current && randnum > previous) {
 457                                 nullmatrix[i][j] = 1;
 458                                 count++;
 459                                 if (count > n) break;
 460                                 else
 461                                     goto nextnum2;
 462                             }
 463                             previous = current;
 464                         }
 465                     }
 466                 }
 467             }
 468
 469             else if (matrix == "sim2") {
 470                 for(int i=0;i<nrows;i++) {
 471                     random_shuffle( co_matrix[i].begin(), co_matrix[i].end() );
 472                 }
 473                 //do this for the scoring since those all have nullmatrix as a parameter
 474                 //nullmatrix gets cleared at the begining of each run
 475                 nullmatrix = co_matrix;
 476             }
 477
 478             else if(matrix == "sim4") {
 479                 for(int i=0;i<nrows;i++) {
 480                     count = 0;
 481                     while(count < rowtotal[i]) {
 482                         previous = 0.0;
 483                         if (m->control_pressed) { return 0; }
 484                         randnum = rand() / double(RAND_MAX);
 485                         for(int j=0;j<ncols;j++) {
 486                             current = probabilityMatrix[ncols * i + j];
 487                             if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
 488                                 nullmatrix[i][j] = 1;
 489                                 count++;
 490                                 previous = 0.0;
 491                                 break;
 492                             }
 493                             previous = current;
 494                         }
 495                     }
 496                 }
 497             }
 498
 499             else if(matrix == "sim3" || matrix == "sim5") {
 500                 //columns
 501                 for(int j=0;j<ncols;j++) {
 502                     count = 0;
 503                     while(count < columntotal[j]) {
 504                         if (m->control_pressed) { return 0; }
 505                         randnum = rand() / double(RAND_MAX);
 506                         for(int i=0;i<nrows;i++) {
 507                             current = probabilityMatrix[ncols * i + j];
 508                             if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
 509                                 nullmatrix[i][j] = 1;
 510                                 count++;
 511                                 previous = 0.0;
 512                                 break;
 513                             }
 514                             previous = current;
 515                         }
 516                     }
 517                 }
 518             }
 519
 520             //swap_checkerboards takes the original matrix and swaps checkerboards
 521             else if(matrix == "sim9") {
 522                 trial.swap_checkerboards (co_matrix, ncols, nrows);
 523                 nullmatrix = co_matrix;
 524             }
 525             else {
 526                 m->mothurOut("[ERROR]: No null model selected!\n\n"); m->control_pressed = true;
 527                 return 1;
 528             }
 529
 530             //run metric on null matrix and add score to the stats vector
 531             if (metric == "cscore"){
 532                 stats.push_back(trial.calc_c_score(nullmatrix, rowtotal, ncols, nrows));
 533             }
 534             else if (metric == "checker") {
 535                 stats.push_back(trial.calc_checker(nullmatrix, rowtotal, ncols, nrows));
 536             }
 537             else if (metric == "vratio") {
 538                 stats.push_back(trial.calc_vratio(nrows, ncols, rowtotal, columntotal));
 539             }
 540             else if (metric == "combo") {
 541                 stats.push_back(trial.calc_combo(nrows, ncols, nullmatrix));
 542             }
 543             else {
 544                 m->mothurOut("[ERROR]: No metric selected!\n\n"); m->control_pressed = true;
 545                 return 1;
 546             }
 547
 548         }
 549
 550
 551
 552         double total = 0.0;
 553         for (int i=0; i<stats.size();i++) { total+=stats[i]; }
 554
 555         double nullMean = double (total/(double)stats.size());
 556
 557         m->mothurOutEndLine(); m->mothurOut("average metric score: " + toString(nullMean)); m->mothurOutEndLine();
 558
 559         //calc_p_value is not a statistical p-value, it's just the average that are either > or < the initscore.
 560         //All it does is show what is expected in a competitively structured community
 561         //zscore is output so p-value can be looked up in a ztable
 562         double pvalue = 0.0;
 563         if (metric == "cscore" || metric == "checker") { pvalue = trial.calc_pvalue_greaterthan (stats, initscore); }
 564         else{ pvalue = trial.calc_pvalue_lessthan (stats, initscore); }
 565
 566         double sd = trial.getSD(runs, stats, nullMean);
 567
 568         double zscore = trial.get_zscore(sd, nullMean, initscore);
 569
 570         m->mothurOut("zscore: " + toString(zscore)); m->mothurOutEndLine();
 571         m->mothurOut("standard deviation: " + toString(sd)); m->mothurOutEndLine();
 572         out << metric << '\t' << thisLookUp[0]->getLabel() << '\t' << nullMean << '\t' << zscore << '\t' << sd << endl;
 573
 574         return 0;
 575     }
 576     catch(exception& e) {
 577         m->errorOut(e, "CooccurrenceCommand", "Cooccurrence");
 578         exit(1);
 579     }
 580 }
 581 //**********************************************************************************************************************
 582
 583