metastatscommand.cpp

   1 /*
   2  *  metastatscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/16/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "metastatscommand.h"
  11 #include "sharedutilities.h"
  12
  13
  14 //**********************************************************************************************************************
  15 vector<string> MetaStatsCommand::setParameters(){
  16         try {
  17                 CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","metastats",false,true,true); parameters.push_back(pshared);
  18                 CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(pdesign);
  19                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
  20                 CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
  21                 CommandParameter pthreshold("threshold", "Number", "", "0.05", "", "", "","",false,false); parameters.push_back(pthreshold);
  22                 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
  23                 CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
  24                 CommandParameter psets("sets", "String", "", "", "", "", "","",false,false); parameters.push_back(psets);
  25                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  26                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  27
  28                 vector<string> myArray;
  29                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  30                 return myArray;
  31         }
  32         catch(exception& e) {
  33                 m->errorOut(e, "MetaStatsCommand", "setParameters");
  34                 exit(1);
  35         }
  36 }
  37 //**********************************************************************************************************************
  38 string MetaStatsCommand::getHelpString(){
  39         try {
  40                 string helpString = "";
  41                 helpString += "This command is based on the Metastats program, White, J.R., Nagarajan, N. & Pop, M. Statistical methods for detecting differentially abundant features in clinical metagenomic samples. PLoS Comput Biol 5, e1000352 (2009).\n";
  42                 helpString += "The metastats command outputs a .metastats file. \n";
  43                 helpString += "The metastats command parameters are shared, iters, threshold, groups, label, design, sets and processors.  The shared and design parameters are required, unless you have valid current files.\n";
  44                 helpString += "The design parameter allows you to assign your groups to sets when you are running metastat. mothur will run all pairwise comparisons of the sets. It is required. \n";
  45                 helpString += "The design file looks like the group file.  It is a 2 column tab delimited file, where the first column is the group name and the second column is the set the group belongs to.\n";
  46                 helpString += "The sets parameter allows you to specify which of the sets in your designfile you would like to analyze. The set names are separated by dashes. THe default is all sets in the designfile.\n";
  47                 helpString += "The iters parameter allows you to set number of bootstrap permutations for estimating null distribution of t statistic.  The default is 1000. \n";
  48                 helpString += "The threshold parameter allows you to set the significance level to reject null hypotheses (default 0.05).\n";
  49                 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like included. The group names are separated by dashes.\n";
  50                 helpString += "The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n";
  51                 helpString += "The processors parameter allows you to specify how many processors you would like to use.  The default is 1. \n";
  52                 helpString += "The metastats command should be in the following format: metastats(design=yourDesignFile).\n";
  53                 helpString += "Example metastats(design=temp.design, groups=A-B-C).\n";
  54                 helpString += "The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n";
  55                 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
  56                 return helpString;
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "MetaStatsCommand", "getHelpString");
  60                 exit(1);
  61         }
  62 }
  63 //**********************************************************************************************************************
  64 string MetaStatsCommand::getOutputPattern(string type) {
  65     try {
  66         string pattern = "";
  67
  68         if (type == "metastats") {  pattern = "[filename],[distance],[group],metastats"; }
  69         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  70
  71         return pattern;
  72     }
  73     catch(exception& e) {
  74         m->errorOut(e, "MetaStatsCommand", "getOutputPattern");
  75         exit(1);
  76     }
  77 }
  78 //**********************************************************************************************************************
  79 MetaStatsCommand::MetaStatsCommand(){
  80         try {
  81                 abort = true; calledHelp = true;
  82                 setParameters();
  83                 vector<string> tempOutNames;
  84                 outputTypes["metastats"] = tempOutNames;
  85         }
  86         catch(exception& e) {
  87                 m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
  88                 exit(1);
  89         }
  90 }
  91 //**********************************************************************************************************************
  92
  93 MetaStatsCommand::MetaStatsCommand(string option) {
  94         try {
  95                 abort = false; calledHelp = false;
  96                 allLines = 1;
  97
  98
  99                 //allow user to run help
 100                 if(option == "help") { help(); abort = true; calledHelp = true; }
 101                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 102
 103                 else {
 104                         vector<string> myArray = setParameters();
 105
 106                         OptionParser parser(option);
 107                         map<string,string> parameters = parser.getParameters();
 108
 109                         ValidParameters validParameter;
 110
 111                         //check to make sure all parameters are valid for command
 112                         map<string,string>::iterator it;
 113                         for (it = parameters.begin(); it != parameters.end(); it++) {
 114                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 115                         }
 116
 117                         //initialize outputTypes
 118                         vector<string> tempOutNames;
 119                         outputTypes["metastats"] = tempOutNames;
 120
 121
 122                         //if the user changes the input directory command factory will send this info to us in the output parameter
 123                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 124                         if (inputDir == "not found"){   inputDir = "";          }
 125                         else {
 126                                 string path;
 127                                 it = parameters.find("design");
 128                                 //user has given a template file
 129                                 if(it != parameters.end()){
 130                                         path = m->hasPath(it->second);
 131                                         //if the user has not given a path then, add inputdir. else leave path alone.
 132                                         if (path == "") {       parameters["design"] = inputDir + it->second;           }
 133                                 }
 134
 135                                 it = parameters.find("shared");
 136                                 //user has given a template file
 137                                 if(it != parameters.end()){
 138                                         path = m->hasPath(it->second);
 139                                         //if the user has not given a path then, add inputdir. else leave path alone.
 140                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 141                                 }
 142
 143                         }
 144
 145                         //check for required parameters
 146                         sharedfile = validParameter.validFile(parameters, "shared", true);
 147                         if (sharedfile == "not open") { abort = true; }
 148                         else if (sharedfile == "not found") {                           //if there is a current shared file, use it
 149                                 sharedfile = m->getSharedFile();
 150                                 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 151                                 else {  m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
 152                         }else { m->setSharedFile(sharedfile); }
 153
 154                         //check for required parameters
 155                         designfile = validParameter.validFile(parameters, "design", true);
 156                         if (designfile == "not open") { abort = true; }
 157                         else if (designfile == "not found") {
 158                                 //if there is a current design file, use it
 159                                 designfile = m->getDesignFile();
 160                                 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
 161                                 else {  m->mothurOut("You have no current designfile and the design parameter is required."); m->mothurOutEndLine(); abort = true; }
 162                         }else { m->setDesignFile(designfile); }
 163
 164                         //if the user changes the output directory command factory will send this info to us in the output parameter
 165                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 166                                 outputDir = "";
 167                                 outputDir += m->hasPath(sharedfile); //if user entered a file with a path then preserve it
 168                         }
 169
 170                         //check for optional parameter and set defaults
 171                         // ...at some point should added some additional type checking...
 172                         label = validParameter.validFile(parameters, "label", false);
 173                         if (label == "not found") { label = ""; }
 174                         else {
 175                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 176                                 else { allLines = 1;  }
 177                         }
 178
 179                         groups = validParameter.validFile(parameters, "groups", false);
 180                         if (groups == "not found") { groups = ""; pickedGroups = false; }
 181                         else {
 182                                 pickedGroups = true;
 183                                 m->splitAtDash(groups, Groups);
 184                                 m->setGroups(Groups);
 185                         }
 186
 187                         sets = validParameter.validFile(parameters, "sets", false);
 188                         if (sets == "not found") { sets = ""; }
 189                         else {
 190                                 m->splitAtDash(sets, Sets);
 191                         }
 192
 193
 194                         string temp = validParameter.validFile(parameters, "iters", false);                     if (temp == "not found") { temp = "1000"; }
 195                         m->mothurConvert(temp, iters);
 196
 197                         temp = validParameter.validFile(parameters, "threshold", false);                        if (temp == "not found") { temp = "0.05"; }
 198                         m->mothurConvert(temp, threshold);
 199
 200                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 201                         m->setProcessors(temp);
 202                         m->mothurConvert(temp, processors);
 203                 }
 204
 205         }
 206         catch(exception& e) {
 207                 m->errorOut(e, "MetaStatsCommand", "MetaStatsCommand");
 208                 exit(1);
 209         }
 210 }
 211 //**********************************************************************************************************************
 212
 213 int MetaStatsCommand::execute(){
 214         try {
 215
 216                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 217
 218                 designMap = new GroupMap(designfile);
 219                 designMap->readDesignMap();
 220
 221                 input = new InputData(sharedfile, "sharedfile");
 222                 lookup = input->getSharedRAbundVectors();
 223                 string lastLabel = lookup[0]->getLabel();
 224
 225                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 226                 set<string> processedLabels;
 227                 set<string> userLabels = labels;
 228
 229                 //setup the pairwise comparions of sets for metastats
 230                 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
 231                 //make sure sets are all in designMap
 232                 SharedUtil* util = new SharedUtil();
 233                 vector<string> dGroups = designMap->getNamesOfGroups();
 234                 util->setGroups(Sets, dGroups);
 235                 delete util;
 236
 237                 int numGroups = Sets.size();
 238                 for (int a=0; a<numGroups; a++) {
 239                         for (int l = 0; l < a; l++) {
 240                                 vector<string> groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]);
 241                                 namesOfGroupCombos.push_back(groups);
 242                         }
 243                 }
 244
 245
 246                 //only 1 combo
 247                 if (numGroups == 2) { processors = 1; }
 248                 else if (numGroups < 2) { m->mothurOut("Not enough sets, I need at least 2 valid sets. Unable to complete command."); m->mothurOutEndLine(); m->control_pressed = true; }
 249
 250         if(processors != 1){
 251             int numPairs = namesOfGroupCombos.size();
 252             int numPairsPerProcessor = numPairs / processors;
 253
 254             for (int i = 0; i < processors; i++) {
 255                 int startPos = i * numPairsPerProcessor;
 256                 if(i == processors - 1){
 257                     numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
 258                 }
 259                 lines.push_back(linePair(startPos, numPairsPerProcessor));
 260             }
 261         }
 262
 263                 //as long as you are not at the end of the file or done wih the lines you want
 264                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 265
 266                         if (m->control_pressed) {  outputTypes.clear(); for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } m->clearGroups(); delete input; delete designMap;  for (int i = 0; i < outputNames.size(); i++) {       m->mothurRemove(outputNames[i]); } return 0; }
 267
 268                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 269
 270                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 271                                 process(lookup);
 272
 273                                 processedLabels.insert(lookup[0]->getLabel());
 274                                 userLabels.erase(lookup[0]->getLabel());
 275                         }
 276
 277                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 278                                 string saveLabel = lookup[0]->getLabel();
 279
 280                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 281                                 lookup = input->getSharedRAbundVectors(lastLabel);
 282                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 283
 284                                 process(lookup);
 285
 286                                 processedLabels.insert(lookup[0]->getLabel());
 287                                 userLabels.erase(lookup[0]->getLabel());
 288
 289                                 //restore real lastlabel to save below
 290                                 lookup[0]->setLabel(saveLabel);
 291                         }
 292
 293                         lastLabel = lookup[0]->getLabel();
 294                         //prevent memory leak
 295                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 296
 297                         if (m->control_pressed) {  outputTypes.clear(); m->clearGroups(); delete input;  delete designMap;  for (int i = 0; i < outputNames.size(); i++) {      m->mothurRemove(outputNames[i]); } return 0; }
 298
 299                         //get next line to process
 300                         lookup = input->getSharedRAbundVectors();
 301                 }
 302
 303                 if (m->control_pressed) { outputTypes.clear(); m->clearGroups(); delete input; delete designMap;  for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); }  return 0; }
 304
 305                 //output error messages about any remaining user labels
 306                 set<string>::iterator it;
 307                 bool needToRun = false;
 308                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 309                         m->mothurOut("Your file does not include the label " + *it);
 310                         if (processedLabels.count(lastLabel) != 1) {
 311                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 312                                 needToRun = true;
 313                         }else {
 314                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 315                         }
 316                 }
 317
 318                 //run last label if you need to
 319                 if (needToRun == true)  {
 320                         for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 321                         lookup = input->getSharedRAbundVectors(lastLabel);
 322
 323                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 324
 325                         process(lookup);
 326
 327                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 328                 }
 329
 330                 //reset groups parameter
 331                 m->clearGroups();
 332                 delete input;
 333                 delete designMap;
 334
 335                 if (m->control_pressed) { outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) {   m->mothurRemove(outputNames[i]); } return 0;}
 336
 337                 m->mothurOutEndLine();
 338                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 339                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 340                 m->mothurOutEndLine();
 341
 342                 return 0;
 343         }
 344         catch(exception& e) {
 345                 m->errorOut(e, "MetaStatsCommand", "execute");
 346                 exit(1);
 347         }
 348 }
 349 //**********************************************************************************************************************
 350
 351 int MetaStatsCommand::process(vector<SharedRAbundVector*>& thisLookUp){
 352         try {
 353
 354
 355                                 if(processors == 1){
 356                                         driver(0, namesOfGroupCombos.size(), thisLookUp);
 357                                 }else{
 358                                         int process = 1;
 359                                         vector<int> processIDS;
 360                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 361                                         //loop through and create all the processes you want
 362                                         while (process != processors) {
 363                                                 int pid = fork();
 364
 365                                                 if (pid > 0) {
 366                                                         processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 367                                                         process++;
 368                                                 }else if (pid == 0){
 369                                                         driver(lines[process].start, lines[process].num, thisLookUp);
 370                                                         exit(0);
 371                                                 }else {
 372                                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 373                                                         for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 374                                                         exit(0);
 375                                                 }
 376                                         }
 377
 378                                         //do my part
 379                                         driver(lines[0].start, lines[0].num, thisLookUp);
 380
 381                                         //force parent to wait until all the processes are done
 382                                         for (int i=0;i<(processors-1);i++) {
 383                                                 int temp = processIDS[i];
 384                                                 wait(&temp);
 385                                         }
 386         #else
 387
 388                     //////////////////////////////////////////////////////////////////////////////////////////////////////
 389                     //Windows version shared memory, so be careful when passing variables through the summarySharedData struct.
 390                     //Above fork() will clone, so memory is separate, but that's not the case with windows,
 391                     //Taking advantage of shared memory to pass results vectors.
 392                     //////////////////////////////////////////////////////////////////////////////////////////////////////
 393
 394                     vector<metastatsData*> pDataArray;
 395                     DWORD   dwThreadIdArray[processors-1];
 396                     HANDLE  hThreadArray[processors-1];
 397
 398                     //Create processor worker threads.
 399                     for( int i=1; i<processors; i++ ){
 400
 401                         //make copy of lookup so we don't get access violations
 402                         vector<SharedRAbundVector*> newLookup;
 403                         vector<string> designMapGroups;
 404                         for (int k = 0; k < thisLookUp.size(); k++) {
 405                             SharedRAbundVector* temp = new SharedRAbundVector();
 406                             temp->setLabel(thisLookUp[k]->getLabel());
 407                             temp->setGroup(thisLookUp[k]->getGroup());
 408                             newLookup.push_back(temp);
 409                             designMapGroups.push_back(designMap->getGroup(thisLookUp[k]->getGroup()));
 410                         }
 411
 412                         //for each bin
 413                         for (int k = 0; k < thisLookUp[0]->getNumBins(); k++) {
 414                             if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
 415                             for (int j = 0; j < thisLookUp.size(); j++) { newLookup[j]->push_back(thisLookUp[j]->getAbundance(k), thisLookUp[j]->getGroup()); }
 416                         }
 417
 418                         // Allocate memory for thread data.
 419                         metastatsData* tempSum = new metastatsData(sharedfile, outputDir, m, lines[i].start, lines[i].num, namesOfGroupCombos, newLookup, designMapGroups, iters, threshold);
 420                         pDataArray.push_back(tempSum);
 421                         processIDS.push_back(i);
 422
 423                         hThreadArray[i-1] = CreateThread(NULL, 0, MyMetastatsThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
 424                     }
 425
 426                     //do my part
 427                                         driver(lines[0].start, lines[0].num, thisLookUp);
 428
 429                     //Wait until all threads have terminated.
 430                     WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 431
 432                     //Close all thread handles and free memory allocations.
 433                     for(int i=0; i < pDataArray.size(); i++){
 434                         for (int j = 0; j < pDataArray[i]->thisLookUp.size(); j++) {  delete pDataArray[i]->thisLookUp[j];  }
 435                         for (int j = 0; j < pDataArray[i]->outputNames.size(); j++) {
 436                             outputNames.push_back(pDataArray[i]->outputNames[j]);
 437                             outputTypes["metastats"].push_back(pDataArray[i]->outputNames[j]);
 438                         }
 439
 440                         CloseHandle(hThreadArray[i]);
 441                         delete pDataArray[i];
 442                     }
 443         #endif
 444
 445                                 }
 446
 447                 return 0;
 448
 449         }
 450         catch(exception& e) {
 451                 m->errorOut(e, "MetaStatsCommand", "process");
 452                 exit(1);
 453         }
 454 }
 455 //**********************************************************************************************************************
 456 int MetaStatsCommand::driver(int start, int num, vector<SharedRAbundVector*>& thisLookUp) {
 457         try {
 458
 459                 //for each combo
 460                 for (int c = start; c < (start+num); c++) {
 461
 462                         //get set names
 463                         string setA = namesOfGroupCombos[c][0];
 464                         string setB = namesOfGroupCombos[c][1];
 465
 466                         //get filename
 467             map<string, string> variables;
 468             variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
 469             variables["[distance]"] = thisLookUp[0]->getLabel();
 470             variables["[group]"] = setA + "-" + setB;
 471                         string outputFileName = getOutputFileName("metastats",variables);
 472                         outputNames.push_back(outputFileName); outputTypes["metastats"].push_back(outputFileName);
 473                         //int nameLength = outputFileName.length();
 474                         //char * output = new char[nameLength];
 475                         //strcpy(output, outputFileName.c_str());
 476
 477                         //build matrix from shared rabunds
 478                         //double** data;
 479                         //data = new double*[thisLookUp[0]->getNumBins()];
 480
 481                         vector< vector<double> > data2; data2.resize(thisLookUp[0]->getNumBins());
 482
 483                         vector<SharedRAbundVector*> subset;
 484                         int setACount = 0;
 485                         int setBCount = 0;
 486                         for (int i = 0; i < thisLookUp.size(); i++) {
 487                                 string thisGroup = thisLookUp[i]->getGroup();
 488
 489                                 //is this group for a set we want to compare??
 490                                 //sorting the sets by putting setB at the back and setA in the front
 491                                 if ((designMap->getGroup(thisGroup) == setB)) {
 492                                         subset.push_back(thisLookUp[i]);
 493                                         setBCount++;
 494                                 }else if ((designMap->getGroup(thisGroup) == setA)) {
 495                                         subset.insert(subset.begin()+setACount, thisLookUp[i]);
 496                                         setACount++;
 497                                 }
 498                         }
 499
 500                         if ((setACount == 0) || (setBCount == 0))  {
 501                                 m->mothurOut("Missing shared info for " + setA + " or " + setB + ". Skipping comparison."); m->mothurOutEndLine();
 502                                 outputNames.pop_back();
 503                         }else {
 504
 505                                 //fill data
 506                                 for (int j = 0; j < thisLookUp[0]->getNumBins(); j++) {
 507                                         //data[j] = new double[subset.size()];
 508                                         data2[j].resize(subset.size(), 0.0);
 509
 510                                         for (int i = 0; i < subset.size(); i++) {
 511                                                 data2[j][i] = (subset[i]->getAbundance(j));
 512                                         }
 513                                 }
 514
 515                                 m->mothurOut("Comparing " + setA + " and " + setB + "..."); m->mothurOutEndLine();
 516                                 //metastat_main(output, thisLookUp[0]->getNumBins(), subset.size(), threshold, iters, data, setACount);
 517
 518                                 m->mothurOutEndLine();
 519                                 MothurMetastats mothurMeta(threshold, iters);
 520                                 mothurMeta.runMetastats(outputFileName , data2, setACount);
 521                                 m->mothurOutEndLine();
 522                                 m->mothurOutEndLine();
 523                         }
 524
 525                         //free memory
 526                         //delete output;
 527                         //for(int i = 0; i < thisLookUp[0]->getNumBins(); i++)  {  delete[] data[i];  }
 528                         //delete[] data;
 529                 }
 530
 531                 return 0;
 532
 533         }
 534         catch(exception& e) {
 535                 m->errorOut(e, "MetaStatsCommand", "driver");
 536                 exit(1);
 537         }
 538 }
 539 //**********************************************************************************************************************