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                         }else { m->setPhylipFile(phylipfile); }
 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                         }else { m->setGroupFile(groupfile); }
 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->setGroups(Groups);
 169                         }
 170
 171                         string temp;
 172                         temp = validParameter.validFile(parameters, "iters", false);                            if (temp == "not found") { temp = "10000"; }
 173                         m->mothurConvert(temp, iters);
 174
 175                         temp = validParameter.validFile(parameters, "cutoff", false);                           if (temp == "not found") { temp = "1.0"; }
 176                         m->mothurConvert(temp, cutOff);
 177
 178                         temp = validParameter.validFile(parameters, "step", false);                             if (temp == "not found") { temp = "0.01"; }
 179                         m->mothurConvert(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(); m->mothurRemove(newGroupFile); 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; m->mothurRemove(newGroupFile); return 0; }
 242                 }
 243
 244
 245                 setGroups();                                                            //set the groups to be analyzed and sorts them
 246
 247                 if (numGroups < 2) { m->mothurOut("[ERROR]: libshuff requires at least 2 groups, you only have " + toString(numGroups) + ", aborting."); m->mothurOutEndLine(); m->control_pressed = true; }
 248
 249                 if (m->control_pressed) { delete groupMap; delete matrix; return 0; }
 250
 251                 /********************************************************************************************/
 252                 //this is needed because when we read the matrix we sort it into groups in alphabetical order
 253                 //the rest of the command and the classes used in this command assume specific order
 254                 /********************************************************************************************/
 255                 matrix->setGroups(groupMap->getNamesOfGroups());
 256                 vector<int> sizes;
 257                 for (int i = 0; i < (groupMap->getNamesOfGroups()).size(); i++) {   sizes.push_back(groupMap->getNumSeqs((groupMap->getNamesOfGroups())[i]));  }
 258                 matrix->setSizes(sizes);
 259
 260
 261                 if(userform == "discrete"){
 262                         form = new DLibshuff(matrix, iters, step, cutOff);
 263                 }
 264                 else{
 265                         form = new SLibshuff(matrix, iters, cutOff);
 266                 }
 267
 268                 savedDXYValues = form->evaluateAll();
 269                 savedMinValues = form->getSavedMins();
 270
 271                 if (m->control_pressed) {  delete form; m->clearGroups(); delete matrix; delete groupMap; return 0; }
 272
 273                 pValueCounts.resize(numGroups);
 274                 for(int i=0;i<numGroups;i++){
 275                         pValueCounts[i].assign(numGroups, 0);
 276                 }
 277
 278                 if (m->control_pressed) {  outputTypes.clear(); delete form; m->clearGroups(); delete matrix; delete groupMap; return 0; }
 279
 280                 Progress* reading = new Progress();
 281
 282                 for(int i=0;i<numGroups-1;i++) {
 283                         for(int j=i+1;j<numGroups;j++) {
 284
 285                                 if (m->control_pressed) {  outputTypes.clear();  delete form; m->clearGroups(); delete matrix; delete groupMap; delete reading; return 0; }
 286
 287                                 reading->newLine(groupNames[i]+'-'+groupNames[j], iters);
 288                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 289                                 int spotj = groupMap->groupIndex[groupNames[j]];
 290
 291                                 for(int p=0;p<iters;p++) {
 292
 293                                         if (m->control_pressed) {  outputTypes.clear(); delete form; m->clearGroups(); delete matrix; delete groupMap; delete reading; return 0; }
 294
 295                                         form->randomizeGroups(spoti,spotj);
 296                                         if(form->evaluatePair(spoti,spotj) >= savedDXYValues[spoti][spotj])     {       pValueCounts[i][j]++;   }
 297                                         if(form->evaluatePair(spotj,spoti) >= savedDXYValues[spotj][spoti])     {       pValueCounts[j][i]++;   }
 298
 299                                         if (m->control_pressed) {  outputTypes.clear(); delete form; m->clearGroups(); delete matrix; delete groupMap; delete reading; return 0; }
 300
 301                                         reading->update(p);
 302                                 }
 303                                 form->resetGroup(spoti);
 304                                 form->resetGroup(spotj);
 305                         }
 306                 }
 307
 308                 if (m->control_pressed) { outputTypes.clear();  delete form; m->clearGroups(); delete matrix; delete groupMap; delete reading; return 0; }
 309
 310                 reading->finish();
 311                 delete reading;
 312
 313                 m->mothurOutEndLine();
 314                 printSummaryFile();
 315                 printCoverageFile();
 316
 317                 //clear out users groups
 318                 m->clearGroups();
 319                 delete form;
 320
 321                 delete matrix; delete groupMap;
 322
 323                 if (m->control_pressed) {  outputTypes.clear(); for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]); } return 0; }
 324
 325                 m->mothurOutEndLine();
 326                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 327                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 328                 m->mothurOutEndLine();
 329
 330                 return 0;
 331         }
 332         catch(exception& e) {
 333                 m->errorOut(e, "LibShuffCommand", "execute");
 334                 exit(1);
 335         }
 336 }
 337
 338 //**********************************************************************************************************************
 339
 340 int LibShuffCommand::printCoverageFile() {
 341         try {
 342
 343                 ofstream outCov;
 344                 summaryFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "libshuff.coverage";
 345                 m->openOutputFile(summaryFile, outCov);
 346                 outputNames.push_back(summaryFile); outputTypes["coverage"].push_back(summaryFile);
 347                 outCov.setf(ios::fixed, ios::floatfield); outCov.setf(ios::showpoint);
 348                 //cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 349
 350                 map<double,vector<int> > allDistances;
 351                 map<double,vector<int> >::iterator it;
 352
 353                 vector<vector<int> > indices(numGroups);
 354                 int numIndices = numGroups * numGroups;
 355
 356                 int index = 0;
 357                 for(int i=0;i<numGroups;i++){
 358                         indices[i].assign(numGroups,0);
 359                         for(int j=0;j<numGroups;j++){
 360                                 indices[i][j] = index++;
 361
 362                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 363                                 int spotj = groupMap->groupIndex[groupNames[j]];
 364
 365                                 for(int k=0;k<savedMinValues[spoti][spotj].size();k++){
 366
 367                                         if(m->control_pressed)  { outCov.close(); return 0; }
 368
 369                                         if(allDistances[savedMinValues[spoti][spotj][k]].size() != 0){
 370                                                 allDistances[savedMinValues[spoti][spotj][k]][indices[i][j]]++;
 371                                         }
 372                                         else{
 373                                                 allDistances[savedMinValues[spoti][spotj][k]].assign(numIndices, 0);
 374                                                 allDistances[savedMinValues[spoti][spotj][k]][indices[i][j]] = 1;
 375                                         }
 376                                 }
 377                         }
 378                 }
 379                 it=allDistances.begin();
 380
 381                 //cout << setprecision(8);
 382
 383                 vector<int> prevRow = it->second;
 384                 it++;
 385
 386                 for(;it!=allDistances.end();it++){
 387                         for(int i=0;i<it->second.size();i++){
 388                                 it->second[i] += prevRow[i];
 389                         }
 390                         prevRow = it->second;
 391                 }
 392
 393                 vector<int> lastRow = allDistances.rbegin()->second;
 394                 outCov << setprecision(8);
 395
 396                 outCov << "dist";
 397                 for (int i = 0; i < numGroups; i++){
 398                         outCov << '\t' << groupNames[i];
 399                 }
 400                 for (int i=0;i<numGroups;i++){
 401                         for(int j=i+1;j<numGroups;j++){
 402                                 if(m->control_pressed)  { outCov.close(); return 0; }
 403                                 outCov << '\t' << groupNames[i] << '-' << groupNames[j] << '\t';
 404                                 outCov << groupNames[j] << '-' << groupNames[i];
 405                         }
 406                 }
 407                 outCov << endl;
 408
 409                 for(it=allDistances.begin();it!=allDistances.end();it++){
 410                         outCov << it->first << '\t';
 411                         for(int i=0;i<numGroups;i++){
 412                                 outCov << it->second[indices[i][i]]/(float)lastRow[indices[i][i]] << '\t';
 413                         }
 414                         for(int i=0;i<numGroups;i++){
 415                                 for(int j=i+1;j<numGroups;j++){
 416                                         if(m->control_pressed)  { outCov.close(); return 0; }
 417
 418                                         outCov << it->second[indices[i][j]]/(float)lastRow[indices[i][j]] << '\t';
 419                                         outCov << it->second[indices[j][i]]/(float)lastRow[indices[j][i]] << '\t';
 420                                 }
 421                         }
 422                         outCov << endl;
 423                 }
 424                 outCov.close();
 425
 426                 return 0;
 427         }
 428         catch(exception& e) {
 429                 m->errorOut(e, "LibShuffCommand", "printCoverageFile");
 430                 exit(1);
 431         }
 432 }
 433
 434 //**********************************************************************************************************************
 435
 436 int LibShuffCommand::printSummaryFile() {
 437         try {
 438
 439                 ofstream outSum;
 440                 summaryFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "libshuff.summary";
 441                 m->openOutputFile(summaryFile, outSum);
 442                 outputNames.push_back(summaryFile); outputTypes["libshuffsummary"].push_back(summaryFile);
 443
 444                 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
 445                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 446
 447                 cout << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl;
 448                 m->mothurOutJustToLog("Comparison\tdCXYScore\tSignificance"); m->mothurOutEndLine();
 449                 outSum << setw(20) << left << "Comparison" << '\t' << setprecision(8) << "dCXYScore" << '\t' << "Significance" << endl;
 450
 451                 int precision = (int)log10(iters);
 452                 for(int i=0;i<numGroups;i++){
 453                         for(int j=i+1;j<numGroups;j++){
 454                                 if(m->control_pressed)  { outSum.close(); return 0; }
 455
 456                                 int spoti = groupMap->groupIndex[groupNames[i]]; //neccessary in case user selects groups so you know where they are in the matrix
 457                                 int spotj = groupMap->groupIndex[groupNames[j]];
 458
 459                                 if(pValueCounts[i][j]){
 460                                         cout << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << setprecision(precision) << pValueCounts[i][j]/(float)iters << endl;
 461                                         m->mothurOutJustToLog(groupNames[i]+"-"+groupNames[j] + "\t" + toString(savedDXYValues[spoti][spotj]) + "\t" + toString((pValueCounts[i][j]/(float)iters))); m->mothurOutEndLine();
 462                                         outSum << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << setprecision(precision) << pValueCounts[i][j]/(float)iters << endl;
 463                                 }
 464                                 else{
 465                                         cout << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << '<' <<setprecision(precision) << 1/(float)iters << endl;
 466                                         m->mothurOutJustToLog(groupNames[i]+"-"+groupNames[j] + "\t" + toString(savedDXYValues[spoti][spotj]) + "\t" + toString((1/(float)iters))); m->mothurOutEndLine();
 467                                         outSum << setw(20) << left << groupNames[i]+'-'+groupNames[j] << '\t' << setprecision(8) << savedDXYValues[spoti][spotj] << '\t' << '<' <<setprecision(precision) << 1/(float)iters << endl;
 468                                 }
 469                                 if(pValueCounts[j][i]){
 470                                         cout << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << setprecision (precision) << pValueCounts[j][i]/(float)iters << endl;
 471                                         m->mothurOutJustToLog(groupNames[j]+"-"+groupNames[i] + "\t" + toString(savedDXYValues[spotj][spoti]) + "\t" + toString((pValueCounts[j][i]/(float)iters))); m->mothurOutEndLine();
 472                                         outSum << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << setprecision (precision) << pValueCounts[j][i]/(float)iters << endl;
 473                                 }
 474                                 else{
 475                                         cout << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << '<' <<setprecision (precision) << 1/(float)iters << endl;
 476                                         m->mothurOutJustToLog(groupNames[j]+"-"+groupNames[i] + "\t" + toString(savedDXYValues[spotj][spoti]) + "\t" + toString((1/(float)iters))); m->mothurOutEndLine();
 477                                         outSum << setw(20) << left << groupNames[j]+'-'+groupNames[i] << '\t' << setprecision(8) << savedDXYValues[spotj][spoti] << '\t' << '<' <<setprecision (precision) << 1/(float)iters << endl;
 478                                 }
 479                         }
 480                 }
 481
 482                 outSum.close();
 483                 return 0;
 484         }
 485         catch(exception& e) {
 486                 m->errorOut(e, "LibShuffCommand", "printSummaryFile");
 487                 exit(1);
 488         }
 489 }
 490
 491 //**********************************************************************************************************************
 492
 493 void LibShuffCommand::setGroups() {
 494         try {
 495                 vector<string> myGroups = m->getGroups();
 496                 //if the user has not entered specific groups to analyze then do them all
 497                 if (m->getNumGroups() == 0) {
 498                         numGroups = groupMap->getNumGroups();
 499                         for (int i=0; i < numGroups; i++) {
 500                                 myGroups.push_back((groupMap->getNamesOfGroups())[i]);
 501                         }
 502                 } else {
 503                         if (savegroups != "all") {
 504                                 //check that groups are valid
 505                                 for (int i = 0; i < myGroups.size(); i++) {
 506                                         if (groupMap->isValidGroup(myGroups[i]) != true) {
 507                                                 m->mothurOut(myGroups[i] + " is not a valid group, and will be disregarded."); m->mothurOutEndLine();
 508                                                 // erase the invalid group from globaldata->Groups
 509                                                 myGroups.erase(myGroups.begin()+i);
 510                                         }
 511                                 }
 512
 513                                 //if the user only entered invalid groups
 514                                 if ((myGroups.size() == 0) || (myGroups.size() == 1)) {
 515                                         numGroups = groupMap->getNumGroups();
 516                                         for (int i=0; i < numGroups; i++) {
 517                                                 myGroups.push_back((groupMap->getNamesOfGroups())[i]);
 518                                         }
 519                                         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();
 520                                 } else { numGroups = myGroups.size(); }
 521                         } else { //users wants all groups
 522                                 numGroups = groupMap->getNumGroups();
 523                                 myGroups.clear();
 524                                 for (int i=0; i < numGroups; i++) {
 525                                         myGroups.push_back((groupMap->getNamesOfGroups())[i]);
 526                                 }
 527                         }
 528                 }
 529
 530                 //sort so labels match
 531                 sort(myGroups.begin(), myGroups.end());
 532
 533                 //sort
 534                 //sort(groupMap->namesOfGroups.begin(), groupMap->namesOfGroups.end());
 535
 536                 for (int i = 0; i < (groupMap->getNamesOfGroups()).size(); i++) {  groupMap->groupIndex[(groupMap->getNamesOfGroups())[i]] = i;  }
 537
 538                 groupNames = myGroups;
 539                 m->setGroups(myGroups);
 540
 541         }
 542         catch(exception& e) {
 543                 m->errorOut(e, "LibShuffCommand", "setGroups");
 544                 exit(1);
 545         }
 546 }
 547
 548 /***********************************************************/