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\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
  90                 else {
  91                         vector<string> myArray = setParameters();
  92
  93                         OptionParser parser(option);
  94                         map<string, string> parameters = parser.getParameters();
  95                         map<string,string>::iterator it;
  96
  97                         ValidParameters validParameter;
  98
  99                         //check to make sure all parameters are valid for command
 100                         for (it = parameters.begin(); it != parameters.end(); it++) {
 101                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 102                         }
 103
 104                         //initialize outputTypes
 105                         vector<string> tempOutNames;
 106                         outputTypes["coverage"] = tempOutNames;
 107                         outputTypes["libshuffsummary"] = tempOutNames;
 108
 109                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 110                         if (inputDir == "not found"){   inputDir = "";          }
 111                         else {
 112                                 string path;
 113                                 it = parameters.find("phylip");
 114                                 //user has given a template file
 115                                 if(it != parameters.end()){
 116                                         path = m->hasPath(it->second);
 117                                         //if the user has not given a path then, add inputdir. else leave path alone.
 118                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
 119                                 }
 120
 121                                 it = parameters.find("group");
 122                                 //user has given a template file
 123                                 if(it != parameters.end()){
 124                                         path = m->hasPath(it->second);
 125                                         //if the user has not given a path then, add inputdir. else leave path alone.
 126                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 127                                 }
 128                         }
 129
 130                         //check for required parameters
 131                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 132                         if (phylipfile == "not open") { phylipfile = ""; abort = true; }
 133                         else if (phylipfile == "not found") {
 134                                 phylipfile = m->getPhylipFile();
 135                                 if (phylipfile != "") {  m->mothurOut("Using " + phylipfile + " as input file for the phylip parameter."); m->mothurOutEndLine(); }
 136                                 else {
 137                                         m->mothurOut("You must provide a phylip file."); m->mothurOutEndLine();
 138                                         abort = true;
 139                                 }
 140                         }
 141
 142                         //check for required parameters
 143                         groupfile = validParameter.validFile(parameters, "group", true);
 144                         if (groupfile == "not open") { groupfile = ""; abort = true; }
 145                         else if (groupfile == "not found") {
 146                                 groupfile = m->getGroupFile();
 147                                 if (groupfile != "") {  m->mothurOut("Using " + groupfile + " as input file for the group parameter."); m->mothurOutEndLine(); }
 148                                 else {
 149                                         m->mothurOut("You must provide a group file."); m->mothurOutEndLine();
 150                                         abort = true;
 151                                 }
 152                         }
 153
 154                         //if the user changes the output directory command factory will send this info to us in the output parameter
 155                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 156                                 outputDir = "";
 157                                 outputDir += m->hasPath(phylipfile); //if user entered a file with a path then preserve it
 158                         }
 159
 160                         //check for optional parameter and set defaults
 161                         // ...at some point should added some additional type checking...
 162                         groups = validParameter.validFile(parameters, "groups", false);
 163                         if (groups == "not found") { groups = ""; savegroups = groups; }
 164                         else {
 165                                 savegroups = groups;
 166                                 m->splitAtDash(groups, Groups);
 167                                 m->Groups = Groups;
 168                         }
 169
 170                         string temp;
 171                         temp = validParameter.validFile(parameters, "iters", false);                            if (temp == "not found") { temp = "10000"; }
 172                         convert(temp, iters);
 173
 174                         temp = validParameter.validFile(parameters, "cutoff", false);                           if (temp == "not found") { temp = "1.0"; }
 175                         convert(temp, cutOff);
 176
 177                         temp = validParameter.validFile(parameters, "step", false);                             if (temp == "not found") { temp = "0.01"; }
 178                         convert(temp, step);
 179
 180                         temp = validParameter.validFile(parameters, "sim", false);                              if (temp == "not found") { temp = "F"; }
 181                         sim = m->isTrue(temp);
 182
 183                         userform = validParameter.validFile(parameters, "form", false);                 if (userform == "not found") { userform = "integral"; }
 184
 185                 }
 186
 187         }
 188         catch(exception& e) {
 189                 m->errorOut(e, "LibShuffCommand", "LibShuffCommand");
 190                 exit(1);
 191         }
 192 }
 193 //**********************************************************************************************************************
 194
 195 int LibShuffCommand::execute(){
 196         try {
 197
 198                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 199
 200                 //read files
 201                 groupMap = new GroupMap(groupfile);
 202                 int error = groupMap->readMap();
 203                 if (error == 1) { delete groupMap; return 0; }
 204
 205                 ifstream in;
 206                 m->openInputFile(phylipfile, in);
 207                 matrix = new FullMatrix(in, groupMap, sim); //reads the matrix file
 208                 in.close();
 209
 210                 if (m->control_pressed) { delete groupMap; delete matrix; return 0; }
 211
 212                 //if files don't match...
 213                 if (matrix->getNumSeqs() < groupMap->getNumSeqs()) {
 214                         m->mothurOut("Your distance file contains " + toString(matrix->getNumSeqs()) + " sequences, and your group file contains " + toString(groupMap->getNumSeqs()) + " sequences.");  m->mothurOutEndLine();
 215                         //create new group file
 216                         if(outputDir == "") { outputDir += m->hasPath(groupfile); }
 217
 218                         string newGroupFile = outputDir + m->getRootName(m->getSimpleName(groupfile)) + "editted.groups";
 219                         outputNames.push_back(newGroupFile);
 220                         ofstream outGroups;
 221                         m->openOutputFile(newGroupFile, outGroups);
 222
 223                         for (int i = 0; i < matrix->getNumSeqs(); i++) {
 224                                 if (m->control_pressed) { delete groupMap; delete matrix; outGroups.close(); remove(newGroupFile.c_str()); return 0; }
 225
 226                                 Names temp = matrix->getRowInfo(i);
 227                                 outGroups << temp.seqName << '\t' << temp.groupName << endl;
 228                         }
 229                         outGroups.close();
 230
 231                         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();
 232
 233                         //read new groupfile
 234                         delete groupMap;
 235                         groupfile = newGroupFile;
 236
 237                         groupMap = new GroupMap(groupfile);
 238                         groupMap->readMap();
 239
 240                         if (m->control_pressed) { delete groupMap; delete matrix; remove(newGroupFile.c_str()); return 0; }
 241                 }
 242
 243
 244                 setGroups();                                                            //set the groups to be analyzed and sorts them
 245
 246                 /********************************************************************************************/
 247                 //this is needed because when we read the matrix we sort it into groups in alphabetical order
 248                 //the rest of the command and the classes used in this command assume specific order
 249                 /********************************************************************************************/
 250                 matrix->setGroups(groupMap->namesOfGroups);
 251                 vector<int> sizes;
 252                 for (int i = 0; i < groupMap->namesOfGroups.size(); i++) {   sizes.push_back(groupMap->getNumSeqs(groupMap->namesOfGroups[i]));  }
 253                 matrix->setSizes(sizes);
 254
 255
 256                 if(userform == "discrete"){
 257                         form = new DLibshuff(matrix, iters, step, cutOff);
 258                 }
 259                 else{
 260                         form = new SLibshuff(matrix, iters, cutOff);
 261                 }
 262
 263                 savedDXYValues = form->evaluateAll();
 264                 savedMinValues = form->getSavedMins();
 265
 266                 if (m->control_pressed) {  delete form; m->Groups.clear(); delete matrix; delete groupMap; return 0; }
 267
 268                 pValueCounts.resize(numGroups);
 269                 for(int i=0;i<numGroups;i++){
 270                         pValueCounts[i].assign(numGroups, 0);
 271                 }
 272
 273                 if (m->control_pressed) {  outputTypes.clear(); delete form; m->Groups.clear(); delete matrix; delete groupMap; return 0; }
 274
 275                 Progress* reading = new Progress();
 276
 277                 for(int i=0;i<numGroups-1;i++) {
 278                         for(int j=i+1;j<numGroups;j++) {
 279
 280                                 if (m->control_pressed) {  outputTypes.clear();  delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 281
 282                                 reading->newLine(groupNames[i]+'-'+groupNames[j], iters);
 283                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 284                                 int spotj = groupMap->groupIndex[groupNames[j]];
 285
 286                                 for(int p=0;p<iters;p++) {
 287
 288                                         if (m->control_pressed) {  outputTypes.clear(); delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 289
 290                                         form->randomizeGroups(spoti,spotj);
 291                                         if(form->evaluatePair(spoti,spotj) >= savedDXYValues[spoti][spotj])     {       pValueCounts[i][j]++;   }
 292                                         if(form->evaluatePair(spotj,spoti) >= savedDXYValues[spotj][spoti])     {       pValueCounts[j][i]++;   }
 293
 294                                         if (m->control_pressed) {  outputTypes.clear(); delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 295
 296                                         reading->update(p);
 297                                 }
 298                                 form->resetGroup(spoti);
 299                                 form->resetGroup(spotj);
 300                         }
 301                 }
 302
 303                 if (m->control_pressed) { outputTypes.clear();  delete form; m->Groups.clear(); delete matrix; delete groupMap; delete reading; return 0; }
 304
 305                 reading->finish();
 306                 delete reading;
 307
 308                 m->mothurOutEndLine();
 309                 printSummaryFile();
 310                 printCoverageFile();
 311
 312                 //clear out users groups
 313                 m->Groups.clear();
 314                 delete form;
 315
 316                 delete matrix; delete groupMap;
 317
 318                 if (m->control_pressed) {  outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } return 0; }
 319
 320                 m->mothurOutEndLine();
 321                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 322                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 323                 m->mothurOutEndLine();
 324
 325                 return 0;
 326         }
 327         catch(exception& e) {
 328                 m->errorOut(e, "LibShuffCommand", "execute");
 329                 exit(1);
 330         }
 331 }
 332
 333 //**********************************************************************************************************************
 334
 335 int LibShuffCommand::printCoverageFile() {
 336         try {
 337
 338                 ofstream outCov;
 339                 summaryFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "libshuff.coverage";
 340                 m->openOutputFile(summaryFile, outCov);
 341                 outputNames.push_back(summaryFile); outputTypes["coverage"].push_back(summaryFile);
 342                 outCov.setf(ios::fixed, ios::floatfield); outCov.setf(ios::showpoint);
 343                 //cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 344
 345                 map<double,vector<int> > allDistances;
 346                 map<double,vector<int> >::iterator it;
 347
 348                 vector<vector<int> > indices(numGroups);
 349                 int numIndices = numGroups * numGroups;
 350
 351                 int index = 0;
 352                 for(int i=0;i<numGroups;i++){
 353                         indices[i].assign(numGroups,0);
 354                         for(int j=0;j<numGroups;j++){
 355                                 indices[i][j] = index++;
 356
 357                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 358                                 int spotj = groupMap->groupIndex[groupNames[j]];
 359
 360                                 for(int k=0;k<savedMinValues[spoti][spotj].size();k++){
 361
 362                                         if(m->control_pressed)  { outCov.close(); return 0; }
 363
 364                                         if(allDistances[savedMinValues[spoti][spotj][k]].size() != 0){
 365                                                 allDistances[savedMinValues[spoti][spotj][k]][indices[i][j]]++;
 366                                         }
 367                                         else{
 368                                                 allDistances[savedMinValues[spoti][spotj][k]].assign(numIndices, 0);
 369                                                 allDistances[savedMinValues[spoti][spotj][k]][indices[i][j]] = 1;
 370                                         }
 371                                 }
 372                         }
 373                 }
 374                 it=allDistances.begin();
 375
 376                 //cout << setprecision(8);
 377
 378                 vector<int> prevRow = it->second;
 379                 it++;
 380
 381                 for(;it!=allDistances.end();it++){
 382                         for(int i=0;i<it->second.size();i++){
 383                                 it->second[i] += prevRow[i];
 384                         }
 385                         prevRow = it->second;
 386                 }
 387
 388                 vector<int> lastRow = allDistances.rbegin()->second;
 389                 outCov << setprecision(8);
 390
 391                 outCov << "dist";
 392                 for (int i = 0; i < numGroups; i++){
 393                         outCov << '\t' << groupNames[i];
 394                 }
 395                 for (int i=0;i<numGroups;i++){
 396                         for(int j=i+1;j<numGroups;j++){
 397                                 if(m->control_pressed)  { outCov.close(); return 0; }
 398                                 outCov << '\t' << groupNames[i] << '-' << groupNames[j] << '\t';
 399                                 outCov << groupNames[j] << '-' << groupNames[i];
 400                         }
 401                 }
 402                 outCov << endl;
 403
 404                 for(it=allDistances.begin();it!=allDistances.end();it++){
 405                         outCov << it->first << '\t';
 406                         for(int i=0;i<numGroups;i++){
 407                                 outCov << it->second[indices[i][i]]/(float)lastRow[indices[i][i]] << '\t';
 408                         }
 409                         for(int i=0;i<numGroups;i++){
 410                                 for(int j=i+1;j<numGroups;j++){
 411                                         if(m->control_pressed)  { outCov.close(); return 0; }
 412
 413                                         outCov << it->second[indices[i][j]]/(float)lastRow[indices[i][j]] << '\t';
 414                                         outCov << it->second[indices[j][i]]/(float)lastRow[indices[j][i]] << '\t';
 415                                 }
 416                         }
 417                         outCov << endl;
 418                 }
 419                 outCov.close();
 420
 421                 return 0;
 422         }
 423         catch(exception& e) {
 424                 m->errorOut(e, "LibShuffCommand", "printCoverageFile");
 425                 exit(1);
 426         }
 427 }
 428
 429 //**********************************************************************************************************************
 430
 431 int LibShuffCommand::printSummaryFile() {
 432         try {
 433
 434                 ofstream outSum;
 435                 summaryFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "libshuff.summary";
 436                 m->openOutputFile(summaryFile, outSum);
 437                 outputNames.push_back(summaryFile); outputTypes["libshuffsummary"].push_back(summaryFile);
 438
 439                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 440                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 441
 442                 cout << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl;
 443                 m->mothurOutJustToLog("Comparison\tdCXYScore\tSignificance"); m->mothurOutEndLine();
 444                 outSum << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl;
 445
 446                 int precision = (int)log10(iters);
 447                 for(int i=0;i<numGroups;i++){
 448                         for(int j=i+1;j<numGroups;j++){
 449                                 if(m->control_pressed)  { outSum.close(); return 0; }
 450
 451                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 452                                 int spotj = groupMap->groupIndex[groupNames[j]];
 453
 454                                 if(pValueCounts[i][j]){
 455                                         cout << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << setprecision(precision) << pValueCounts[i][j]/(float)iters << endl;
 456                                         m->mothurOutJustToLog(groupNames[i]+"-"+groupNames[j] + "\t" + toString(savedDXYValues[spoti][spotj]) + "\t" + toString((pValueCounts[i][j]/(float)iters))); m->mothurOutEndLine();
 457                                         outSum << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << setprecision(precision) << pValueCounts[i][j]/(float)iters << endl;
 458                                 }
 459                                 else{
 460                                         cout << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << '<' <<setprecision(precision) << 1/(float)iters << endl;
 461                                         m->mothurOutJustToLog(groupNames[i]+"-"+groupNames[j] + "\t" + toString(savedDXYValues[spoti][spotj]) + "\t" + toString((1/(float)iters))); m->mothurOutEndLine();
 462                                         outSum << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << '<' <<setprecision(precision) << 1/(float)iters << endl;
 463                                 }
 464                                 if(pValueCounts[j][i]){
 465                                         cout << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << setprecision (precision) << pValueCounts[j][i]/(float)iters << endl;
 466                                         m->mothurOutJustToLog(groupNames[j]+"-"+groupNames[i] + "\t" + toString(savedDXYValues[spotj][spoti]) + "\t" + toString((pValueCounts[j][i]/(float)iters))); m->mothurOutEndLine();
 467                                         outSum << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << setprecision (precision) << pValueCounts[j][i]/(float)iters << endl;
 468                                 }
 469                                 else{
 470                                         cout << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << '<' <<setprecision (precision) << 1/(float)iters << endl;
 471                                         m->mothurOutJustToLog(groupNames[j]+"-"+groupNames[i] + "\t" + toString(savedDXYValues[spotj][spoti]) + "\t" + toString((1/(float)iters))); m->mothurOutEndLine();
 472                                         outSum << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << '<' <<setprecision (precision) << 1/(float)iters << endl;
 473                                 }
 474                         }
 475                 }
 476
 477                 outSum.close();
 478                 return 0;
 479         }
 480         catch(exception& e) {
 481                 m->errorOut(e, "LibShuffCommand", "printSummaryFile");
 482                 exit(1);
 483         }
 484 }
 485
 486 //**********************************************************************************************************************
 487
 488 void LibShuffCommand::setGroups() {
 489         try {
 490                 //if the user has not entered specific groups to analyze then do them all
 491                 if (m->Groups.size() == 0) {
 492                         numGroups = groupMap->getNumGroups();
 493                         for (int i=0; i < numGroups; i++) {
 494                                 m->Groups.push_back(groupMap->namesOfGroups[i]);
 495                         }
 496                 } else {
 497                         if (savegroups != "all") {
 498                                 //check that groups are valid
 499                                 for (int i = 0; i < m->Groups.size(); i++) {
 500                                         if (groupMap->isValidGroup(m->Groups[i]) != true) {
 501                                                 m->mothurOut(m->Groups[i] + " is not a valid group, and will be disregarded."); m->mothurOutEndLine();
 502                                                 // erase the invalid group from globaldata->Groups
 503                                                 m->Groups.erase(m->Groups.begin()+i);
 504                                         }
 505                                 }
 506
 507                                 //if the user only entered invalid groups
 508                                 if ((m->Groups.size() == 0) || (m->Groups.size() == 1)) {
 509                                         numGroups = groupMap->getNumGroups();
 510                                         for (int i=0; i < numGroups; i++) {
 511                                                 m->Groups.push_back(groupMap->namesOfGroups[i]);
 512                                         }
 513                                         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();
 514                                 } else { numGroups = m->Groups.size(); }
 515                         } else { //users wants all groups
 516                                 numGroups = groupMap->getNumGroups();
 517                                 m->Groups.clear();
 518                                 for (int i=0; i < numGroups; i++) {
 519                                         m->Groups.push_back(groupMap->namesOfGroups[i]);
 520                                 }
 521                         }
 522                 }
 523
 524                 //sort so labels match
 525                 sort(m->Groups.begin(), m->Groups.end());
 526
 527                 //sort
 528                 sort(groupMap->namesOfGroups.begin(), groupMap->namesOfGroups.end());
 529
 530                 for (int i = 0; i < groupMap->namesOfGroups.size(); i++) {  groupMap->groupIndex[groupMap->namesOfGroups[i]] = i;  }
 531
 532                 groupNames = m->Groups;
 533
 534         }
 535         catch(exception& e) {
 536                 m->errorOut(e, "LibShuffCommand", "setGroups");
 537                 exit(1);
 538         }
 539 }
 540
 541 /***********************************************************/