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