libshuffcommand.cpp

   1 /*
   2  *  libshuffcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 3/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 /* This class is designed to implement an integral form of the Cramer-von Mises statistic.
  11         you may refer to the "Integration of Microbial Ecology and Statistics: A Test To Compare Gene Libraries"
  12         paper in Applied and Environmental Microbiology, Sept. 2004, p. 5485-5492 0099-2240/04/$8.00+0
  13         DOI: 10.1128/AEM.70.9.5485-5492.2004 Copyright 2004 American Society for Microbiology for more information. */
  14
  15
  16 #include "libshuffcommand.h"
  17 #include "libshuff.h"
  18 #include "slibshuff.h"
  19 #include "dlibshuff.h"
  20
  21
  22
  23 //**********************************************************************************************************************
  24 vector<string> LibShuffCommand::setParameters(){
  25         try {
  26                 CommandParameter pphylip("phylip", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pphylip);
  27                 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup);
  28                 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
  29                 CommandParameter piters("iters", "Number", "", "10000", "", "", "",false,false); parameters.push_back(piters);
  30                 CommandParameter pstep("step", "Number", "", "0.01", "", "", "",false,false); parameters.push_back(pstep);
  31                 CommandParameter pcutoff("cutoff", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pcutoff);
  32                 CommandParameter pform("form", "Multiple", "discrete-integral", "integral", "", "", "",false,false); parameters.push_back(pform);
  33                 CommandParameter psim("sim", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(psim);
  34                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  35                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  36
  37                 vector<string> myArray;
  38                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  39                 return myArray;
  40         }
  41         catch(exception& e) {
  42                 m->errorOut(e, "LibShuffCommand", "setParameters");
  43                 exit(1);
  44         }
  45 }
  46 //**********************************************************************************************************************
  47 string LibShuffCommand::getHelpString(){
  48         try {
  49                 string helpString = "";
  50                 helpString += "The libshuff command parameters are phylip, group, sim, groups, iters, step, form and cutoff.  phylip and group parameters are required, unless you have valid current files.\n";
  51                 helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed.  You must enter at least 2 valid groups.\n";
  52                 helpString += "The group names are separated by dashes.  The iters parameter allows you to specify how many random matrices you would like compared to your matrix.\n";
  53                 helpString += "The step parameter allows you to specify change in distance you would like between each output if you are using the discrete form.\n";
  54                 helpString += "The form parameter allows you to specify if you would like to analyze your matrix using the discrete or integral form. Your options are integral or discrete.\n";
  55                 helpString += "The libshuff command should be in the following format: libshuff(groups=yourGroups, iters=yourIters, cutOff=yourCutOff, form=yourForm, step=yourStep).\n";
  56                 helpString += "Example libshuff(groups=A-B-C, iters=500, form=discrete, step=0.01, cutOff=2.0).\n";
  57                 helpString += "The default value for groups is all the groups in your groupfile, iters is 10000, cutoff is 1.0, form is integral and step is 0.01.\n";
  58                 helpString += "The libshuff command output two files: .coverage and .slsummary their descriptions are in the manual.\n";
  59                 helpString += "Note: No spaces between parameter labels (i.e. iters), '=' and parameters (i.e.yourIters).\n";
  60                 return helpString;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "LibShuffCommand", "getHelpString");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 LibShuffCommand::LibShuffCommand(){
  69         try {
  70                 abort = true; calledHelp = true;
  71                 setParameters();
  72                 vector<string> tempOutNames;
  73                 outputTypes["coverage"] = tempOutNames;
  74                 outputTypes["libshuffsummary"] = tempOutNames;
  75         }
  76         catch(exception& e) {
  77                 m->errorOut(e, "LibShuffCommand", "LibShuffCommand");
  78                 exit(1);
  79         }
  80 }
  81
  82 //**********************************************************************************************************************
  83 LibShuffCommand::LibShuffCommand(string option)  {
  84         try {
  85                 abort = false; calledHelp = false;
  86
  87                 //allow user to run help
  88                 if(option == "help") { help(); abort = true; calledHelp = true; }
  89                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  90
  91                 else {
  92                         vector<string> myArray = setParameters();
  93
  94                         OptionParser parser(option);
  95                         map<string, string> parameters = parser.getParameters();
  96                         map<string,string>::iterator it;
  97
  98                         ValidParameters validParameter;
  99
 100                         //check to make sure all parameters are valid for command
 101                         for (it = parameters.begin(); it != parameters.end(); it++) {
 102                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 103                         }
 104
 105                         //initialize outputTypes
 106                         vector<string> tempOutNames;
 107                         outputTypes["coverage"] = tempOutNames;
 108                         outputTypes["libshuffsummary"] = tempOutNames;
 109
 110                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 111                         if (inputDir == "not found"){   inputDir = "";          }
 112                         else {
 113                                 string path;
 114                                 it = parameters.find("phylip");
 115                                 //user has given a template file
 116                                 if(it != parameters.end()){
 117                                         path = m->hasPath(it->second);
 118                                         //if the user has not given a path then, add inputdir. else leave path alone.
 119                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
 120                                 }
 121
 122                                 it = parameters.find("group");
 123                                 //user has given a template file
 124                                 if(it != parameters.end()){
 125                                         path = m->hasPath(it->second);
 126                                         //if the user has not given a path then, add inputdir. else leave path alone.
 127                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 128                                 }
 129                         }
 130
 131                         //check for required parameters
 132                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 133                         if (phylipfile == "not open") { phylipfile = ""; abort = true; }
 134                         else if (phylipfile == "not found") {
 135                                 phylipfile = m->getPhylipFile();
 136                                 if (phylipfile != "") {  m->mothurOut("Using " + phylipfile + " as input file for the phylip parameter."); m->mothurOutEndLine(); }
 137                                 else {
 138                                         m->mothurOut("You must provide a phylip file."); m->mothurOutEndLine();
 139                                         abort = true;
 140                                 }
 141                         }
 142
 143                         //check for required parameters
 144                         groupfile = validParameter.validFile(parameters, "group", true);
 145                         if (groupfile == "not open") { groupfile = ""; abort = true; }
 146                         else if (groupfile == "not found") {
 147                                 groupfile = m->getGroupFile();
 148                                 if (groupfile != "") {  m->mothurOut("Using " + groupfile + " as input file for the group parameter."); m->mothurOutEndLine(); }
 149                                 else {
 150                                         m->mothurOut("You must provide a group file."); m->mothurOutEndLine();
 151                                         abort = true;
 152                                 }
 153                         }
 154
 155                         //if the user changes the output directory command factory will send this info to us in the output parameter
 156                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 157                                 outputDir = "";
 158                                 outputDir += m->hasPath(phylipfile); //if user entered a file with a path then preserve it
 159                         }
 160
 161                         //check for optional parameter and set defaults
 162                         // ...at some point should added some additional type checking...
 163                         groups = validParameter.validFile(parameters, "groups", false);
 164                         if (groups == "not found") { groups = ""; savegroups = groups; }
 165                         else {
 166                                 savegroups = groups;
 167                                 m->splitAtDash(groups, Groups);
 168                                 m->Groups = Groups;
 169                         }
 170
 171                         string temp;
 172                         temp = validParameter.validFile(parameters, "iters", false);                            if (temp == "not found") { temp = "10000"; }
 173                         convert(temp, iters);
 174
 175                         temp = validParameter.validFile(parameters, "cutoff", false);                           if (temp == "not found") { temp = "1.0"; }
 176                         convert(temp, cutOff);
 177
 178                         temp = validParameter.validFile(parameters, "step", false);                             if (temp == "not found") { temp = "0.01"; }
 179                         convert(temp, step);
 180
 181                         temp = validParameter.validFile(parameters, "sim", false);                              if (temp == "not found") { temp = "F"; }
 182                         sim = m->isTrue(temp);
 183
 184                         userform = validParameter.validFile(parameters, "form", false);                 if (userform == "not found") { userform = "integral"; }
 185
 186                 }
 187
 188         }
 189         catch(exception& e) {
 190                 m->errorOut(e, "LibShuffCommand", "LibShuffCommand");
 191                 exit(1);
 192         }
 193 }
 194 //**********************************************************************************************************************
 195
 196 int LibShuffCommand::execute(){
 197         try {
 198
 199                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 200
 201                 //read files
 202                 groupMap = new GroupMap(groupfile);
 203                 int error = groupMap->readMap();
 204                 if (error == 1) { delete groupMap; return 0; }
 205
 206                 ifstream in;
 207                 m->openInputFile(phylipfile, in);
 208                 matrix = new FullMatrix(in, groupMap, sim); //reads the matrix file
 209                 in.close();
 210
 211                 if (m->control_pressed) { delete groupMap; delete matrix; return 0; }
 212
 213                 //if files don't match...
 214                 if (matrix->getNumSeqs() < groupMap->getNumSeqs()) {
 215                         m->mothurOut("Your distance file contains " + toString(matrix->getNumSeqs()) + " sequences, and your group file contains " + toString(groupMap->getNumSeqs()) + " sequences.");  m->mothurOutEndLine();
 216                         //create new group file
 217                         if(outputDir == "") { outputDir += m->hasPath(groupfile); }
 218
 219                         string newGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "editted.groups";
 220                         outputNames.push_back(newGroupFile);
 221                         ofstream outGroups;
 222                         m->openOutputFile(newGroupFile, outGroups);
 223
 224                         for (int i = 0; i < matrix->getNumSeqs(); i++) {
 225                                 if (m->control_pressed) { delete groupMap; delete matrix; outGroups.close(); remove(newGroupFile.c_str()); return 0; }
 226
 227                                 Names temp = matrix->getRowInfo(i);
 228                                 outGroups << temp.seqName << '\t' << temp.groupName << endl;
 229                         }
 230                         outGroups.close();
 231
 232                         m->mothurOut(newGroupFile + " is a new group file containing only the sequence that are in your distance file. I will read this file instead."); m->mothurOutEndLine();
 233
 234                         //read new groupfile
 235                         delete groupMap;
 236                         groupfile = newGroupFile;
 237
 238                         groupMap = new GroupMap(groupfile);
 239                         groupMap->readMap();
 240
 241                         if (m->control_pressed) { delete groupMap; delete matrix; remove(newGroupFile.c_str()); return 0; }
 242                 }
 243
 244
 245                 setGroups();                                                            //set the groups to be analyzed and sorts them
 246
 247                 /********************************************************************************************/
 248                 //this is needed because when we read the matrix we sort it into groups in alphabetical order
 249                 //the rest of the command and the classes used in this command assume specific order
 250                 /********************************************************************************************/
 251                 matrix->setGroups(groupMap->namesOfGroups);
 252                 vector<int> sizes;
 253                 for (int i = 0; i < groupMap->namesOfGroups.size(); i++) {   sizes.push_back(groupMap->getNumSeqs(groupMap->namesOfGroups[i]));  }
 254                 matrix->setSizes(sizes);
 255
 256
 257                 if(userform == "discrete"){
 258                         form = new DLibshuff(matrix, iters, step, cutOff);
 259                 }
 260                 else{
 261                         form = new SLibshuff(matrix, iters, cutOff);
 262                 }
 263
 264                 savedDXYValues = form->evaluateAll();
 265                 savedMinValues = form->getSavedMins();
 266
 267                 if (m->control_pressed) {  delete form; m->Groups.clear(); delete matrix; delete groupMap; return 0; }
 268
 269                 pValueCounts.resize(numGroups);
 270                 for(int i=0;i<numGroups;i++){
 271                         pValueCounts[i].assign(numGroups, 0);
 272                 }
 273
 274                 if (m->control_pressed) {  outputTypes.clear(); delete form; m->Groups.clear(); delete matrix; delete groupMap; return 0; }
 275
 276                 Progress* reading = new Progress();
 277
 278                 for(int i=0;i<numGroups-1;i++) {
 279                         for(int j=i+1;j<numGroups;j++) {
 280
 281                                 if (m->control_pressed) {  outputTypes.clear();  delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 282
 283                                 reading->newLine(groupNames[i]+'-'+groupNames[j], iters);
 284                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 285                                 int spotj = groupMap->groupIndex[groupNames[j]];
 286
 287                                 for(int p=0;p<iters;p++) {
 288
 289                                         if (m->control_pressed) {  outputTypes.clear(); delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 290
 291                                         form->randomizeGroups(spoti,spotj);
 292                                         if(form->evaluatePair(spoti,spotj) >= savedDXYValues[spoti][spotj])     {       pValueCounts[i][j]++;   }
 293                                         if(form->evaluatePair(spotj,spoti) >= savedDXYValues[spotj][spoti])     {       pValueCounts[j][i]++;   }
 294
 295                                         if (m->control_pressed) {  outputTypes.clear(); delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 296
 297                                         reading->update(p);
 298                                 }
 299                                 form->resetGroup(spoti);
 300                                 form->resetGroup(spotj);
 301                         }
 302                 }
 303
 304                 if (m->control_pressed) { outputTypes.clear();  delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 305
 306                 reading->finish();
 307                 delete reading;
 308
 309                 m->mothurOutEndLine();
 310                 printSummaryFile();
 311                 printCoverageFile();
 312
 313                 //clear out users groups
 314                 m->Groups.clear();
 315                 delete form;
 316
 317                 delete matrix; delete groupMap;
 318
 319                 if (m->control_pressed) {  outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } return 0; }
 320
 321                 m->mothurOutEndLine();
 322                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 323                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 324                 m->mothurOutEndLine();
 325
 326                 return 0;
 327         }
 328         catch(exception& e) {
 329                 m->errorOut(e, "LibShuffCommand", "execute");
 330                 exit(1);
 331         }
 332 }
 333
 334 //**********************************************************************************************************************
 335
 336 int LibShuffCommand::printCoverageFile() {
 337         try {
 338
 339                 ofstream outCov;
 340                 summaryFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "libshuff.coverage";
 341                 m->openOutputFile(summaryFile, outCov);
 342                 outputNames.push_back(summaryFile); outputTypes["coverage"].push_back(summaryFile);
 343                 outCov.setf(ios::fixed, ios::floatfield); outCov.setf(ios::showpoint);
 344                 //cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 345
 346                 map<double,vector<int> > allDistances;
 347                 map<double,vector<int> >::iterator it;
 348
 349                 vector<vector<int> > indices(numGroups);
 350                 int numIndices = numGroups * numGroups;
 351
 352                 int index = 0;
 353                 for(int i=0;i<numGroups;i++){
 354                         indices[i].assign(numGroups,0);
 355                         for(int j=0;j<numGroups;j++){
 356                                 indices[i][j] = index++;
 357
 358                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 359                                 int spotj = groupMap->groupIndex[groupNames[j]];
 360
 361                                 for(int k=0;k<savedMinValues[spoti][spotj].size();k++){
 362
 363                                         if(m->control_pressed)  { outCov.close(); return 0; }
 364
 365                                         if(allDistances[savedMinValues[spoti][spotj][k]].size() != 0){
 366                                                 allDistances[savedMinValues[spoti][spotj][k]][indices[i][j]]++;
 367                                         }
 368                                         else{
 369                                                 allDistances[savedMinValues[spoti][spotj][k]].assign(numIndices, 0);
 370                                                 allDistances[savedMinValues[spoti][spotj][k]][indices[i][j]] = 1;
 371                                         }
 372                                 }
 373                         }
 374                 }
 375                 it=allDistances.begin();
 376
 377                 //cout << setprecision(8);
 378
 379                 vector<int> prevRow = it->second;
 380                 it++;
 381
 382                 for(;it!=allDistances.end();it++){
 383                         for(int i=0;i<it->second.size();i++){
 384                                 it->second[i] += prevRow[i];
 385                         }
 386                         prevRow = it->second;
 387                 }
 388
 389                 vector<int> lastRow = allDistances.rbegin()->second;
 390                 outCov << setprecision(8);
 391
 392                 outCov << "dist";
 393                 for (int i = 0; i < numGroups; i++){
 394                         outCov << '\t' << groupNames[i];
 395                 }
 396                 for (int i=0;i<numGroups;i++){
 397                         for(int j=i+1;j<numGroups;j++){
 398                                 if(m->control_pressed)  { outCov.close(); return 0; }
 399                                 outCov << '\t' << groupNames[i] << '-' << groupNames[j] << '\t';
 400                                 outCov << groupNames[j] << '-' << groupNames[i];
 401                         }
 402                 }
 403                 outCov << endl;
 404
 405                 for(it=allDistances.begin();it!=allDistances.end();it++){
 406                         outCov << it->first << '\t';
 407                         for(int i=0;i<numGroups;i++){
 408                                 outCov << it->second[indices[i][i]]/(float)lastRow[indices[i][i]] << '\t';
 409                         }
 410                         for(int i=0;i<numGroups;i++){
 411                                 for(int j=i+1;j<numGroups;j++){
 412                                         if(m->control_pressed)  { outCov.close(); return 0; }
 413
 414                                         outCov << it->second[indices[i][j]]/(float)lastRow[indices[i][j]] << '\t';
 415                                         outCov << it->second[indices[j][i]]/(float)lastRow[indices[j][i]] << '\t';
 416                                 }
 417                         }
 418                         outCov << endl;
 419                 }
 420                 outCov.close();
 421
 422                 return 0;
 423         }
 424         catch(exception& e) {
 425                 m->errorOut(e, "LibShuffCommand", "printCoverageFile");
 426                 exit(1);
 427         }
 428 }
 429
 430 //**********************************************************************************************************************
 431
 432 int LibShuffCommand::printSummaryFile() {
 433         try {
 434
 435                 ofstream outSum;
 436                 summaryFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "libshuff.summary";
 437                 m->openOutputFile(summaryFile, outSum);
 438                 outputNames.push_back(summaryFile); outputTypes["libshuffsummary"].push_back(summaryFile);
 439
 440                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 441                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 442
 443                 cout << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl;
 444                 m->mothurOutJustToLog("Comparison\tdCXYScore\tSignificance"); m->mothurOutEndLine();
 445                 outSum << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl;
 446
 447                 int precision = (int)log10(iters);
 448                 for(int i=0;i<numGroups;i++){
 449                         for(int j=i+1;j<numGroups;j++){
 450                                 if(m->control_pressed)  { outSum.close(); return 0; }
 451
 452                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 453                                 int spotj = groupMap->groupIndex[groupNames[j]];
 454
 455                                 if(pValueCounts[i][j]){
 456                                         cout << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << setprecision(precision) << pValueCounts[i][j]/(float)iters << endl;
 457                                         m->mothurOutJustToLog(groupNames[i]+"-"+groupNames[j] + "\t" + toString(savedDXYValues[spoti][spotj]) + "\t" + toString((pValueCounts[i][j]/(float)iters))); m->mothurOutEndLine();
 458                                         outSum << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << setprecision(precision) << pValueCounts[i][j]/(float)iters << endl;
 459                                 }
 460                                 else{
 461                                         cout << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << '<' <<setprecision(precision) << 1/(float)iters << endl;
 462                                         m->mothurOutJustToLog(groupNames[i]+"-"+groupNames[j] + "\t" + toString(savedDXYValues[spoti][spotj]) + "\t" + toString((1/(float)iters))); m->mothurOutEndLine();
 463                                         outSum << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << '<' <<setprecision(precision) << 1/(float)iters << endl;
 464                                 }
 465                                 if(pValueCounts[j][i]){
 466                                         cout << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << setprecision (precision) << pValueCounts[j][i]/(float)iters << endl;
 467                                         m->mothurOutJustToLog(groupNames[j]+"-"+groupNames[i] + "\t" + toString(savedDXYValues[spotj][spoti]) + "\t" + toString((pValueCounts[j][i]/(float)iters))); m->mothurOutEndLine();
 468                                         outSum << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << setprecision (precision) << pValueCounts[j][i]/(float)iters << endl;
 469                                 }
 470                                 else{
 471                                         cout << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << '<' <<setprecision (precision) << 1/(float)iters << endl;
 472                                         m->mothurOutJustToLog(groupNames[j]+"-"+groupNames[i] + "\t" + toString(savedDXYValues[spotj][spoti]) + "\t" + toString((1/(float)iters))); m->mothurOutEndLine();
 473                                         outSum << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << '<' <<setprecision (precision) << 1/(float)iters << endl;
 474                                 }
 475                         }
 476                 }
 477
 478                 outSum.close();
 479                 return 0;
 480         }
 481         catch(exception& e) {
 482                 m->errorOut(e, "LibShuffCommand", "printSummaryFile");
 483                 exit(1);
 484         }
 485 }
 486
 487 //**********************************************************************************************************************
 488
 489 void LibShuffCommand::setGroups() {
 490         try {
 491                 //if the user has not entered specific groups to analyze then do them all
 492                 if (m->Groups.size() == 0) {
 493                         numGroups = groupMap->getNumGroups();
 494                         for (int i=0; i < numGroups; i++) {
 495                                 m->Groups.push_back(groupMap->namesOfGroups[i]);
 496                         }
 497                 } else {
 498                         if (savegroups != "all") {
 499                                 //check that groups are valid
 500                                 for (int i = 0; i < m->Groups.size(); i++) {
 501                                         if (groupMap->isValidGroup(m->Groups[i]) != true) {
 502                                                 m->mothurOut(m->Groups[i] + " is not a valid group, and will be disregarded."); m->mothurOutEndLine();
 503                                                 // erase the invalid group from globaldata->Groups
 504                                                 m->Groups.erase(m->Groups.begin()+i);
 505                                         }
 506                                 }
 507
 508                                 //if the user only entered invalid groups
 509                                 if ((m->Groups.size() == 0) || (m->Groups.size() == 1)) {
 510                                         numGroups = groupMap->getNumGroups();
 511                                         for (int i=0; i < numGroups; i++) {
 512                                                 m->Groups.push_back(groupMap->namesOfGroups[i]);
 513                                         }
 514                                         m->mothurOut("When using the groups parameter you must have at least 2 valid groups. I will run the command using all the groups in your groupfile."); m->mothurOutEndLine();
 515                                 } else { numGroups = m->Groups.size(); }
 516                         } else { //users wants all groups
 517                                 numGroups = groupMap->getNumGroups();
 518                                 m->Groups.clear();
 519                                 for (int i=0; i < numGroups; i++) {
 520                                         m->Groups.push_back(groupMap->namesOfGroups[i]);
 521                                 }
 522                         }
 523                 }
 524
 525                 //sort so labels match
 526                 sort(m->Groups.begin(), m->Groups.end());
 527
 528                 //sort
 529                 sort(groupMap->namesOfGroups.begin(), groupMap->namesOfGroups.end());
 530
 531                 for (int i = 0; i < groupMap->namesOfGroups.size(); i++) {  groupMap->groupIndex[groupMap->namesOfGroups[i]] = i;  }
 532
 533                 groupNames = m->Groups;
 534
 535         }
 536         catch(exception& e) {
 537                 m->errorOut(e, "LibShuffCommand", "setGroups");
 538                 exit(1);
 539         }
 540 }
 541
 542 /***********************************************************/