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::setParameters(){
  18         try {
  19                 CommandParameter plist("list", "InputTypes", "", "", "none", "atleast", "none","list",false,false,true); parameters.push_back(plist);
  20                 CommandParameter prabund("rabund", "InputTypes", "", "", "none", "atleast", "none","rabund",false,false,true); parameters.push_back(prabund);
  21                 CommandParameter psabund("sabund", "InputTypes", "", "", "none", "atleast", "none","sabund",false,false,true); parameters.push_back(psabund);
  22                 CommandParameter pshared("shared", "InputTypes", "", "", "none", "atleast", "none","shared",false,false,true); parameters.push_back(pshared);
  23         CommandParameter pcount("count", "InputTypes", "", "", "CountGroup", "none", "none","count",false,false); parameters.push_back(pcount);
  24                 CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","group",false,false); parameters.push_back(pgroup);
  25                 CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups);
  26                 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
  27                 CommandParameter pnseqs("nseqs", "Number", "", "0", "", "", "","",false,true,true); parameters.push_back(pnseqs);
  28                 CommandParameter pbygroup("bygroup", "Boolean", "", "f", "", "", "","",false,false); parameters.push_back(pbygroup);
  29                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  30                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  31
  32                 vector<string> myArray;
  33                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  34                 return myArray;
  35         }
  36         catch(exception& e) {
  37                 m->errorOut(e, "RemoveRareCommand", "setParameters");
  38                 exit(1);
  39         }
  40 }
  41 //**********************************************************************************************************************
  42 string RemoveRareCommand::getHelpString(){
  43         try {
  44                 string helpString = "";
  45                 helpString += "The remove.rare command parameters are list, rabund, sabund, shared, group, count, label, groups, bygroup and nseqs.\n";
  46                 helpString += "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";
  47                 helpString += "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";
  48                 helpString += "The label parameter is used to analyze specific labels in your input. default=all. You may separate label names with dashes.\n";
  49                 helpString += "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";
  50                 helpString += "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";
  51                 helpString += "The nseqs parameter allows you to set the cutoff for an otu to be deemed rare. It is required.\n";
  52                 helpString += "The remove.rare command should be in the following format: remove.rare(shared=yourSharedFile, nseqs=yourRareCutoff).\n";
  53                 helpString += "Example remove.rare(shared=amazon.fn.shared, nseqs=2).\n";
  54                 helpString += "Note: No spaces between parameter labels (i.e. shared), '=' and parameters (i.e.yourSharedFile).\n";
  55                 return helpString;
  56         }
  57         catch(exception& e) {
  58                 m->errorOut(e, "RemoveRareCommand", "getHelpString");
  59                 exit(1);
  60         }
  61 }
  62 //**********************************************************************************************************************
  63 string RemoveRareCommand::getOutputPattern(string type) {
  64     try {
  65         string pattern = "";
  66
  67         if (type == "rabund")            {   pattern = "[filename],pick,[extension]";    }
  68         else if (type == "sabund")    {   pattern = "[filename],pick,[extension]";    }
  69         else if (type == "group")       {   pattern = "[filename],pick,[extension]";    }
  70         else if (type == "count")       {   pattern = "[filename],pick,[extension]";    }
  71         else if (type == "list")        {   pattern = "[filename],[tag],pick,[extension]";    }
  72         else if (type == "shared")      {   pattern = "[filename],[tag],pick,[extension]";    }
  73         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  74
  75         return pattern;
  76     }
  77     catch(exception& e) {
  78         m->errorOut(e, "RemoveRareCommand", "getOutputPattern");
  79         exit(1);
  80     }
  81 }
  82 //**********************************************************************************************************************
  83 RemoveRareCommand::RemoveRareCommand(){
  84         try {
  85                 abort = true; calledHelp = true;
  86                 setParameters();
  87                 vector<string> tempOutNames;
  88                 outputTypes["rabund"] = tempOutNames;
  89                 outputTypes["sabund"] = tempOutNames;
  90                 outputTypes["list"] = tempOutNames;
  91                 outputTypes["group"] = tempOutNames;
  92         outputTypes["count"] = tempOutNames;
  93                 outputTypes["shared"] = tempOutNames;
  94         }
  95         catch(exception& e) {
  96                 m->errorOut(e, "RemoveRareCommand", "RemoveRareCommand");
  97                 exit(1);
  98         }
  99 }
 100 //**********************************************************************************************************************
 101 RemoveRareCommand::RemoveRareCommand(string option)  {
 102         try {
 103                 abort = false; calledHelp = false;
 104                 allLines = 1;
 105
 106                 //allow user to run help
 107                 if(option == "help") { help(); abort = true; calledHelp = true; }
 108                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 109
 110                 else {
 111                         vector<string> myArray = setParameters();
 112
 113                         OptionParser parser(option);
 114                         map<string,string> parameters = parser.getParameters();
 115
 116                         ValidParameters validParameter;
 117                         map<string,string>::iterator it;
 118
 119                         //check to make sure all parameters are valid for command
 120                         for (it = parameters.begin(); it != parameters.end(); it++) {
 121                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 122                         }
 123
 124                         //initialize outputTypes
 125                         vector<string> tempOutNames;
 126                         outputTypes["rabund"] = tempOutNames;
 127                         outputTypes["sabund"] = tempOutNames;
 128                         outputTypes["list"] = tempOutNames;
 129                         outputTypes["group"] = tempOutNames;
 130                         outputTypes["shared"] = tempOutNames;
 131             outputTypes["count"] = tempOutNames;
 132
 133                         //if the user changes the output directory command factory will send this info to us in the output parameter
 134                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
 135
 136                         //if the user changes the input directory command factory will send this info to us in the output parameter
 137                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 138                         if (inputDir == "not found"){   inputDir = "";          }
 139                         else {
 140                                 string path;
 141                                 it = parameters.find("list");
 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["list"] = inputDir + it->second;             }
 147                                 }
 148
 149                                 it = parameters.find("group");
 150                                 //user has given a template file
 151                                 if(it != parameters.end()){
 152                                         path = m->hasPath(it->second);
 153                                         //if the user has not given a path then, add inputdir. else leave path alone.
 154                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 155                                 }
 156
 157                                 it = parameters.find("sabund");
 158                                 //user has given a template file
 159                                 if(it != parameters.end()){
 160                                         path = m->hasPath(it->second);
 161                                         //if the user has not given a path then, add inputdir. else leave path alone.
 162                                         if (path == "") {       parameters["sabund"] = inputDir + it->second;           }
 163                                 }
 164
 165                                 it = parameters.find("rabund");
 166                                 //user has given a template file
 167                                 if(it != parameters.end()){
 168                                         path = m->hasPath(it->second);
 169                                         //if the user has not given a path then, add inputdir. else leave path alone.
 170                                         if (path == "") {       parameters["rabund"] = inputDir + it->second;           }
 171                                 }
 172
 173                                 it = parameters.find("shared");
 174                                 //user has given a template file
 175                                 if(it != parameters.end()){
 176                                         path = m->hasPath(it->second);
 177                                         //if the user has not given a path then, add inputdir. else leave path alone.
 178                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 179                                 }
 180
 181                 it = parameters.find("count");
 182                                 //user has given a template file
 183                                 if(it != parameters.end()){
 184                                         path = m->hasPath(it->second);
 185                                         //if the user has not given a path then, add inputdir. else leave path alone.
 186                                         if (path == "") {       parameters["count"] = inputDir + it->second;            }
 187                                 }
 188                         }
 189
 190
 191                         //check for file parameters
 192                         listfile = validParameter.validFile(parameters, "list", true);
 193                         if (listfile == "not open") { abort = true; }
 194                         else if (listfile == "not found") {  listfile = "";  }
 195                         else { m->setListFile(listfile); }
 196
 197                         sabundfile = validParameter.validFile(parameters, "sabund", true);
 198                         if (sabundfile == "not open") { abort = true; }
 199                         else if (sabundfile == "not found") {  sabundfile = "";  }
 200                         else { m->setSabundFile(sabundfile); }
 201
 202                         rabundfile = validParameter.validFile(parameters, "rabund", true);
 203                         if (rabundfile == "not open") { abort = true; }
 204                         else if (rabundfile == "not found") {  rabundfile = "";  }
 205                         else { m->setRabundFile(rabundfile); }
 206
 207                         groupfile = validParameter.validFile(parameters, "group", true);
 208                         if (groupfile == "not open") { groupfile = ""; abort = true; }
 209                         else if (groupfile == "not found") {  groupfile = "";  }
 210                         else { m->setGroupFile(groupfile); }
 211
 212                         sharedfile = validParameter.validFile(parameters, "shared", true);
 213                         if (sharedfile == "not open") { sharedfile = "";  abort = true; }
 214                         else if (sharedfile == "not found") {  sharedfile = "";  }
 215                         else { m->setSharedFile(sharedfile); }
 216
 217             countfile = validParameter.validFile(parameters, "count", true);
 218                         if (countfile == "not open") { countfile = ""; abort = true; }
 219                         else if (countfile == "not found") { countfile = "";  }
 220                         else { m->setCountTableFile(countfile); }
 221
 222             if ((groupfile != "") && (countfile != "")) {
 223                 m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true;
 224             }
 225
 226                         if ((sharedfile == "") && (listfile == "") && (rabundfile == "") && (sabundfile == "")) {
 227                                 //is there are current file available for any of these?
 228                                 //give priority to shared, then list, then rabund, then sabund
 229                                 //if there is a current shared file, use it
 230                                 sharedfile = m->getSharedFile();
 231                                 if (sharedfile != "") {  m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 232                                 else {
 233                                         listfile = m->getListFile();
 234                                         if (listfile != "") {  m->mothurOut("Using " + listfile + " as input file for the list parameter."); m->mothurOutEndLine(); }
 235                                         else {
 236                                                 rabundfile = m->getRabundFile();
 237                                                 if (rabundfile != "") {  m->mothurOut("Using " + rabundfile + " as input file for the rabund parameter."); m->mothurOutEndLine(); }
 238                                                 else {
 239                                                         sabundfile = m->getSabundFile();
 240                                                         if (sabundfile != "") {  m->mothurOut("Using " + sabundfile + " as input file for the sabund parameter."); m->mothurOutEndLine(); }
 241                                                         else {
 242                                                                 m->mothurOut("No valid current files. You must provide a list, sabund, rabund or shared file."); m->mothurOutEndLine();
 243                                                                 abort = true;
 244                                                         }
 245                                                 }
 246                                         }
 247                                 }
 248                         }
 249
 250                         groups = validParameter.validFile(parameters, "groups", false);
 251                         if (groups == "not found") { groups = "all"; }
 252                         m->splitAtDash(groups, Groups);
 253
 254                         label = validParameter.validFile(parameters, "label", false);
 255                         if (label == "not found") { label = ""; }
 256                         else {
 257                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 258                                 else { allLines = 1;  }
 259                         }
 260
 261                         string temp = validParameter.validFile(parameters, "nseqs", false);
 262                         if (temp == "not found") { m->mothurOut("nseqs is a required parameter."); m->mothurOutEndLine(); abort = true; }
 263                         else { m->mothurConvert(temp, nseqs); }
 264
 265                         temp = validParameter.validFile(parameters, "bygroup", false);   if (temp == "not found") { temp = "f"; }
 266                         byGroup = m->isTrue(temp);
 267
 268                         if (byGroup && (sharedfile == "")) { m->mothurOut("The byGroup parameter is only valid with a shared file."); m->mothurOutEndLine(); }
 269
 270                         if (((groupfile != "") || (countfile != "")) && (listfile == "")) { m->mothurOut("A group or count file is only valid with a list file."); m->mothurOutEndLine(); groupfile = ""; countfile = ""; }
 271                 }
 272
 273         }
 274         catch(exception& e) {
 275                 m->errorOut(e, "RemoveRareCommand", "RemoveRareCommand");
 276                 exit(1);
 277         }
 278 }
 279 //**********************************************************************************************************************
 280
 281 int RemoveRareCommand::execute(){
 282         try {
 283
 284                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 285
 286                 if (m->control_pressed) { return 0; }
 287
 288                 //read through the correct file and output lines you want to keep
 289                 if (sabundfile != "")           {               processSabund();        }
 290                 if (rabundfile != "")           {               processRabund();        }
 291                 if (listfile != "")                     {               processList();          }
 292                 if (sharedfile != "")           {               processShared();        }
 293
 294                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
 295
 296                 if (outputNames.size() != 0) {
 297                         m->mothurOutEndLine();
 298                         m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 299                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 300                         m->mothurOutEndLine();
 301
 302                         //set rabund file as new current rabundfile
 303                         string current = "";
 304                         itTypes = outputTypes.find("rabund");
 305                         if (itTypes != outputTypes.end()) {
 306                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRabundFile(current); }
 307                         }
 308
 309                         itTypes = outputTypes.find("sabund");
 310                         if (itTypes != outputTypes.end()) {
 311                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSabundFile(current); }
 312                         }
 313
 314                         itTypes = outputTypes.find("group");
 315                         if (itTypes != outputTypes.end()) {
 316                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
 317                         }
 318
 319                         itTypes = outputTypes.find("list");
 320                         if (itTypes != outputTypes.end()) {
 321                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
 322                         }
 323
 324                         itTypes = outputTypes.find("shared");
 325                         if (itTypes != outputTypes.end()) {
 326                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSharedFile(current); }
 327                         }
 328
 329             itTypes = outputTypes.find("count");
 330                         if (itTypes != outputTypes.end()) {
 331                                 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
 332                         }
 333                 }
 334
 335                 return 0;
 336         }
 337
 338         catch(exception& e) {
 339                 m->errorOut(e, "RemoveRareCommand", "execute");
 340                 exit(1);
 341         }
 342 }
 343
 344 //**********************************************************************************************************************
 345 int RemoveRareCommand::processList(){
 346         try {
 347
 348                 //you must provide a label because the names in the listfile need to be consistent
 349                 string thisLabel = "";
 350                 if (allLines) { m->mothurOut("For the listfile you must select one label, using first label in your listfile."); m->mothurOutEndLine(); }
 351                 else if (labels.size() > 1) { m->mothurOut("For the listfile you must select one label, using " + (*labels.begin()) + "."); m->mothurOutEndLine(); thisLabel = *labels.begin(); }
 352                 else { thisLabel = *labels.begin(); }
 353
 354                 InputData input(listfile, "list");
 355                 ListVector* list = input.getListVector();
 356
 357                 //get first one or the one we want
 358                 if (thisLabel != "") {
 359                         //use smart distancing
 360                         set<string> userLabels; userLabels.insert(thisLabel);
 361                         set<string> processedLabels;
 362                         string lastLabel = list->getLabel();
 363                         while((list != NULL) && (userLabels.size() != 0)) {
 364                                 if(userLabels.count(list->getLabel()) == 1){
 365                                         processedLabels.insert(list->getLabel());
 366                                         userLabels.erase(list->getLabel());
 367                                         break;
 368                                 }
 369
 370                                 if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 371                                         processedLabels.insert(list->getLabel());
 372                                         userLabels.erase(list->getLabel());
 373                                         delete list;
 374                                         list = input.getListVector(lastLabel);
 375                                         break;
 376                                 }
 377                                 lastLabel = list->getLabel();
 378                                 delete list;
 379                                 list = input.getListVector();
 380                         }
 381                         if (userLabels.size() != 0) {
 382                                 m->mothurOut("Your file does not include the label " + thisLabel + ". I will use " + lastLabel + ".");  m->mothurOutEndLine();
 383                                 list = input.getListVector(lastLabel);
 384                         }
 385                 }
 386
 387         string thisOutputDir = outputDir;
 388                 if (outputDir == "") {  thisOutputDir += m->hasPath(listfile);  }
 389         map<string, string> variables;
 390         variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(listfile));
 391         variables["[extension]"] = m->getExtension(listfile);
 392         variables["[tag]"] = list->getLabel();
 393                 string outputFileName = getOutputFileName("list", variables);
 394         variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(groupfile));
 395         variables["[extension]"] = m->getExtension(groupfile);
 396                 string outputGroupFileName = getOutputFileName("group", variables);
 397         variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(countfile));
 398         variables["[extension]"] = m->getExtension(countfile);
 399         string outputCountFileName = getOutputFileName("count", variables);
 400
 401                 ofstream out, outGroup;
 402                 m->openOutputFile(outputFileName, out);
 403
 404                 bool wroteSomething = false;
 405
 406
 407                 //if groupfile is given then use it
 408                 GroupMap* groupMap;
 409         CountTable ct;
 410                 if (groupfile != "") {
 411                         groupMap = new GroupMap(groupfile); groupMap->readMap();
 412                         SharedUtil util;
 413                         vector<string> namesGroups = groupMap->getNamesOfGroups();
 414                         util.setGroups(Groups, namesGroups);
 415                         m->openOutputFile(outputGroupFileName, outGroup);
 416                 }else if (countfile != "") {
 417             ct.readTable(countfile, true, false);
 418             if (ct.hasGroupInfo()) {
 419                 vector<string> namesGroups = ct.getNamesOfGroups();
 420                 SharedUtil util;
 421                 util.setGroups(Groups, namesGroups);
 422             }
 423         }
 424
 425
 426                 if (list != NULL) {
 427
 428             vector<string> binLabels = list->getLabels();
 429             vector<string> newLabels;
 430
 431                         //make a new list vector
 432                         ListVector newList;
 433                         newList.setLabel(list->getLabel());
 434
 435                         //for each bin
 436                         for (int i = 0; i < list->getNumBins(); i++) {
 437                                 if (m->control_pressed) {  if (groupfile != "") { delete groupMap; outGroup.close(); m->mothurRemove(outputGroupFileName); } out.close();  m->mothurRemove(outputFileName);  return 0; }
 438
 439                                 //parse out names that are in accnos file
 440                                 string binnames = list->get(i);
 441                                 vector<string> names;
 442                                 string saveBinNames = binnames;
 443                                 m->splitAtComma(binnames, names);
 444                 int binsize = names.size();
 445
 446                                 vector<string> newGroupFile;
 447                                 if (groupfile != "") {
 448                                         vector<string> newNames;
 449                                         saveBinNames = "";
 450                                         for(int k = 0; k < names.size(); k++) {
 451                                                 string group = groupMap->getGroup(names[k]);
 452
 453                                                 if (m->inUsersGroups(group, Groups)) {
 454                                                         newGroupFile.push_back(names[k] + "\t" + group);
 455
 456                                                         newNames.push_back(names[k]);
 457                                                         saveBinNames += names[k] + ",";
 458                                                 }
 459                                         }
 460                                         names = newNames; binsize = names.size();
 461                                         saveBinNames = saveBinNames.substr(0, saveBinNames.length()-1);
 462                                 }else if (countfile != "") {
 463                                         saveBinNames = "";
 464                     binsize = 0;
 465                                         for(int k = 0; k < names.size(); k++) {
 466                         if (ct.hasGroupInfo()) {
 467                             vector<string> thisSeqsGroups = ct.getGroups(names[k]);
 468
 469                             int thisSeqsCount = 0;
 470                             for (int n = 0; n < thisSeqsGroups.size(); n++) {
 471                                 if (m->inUsersGroups(thisSeqsGroups[n], Groups)) {
 472                                     thisSeqsCount += ct.getGroupCount(names[k], thisSeqsGroups[n]);
 473                                 }
 474                             }
 475                             binsize += thisSeqsCount;
 476                             //if you don't have any seqs from the groups the user wants, then remove you.
 477                             if (thisSeqsCount == 0) { newGroupFile.push_back(names[k]); }
 478                             else { saveBinNames += names[k] + ","; }
 479                         }else {
 480                             binsize += ct.getNumSeqs(names[k]);
 481                             saveBinNames += names[k] + ",";
 482                         }
 483                                         }
 484                                         saveBinNames = saveBinNames.substr(0, saveBinNames.length()-1);
 485                 }
 486
 487                                 if (binsize > nseqs) { //keep bin
 488                                         newList.push_back(saveBinNames);
 489                     newLabels.push_back(binLabels[i]);
 490                                         if (groupfile != "") {  for(int k = 0; k < newGroupFile.size(); k++) { outGroup << newGroupFile[k] << endl; }  }
 491                     else if (countfile != "") { for(int k = 0; k < newGroupFile.size(); k++) {  ct.remove(newGroupFile[k]); } }
 492                                 }else {  if (countfile != "") {  for(int k = 0; k < names.size(); k++) {  ct.remove(names[k]); } }  }
 493                         }
 494
 495                         //print new listvector
 496                         if (newList.getNumBins() != 0) {
 497                                 wroteSomething = true;
 498                                 newList.setLabels(newLabels);
 499                 newList.printHeaders(out);
 500                 newList.print(out);
 501                         }
 502                 }
 503
 504                 out.close();
 505                 if (groupfile != "") { outGroup.close(); outputTypes["group"].push_back(outputGroupFileName); outputNames.push_back(outputGroupFileName); }
 506         if (countfile != "") {
 507             if (ct.hasGroupInfo()) {
 508                 vector<string> allGroups = ct.getNamesOfGroups();
 509                 for (int i = 0; i < allGroups.size(); i++) {
 510                     if (!m->inUsersGroups(allGroups[i], Groups)) { ct.removeGroup(allGroups[i]); }
 511                 }
 512
 513             }
 514             ct.printTable(outputCountFileName);
 515             outputTypes["count"].push_back(outputCountFileName); outputNames.push_back(outputCountFileName);
 516         }
 517
 518                 if (wroteSomething == false) {  m->mothurOut("Your file contains only rare sequences."); m->mothurOutEndLine();  }
 519                 outputTypes["list"].push_back(outputFileName); outputNames.push_back(outputFileName);
 520
 521                 return 0;
 522         }
 523         catch(exception& e) {
 524                 m->errorOut(e, "RemoveRareCommand", "processList");
 525                 exit(1);
 526         }
 527 }
 528 //**********************************************************************************************************************
 529 int RemoveRareCommand::processSabund(){
 530         try {
 531                 string thisOutputDir = outputDir;
 532                 if (outputDir == "") {  thisOutputDir += m->hasPath(sabundfile);  }
 533         map<string, string> variables;
 534         variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(sabundfile));
 535         variables["[extension]"] = m->getExtension(sabundfile);
 536                 string outputFileName = getOutputFileName("sabund", variables);
 537                 outputTypes["sabund"].push_back(outputFileName); outputNames.push_back(outputFileName);
 538
 539                 ofstream out;
 540                 m->openOutputFile(outputFileName, out);
 541
 542                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 543                 InputData input(sabundfile, "sabund");
 544                 SAbundVector* sabund = input.getSAbundVector();
 545                 string lastLabel = sabund->getLabel();
 546                 set<string> processedLabels;
 547                 set<string> userLabels = labels;
 548
 549                 while((sabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 550
 551                         if (m->control_pressed) { delete sabund; out.close(); return 0; }
 552
 553                         if(allLines == 1 || labels.count(sabund->getLabel()) == 1){
 554
 555                                 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
 556                                 processedLabels.insert(sabund->getLabel());
 557                                 userLabels.erase(sabund->getLabel());
 558
 559                                 if (sabund->getMaxRank() > nseqs) {
 560                                         for(int i = 1; i <=nseqs; i++) {  sabund->set(i, 0); }
 561                                 }else { sabund->clear(); }
 562
 563                                 if (sabund->getNumBins() > 0) { sabund->print(out); }
 564                         }
 565
 566                         if ((m->anyLabelsToProcess(sabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 567                                 string saveLabel = sabund->getLabel();
 568
 569                                 delete sabund;
 570                                 sabund = input.getSAbundVector(lastLabel);
 571
 572                                 m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
 573                                 processedLabels.insert(sabund->getLabel());
 574                                 userLabels.erase(sabund->getLabel());
 575
 576                                 if (sabund->getMaxRank() > nseqs) {
 577                                         for(int i = 1; i <=nseqs; i++) {  sabund->set(i, 0); }
 578                                 }else { sabund->clear(); }
 579
 580                                 if (sabund->getNumBins() > 0) { sabund->print(out); }
 581
 582                                 //restore real lastlabel to save below
 583                                 sabund->setLabel(saveLabel);
 584                         }
 585
 586                         lastLabel = sabund->getLabel();
 587
 588                         delete sabund;
 589                         sabund = input.getSAbundVector();
 590                 }
 591
 592                 if (m->control_pressed) {  out.close(); return 0; }
 593
 594                 //output error messages about any remaining user labels
 595                 set<string>::iterator it;
 596                 bool needToRun = false;
 597                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 598                         m->mothurOut("Your file does not include the label " + *it);
 599                         if (processedLabels.count(lastLabel) != 1) {
 600                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 601                                 needToRun = true;
 602                         }else {
 603                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 604                         }
 605                 }
 606
 607                 //run last label if you need to
 608                 if (needToRun == true)  {
 609                         if (sabund != NULL) {   delete sabund;  }
 610                         sabund = input.getSAbundVector(lastLabel);
 611
 612                         m->mothurOut(sabund->getLabel()); m->mothurOutEndLine();
 613
 614                         if (sabund->getMaxRank() > nseqs) {
 615                                 for(int i = 1; i <=nseqs; i++) {  sabund->set(i, 0); }
 616                         }else { sabund->clear(); }
 617
 618                         if (sabund->getNumBins() > 0) { sabund->print(out); }
 619
 620                         delete sabund;
 621                 }
 622
 623                 return 0;
 624         }
 625         catch(exception& e) {
 626                 m->errorOut(e, "RemoveRareCommand", "processSabund");
 627                 exit(1);
 628         }
 629 }
 630 //**********************************************************************************************************************
 631 int RemoveRareCommand::processRabund(){
 632         try {
 633                 string thisOutputDir = outputDir;
 634                 if (outputDir == "") {  thisOutputDir += m->hasPath(rabundfile);  }
 635         map<string, string> variables;
 636         variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(rabundfile));
 637         variables["[extension]"] = m->getExtension(rabundfile);
 638                 string outputFileName = getOutputFileName("rabund", variables);
 639                 outputTypes["rabund"].push_back(outputFileName); outputNames.push_back(outputFileName);
 640
 641                 ofstream out;
 642                 m->openOutputFile(outputFileName, out);
 643
 644                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 645                 InputData input(rabundfile, "rabund");
 646                 RAbundVector* rabund = input.getRAbundVector();
 647                 string lastLabel = rabund->getLabel();
 648                 set<string> processedLabels;
 649                 set<string> userLabels = labels;
 650
 651                 while((rabund != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 652
 653                         if (m->control_pressed) { delete rabund; out.close(); return 0; }
 654
 655                         if(allLines == 1 || labels.count(rabund->getLabel()) == 1){
 656
 657                                 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
 658                                 processedLabels.insert(rabund->getLabel());
 659                                 userLabels.erase(rabund->getLabel());
 660
 661                                 RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
 662                                 for (int i = 0; i < rabund->getNumBins(); i++) {
 663                                         if (rabund->get(i) > nseqs) {
 664                                                 newRabund.push_back(rabund->get(i));
 665                                         }
 666                                 }
 667                                 if (newRabund.getNumBins() > 0) { newRabund.print(out); }
 668                         }
 669
 670                         if ((m->anyLabelsToProcess(rabund->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 671                                 string saveLabel = rabund->getLabel();
 672
 673                                 delete rabund;
 674                                 rabund = input.getRAbundVector(lastLabel);
 675
 676                                 m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
 677                                 processedLabels.insert(rabund->getLabel());
 678                                 userLabels.erase(rabund->getLabel());
 679
 680                                 RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
 681                                 for (int i = 0; i < rabund->getNumBins(); i++) {
 682                                         if (rabund->get(i) > nseqs) {
 683                                                 newRabund.push_back(rabund->get(i));
 684                                         }
 685                                 }
 686                                 if (newRabund.getNumBins() > 0) { newRabund.print(out); }
 687
 688                                 //restore real lastlabel to save below
 689                                 rabund->setLabel(saveLabel);
 690                         }
 691
 692                         lastLabel = rabund->getLabel();
 693
 694                         delete rabund;
 695                         rabund = input.getRAbundVector();
 696                 }
 697
 698                 if (m->control_pressed) {  out.close(); return 0; }
 699
 700                 //output error messages about any remaining user labels
 701                 set<string>::iterator it;
 702                 bool needToRun = false;
 703                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 704                         m->mothurOut("Your file does not include the label " + *it);
 705                         if (processedLabels.count(lastLabel) != 1) {
 706                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 707                                 needToRun = true;
 708                         }else {
 709                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 710                         }
 711                 }
 712
 713                 //run last label if you need to
 714                 if (needToRun == true)  {
 715                         if (rabund != NULL) {   delete rabund;  }
 716                         rabund = input.getRAbundVector(lastLabel);
 717
 718                         m->mothurOut(rabund->getLabel()); m->mothurOutEndLine();
 719
 720                         RAbundVector newRabund; newRabund.setLabel(rabund->getLabel());
 721                         for (int i = 0; i < rabund->getNumBins(); i++) {
 722                                 if (rabund->get(i) > nseqs) {
 723                                         newRabund.push_back(rabund->get(i));
 724                                 }
 725                         }
 726                         if (newRabund.getNumBins() > 0) { newRabund.print(out); }
 727
 728                         delete rabund;
 729                 }
 730
 731                 return 0;
 732         }
 733         catch(exception& e) {
 734                 m->errorOut(e, "RemoveRareCommand", "processRabund");
 735                 exit(1);
 736         }
 737 }
 738 //**********************************************************************************************************************
 739 int RemoveRareCommand::processShared(){
 740         try {
 741                 m->setGroups(Groups);
 742
 743                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 744                 InputData input(sharedfile, "sharedfile");
 745                 vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
 746                 string lastLabel = lookup[0]->getLabel();
 747                 set<string> processedLabels;
 748                 set<string> userLabels = labels;
 749
 750                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 751
 752                         if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }   return 0; }
 753
 754                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 755
 756                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 757                                 processedLabels.insert(lookup[0]->getLabel());
 758                                 userLabels.erase(lookup[0]->getLabel());
 759
 760                                 processLookup(lookup);
 761                         }
 762
 763                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 764                                 string saveLabel = lookup[0]->getLabel();
 765
 766                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 767                                 lookup = input.getSharedRAbundVectors(lastLabel);
 768
 769                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 770                                 processedLabels.insert(lookup[0]->getLabel());
 771                                 userLabels.erase(lookup[0]->getLabel());
 772
 773                                 processLookup(lookup);
 774
 775                                 //restore real lastlabel to save below
 776                                 lookup[0]->setLabel(saveLabel);
 777                         }
 778
 779                         lastLabel = lookup[0]->getLabel();
 780
 781                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 782                         lookup = input.getSharedRAbundVectors();
 783                 }
 784
 785                 if (m->control_pressed) {  return 0; }
 786
 787                 //output error messages about any remaining user labels
 788                 set<string>::iterator it;
 789                 bool needToRun = false;
 790                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 791                         m->mothurOut("Your file does not include the label " + *it);
 792                         if (processedLabels.count(lastLabel) != 1) {
 793                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 794                                 needToRun = true;
 795                         }else {
 796                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 797                         }
 798                 }
 799
 800                 //run last label if you need to
 801                 if (needToRun == true)  {
 802                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       }  }
 803                         lookup = input.getSharedRAbundVectors(lastLabel);
 804
 805                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 806                         processLookup(lookup);
 807
 808                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 809                 }
 810
 811                 return 0;
 812         }
 813         catch(exception& e) {
 814                 m->errorOut(e, "RemoveRareCommand", "processSabund");
 815                 exit(1);
 816         }
 817 }
 818 //**********************************************************************************************************************
 819 int RemoveRareCommand::processLookup(vector<SharedRAbundVector*>& lookup){
 820         try {
 821
 822         string thisOutputDir = outputDir;
 823                 if (outputDir == "") {  thisOutputDir += m->hasPath(sharedfile);  }
 824         map<string, string> variables;
 825         variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(sharedfile));
 826         variables["[extension]"] = m->getExtension(sharedfile);
 827         variables["[tag]"] = lookup[0]->getLabel();
 828                 string outputFileName = getOutputFileName("shared", variables);
 829                 outputTypes["shared"].push_back(outputFileName); outputNames.push_back(outputFileName);
 830
 831                 ofstream out;
 832                 m->openOutputFile(outputFileName, out);
 833
 834                 vector<SharedRAbundVector> newRabunds;  newRabunds.resize(lookup.size());
 835         vector<string> headers;
 836                 for (int i = 0; i < lookup.size(); i++) {
 837                         newRabunds[i].setGroup(lookup[i]->getGroup());
 838                         newRabunds[i].setLabel(lookup[i]->getLabel());
 839                 }
 840
 841                 if (byGroup) {
 842
 843                         //for each otu
 844                         for (int i = 0; i < lookup[0]->getNumBins(); i++) {
 845                                 bool allZero = true;
 846
 847                                 if (m->control_pressed) { out.close(); return 0; }
 848
 849                                 //for each group
 850                                 for (int j = 0; j < lookup.size(); j++) {
 851
 852                                         //are you rare?
 853                                         if (lookup[j]->getAbundance(i) > nseqs) {
 854                                                 newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
 855                                                 allZero = false;
 856                                         }else {
 857                                                 newRabunds[j].push_back(0, newRabunds[j].getGroup());
 858                                         }
 859                                 }
 860
 861                                 //eliminates zero otus
 862                                 if (allZero) { for (int j = 0; j < newRabunds.size(); j++) {  newRabunds[j].pop_back(); } }
 863                 else { headers.push_back(m->currentSharedBinLabels[i]); }
 864                         }
 865                 }else {
 866                         //for each otu
 867                         for (int i = 0; i < lookup[0]->getNumBins(); i++) {
 868
 869                                 if (m->control_pressed) { out.close(); return 0; }
 870
 871                                 int totalAbund = 0;
 872                                 //get total otu abundance
 873                                 for (int j = 0; j < lookup.size(); j++) {
 874                                         newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
 875                                         totalAbund += lookup[j]->getAbundance(i);
 876                                 }
 877
 878                                 //eliminates otus below rare cutoff
 879                                 if (totalAbund <= nseqs) { for (int j = 0; j < newRabunds.size(); j++) {  newRabunds[j].pop_back(); } }
 880                 else { headers.push_back(m->currentSharedBinLabels[i]); }
 881                         }
 882                 }
 883
 884                 //do we have any otus above the rare cutoff
 885                 if (newRabunds[0].getNumBins() != 0) {
 886             out << "label\tGroup\tnumOtus\t";
 887             for (int j = 0; j < headers.size(); j++) { out << headers[j] << '\t'; }
 888             out << endl;
 889                         for (int j = 0; j < newRabunds.size(); j++) {
 890                                 out << newRabunds[j].getLabel() << '\t' << newRabunds[j].getGroup() << '\t';
 891                                 newRabunds[j].print(out);
 892                         }
 893                 }
 894
 895         out.close();
 896
 897                 return 0;
 898         }
 899         catch(exception& e) {
 900                 m->errorOut(e, "RemoveRareCommand", "processLookup");
 901                 exit(1);
 902         }
 903 }
 904 //**********************************************************************************************************************
 905
 906
 907
 908