removerarecommand.cpp

   1 /*
   2  *  removerarecommand.cpp
   3  *  mothur
   4  *
   5  *  Created by westcott on 1/21/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "removerarecommand.h"
  11 #include "sequence.hpp"
  12 #include "groupmap.h"
  13 #include "sharedutilities.h"
  14 #include "inputdata.h"
  15
  16 //**********************************************************************************************************************
  17 vector<string> RemoveRareCommand::getValidParameters(){
  18         try {
  19                 string Array[] =  {"rabund","sabund", "group", "list", "shared","bygroup","nseqs","groups","label","outputdir","inputdir"};
  20                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  21                 return myArray;
  22         }
  23         catch(exception& e) {
  24                 m->errorOut(e, "RemoveRareCommand", "getValidParameters");
  25                 exit(1);
  26         }
  27 }
  28 //**********************************************************************************************************************
  29 RemoveRareCommand::RemoveRareCommand(){
  30         try {
  31                 abort = true; calledHelp = true;
  32                 vector<string> tempOutNames;
  33                 outputTypes["rabund"] = tempOutNames;
  34                 outputTypes["sabund"] = tempOutNames;
  35                 outputTypes["list"] = tempOutNames;
  36                 outputTypes["group"] = tempOutNames;
  37                 outputTypes["shared"] = tempOutNames;
  38         }
  39         catch(exception& e) {
  40                 m->errorOut(e, "RemoveRareCommand", "RemoveRareCommand");
  41                 exit(1);
  42         }
  43 }
  44 //**********************************************************************************************************************
  45 vector<string> RemoveRareCommand::getRequiredParameters(){
  46         try {
  47                 string Array[] =  {"nseqs"};
  48                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  49                 return myArray;
  50         }
  51         catch(exception& e) {
  52                 m->errorOut(e, "RemoveRareCommand", "getRequiredParameters");
  53                 exit(1);
  54         }
  55 }
  56 //**********************************************************************************************************************
  57 vector<string> RemoveRareCommand::getRequiredFiles(){
  58         try {
  59                 vector<string> myArray;
  60                 return myArray;
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "RemoveRareCommand", "getRequiredFiles");
  64                 exit(1);
  65         }
  66 }
  67 //**********************************************************************************************************************
  68 RemoveRareCommand::RemoveRareCommand(string option)  {
  69         try {
  70                 globaldata = GlobalData::getInstance();
  71                 abort = false; calledHelp = false;
  72                 allLines = 1;
  73
  74                 //allow user to run help
  75                 if(option == "help") { help(); abort = true; calledHelp = true; }
  76
  77                 else {
  78                         //valid paramters for this command
  79                         string Array[] =  {"rabund","sabund", "group", "list", "shared", "bygroup", "nseqs","groups","label","outputdir","inputdir"};
  80                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  81
  82                         OptionParser parser(option);
  83                         map<string,string> parameters = parser.getParameters();
  84
  85                         ValidParameters validParameter;
  86                         map<string,string>::iterator it;
  87
  88                         //check to make sure all parameters are valid for command
  89                         for (it = parameters.begin(); it != parameters.end(); it++) {
  90                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  91                         }
  92
  93                         //initialize outputTypes
  94                         vector<string> tempOutNames;
  95                         outputTypes["rabund"] = tempOutNames;
  96                         outputTypes["sabund"] = tempOutNames;
  97                         outputTypes["list"] = tempOutNames;
  98                         outputTypes["group"] = tempOutNames;
  99                         outputTypes["shared"] = tempOutNames;
 100
 101                         //if the user changes the output directory command factory will send this info to us in the output parameter
 102                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
 103
 104                         //if the user changes the input directory command factory will send this info to us in the output parameter
 105                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 106                         if (inputDir == "not found"){   inputDir = "";          }
 107                         else {
 108                                 string path;
 109                                 it = parameters.find("list");
 110                                 //user has given a template file
 111                                 if(it != parameters.end()){
 112                                         path = m->hasPath(it->second);
 113                                         //if the user has not given a path then, add inputdir. else leave path alone.
 114                                         if (path == "") {       parameters["list"] = inputDir + it->second;             }
 115                                 }
 116
 117                                 it = parameters.find("group");
 118                                 //user has given a template file
 119                                 if(it != parameters.end()){
 120                                         path = m->hasPath(it->second);
 121                                         //if the user has not given a path then, add inputdir. else leave path alone.
 122                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 123                                 }
 124
 125                                 it = parameters.find("sabund");
 126                                 //user has given a template file
 127                                 if(it != parameters.end()){
 128                                         path = m->hasPath(it->second);
 129                                         //if the user has not given a path then, add inputdir. else leave path alone.
 130                                         if (path == "") {       parameters["sabund"] = inputDir + it->second;           }
 131                                 }
 132
 133                                 it = parameters.find("rabund");
 134                                 //user has given a template file
 135                                 if(it != parameters.end()){
 136                                         path = m->hasPath(it->second);
 137                                         //if the user has not given a path then, add inputdir. else leave path alone.
 138                                         if (path == "") {       parameters["rabund"] = inputDir + it->second;           }
 139                                 }
 140
 141                                 it = parameters.find("shared");
 142                                 //user has given a template file
 143                                 if(it != parameters.end()){
 144                                         path = m->hasPath(it->second);
 145                                         //if the user has not given a path then, add inputdir. else leave path alone.
 146                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 147                                 }
 148                         }
 149
 150
 151                         //check for file parameters
 152                         listfile = validParameter.validFile(parameters, "list", true);
 153                         if (listfile == "not open") { abort = true; }
 154                         else if (listfile == "not found") {  listfile = "";  }
 155
 156                         sabundfile = validParameter.validFile(parameters, "sabund", true);
 157                         if (sabundfile == "not open") { abort = true; }
 158                         else if (sabundfile == "not found") {  sabundfile = "";  }
 159
 160                         rabundfile = validParameter.validFile(parameters, "rabund", true);
 161                         if (rabundfile == "not open") { abort = true; }
 162                         else if (rabundfile == "not found") {  rabundfile = "";  }
 163
 164                         groupfile = validParameter.validFile(parameters, "group", true);
 165                         if (groupfile == "not open") { groupfile = ""; abort = true; }
 166                         else if (groupfile == "not found") {  groupfile = "";  }
 167
 168                         sharedfile = validParameter.validFile(parameters, "shared", true);
 169                         if (sharedfile == "not open") { sharedfile = "";  abort = true; }
 170                         else if (sharedfile == "not found") {  sharedfile = "";  }
 171
 172                         if ((sabundfile == "") && (rabundfile == "")  && (sharedfile == "") && (listfile == ""))  { m->mothurOut("You must provide at least one of the following: rabund, sabund, shared or list."); m->mothurOutEndLine(); abort = true; }
 173
 174                         groups = validParameter.validFile(parameters, "groups", false);
 175                         if (groups == "not found") { groups = "all"; }
 176                         m->splitAtDash(groups, Groups);
 177
 178                         label = validParameter.validFile(parameters, "label", false);
 179                         if (label == "not found") { label = ""; }
 180                         else {
 181                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 182                                 else { allLines = 1;  }
 183                         }
 184
 185                         string temp = validParameter.validFile(parameters, "nseqs", false);
 186                         if (temp == "not found") { m->mothurOut("nseqs is a required parameter."); m->mothurOutEndLine(); abort = true; }
 187                         else { convert(temp, nseqs); }
 188
 189                         temp = validParameter.validFile(parameters, "bygroup", false);   if (temp == "not found") { temp = "f"; }
 190                         byGroup = m->isTrue(temp);
 191
 192                         if (byGroup && (sharedfile == "")) { m->mothurOut("The byGroup parameter is only valid with a shared file."); m->mothurOutEndLine(); }
 193
 194                         if ((groupfile != "") && (listfile == "")) { m->mothurOut("A groupfile is only valid with a list file."); m->mothurOutEndLine(); groupfile = ""; }
 195                 }
 196
 197         }
 198         catch(exception& e) {
 199                 m->errorOut(e, "RemoveRareCommand", "RemoveRareCommand");
 200                 exit(1);
 201         }
 202 }
 203 //**********************************************************************************************************************
 204
 205 void RemoveRareCommand::help(){
 206         try {
 207                 m->mothurOut("The remove.rare command parameters are list, rabund, sabund, shared, group, label, groups, bygroup and nseqs.\n");
 208                 m->mothurOut("The remove.rare command reads one of the following file types: list, rabund, sabund or shared file. It outputs a new file after removing the rare otus.\n");
 209                 m->mothurOut("The groups parameter allows you to specify which of the groups you would like analyzed.  Default=all. You may separate group names with dashes.\n");
 210                 m->mothurOut("The label parameter is used to analyze specific labels in your input. default=all. You may separate label names with dashes.\n");
 211                 m->mothurOut("The bygroup parameter is only valid with the shared file. default=f, meaning remove any OTU that has nseqs or fewer sequences across all groups.\n");
 212                 m->mothurOut("bygroups=T means remove any OTU that has nseqs or fewer sequences in each group (if groupA has 1 sequence and group B has 100 sequences in OTUZ and nseqs=1, then set the groupA count for OTUZ to 0 and keep groupB's count at 100.) \n");
 213                 m->mothurOut("The nseqs parameter allows you to set the cutoff for an otu to be deemed rare. It is required.\n");
 214                 m->mothurOut("The remove.rare command should be in the following format: remove.rare(shared=yourSharedFile, nseqs=yourRareCutoff).\n");
 215                 m->mothurOut("Example remove.rare(shared=amazon.fn.shared, nseqs=2).\n");
 216                 m->mothurOut("Note: No spaces between parameter labels (i.e. shared), '=' and parameters (i.e.yourSharedFile).\n\n");
 217         }
 218         catch(exception& e) {
 219                 m->errorOut(e, "RemoveRareCommand", "help");
 220                 exit(1);
 221         }
 222 }
 223
 224 //**********************************************************************************************************************
 225
 226 int RemoveRareCommand::execute(){
 227         try {
 228
 229                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 230
 231                 if (m->control_pressed) { return 0; }
 232
 233                 //read through the correct file and output lines you want to keep
 234                 if (sabundfile != "")           {               processSabund();        }
 235                 if (rabundfile != "")           {               processRabund();        }
 236                 if (listfile != "")                     {               processList();          }
 237                 if (sharedfile != "")           {               processShared();        }
 238
 239                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str()); } return 0; }
 240
 241                 if (outputNames.size() != 0) {
 242                         m->mothurOutEndLine();
 243                         m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 244                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 245                         m->mothurOutEndLine();
 246                 }
 247
 248                 return 0;
 249         }
 250
 251         catch(exception& e) {
 252                 m->errorOut(e, "RemoveRareCommand", "execute");
 253                 exit(1);
 254         }
 255 }
 256
 257 //**********************************************************************************************************************
 258 int RemoveRareCommand::processList(){
 259         try {
 260                 string thisOutputDir = outputDir;
 261                 if (outputDir == "") {  thisOutputDir += m->hasPath(listfile);  }
 262                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(listfile)) + "pick" +  m->getExtension(listfile);
 263                 string outputGroupFileName = thisOutputDir + m->getRootName(m->getSimpleName(groupfile)) + "pick" +  m->getExtension(groupfile);
 264
 265                 ofstream out, outGroup;
 266                 m->openOutputFile(outputFileName, out);
 267
 268                 bool wroteSomething = false;
 269
 270                 //you must provide a label because the names in the listfile need to be consistent
 271                 string thisLabel = "";
 272                 if (allLines) { m->mothurOut("For the listfile you must select one label, using first label in your listfile."); m->mothurOutEndLine(); }
 273                 else if (labels.size() > 1) { m->mothurOut("For the listfile you must select one label, using " + (*labels.begin()) + "."); m->mothurOutEndLine(); thisLabel = *labels.begin(); }
 274                 else { thisLabel = *labels.begin(); }
 275
 276                 InputData input(listfile, "list");
 277                 ListVector* list = input.getListVector();
 278
 279                 //get first one or the one we want
 280                 if (thisLabel != "") {
 281                         //use smart distancing
 282                         set<string> userLabels; userLabels.insert(thisLabel);
 283                         set<string> processedLabels;
 284                         string lastLabel = list->getLabel();
 285                         while((list != NULL) && (userLabels.size() != 0)) {
 286                                 if(userLabels.count(list->getLabel()) == 1){
 287                                         processedLabels.insert(list->getLabel());
 288                                         userLabels.erase(list->getLabel());
 289                                         break;
 290                                 }
 291
 292                                 if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 293                                         processedLabels.insert(list->getLabel());
 294                                         userLabels.erase(list->getLabel());
 295                                         delete list;
 296                                         list = input.getListVector(lastLabel);
 297                                         break;
 298                                 }
 299                                 lastLabel = list->getLabel();
 300                                 delete list;
 301                                 list = input.getListVector();
 302                         }
 303                         if (userLabels.size() != 0) {
 304                                 m->mothurOut("Your file does not include the label " + thisLabel + ". I will use " + lastLabel + ".");  m->mothurOutEndLine();
 305                                 list = input.getListVector(lastLabel);
 306                         }
 307                 }
 308
 309                 //if groupfile is given then use it
 310                 GroupMap* groupMap;
 311                 if (groupfile != "") {
 312                         groupMap = new GroupMap(groupfile); groupMap->readMap();
 313                         SharedUtil util;
 314                         util.setGroups(Groups, groupMap->namesOfGroups);
 315                         m->openOutputFile(outputGroupFileName, outGroup);
 316                 }
 317
 318
 319                 if (list != NULL) {
 320                         //make a new list vector
 321                         ListVector newList;
 322                         newList.setLabel(list->getLabel());
 323
 324                         //for each bin
 325                         for (int i = 0; i < list->getNumBins(); i++) {
 326                                 if (m->control_pressed) {  if (groupfile != "") { delete groupMap; outGroup.close(); remove(outputGroupFileName.c_str()); } out.close();  remove(outputFileName.c_str());  return 0; }
 327
 328                                 //parse out names that are in accnos file
 329                                 string binnames = list->get(i);
 330                                 vector<string> names;
 331                                 string saveBinNames = binnames;
 332                                 m->splitAtComma(binnames, names);
 333
 334                                 vector<string> newGroupFile;
 335                                 if (groupfile != "") {
 336                                         vector<string> newNames;
 337                                         saveBinNames = "";
 338                                         for(int k = 0; k < names.size(); k++) {
 339                                                 string group = groupMap->getGroup(names[k]);
 340
 341                                                 if (m->inUsersGroups(group, Groups)) {
 342                                                         newGroupFile.push_back(names[k] + "\t" + group);
 343
 344                                                         newNames.push_back(names[k]);
 345                                                         saveBinNames += names[k] + ",";
 346                                                 }
 347                                         }
 348                                         names = newNames;
 349                                         saveBinNames = saveBinNames.substr(0, saveBinNames.length()-1);
 350                                 }
 351
 352                                 if (names.size() > nseqs) { //keep bin
 353                                         newList.push_back(saveBinNames);
 354                                         for(int k = 0; k < newGroupFile.size(); k++) { outGroup << newGroupFile[k] << endl; }
 355                                 }
 356                         }
 357
 358                         //print new listvector
 359                         if (newList.getNumBins() != 0) {
 360                                 wroteSomething = true;
 361                                 newList.print(out);
 362                         }
 363                 }
 364
 365                 out.close();
 366                 if (groupfile != "") { outGroup.close(); outputTypes["group"].push_back(outputGroupFileName); outputNames.push_back(outputGroupFileName); }
 367
 368                 if (wroteSomething == false) {  m->mothurOut("Your file contains only rare sequences."); m->mothurOutEndLine();  }
 369                 outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName);
 370
 371                 return 0;
 372         }
 373         catch(exception& e) {
 374                 m->errorOut(e, "RemoveRareCommand", "processList");
 375                 exit(1);
 376         }
 377 }
 378 //**********************************************************************************************************************
 379 int RemoveRareCommand::processSabund(){
 380         try {
 381                 string thisOutputDir = outputDir;
 382                 if (outputDir == "") {  thisOutputDir += m->hasPath(sabundfile);  }
 383                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(sabundfile)) + "pick" +  m->getExtension(sabundfile);
 384                 outputTypes["sabund"].push_back(outputFileName); outputNames.push_back(outputFileName);
 385
 386                 ofstream out;
 387                 m->openOutputFile(outputFileName, out);
 388
 389                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 390                 InputData input(sabundfile, "sabund");
 391                 SAbundVector* sabund = input.getSAbundVector();
 392                 string lastLabel = sabund->getLabel();
 393                 set<string> processedLabels;
 394                 set<string> userLabels = labels;
 395
 396                 while((sabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 397
 398                         if (m->control_pressed) { delete sabund; out.close(); return 0; }
 399
 400                         if(allLines == 1 || labels.count(sabund->getLabel()) == 1){
 401
 402                                 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
 403                                 processedLabels.insert(sabund->getLabel());
 404                                 userLabels.erase(sabund->getLabel());
 405
 406                                 if (sabund->getMaxRank() > nseqs) {
 407                                         for(int i = 1; i <=nseqs; i++) {  sabund->set(i, 0); }
 408                                 }else { sabund->clear(); }
 409
 410                                 if (sabund->getNumBins() > 0) { sabund->print(out); }
 411                         }
 412
 413                         if ((m->anyLabelsToProcess(sabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 414                                 string saveLabel = sabund->getLabel();
 415
 416                                 delete sabund;
 417                                 sabund = input.getSAbundVector(lastLabel);
 418
 419                                 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
 420                                 processedLabels.insert(sabund->getLabel());
 421                                 userLabels.erase(sabund->getLabel());
 422
 423                                 if (sabund->getMaxRank() > nseqs) {
 424                                         for(int i = 1; i <=nseqs; i++) {  sabund->set(i, 0); }
 425                                 }else { sabund->clear(); }
 426
 427                                 if (sabund->getNumBins() > 0) { sabund->print(out); }
 428
 429                                 //restore real lastlabel to save below
 430                                 sabund->setLabel(saveLabel);
 431                         }
 432
 433                         lastLabel = sabund->getLabel();
 434
 435                         delete sabund;
 436                         sabund = input.getSAbundVector();
 437                 }
 438
 439                 if (m->control_pressed) {  out.close(); return 0; }
 440
 441                 //output error messages about any remaining user labels
 442                 set<string>::iterator it;
 443                 bool needToRun = false;
 444                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 445                         m->mothurOut("Your file does not include the label " + *it);
 446                         if (processedLabels.count(lastLabel) != 1) {
 447                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 448                                 needToRun = true;
 449                         }else {
 450                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 451                         }
 452                 }
 453
 454                 //run last label if you need to
 455                 if (needToRun == true)  {
 456                         if (sabund != NULL) {   delete sabund;  }
 457                         sabund = input.getSAbundVector(lastLabel);
 458
 459                         m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
 460
 461                         if (sabund->getMaxRank() > nseqs) {
 462                                 for(int i = 1; i <=nseqs; i++) {  sabund->set(i, 0); }
 463                         }else { sabund->clear(); }
 464
 465                         if (sabund->getNumBins() > 0) { sabund->print(out); }
 466
 467                         delete sabund;
 468                 }
 469
 470                 return 0;
 471         }
 472         catch(exception& e) {
 473                 m->errorOut(e, "RemoveRareCommand", "processSabund");
 474                 exit(1);
 475         }
 476 }
 477 //**********************************************************************************************************************
 478 int RemoveRareCommand::processRabund(){
 479         try {
 480                 string thisOutputDir = outputDir;
 481                 if (outputDir == "") {  thisOutputDir += m->hasPath(rabundfile);  }
 482                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(rabundfile)) + "pick" +  m->getExtension(rabundfile);
 483                 outputTypes["rabund"].push_back(outputFileName); outputNames.push_back(outputFileName);
 484
 485                 ofstream out;
 486                 m->openOutputFile(outputFileName, out);
 487
 488                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 489                 InputData input(rabundfile, "rabund");
 490                 RAbundVector* rabund = input.getRAbundVector();
 491                 string lastLabel = rabund->getLabel();
 492                 set<string> processedLabels;
 493                 set<string> userLabels = labels;
 494
 495                 while((rabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 496
 497                         if (m->control_pressed) { delete rabund; out.close(); return 0; }
 498
 499                         if(allLines == 1 || labels.count(rabund->getLabel()) == 1){
 500
 501                                 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
 502                                 processedLabels.insert(rabund->getLabel());
 503                                 userLabels.erase(rabund->getLabel());
 504
 505                                 RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
 506                                 for (int i = 0; i < rabund->getNumBins(); i++) {
 507                                         if (rabund->get(i) > nseqs) {
 508                                                 newRabund.push_back(rabund->get(i));
 509                                         }
 510                                 }
 511                                 if (newRabund.getNumBins() > 0) { newRabund.print(out); }
 512                         }
 513
 514                         if ((m->anyLabelsToProcess(rabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 515                                 string saveLabel = rabund->getLabel();
 516
 517                                 delete rabund;
 518                                 rabund = input.getRAbundVector(lastLabel);
 519
 520                                 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
 521                                 processedLabels.insert(rabund->getLabel());
 522                                 userLabels.erase(rabund->getLabel());
 523
 524                                 RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
 525                                 for (int i = 0; i < rabund->getNumBins(); i++) {
 526                                         if (rabund->get(i) > nseqs) {
 527                                                 newRabund.push_back(rabund->get(i));
 528                                         }
 529                                 }
 530                                 if (newRabund.getNumBins() > 0) { newRabund.print(out); }
 531
 532                                 //restore real lastlabel to save below
 533                                 rabund->setLabel(saveLabel);
 534                         }
 535
 536                         lastLabel = rabund->getLabel();
 537
 538                         delete rabund;
 539                         rabund = input.getRAbundVector();
 540                 }
 541
 542                 if (m->control_pressed) {  out.close(); return 0; }
 543
 544                 //output error messages about any remaining user labels
 545                 set<string>::iterator it;
 546                 bool needToRun = false;
 547                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 548                         m->mothurOut("Your file does not include the label " + *it);
 549                         if (processedLabels.count(lastLabel) != 1) {
 550                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 551                                 needToRun = true;
 552                         }else {
 553                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 554                         }
 555                 }
 556
 557                 //run last label if you need to
 558                 if (needToRun == true)  {
 559                         if (rabund != NULL) {   delete rabund;  }
 560                         rabund = input.getRAbundVector(lastLabel);
 561
 562                         m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
 563
 564                         RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
 565                         for (int i = 0; i < rabund->getNumBins(); i++) {
 566                                 if (rabund->get(i) > nseqs) {
 567                                         newRabund.push_back(rabund->get(i));
 568                                 }
 569                         }
 570                         if (newRabund.getNumBins() > 0) { newRabund.print(out); }
 571
 572                         delete rabund;
 573                 }
 574
 575                 return 0;
 576         }
 577         catch(exception& e) {
 578                 m->errorOut(e, "RemoveRareCommand", "processRabund");
 579                 exit(1);
 580         }
 581 }
 582 //**********************************************************************************************************************
 583 int RemoveRareCommand::processShared(){
 584         try {
 585                 globaldata->Groups = Groups;
 586
 587                 string thisOutputDir = outputDir;
 588                 if (outputDir == "") {  thisOutputDir += m->hasPath(sharedfile);  }
 589                 string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(sharedfile)) + "pick" +  m->getExtension(sharedfile);
 590                 outputTypes["shared"].push_back(outputFileName); outputNames.push_back(outputFileName);
 591
 592                 ofstream out;
 593                 m->openOutputFile(outputFileName, out);
 594
 595                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 596                 InputData input(sharedfile, "sharedfile");
 597                 vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
 598                 string lastLabel = lookup[0]->getLabel();
 599                 set<string> processedLabels;
 600                 set<string> userLabels = labels;
 601
 602                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 603
 604                         if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }  out.close(); return 0; }
 605
 606                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 607
 608                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 609                                 processedLabels.insert(lookup[0]->getLabel());
 610                                 userLabels.erase(lookup[0]->getLabel());
 611
 612                                 processLookup(lookup, out);
 613                         }
 614
 615                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 616                                 string saveLabel = lookup[0]->getLabel();
 617
 618                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 619                                 lookup = input.getSharedRAbundVectors(lastLabel);
 620
 621                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 622                                 processedLabels.insert(lookup[0]->getLabel());
 623                                 userLabels.erase(lookup[0]->getLabel());
 624
 625                                 processLookup(lookup, out);
 626
 627                                 //restore real lastlabel to save below
 628                                 lookup[0]->setLabel(saveLabel);
 629                         }
 630
 631                         lastLabel = lookup[0]->getLabel();
 632
 633                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 634                         lookup = input.getSharedRAbundVectors();
 635                 }
 636
 637                 if (m->control_pressed) {  out.close(); return 0; }
 638
 639                 //output error messages about any remaining user labels
 640                 set<string>::iterator it;
 641                 bool needToRun = false;
 642                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 643                         m->mothurOut("Your file does not include the label " + *it);
 644                         if (processedLabels.count(lastLabel) != 1) {
 645                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 646                                 needToRun = true;
 647                         }else {
 648                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 649                         }
 650                 }
 651
 652                 //run last label if you need to
 653                 if (needToRun == true)  {
 654                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       }  }
 655                         lookup = input.getSharedRAbundVectors(lastLabel);
 656
 657                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 658
 659                         processLookup(lookup, out);
 660
 661                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 662                 }
 663
 664                 return 0;
 665         }
 666         catch(exception& e) {
 667                 m->errorOut(e, "RemoveRareCommand", "processSabund");
 668                 exit(1);
 669         }
 670 }
 671 //**********************************************************************************************************************
 672 int RemoveRareCommand::processLookup(vector<SharedRAbundVector*>& lookup, ofstream& out){
 673         try {
 674
 675                 vector<SharedRAbundVector> newRabunds;  newRabunds.resize(lookup.size());
 676                 for (int i = 0; i < lookup.size(); i++) {
 677                         newRabunds[i].setGroup(lookup[i]->getGroup());
 678                         newRabunds[i].setLabel(lookup[i]->getLabel());
 679                 }
 680
 681                 if (byGroup) {
 682
 683                         //for each otu
 684                         for (int i = 0; i < lookup[0]->getNumBins(); i++) {
 685                                 bool allZero = true;
 686
 687                                 if (m->control_pressed) { return 0; }
 688
 689                                 //for each group
 690                                 for (int j = 0; j < lookup.size(); j++) {
 691
 692                                         //are you rare?
 693                                         if (lookup[j]->getAbundance(i) > nseqs) {
 694                                                 newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
 695                                                 allZero = false;
 696                                         }else {
 697                                                 newRabunds[j].push_back(0, newRabunds[j].getGroup());
 698                                         }
 699                                 }
 700
 701                                 //eliminates zero otus
 702                                 if (allZero) { for (int j = 0; j < newRabunds.size(); j++) {  newRabunds[j].pop_back(); } }
 703                         }
 704                 }else {
 705                         //for each otu
 706                         for (int i = 0; i < lookup[0]->getNumBins(); i++) {
 707
 708                                 if (m->control_pressed) { return 0; }
 709
 710                                 int totalAbund = 0;
 711                                 //get total otu abundance
 712                                 for (int j = 0; j < lookup.size(); j++) {
 713                                         newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
 714                                         totalAbund += lookup[j]->getAbundance(i);
 715                                 }
 716
 717                                 //eliminates otus below rare cutoff
 718                                 if (totalAbund <= nseqs) { for (int j = 0; j < newRabunds.size(); j++) {  newRabunds[j].pop_back(); } }
 719                         }
 720                 }
 721
 722                 //do we have any otus above the rare cutoff
 723                 if (newRabunds[0].getNumBins() != 0) {
 724                         for (int j = 0; j < newRabunds.size(); j++) {
 725                                 out << newRabunds[j].getLabel() << '\t' << newRabunds[j].getGroup() << '\t';
 726                                 newRabunds[j].print(out);
 727                         }
 728                 }
 729
 730                 return 0;
 731         }
 732         catch(exception& e) {
 733                 m->errorOut(e, "RemoveRareCommand", "processLookup");
 734                 exit(1);
 735         }
 736 }
 737 //**********************************************************************************************************************
 738
 739
 740
 741