subsamplecommand.cpp

   1 /*
   2  *  subsamplecommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 10/27/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "subsamplecommand.h"
  11 #include "sharedutilities.h"
  12
  13 //**********************************************************************************************************************
  14 vector<string> SubSampleCommand::getValidParameters(){
  15         try {
  16                 string Array[] =  {"fasta", "group", "list","shared","rabund","persample", "name","sabund","size","groups","label","outputdir","inputdir"};
  17                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  18                 return myArray;
  19         }
  20         catch(exception& e) {
  21                 m->errorOut(e, "SubSampleCommand", "getValidParameters");
  22                 exit(1);
  23         }
  24 }
  25 //**********************************************************************************************************************
  26 SubSampleCommand::SubSampleCommand(){
  27         try {
  28                 abort = true;
  29                 //initialize outputTypes
  30                 vector<string> tempOutNames;
  31                 outputTypes["shared"] = tempOutNames;
  32                 outputTypes["list"] = tempOutNames;
  33                 outputTypes["rabund"] = tempOutNames;
  34                 outputTypes["sabund"] = tempOutNames;
  35                 outputTypes["fasta"] = tempOutNames;
  36                 outputTypes["name"] = tempOutNames;
  37                 outputTypes["group"] = tempOutNames;
  38         }
  39         catch(exception& e) {
  40                 m->errorOut(e, "SubSampleCommand", "GetRelAbundCommand");
  41                 exit(1);
  42         }
  43 }
  44 //**********************************************************************************************************************
  45 vector<string> SubSampleCommand::getRequiredParameters(){
  46         try {
  47                 string Array[] =  {"fasta","list","shared","rabund", "sabund","or"};
  48                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  49                 return myArray;
  50         }
  51         catch(exception& e) {
  52                 m->errorOut(e, "SubSampleCommand", "getRequiredParameters");
  53                 exit(1);
  54         }
  55 }
  56 //**********************************************************************************************************************
  57 vector<string> SubSampleCommand::getRequiredFiles(){
  58         try {
  59                 vector<string> myArray;
  60                 return myArray;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "SubSampleCommand", "getRequiredFiles");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 SubSampleCommand::SubSampleCommand(string option) {
  69         try {
  70                 globaldata = GlobalData::getInstance();
  71                 abort = false;
  72                 allLines = 1;
  73                 labels.clear();
  74
  75                 //allow user to run help
  76                 if(option == "help") { help(); abort = true; }
  77
  78                 else {
  79                         //valid paramters for this command
  80                         string Array[] =  {"fasta", "group", "list","shared","rabund","persample", "sabund","name","size","groups","label","outputdir","inputdir"};
  81                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  82
  83                         OptionParser parser(option);
  84                         map<string,string> parameters = parser.getParameters();
  85
  86                         ValidParameters validParameter;
  87
  88                         //check to make sure all parameters are valid for command
  89                         map<string,string>::iterator it;
  90                         for (it = parameters.begin(); it != parameters.end(); it++) {
  91                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  92                         }
  93
  94                         //initialize outputTypes
  95                         vector<string> tempOutNames;
  96                         outputTypes["shared"] = tempOutNames;
  97                         outputTypes["list"] = tempOutNames;
  98                         outputTypes["rabund"] = tempOutNames;
  99                         outputTypes["sabund"] = tempOutNames;
 100                         outputTypes["fasta"] = tempOutNames;
 101                         outputTypes["name"] = tempOutNames;
 102                         outputTypes["group"] = tempOutNames;
 103
 104                         //if the user changes the output directory command factory will send this info to us in the output parameter
 105                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 106
 107                         //if the user changes the input directory command factory will send this info to us in the output parameter
 108                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 109                         if (inputDir == "not found"){   inputDir = "";          }
 110                         else {
 111                                 string path;
 112                                 it = parameters.find("list");
 113                                 //user has given a template file
 114                                 if(it != parameters.end()){
 115                                         path = m->hasPath(it->second);
 116                                         //if the user has not given a path then, add inputdir. else leave path alone.
 117                                         if (path == "") {       parameters["list"] = inputDir + it->second;             }
 118                                 }
 119
 120                                 it = parameters.find("fasta");
 121                                 //user has given a template file
 122                                 if(it != parameters.end()){
 123                                         path = m->hasPath(it->second);
 124                                         //if the user has not given a path then, add inputdir. else leave path alone.
 125                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
 126                                 }
 127
 128                                 it = parameters.find("shared");
 129                                 //user has given a template file
 130                                 if(it != parameters.end()){
 131                                         path = m->hasPath(it->second);
 132                                         //if the user has not given a path then, add inputdir. else leave path alone.
 133                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 134                                 }
 135
 136                                 it = parameters.find("group");
 137                                 //user has given a template file
 138                                 if(it != parameters.end()){
 139                                         path = m->hasPath(it->second);
 140                                         //if the user has not given a path then, add inputdir. else leave path alone.
 141                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 142                                 }
 143
 144                                 it = parameters.find("sabund");
 145                                 //user has given a template file
 146                                 if(it != parameters.end()){
 147                                         path = m->hasPath(it->second);
 148                                         //if the user has not given a path then, add inputdir. else leave path alone.
 149                                         if (path == "") {       parameters["sabund"] = inputDir + it->second;           }
 150                                 }
 151
 152                                 it = parameters.find("rabund");
 153                                 //user has given a template file
 154                                 if(it != parameters.end()){
 155                                         path = m->hasPath(it->second);
 156                                         //if the user has not given a path then, add inputdir. else leave path alone.
 157                                         if (path == "") {       parameters["rabund"] = inputDir + it->second;           }
 158                                 }
 159
 160                                 it = parameters.find("name");
 161                                 //user has given a template file
 162                                 if(it != parameters.end()){
 163                                         path = m->hasPath(it->second);
 164                                         //if the user has not given a path then, add inputdir. else leave path alone.
 165                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 166                                 }
 167                         }
 168
 169                         //check for required parameters
 170                         listfile = validParameter.validFile(parameters, "list", true);
 171                         if (listfile == "not open") { listfile = ""; abort = true; }
 172                         else if (listfile == "not found") { listfile = ""; }
 173
 174                         sabundfile = validParameter.validFile(parameters, "sabund", true);
 175                         if (sabundfile == "not open") { sabundfile = ""; abort = true; }
 176                         else if (sabundfile == "not found") { sabundfile = ""; }
 177
 178                         rabundfile = validParameter.validFile(parameters, "rabund", true);
 179                         if (rabundfile == "not open") { rabundfile = ""; abort = true; }
 180                         else if (rabundfile == "not found") { rabundfile = ""; }
 181
 182                         fastafile = validParameter.validFile(parameters, "fasta", true);
 183                         if (fastafile == "not open") { fastafile = ""; abort = true; }
 184                         else if (fastafile == "not found") { fastafile = ""; }
 185
 186                         sharedfile = validParameter.validFile(parameters, "shared", true);
 187                         if (sharedfile == "not open") { sharedfile = ""; abort = true; }
 188                         else if (sharedfile == "not found") { sharedfile = ""; }
 189
 190                         namefile = validParameter.validFile(parameters, "name", true);
 191                         if (namefile == "not open") { namefile = ""; abort = true; }
 192                         else if (namefile == "not found") { namefile = ""; }
 193
 194                         groupfile = validParameter.validFile(parameters, "group", true);
 195                         if (groupfile == "not open") { groupfile = ""; abort = true; }
 196                         else if (groupfile == "not found") { groupfile = ""; }
 197
 198
 199                         //check for optional parameter and set defaults
 200                         // ...at some point should added some additional type checking...
 201                         label = validParameter.validFile(parameters, "label", false);
 202                         if (label == "not found") { label = ""; }
 203                         else {
 204                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 205                                 else { allLines = 1;  }
 206                         }
 207
 208                         groups = validParameter.validFile(parameters, "groups", false);
 209                         if (groups == "not found") { groups = ""; pickedGroups = false; }
 210                         else {
 211                                 pickedGroups = true;
 212                                 m->splitAtDash(groups, Groups);
 213                                 globaldata->Groups = Groups;
 214                         }
 215
 216                         string temp = validParameter.validFile(parameters, "size", false);              if (temp == "not found"){       temp = "0";             }
 217                         convert(temp, size);
 218
 219                         temp = validParameter.validFile(parameters, "persample", false);                if (temp == "not found"){       temp = "f";             }
 220                         persample = m->isTrue(temp);
 221
 222                         if (groupfile == "") { persample = false; }
 223
 224                         if ((namefile != "") && (fastafile == "")) { m->mothurOut("You may only use a namefile with a fastafile."); m->mothurOutEndLine(); abort = true; }
 225
 226                         if ((fastafile == "") && (listfile == "") && (sabundfile == "") && (rabundfile == "") && (sharedfile == "")) {
 227                                 m->mothurOut("You must provide a fasta, list, sabund, rabund or shared file as an input file."); m->mothurOutEndLine(); abort = true; }
 228
 229                         if (pickedGroups && ((groupfile == "") && (sharedfile == ""))) {
 230                                 m->mothurOut("You cannot pick groups without a valid group file or shared file."); m->mothurOutEndLine(); abort = true; }
 231
 232                         if ((groupfile != "") && ((fastafile == "") && (listfile == ""))) {
 233                                 m->mothurOut("Group file only valid with listfile or fastafile."); m->mothurOutEndLine(); abort = true; }
 234
 235                         if ((groupfile != "") && ((fastafile != "") && (listfile != ""))) {
 236                                 m->mothurOut("A new group file can only be made from the subsample of a listfile or fastafile, not both. Please correct."); m->mothurOutEndLine(); abort = true; }
 237
 238                 }
 239
 240         }
 241         catch(exception& e) {
 242                 m->errorOut(e, "SubSampleCommand", "SubSampleCommand");
 243                 exit(1);
 244         }
 245 }
 246
 247 //**********************************************************************************************************************
 248
 249 void SubSampleCommand::help(){
 250         try {
 251                 m->mothurOut("The sub.sample command is designed to be used as a way to normalize your data, or create a smaller set from your original set.\n");
 252                 m->mothurOut("The sub.sample command parameters are fasta, name, list, group, rabund, sabund, shared, groups, size, persample and label.  You must provide a fasta, list, sabund, rabund or shared file as an input file.\n");
 253                 m->mothurOut("The namefile is only used with the fasta file, not with the listfile, because the list file should contain all sequences.\n");
 254                 m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included. The group names are separated by dashes.\n");
 255                 m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n");
 256                 m->mothurOut("The size parameter allows you indicate the size of your subsample.\n");
 257                 m->mothurOut("The persample parameter allows you indicate you want to select subsample of the same size from each of your groups, default=false. It is only used with the list and fasta files if a groupfile is given.\n");
 258                 m->mothurOut("persample=false will select a random set of sequences of the size you select, but the number of seqs from each group may differ.\n");
 259                 m->mothurOut("The size parameter is not set: with shared file size=number of seqs in smallest sample, with all other files if a groupfile is given and persample=true, then size=number of seqs in smallest sample, otherwise size=10% of number of seqs.\n");
 260                 m->mothurOut("The sub.sample command should be in the following format: sub.sample(list=yourListFile, group=yourGroupFile, groups=yourGroups, label=yourLabels).\n");
 261                 m->mothurOut("Example sub.sample(list=abrecovery.fn.list, group=abrecovery.groups, groups=B-C, size=20).\n");
 262                 m->mothurOut("The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n");
 263                 m->mothurOut("The sub.sample command outputs a .subsample file.\n");
 264                 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
 265
 266         }
 267         catch(exception& e) {
 268                 m->errorOut(e, "SubSampleCommand", "help");
 269                 exit(1);
 270         }
 271 }
 272
 273 //**********************************************************************************************************************
 274
 275 SubSampleCommand::~SubSampleCommand(){}
 276
 277 //**********************************************************************************************************************
 278
 279 int SubSampleCommand::execute(){
 280         try {
 281
 282                 if (abort == true) { return 0; }
 283
 284                 if (sharedfile != "")   {   getSubSampleShared();       }
 285                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); return 0; } }
 286
 287                 if (listfile != "")             {   getSubSampleList();         }
 288                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); return 0; } }
 289
 290                 if (rabundfile != "")   {   getSubSampleRabund();       }
 291                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); return 0; } }
 292
 293                 if (sabundfile != "")   {   getSubSampleSabund();       }
 294                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); return 0; } }
 295
 296                 if (fastafile != "")    {   getSubSampleFasta();        }
 297                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); return 0; } }
 298
 299
 300                 m->mothurOutEndLine();
 301                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 302                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 303                 m->mothurOutEndLine();
 304
 305                 return 0;
 306         }
 307         catch(exception& e) {
 308                 m->errorOut(e, "SubSampleCommand", "execute");
 309                 exit(1);
 310         }
 311 }
 312 //**********************************************************************************************************************
 313 int SubSampleCommand::getSubSampleFasta() {
 314         try {
 315
 316                 if (namefile != "") { readNames(); }    //fills names with all names in namefile.
 317                 else { getNames(); }//no name file, so get list of names to pick from
 318
 319                 GroupMap* groupMap;
 320                 if (groupfile != "") {
 321
 322                         groupMap = new GroupMap(groupfile);
 323                         groupMap->readMap();
 324
 325                         //takes care of user setting groupNames that are invalid or setting groups=all
 326                         SharedUtil* util = new SharedUtil();
 327                         util->setGroups(Groups, groupMap->namesOfGroups);
 328                         delete util;
 329
 330                         //file mismatch quit
 331                         if (names.size() != groupMap->getNumSeqs()) {
 332                                 m->mothurOut("[ERROR]: your fasta file contains " + toString(names.size()) + " sequences, and your groupfile contains " + toString(groupMap->getNumSeqs()) + ", please correct.");
 333                                 m->mothurOutEndLine();
 334                                 delete groupMap;
 335                                 return 0;
 336                         }
 337                 }
 338
 339                 if (m->control_pressed) { return 0; }
 340
 341                 string thisOutputDir = outputDir;
 342                 if (outputDir == "") {  thisOutputDir += m->hasPath(fastafile);  }
 343                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(fastafile)) + "subsample" + m->getExtension(fastafile);
 344
 345                 ofstream out;
 346                 m->openOutputFile(outputFileName, out);
 347                 outputTypes["fasta"].push_back(outputFileName);  outputNames.push_back(outputFileName);
 348
 349                 //make sure that if your picked groups size is not too big
 350                 int thisSize = names.size();
 351                 if (persample) {
 352                         if (size == 0) { //user has not set size, set size = smallest samples size
 353                                 size = groupMap->getNumSeqs(Groups[0]);
 354                                 for (int i = 1; i < Groups.size(); i++) {
 355                                         int thisSize = groupMap->getNumSeqs(Groups[i]);
 356
 357                                         if (thisSize < size) {  size = thisSize;        }
 358                                 }
 359                         }else { //make sure size is not too large
 360                                 int smallestSize = groupMap->getNumSeqs(Groups[0]);
 361                                 for (int i = 1; i < Groups.size(); i++) {
 362                                         int thisSize = groupMap->getNumSeqs(Groups[i]);
 363
 364                                         if (thisSize < smallestSize) {  smallestSize = thisSize;        }
 365                                 }
 366                                 if (smallestSize < size) { size = smallestSize; m->mothurOut("You have selected a size that is larger than your smallest sample, using your samllest sample size, " + toString(smallestSize) + "."); m->mothurOutEndLine(); }
 367                         }
 368
 369                         m->mothurOut("Sampling " + toString(size) + " from each group."); m->mothurOutEndLine();
 370                 }else {
 371                         if (pickedGroups) {
 372                                 int total = 0;
 373                                 for(int i = 0; i < Groups.size(); i++) {
 374                                         total += groupMap->getNumSeqs(Groups[i]);
 375                                 }
 376
 377                                 if (size == 0) { //user has not set size, set size = 10% samples size
 378                                         size = int (total * 0.10);
 379                                 }
 380
 381                                 if (total < size) {
 382                                         if (size != 0) {
 383                                                 m->mothurOut("Your size is too large for the number of groups you selected. Adjusting to " + toString(int (total * 0.10)) + "."); m->mothurOutEndLine();
 384                                         }
 385                                         size = int (total * 0.10);
 386                                 }
 387
 388                                 m->mothurOut("Sampling " + toString(size) + " from " + toString(total) + "."); m->mothurOutEndLine();
 389                         }
 390
 391                         if (size == 0) { //user has not set size, set size = 10% samples size
 392                                 size = int (names.size() * 0.10);
 393                         }
 394
 395                         if (size > thisSize) { m->mothurOut("Your fasta file only contains " + toString(thisSize) + " sequences. Setting size to " + toString(thisSize) + "."); m->mothurOutEndLine();
 396                                 size = thisSize;
 397                         }
 398
 399                         if (!pickedGroups) { m->mothurOut("Sampling " + toString(size) + " from " + toString(thisSize) + "."); m->mothurOutEndLine(); }
 400
 401                 }
 402                 random_shuffle(names.begin(), names.end());
 403
 404                 set<string> subset; //dont want repeat sequence names added
 405                 if (persample) {
 406                         for (int i = 0; i < Groups.size(); i++) {
 407
 408                                 //randomly select a subset of those names from this group to include in the subsample
 409                                 for (int j = 0; j < size; j++) {
 410
 411                                         if (m->control_pressed) { return 0; }
 412
 413                                         //get random sequence to add, making sure we have not already added it
 414                                         bool done = false;
 415                                         int myrand;
 416                                         while (!done) {
 417                                                 myrand = (int)((float)(rand()) / (RAND_MAX / (thisSize-1) + 1));
 418
 419                                                 if (subset.count(names[myrand]) == 0)  {
 420
 421                                                         string group = groupMap->getGroup(names[myrand]);
 422                                                         if (group == "not found") { m->mothurOut("[ERROR]: " + names[myrand] + " is not in your groupfile. please correct."); m->mothurOutEndLine(); group = "NOTFOUND"; }
 423
 424                                                         if (group == Groups[i]) {       subset.insert(names[myrand]); break;    }
 425                                                 }
 426                                         }
 427                                 }
 428                         }
 429                 }else {
 430                         //randomly select a subset of those names to include in the subsample
 431                         for (int j = 0; j < size; j++) {
 432
 433                                 if (m->control_pressed) { return 0; }
 434
 435                                 //get random sequence to add, making sure we have not already added it
 436                                 bool done = false;
 437                                 int myrand;
 438                                 while (!done) {
 439                                         myrand = (int)((float)(rand()) / (RAND_MAX / (thisSize-1) + 1));
 440
 441                                         if (subset.count(names[myrand]) == 0)  {
 442
 443                                                 if (groupfile != "") { //if there is a groupfile given fill in group info
 444                                                         string group = groupMap->getGroup(names[myrand]);
 445                                                         if (group == "not found") { m->mothurOut("[ERROR]: " + names[myrand] + " is not in your groupfile. please correct."); m->mothurOutEndLine(); group = "NOTFOUND"; }
 446
 447                                                         if (pickedGroups) { //if hte user picked groups, we only want to keep the names of sequences from those groups
 448                                                                 if (m->inUsersGroups(group, Groups)) {
 449                                                                         subset.insert(names[myrand]); break;
 450                                                                 }
 451                                                         }else{
 452                                                                 subset.insert(names[myrand]); break;
 453                                                         }
 454                                                 }else{ //save everyone, group
 455                                                         subset.insert(names[myrand]); break;
 456                                                 }
 457                                         }
 458                                 }
 459                         }
 460                 }
 461                 //read through fasta file outputting only the names on the subsample list
 462                 ifstream in;
 463                 m->openInputFile(fastafile, in);
 464
 465                 string thisname;
 466                 int count = 0;
 467                 map<string, vector<string> >::iterator itNameMap;
 468
 469                 while(!in.eof()){
 470
 471                         if (m->control_pressed) { in.close(); out.close();  return 0; }
 472
 473                         Sequence currSeq(in);
 474                         thisname = currSeq.getName();
 475
 476                         if (thisname != "") {
 477
 478                                 //does the subset contain a sequence that this sequence represents
 479                                 itNameMap = nameMap.find(thisname);
 480                                 if (itNameMap != nameMap.end()) {
 481                                         vector<string> nameRepresents = itNameMap->second;
 482
 483                                         for (int i = 0; i < nameRepresents.size(); i++){
 484                                                 if (subset.count(nameRepresents[i]) != 0) {
 485                                                         out << ">" << nameRepresents[i] << endl << currSeq.getAligned() << endl;
 486                                                         count++;
 487                                                 }
 488                                         }
 489                                 }else{
 490                                         m->mothurOut("[ERROR]: " + thisname + " is not in your namefile, please correct."); m->mothurOutEndLine();
 491                                 }
 492                         }
 493                         m->gobble(in);
 494                 }
 495                 in.close();
 496                 out.close();
 497
 498                 if (count != subset.size()) {
 499                         m->mothurOut("[ERROR]: The subset selected contained " + toString(subset.size()) + " sequences, but I only found " + toString(count) + " of those in the fastafile."); m->mothurOutEndLine();
 500                 }
 501
 502                 //if a groupfile is provided read through the group file only outputting the names on the subsample list
 503                 if (groupfile != "") {
 504
 505                         string groupOutputDir = outputDir;
 506                         if (outputDir == "") {  groupOutputDir += m->hasPath(groupfile);  }
 507                         string groupOutputFileName = groupOutputDir + m->getRootName(m->getSimpleName(groupfile)) + "subsample" + m->getExtension(groupfile);
 508
 509                         ofstream outGroup;
 510                         m->openOutputFile(groupOutputFileName, outGroup);
 511                         outputTypes["group"].push_back(groupOutputFileName);  outputNames.push_back(groupOutputFileName);
 512
 513                         ifstream inGroup;
 514                         m->openInputFile(groupfile, inGroup);
 515                         string name, group;
 516
 517                         while(!inGroup.eof()){
 518
 519                                 if (m->control_pressed) { inGroup.close(); outGroup.close(); return 0; }
 520
 521                                 inGroup >> name;        m->gobble(inGroup);                     //read from first column
 522                                 inGroup >> group;                       //read from second column
 523
 524                                 //if this name is in the accnos file
 525                                 if (subset.count(name) != 0) {
 526                                         outGroup << name << '\t' << group << endl;
 527                                         subset.erase(name);
 528                                 }
 529
 530                                 m->gobble(inGroup);
 531                         }
 532                         inGroup.close();
 533                         outGroup.close();
 534
 535                         //sanity check
 536                         if (subset.size() != 0) {
 537                                 m->mothurOut("Your groupfile does not match your fasta file."); m->mothurOutEndLine();
 538                                 for (set<string>::iterator it = subset.begin(); it != subset.end(); it++) {
 539                                         m->mothurOut("[ERROR]: " + *it + " is missing from your groupfile."); m->mothurOutEndLine();
 540                                 }
 541                         }
 542                 }
 543
 544                 return 0;
 545
 546         }
 547         catch(exception& e) {
 548                 m->errorOut(e, "SubSampleCommand", "getSubSampleFasta");
 549                 exit(1);
 550         }
 551 }
 552 //**********************************************************************************************************************
 553 int SubSampleCommand::getNames() {
 554         try {
 555
 556                 ifstream in;
 557                 m->openInputFile(fastafile, in);
 558
 559                 string thisname;
 560                 while(!in.eof()){
 561
 562                         if (m->control_pressed) { in.close(); return 0; }
 563
 564                         Sequence currSeq(in);
 565                         thisname = currSeq.getName();
 566
 567                         if (thisname != "") {
 568                                 vector<string> temp; temp.push_back(thisname);
 569                                 nameMap[thisname] = temp;
 570                                 names.push_back(thisname);
 571                         }
 572                         m->gobble(in);
 573                 }
 574                 in.close();
 575
 576                 return 0;
 577
 578         }
 579         catch(exception& e) {
 580                 m->errorOut(e, "SubSampleCommand", "getNames");
 581                 exit(1);
 582         }
 583 }
 584 //**********************************************************************************************************************
 585 int SubSampleCommand::readNames() {
 586         try {
 587
 588                 ifstream in;
 589                 m->openInputFile(namefile, in);
 590
 591                 string thisname, repnames;
 592                 map<string, vector<string> >::iterator it;
 593
 594                 while(!in.eof()){
 595
 596                         if (m->control_pressed) { in.close(); return 0; }
 597
 598                         in >> thisname;         m->gobble(in);          //read from first column
 599                         in >> repnames;                 //read from second column
 600
 601                         it = nameMap.find(thisname);
 602                         if (it == nameMap.end()) {
 603
 604                                 vector<string> splitRepNames;
 605                                 m->splitAtComma(repnames, splitRepNames);
 606
 607                                 nameMap[thisname] = splitRepNames;
 608                                 for (int i = 0; i < splitRepNames.size(); i++) { names.push_back(splitRepNames[i]); }
 609
 610                         }else{  m->mothurOut(thisname + " is already in namesfile. I will use first definition."); m->mothurOutEndLine();  }
 611
 612                         m->gobble(in);
 613                 }
 614                 in.close();
 615
 616                 return 0;
 617
 618         }
 619         catch(exception& e) {
 620                 m->errorOut(e, "SubSampleCommand", "readNames");
 621                 exit(1);
 622         }
 623 }
 624 //**********************************************************************************************************************
 625 int SubSampleCommand::getSubSampleShared() {
 626         try {
 627
 628                 string thisOutputDir = outputDir;
 629                 if (outputDir == "") {  thisOutputDir += m->hasPath(sharedfile);  }
 630                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(sharedfile)) + "subsample" + m->getExtension(sharedfile);
 631
 632                 ofstream out;
 633                 m->openOutputFile(outputFileName, out);
 634                 outputTypes["shared"].push_back(outputFileName);  outputNames.push_back(outputFileName);
 635
 636                 InputData* input = new InputData(sharedfile, "sharedfile");
 637                 vector<SharedRAbundVector*> lookup = input->getSharedRAbundVectors();
 638                 string lastLabel = lookup[0]->getLabel();
 639
 640                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 641                 set<string> processedLabels;
 642                 set<string> userLabels = labels;
 643
 644                 if (size == 0) { //user has not set size, set size = smallest samples size
 645                         size = lookup[0]->getNumSeqs();
 646                         for (int i = 1; i < lookup.size(); i++) {
 647                                 int thisSize = lookup[i]->getNumSeqs();
 648
 649                                 if (thisSize < size) {  size = thisSize;        }
 650                         }
 651                 }
 652
 653                 m->mothurOut("Sampling " + toString(size) + " from " + toString(lookup[0]->getNumSeqs()) + "."); m->mothurOutEndLine();
 654
 655                 //as long as you are not at the end of the file or done wih the lines you want
 656                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 657                         if (m->control_pressed) {  delete input; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; } out.close(); return 0;  }
 658
 659                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 660
 661                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 662
 663                                 processShared(lookup, out);
 664
 665                                 processedLabels.insert(lookup[0]->getLabel());
 666                                 userLabels.erase(lookup[0]->getLabel());
 667                         }
 668
 669                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 670                                 string saveLabel = lookup[0]->getLabel();
 671
 672                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 673
 674                                 lookup = input->getSharedRAbundVectors(lastLabel);
 675                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 676
 677                                 processShared(lookup, out);
 678
 679                                 processedLabels.insert(lookup[0]->getLabel());
 680                                 userLabels.erase(lookup[0]->getLabel());
 681
 682                                 //restore real lastlabel to save below
 683                                 lookup[0]->setLabel(saveLabel);
 684                         }
 685
 686                         lastLabel = lookup[0]->getLabel();
 687                         //prevent memory leak
 688                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 689
 690                         //get next line to process
 691                         lookup = input->getSharedRAbundVectors();
 692                 }
 693
 694
 695                 if (m->control_pressed) {  out.close(); return 0;  }
 696
 697                 //output error messages about any remaining user labels
 698                 set<string>::iterator it;
 699                 bool needToRun = false;
 700                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 701                         m->mothurOut("Your file does not include the label " + *it);
 702                         if (processedLabels.count(lastLabel) != 1) {
 703                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 704                                 needToRun = true;
 705                         }else {
 706                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 707                         }
 708                 }
 709
 710                 //run last label if you need to
 711                 if (needToRun == true)  {
 712                         for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 713                         lookup = input->getSharedRAbundVectors(lastLabel);
 714
 715                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 716
 717                         processShared(lookup, out);
 718
 719                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 720                 }
 721
 722                 delete input;
 723                 out.close();
 724
 725                 return 0;
 726
 727         }
 728         catch(exception& e) {
 729                 m->errorOut(e, "SubSampleCommand", "getSubSampleShared");
 730                 exit(1);
 731         }
 732 }
 733 //**********************************************************************************************************************
 734 int SubSampleCommand::processShared(vector<SharedRAbundVector*>& thislookup, ofstream& out) {
 735         try {
 736
 737                 int numBins = thislookup[0]->getNumBins();
 738                 for (int i = 0; i < thislookup.size(); i++) {
 739                         int thisSize = thislookup[i]->getNumSeqs();
 740
 741                         if (thisSize != size) {
 742
 743                                 string thisgroup = thislookup[i]->getGroup();
 744
 745                                 OrderVector* order = new OrderVector();
 746                                 for(int p=0;p<numBins;p++){
 747                                         for(int j=0;j<thislookup[i]->getAbundance(p);j++){
 748                                                 order->push_back(p);
 749                                         }
 750                                 }
 751                                 random_shuffle(order->begin(), order->end());
 752
 753                                 SharedRAbundVector* temp = new SharedRAbundVector(numBins);
 754                                 temp->setLabel(thislookup[i]->getLabel());
 755                                 temp->setGroup(thislookup[i]->getGroup());
 756
 757                                 delete thislookup[i];
 758                                 thislookup[i] = temp;
 759
 760
 761                                 for (int j = 0; j < size; j++) {
 762
 763                                         if (m->control_pressed) { delete order; return 0; }
 764
 765                                         //get random number to sample from order between 0 and thisSize-1.
 766                                         int myrand = (int)((float)(rand()) / (RAND_MAX / (thisSize-1) + 1));
 767
 768                                         int bin = order->get(myrand);
 769
 770                                         int abund = thislookup[i]->getAbundance(bin);
 771                                         thislookup[i]->set(bin, (abund+1), thisgroup);
 772                                 }
 773                                 delete order;
 774                         }
 775                 }
 776
 777                 //subsampling may have created some otus with no sequences in them
 778                 eliminateZeroOTUS(thislookup);
 779
 780                 if (m->control_pressed) { return 0; }
 781
 782                 for (int i = 0; i < thislookup.size(); i++) {
 783                         out << thislookup[i]->getLabel() << '\t' << thislookup[i]->getGroup() << '\t';
 784                         thislookup[i]->print(out);
 785                 }
 786
 787                 return 0;
 788
 789         }
 790         catch(exception& e) {
 791                 m->errorOut(e, "SubSampleCommand", "processShared");
 792                 exit(1);
 793         }
 794 }
 795 //**********************************************************************************************************************
 796 int SubSampleCommand::getSubSampleList() {
 797         try {
 798
 799                 string thisOutputDir = outputDir;
 800                 if (outputDir == "") {  thisOutputDir += m->hasPath(listfile);  }
 801                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(listfile)) + "subsample" + m->getExtension(listfile);
 802
 803                 ofstream out;
 804                 m->openOutputFile(outputFileName, out);
 805                 outputTypes["list"].push_back(outputFileName);  outputNames.push_back(outputFileName);
 806
 807                 InputData* input = new InputData(listfile, "list");
 808                 ListVector* list = input->getListVector();
 809                 string lastLabel = list->getLabel();
 810
 811                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 812                 set<string> processedLabels;
 813                 set<string> userLabels = labels;
 814
 815                 ofstream outGroup;
 816                 GroupMap* groupMap;
 817                 if (groupfile != "") {
 818
 819                         groupMap = new GroupMap(groupfile);
 820                         groupMap->readMap();
 821
 822                         //takes care of user setting groupNames that are invalid or setting groups=all
 823                         SharedUtil* util = new SharedUtil();
 824                         util->setGroups(Groups, groupMap->namesOfGroups);
 825                         delete util;
 826
 827                         //create outputfiles
 828                         string groupOutputDir = outputDir;
 829                         if (outputDir == "") {  groupOutputDir += m->hasPath(groupfile);  }
 830                         string groupOutputFileName = groupOutputDir + m->getRootName(m->getSimpleName(groupfile)) + "subsample" + m->getExtension(groupfile);
 831
 832                         m->openOutputFile(groupOutputFileName, outGroup);
 833                         outputTypes["group"].push_back(groupOutputFileName);  outputNames.push_back(groupOutputFileName);
 834
 835                         //file mismatch quit
 836                         if (list->getNumSeqs() != groupMap->getNumSeqs()) {
 837                                 m->mothurOut("[ERROR]: your list file contains " + toString(list->getNumSeqs()) + " sequences, and your groupfile contains " + toString(groupMap->getNumSeqs()) + ", please correct.");
 838                                 m->mothurOutEndLine();
 839                                 delete groupMap;
 840                                 delete list;
 841                                 delete input;
 842                                 out.close();
 843                                 outGroup.close();
 844                                 return 0;
 845                         }
 846                 }
 847
 848                 //make sure that if your picked groups size is not too big
 849                 if (persample) {
 850                         if (size == 0) { //user has not set size, set size = smallest samples size
 851                                 size = groupMap->getNumSeqs(Groups[0]);
 852                                 for (int i = 1; i < Groups.size(); i++) {
 853                                         int thisSize = groupMap->getNumSeqs(Groups[i]);
 854
 855                                         if (thisSize < size) {  size = thisSize;        }
 856                                 }
 857                         }else { //make sure size is not too large
 858                                 int smallestSize = groupMap->getNumSeqs(Groups[0]);
 859                                 for (int i = 1; i < Groups.size(); i++) {
 860                                         int thisSize = groupMap->getNumSeqs(Groups[i]);
 861
 862                                         if (thisSize < smallestSize) {  smallestSize = thisSize;        }
 863                                 }
 864                                 if (smallestSize < size) { size = smallestSize; m->mothurOut("You have selected a size that is larger than your smallest sample, using your samllest sample size, " + toString(smallestSize) + "."); m->mothurOutEndLine(); }
 865                         }
 866
 867                         m->mothurOut("Sampling " + toString(size) + " from each group."); m->mothurOutEndLine();
 868                 }else{
 869                         if (pickedGroups) {
 870                                 int total = 0;
 871                                 for(int i = 0; i < Groups.size(); i++) {
 872                                         total += groupMap->getNumSeqs(Groups[i]);
 873                                 }
 874
 875                                 if (size == 0) { //user has not set size, set size = 10% samples size
 876                                         size = int (total * 0.10);
 877                                 }
 878
 879                                 if (total < size) {
 880                                         m->mothurOut("Your size is too large for the number of groups you selected. Adjusting to " + toString(int (total * 0.10)) + "."); m->mothurOutEndLine();
 881                                         size = int (total * 0.10);
 882                                 }
 883
 884                                 m->mothurOut("Sampling " + toString(size) + " from " + toString(total) + "."); m->mothurOutEndLine();
 885                         }else{
 886
 887                                 if (size == 0) { //user has not set size, set size = 10% samples size
 888                                         size = int (list->getNumSeqs() * 0.10);
 889                                 }
 890
 891                                 int thisSize = list->getNumSeqs();
 892                                 if (size > thisSize) { m->mothurOut("Your list file only contains " + toString(thisSize) + " sequences. Setting size to " + toString(thisSize) + "."); m->mothurOutEndLine();
 893                                         size = thisSize;
 894                                 }
 895
 896                                 m->mothurOut("Sampling " + toString(size) + " from " + toString(list->getNumSeqs()) + "."); m->mothurOutEndLine();
 897                         }
 898                 }
 899
 900
 901                 //fill names
 902                 for (int i = 0; i < list->getNumBins(); i++) {
 903                         string binnames = list->get(i);
 904
 905                         //parse names
 906                         string individual = "";
 907                         int length = binnames.length();
 908                         for(int j=0;j<length;j++){
 909                                 if(binnames[j] == ','){
 910
 911                                         if (groupfile != "") { //if there is a groupfile given fill in group info
 912                                                 string group = groupMap->getGroup(individual);
 913                                                 if (group == "not found") { m->mothurOut("[ERROR]: " + individual + " is not in your groupfile. please correct."); m->mothurOutEndLine(); group = "NOTFOUND"; }
 914
 915                                                 if (pickedGroups) { //if hte user picked groups, we only want to keep the names of sequences from those groups
 916                                                         if (m->inUsersGroups(group, Groups)) {
 917                                                                 names.push_back(individual);
 918                                                         }
 919                                                 }else{
 920                                                         names.push_back(individual);
 921                                                 }
 922                                         }else{ //save everyone, group
 923                                                 names.push_back(individual);
 924                                         }
 925                                         individual = "";
 926                                 }
 927                                 else{
 928                                         individual += binnames[j];
 929                                 }
 930                         }
 931                         //save last name
 932                         if (groupfile != "") { //if there is a groupfile given fill in group info
 933                                 string group = groupMap->getGroup(individual);
 934                                 if (group == "not found") { m->mothurOut("[ERROR]: " + individual + " is not in your groupfile. please correct."); m->mothurOutEndLine(); group = "NOTFOUND"; }
 935
 936                                 if (pickedGroups) { //if hte user picked groups, we only want to keep the names of sequences from those groups
 937                                         if (m->inUsersGroups(group, Groups)) {
 938                                                 names.push_back(individual);
 939                                         }
 940                                 }else{
 941                                         names.push_back(individual);
 942                                 }
 943                         }else{ //save everyone, group
 944                                 names.push_back(individual);
 945                         }
 946                 }
 947
 948                 random_shuffle(names.begin(), names.end());
 949
 950                 //randomly select a subset of those names to include in the subsample
 951                 set<string> subset; //dont want repeat sequence names added
 952                 if (persample) {
 953                         for (int i = 0; i < Groups.size(); i++) {
 954
 955                                 for (int j = 0; j < size; j++) {
 956
 957                                         if (m->control_pressed) { break; }
 958
 959                                         //get random sequence to add, making sure we have not already added it
 960                                         bool done = false;
 961                                         int myrand;
 962                                         while (!done) {
 963                                                 myrand = (int)((float)(rand()) / (RAND_MAX / (names.size()-1) + 1));
 964
 965                                                 if (subset.count(names[myrand]) == 0) { //you are not already added
 966                                                         if (groupMap->getGroup(names[myrand]) == Groups[i])  { subset.insert(names[myrand]); break;     }
 967                                                 }
 968                                         }
 969                                 }
 970                         }
 971                 }else{
 972                         for (int j = 0; j < size; j++) {
 973
 974                                 if (m->control_pressed) { break; }
 975
 976                                 //get random sequence to add, making sure we have not already added it
 977                                 bool done = false;
 978                                 int myrand;
 979                                 while (!done) {
 980                                         myrand = (int)((float)(rand()) / (RAND_MAX / (names.size()-1) + 1));
 981
 982                                         if (subset.count(names[myrand]) == 0)  { subset.insert(names[myrand]); break;   }
 983                                 }
 984                         }
 985                 }
 986
 987                 if (groupfile != "") {
 988                         //write out new groupfile
 989                         for (set<string>::iterator it = subset.begin(); it != subset.end(); it++) {
 990                                 string group = groupMap->getGroup(*it);
 991                                 if (group == "not found") { group = "NOTFOUND"; }
 992
 993                                 outGroup << *it << '\t' << group << endl;
 994                         }
 995                         outGroup.close(); delete groupMap;
 996                 }
 997
 998
 999                 //as long as you are not at the end of the file or done wih the lines you want
1000                 while((list != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
1001
1002                         if (m->control_pressed) {  delete list; delete input; out.close();  return 0;  }
1003
1004                         if(allLines == 1 || labels.count(list->getLabel()) == 1){
1005
1006                                 m->mothurOut(list->getLabel()); m->mothurOutEndLine();
1007
1008                                 processList(list, out, subset);
1009
1010                                 processedLabels.insert(list->getLabel());
1011                                 userLabels.erase(list->getLabel());
1012                         }
1013
1014                         if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1015                                 string saveLabel = list->getLabel();
1016
1017                                 delete list;
1018
1019                                 list = input->getListVector(lastLabel);
1020                                 m->mothurOut(list->getLabel()); m->mothurOutEndLine();
1021
1022                                 processList(list, out, subset);
1023
1024                                 processedLabels.insert(list->getLabel());
1025                                 userLabels.erase(list->getLabel());
1026
1027                                 //restore real lastlabel to save below
1028                                 list->setLabel(saveLabel);
1029                         }
1030
1031                         lastLabel = list->getLabel();
1032
1033                         delete list; list = NULL;
1034
1035                         //get next line to process
1036                         list = input->getListVector();
1037                 }
1038
1039
1040                 if (m->control_pressed) {  if (list != NULL) { delete list; } delete input; out.close(); return 0;  }
1041
1042                 //output error messages about any remaining user labels
1043                 set<string>::iterator it;
1044                 bool needToRun = false;
1045                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1046                         m->mothurOut("Your file does not include the label " + *it);
1047                         if (processedLabels.count(lastLabel) != 1) {
1048                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1049                                 needToRun = true;
1050                         }else {
1051                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1052                         }
1053                 }
1054
1055                 //run last label if you need to
1056                 if (needToRun == true)  {
1057                         if (list != NULL) { delete list; }
1058
1059                         list = input->getListVector(lastLabel);
1060
1061                         m->mothurOut(list->getLabel()); m->mothurOutEndLine();
1062
1063                         processList(list, out, subset);
1064
1065                         delete list; list = NULL;
1066                 }
1067
1068                 out.close();
1069                 if (list != NULL) { delete list; }
1070                 delete input;
1071
1072                 return 0;
1073
1074         }
1075         catch(exception& e) {
1076                 m->errorOut(e, "SubSampleCommand", "getSubSampleList");
1077                 exit(1);
1078         }
1079 }
1080 //**********************************************************************************************************************
1081 int SubSampleCommand::processList(ListVector*& list, ofstream& out, set<string>& subset) {
1082         try {
1083
1084                 int numBins = list->getNumBins();
1085
1086                 ListVector* temp = new ListVector();
1087                 temp->setLabel(list->getLabel());
1088
1089                 for (int i = 0; i < numBins; i++) {
1090
1091                         if (m->control_pressed) { break; }
1092
1093                         string binnames = list->get(i);
1094
1095                         //parse names
1096                         string individual = "";
1097                         string newNames = "";
1098                         int length = binnames.length();
1099                         for(int j=0;j<length;j++){
1100                                 if(binnames[j] == ','){
1101                                         if (subset.count(individual) != 0) {  newNames += individual + ",";  }
1102                                         individual = "";
1103                                 }else{
1104                                         individual += binnames[j];
1105                                 }
1106                         }
1107                         if (subset.count(individual) != 0) {  newNames += individual;  }
1108
1109
1110                         //if there are names in this bin add to new list
1111                         if (newNames != "") {
1112                                 newNames = newNames.substr(0, newNames.length()-1); //rip off extra comma
1113                                 temp->push_back(newNames);
1114                         }
1115                 }
1116
1117                 delete list;
1118                 list = temp;
1119
1120                 if (m->control_pressed) { return 0; }
1121
1122                 list->print(out);
1123
1124                 return 0;
1125
1126         }
1127         catch(exception& e) {
1128                 m->errorOut(e, "SubSampleCommand", "processList");
1129                 exit(1);
1130         }
1131 }
1132 //**********************************************************************************************************************
1133 int SubSampleCommand::getSubSampleRabund() {
1134         try {
1135
1136                 string thisOutputDir = outputDir;
1137                 if (outputDir == "") {  thisOutputDir += m->hasPath(rabundfile);  }
1138                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(rabundfile)) + "subsample" + m->getExtension(rabundfile);
1139
1140                 ofstream out;
1141                 m->openOutputFile(outputFileName, out);
1142                 outputTypes["rabund"].push_back(outputFileName);  outputNames.push_back(outputFileName);
1143
1144                 InputData* input = new InputData(rabundfile, "rabund");
1145                 RAbundVector* rabund = input->getRAbundVector();
1146                 string lastLabel = rabund->getLabel();
1147
1148                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1149                 set<string> processedLabels;
1150                 set<string> userLabels = labels;
1151
1152                 if (size == 0) { //user has not set size, set size = 10%
1153                         size = int((rabund->getNumSeqs()) * 0.10);
1154                 }
1155
1156                 m->mothurOut("Sampling " + toString(size) + " from " + toString(rabund->getNumSeqs()) + "."); m->mothurOutEndLine();
1157
1158                 //as long as you are not at the end of the file or done wih the lines you want
1159                 while((rabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
1160                         if (m->control_pressed) {  delete input; delete rabund; out.close(); return 0;  }
1161
1162                         if(allLines == 1 || labels.count(rabund->getLabel()) == 1){
1163
1164                                 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
1165
1166                                 processRabund(rabund, out);
1167
1168                                 processedLabels.insert(rabund->getLabel());
1169                                 userLabels.erase(rabund->getLabel());
1170                         }
1171
1172                         if ((m->anyLabelsToProcess(rabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1173                                 string saveLabel = rabund->getLabel();
1174
1175                                 delete rabund;
1176
1177                                 rabund = input->getRAbundVector(lastLabel);
1178                                 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
1179
1180                                 processRabund(rabund, out);
1181
1182                                 processedLabels.insert(rabund->getLabel());
1183                                 userLabels.erase(rabund->getLabel());
1184
1185                                 //restore real lastlabel to save below
1186                                 rabund->setLabel(saveLabel);
1187                         }
1188
1189                         lastLabel = rabund->getLabel();
1190
1191                         //prevent memory leak
1192                         delete rabund; rabund = NULL;
1193
1194                         //get next line to process
1195                         rabund = input->getRAbundVector();
1196                 }
1197
1198
1199                 if (m->control_pressed) {  out.close(); return 0;  }
1200
1201                 //output error messages about any remaining user labels
1202                 set<string>::iterator it;
1203                 bool needToRun = false;
1204                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1205                         m->mothurOut("Your file does not include the label " + *it);
1206                         if (processedLabels.count(lastLabel) != 1) {
1207                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1208                                 needToRun = true;
1209                         }else {
1210                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1211                         }
1212                 }
1213
1214                 //run last label if you need to
1215                 if (needToRun == true)  {
1216                         if (rabund != NULL) { delete rabund; }
1217
1218                         rabund = input->getRAbundVector(lastLabel);
1219
1220                         m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
1221
1222                         processRabund(rabund, out);
1223
1224                         delete rabund;
1225                 }
1226
1227                 delete input;
1228                 out.close();
1229
1230                 return 0;
1231
1232         }
1233         catch(exception& e) {
1234                 m->errorOut(e, "SubSampleCommand", "getSubSampleRabund");
1235                 exit(1);
1236         }
1237 }
1238 //**********************************************************************************************************************
1239 int SubSampleCommand::processRabund(RAbundVector*& rabund, ofstream& out) {
1240         try {
1241
1242                 int numBins = rabund->getNumBins();
1243                 int thisSize = rabund->getNumSeqs();
1244
1245                 if (thisSize != size) {
1246
1247                         OrderVector* order = new OrderVector();
1248                         for(int p=0;p<numBins;p++){
1249                                 for(int j=0;j<rabund->get(p);j++){
1250                                         order->push_back(p);
1251                                 }
1252                         }
1253                         random_shuffle(order->begin(), order->end());
1254
1255                         RAbundVector* temp = new RAbundVector(numBins);
1256                         temp->setLabel(rabund->getLabel());
1257
1258                         delete rabund;
1259                         rabund = temp;
1260
1261                         for (int j = 0; j < size; j++) {
1262
1263                                 if (m->control_pressed) { delete order; return 0; }
1264
1265                                 //get random number to sample from order between 0 and thisSize-1.
1266                                 int myrand = (int)((float)(rand()) / (RAND_MAX / (thisSize-1) + 1));
1267
1268                                 int bin = order->get(myrand);
1269
1270                                 int abund = rabund->get(bin);
1271                                 rabund->set(bin, (abund+1));
1272                         }
1273
1274                         delete order;
1275                 }
1276
1277                 if (m->control_pressed) { return 0; }
1278
1279                 rabund->print(out);
1280
1281                 return 0;
1282
1283         }
1284         catch(exception& e) {
1285                 m->errorOut(e, "SubSampleCommand", "processRabund");
1286                 exit(1);
1287         }
1288 }
1289 //**********************************************************************************************************************
1290 int SubSampleCommand::getSubSampleSabund() {
1291         try {
1292
1293                 string thisOutputDir = outputDir;
1294                 if (outputDir == "") {  thisOutputDir += m->hasPath(sabundfile);  }
1295                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(sabundfile)) + "subsample" + m->getExtension(sabundfile);
1296
1297                 ofstream out;
1298                 m->openOutputFile(outputFileName, out);
1299                 outputTypes["sabund"].push_back(outputFileName);  outputNames.push_back(outputFileName);
1300
1301                 InputData* input = new InputData(sabundfile, "sabund");
1302                 SAbundVector* sabund = input->getSAbundVector();
1303                 string lastLabel = sabund->getLabel();
1304
1305                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1306                 set<string> processedLabels;
1307                 set<string> userLabels = labels;
1308
1309                 if (size == 0) { //user has not set size, set size = 10%
1310                         size = int((sabund->getNumSeqs()) * 0.10);
1311                 }
1312
1313                 m->mothurOut("Sampling " + toString(size) + " from " + toString(sabund->getNumSeqs()) + "."); m->mothurOutEndLine();
1314
1315                 //as long as you are not at the end of the file or done wih the lines you want
1316                 while((sabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
1317                         if (m->control_pressed) {  delete input; delete sabund; out.close(); return 0;  }
1318
1319                         if(allLines == 1 || labels.count(sabund->getLabel()) == 1){
1320
1321                                 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
1322
1323                                 processSabund(sabund, out);
1324
1325                                 processedLabels.insert(sabund->getLabel());
1326                                 userLabels.erase(sabund->getLabel());
1327                         }
1328
1329                         if ((m->anyLabelsToProcess(sabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1330                                 string saveLabel = sabund->getLabel();
1331
1332                                 delete sabund;
1333
1334                                 sabund = input->getSAbundVector(lastLabel);
1335                                 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
1336
1337                                 processSabund(sabund, out);
1338
1339                                 processedLabels.insert(sabund->getLabel());
1340                                 userLabels.erase(sabund->getLabel());
1341
1342                                 //restore real lastlabel to save below
1343                                 sabund->setLabel(saveLabel);
1344                         }
1345
1346                         lastLabel = sabund->getLabel();
1347
1348                         //prevent memory leak
1349                         delete sabund; sabund = NULL;
1350
1351                         //get next line to process
1352                         sabund = input->getSAbundVector();
1353                 }
1354
1355
1356                 if (m->control_pressed) {  out.close(); return 0;  }
1357
1358                 //output error messages about any remaining user labels
1359                 set<string>::iterator it;
1360                 bool needToRun = false;
1361                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1362                         m->mothurOut("Your file does not include the label " + *it);
1363                         if (processedLabels.count(lastLabel) != 1) {
1364                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1365                                 needToRun = true;
1366                         }else {
1367                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1368                         }
1369                 }
1370
1371                 //run last label if you need to
1372                 if (needToRun == true)  {
1373                         if (sabund != NULL) { delete sabund; }
1374
1375                         sabund = input->getSAbundVector(lastLabel);
1376
1377                         m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
1378
1379                         processSabund(sabund, out);
1380
1381                         delete sabund;
1382                 }
1383
1384                 delete input;
1385                 out.close();
1386
1387                 return 0;
1388
1389         }
1390         catch(exception& e) {
1391                 m->errorOut(e, "SubSampleCommand", "getSubSampleSabund");
1392                 exit(1);
1393         }
1394 }
1395 //**********************************************************************************************************************
1396 int SubSampleCommand::processSabund(SAbundVector*& sabund, ofstream& out) {
1397         try {
1398
1399                 RAbundVector* rabund = new RAbundVector();
1400                 *rabund = sabund->getRAbundVector();
1401
1402                 int numBins = rabund->getNumBins();
1403                 int thisSize = rabund->getNumSeqs();
1404
1405                 if (thisSize != size) {
1406
1407                         OrderVector* order = new OrderVector();
1408                         for(int p=0;p<numBins;p++){
1409                                 for(int j=0;j<rabund->get(p);j++){
1410                                         order->push_back(p);
1411                                 }
1412                         }
1413                         random_shuffle(order->begin(), order->end());
1414
1415                         RAbundVector* temp = new RAbundVector(numBins);
1416                         temp->setLabel(rabund->getLabel());
1417
1418                         delete rabund;
1419                         rabund = temp;
1420
1421                         for (int j = 0; j < size; j++) {
1422
1423                                 if (m->control_pressed) { delete order; return 0; }
1424
1425                                 //get random number to sample from order between 0 and thisSize-1.
1426                                 int myrand = (int)((float)(rand()) / (RAND_MAX / (thisSize-1) + 1));
1427
1428                                 int bin = order->get(myrand);
1429
1430                                 int abund = rabund->get(bin);
1431                                 rabund->set(bin, (abund+1));
1432                         }
1433
1434                         delete order;
1435                 }
1436
1437                 if (m->control_pressed) { return 0; }
1438
1439                 delete sabund;
1440                 sabund = new SAbundVector();
1441                 *sabund = rabund->getSAbundVector();
1442                 delete rabund;
1443
1444                 sabund->print(out);
1445
1446                 return 0;
1447
1448         }
1449         catch(exception& e) {
1450                 m->errorOut(e, "SubSampleCommand", "processSabund");
1451                 exit(1);
1452         }
1453 }
1454 //**********************************************************************************************************************
1455 int SubSampleCommand::eliminateZeroOTUS(vector<SharedRAbundVector*>& thislookup) {
1456         try {
1457
1458                 vector<SharedRAbundVector*> newLookup;
1459                 for (int i = 0; i < thislookup.size(); i++) {
1460                         SharedRAbundVector* temp = new SharedRAbundVector();
1461                         temp->setLabel(thislookup[i]->getLabel());
1462                         temp->setGroup(thislookup[i]->getGroup());
1463                         newLookup.push_back(temp);
1464                 }
1465
1466                 //for each bin
1467                 for (int i = 0; i < thislookup[0]->getNumBins(); i++) {
1468                         if (m->control_pressed) { for (int j = 0; j < newLookup.size(); j++) {  delete newLookup[j];  } return 0; }
1469
1470                         //look at each sharedRabund and make sure they are not all zero
1471                         bool allZero = true;
1472                         for (int j = 0; j < thislookup.size(); j++) {
1473                                 if (thislookup[j]->getAbundance(i) != 0) { allZero = false;  break;  }
1474                         }
1475
1476                         //if they are not all zero add this bin
1477                         if (!allZero) {
1478                                 for (int j = 0; j < thislookup.size(); j++) {
1479                                         newLookup[j]->push_back(thislookup[j]->getAbundance(i), thislookup[j]->getGroup());
1480                                 }
1481                         }
1482                 }
1483
1484                 for (int j = 0; j < thislookup.size(); j++) {  delete thislookup[j];  }
1485                 thislookup.clear();
1486
1487                 thislookup = newLookup;
1488
1489                 return 0;
1490
1491         }
1492         catch(exception& e) {
1493                 m->errorOut(e, "SubSampleCommand", "eliminateZeroOTUS");
1494                 exit(1);
1495         }
1496 }
1497
1498 //**********************************************************************************************************************
1499
1500
1501