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