pipelinepdscommand.cpp

   1 /*
   2  *  pipelinepdscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 10/5/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "pipelinepdscommand.h"
  11 #include "sffinfocommand.h"
  12 #include "commandoptionparser.hpp"
  13
  14 //**********************************************************************************************************************
  15 vector<string> PipelineCommand::getValidParameters(){
  16         try {
  17                 string Array[] =  {"sff","oligos","align","chimera","classify","taxonomy","pipeline","processors","outputdir","inputdir"};
  18                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  19                 return myArray;
  20         }
  21         catch(exception& e) {
  22                 m->errorOut(e, "PipelineCommand", "getValidParameters");
  23                 exit(1);
  24         }
  25 }
  26 //**********************************************************************************************************************
  27 vector<string> PipelineCommand::getRequiredParameters(){
  28         try {
  29                 vector<string> myArray;
  30                 return myArray;
  31         }
  32         catch(exception& e) {
  33                 m->errorOut(e, "PipelineCommand", "getRequiredParameters");
  34                 exit(1);
  35         }
  36 }
  37 //**********************************************************************************************************************
  38 vector<string> PipelineCommand::getRequiredFiles(){
  39         try {
  40                 vector<string> myArray;
  41                 return myArray;
  42         }
  43         catch(exception& e) {
  44                 m->errorOut(e, "PipelineCommand", "getRequiredFiles");
  45                 exit(1);
  46         }
  47 }
  48 //**********************************************************************************************************************
  49 PipelineCommand::PipelineCommand(string option) {
  50         try {
  51                 cFactory = CommandFactory::getInstance();
  52                 abort = false;
  53
  54                 //allow user to run help
  55                 if(option == "help") { help(); abort = true; }
  56
  57                 else {
  58
  59                         //valid paramters for this command
  60                         string AlignArray[] =  {"sff","oligos","align","chimera","classify","taxonomy","pipeline","processors","outputdir","inputdir"};
  61                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  62
  63                         OptionParser parser(option);
  64                         map<string, string> parameters = parser.getParameters();
  65
  66                         ValidParameters validParameter;
  67                         map<string, string>::iterator it;
  68
  69                         //check to make sure all parameters are valid for command
  70                         for (it = parameters.begin(); it != parameters.end(); it++) {
  71                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  72                         }
  73
  74                         //if the user changes the input directory command factory will send this info to us in the output parameter
  75                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  76                         if (inputDir == "not found"){   inputDir = "";          }
  77                         else {
  78                                 string path;
  79                                 it = parameters.find("sff");
  80                                 //user has given a template file
  81                                 if(it != parameters.end()){
  82                                         path = m->hasPath(it->second);
  83                                         //if the user has not given a path then, add inputdir. else leave path alone.
  84                                         if (path == "") {       parameters["sff"] = inputDir + it->second;              }
  85                                 }
  86
  87                                 it = parameters.find("oligos");
  88                                 //user has given a template file
  89                                 if(it != parameters.end()){
  90                                         path = m->hasPath(it->second);
  91                                         //if the user has not given a path then, add inputdir. else leave path alone.
  92                                         if (path == "") {       parameters["oligos"] = inputDir + it->second;           }
  93                                 }
  94
  95                                 it = parameters.find("align");
  96                                 //user has given a template file
  97                                 if(it != parameters.end()){
  98                                         path = m->hasPath(it->second);
  99                                         //if the user has not given a path then, add inputdir. else leave path alone.
 100                                         if (path == "") {       parameters["align"] = inputDir + it->second;            }
 101                                 }
 102
 103                                 it = parameters.find("chimera");
 104                                 //user has given a template file
 105                                 if(it != parameters.end()){
 106                                         path = m->hasPath(it->second);
 107                                         //if the user has not given a path then, add inputdir. else leave path alone.
 108                                         if (path == "") {       parameters["chimera"] = inputDir + it->second;          }
 109                                 }
 110
 111                                 it = parameters.find("classify");
 112                                 //user has given a template file
 113                                 if(it != parameters.end()){
 114                                         path = m->hasPath(it->second);
 115                                         //if the user has not given a path then, add inputdir. else leave path alone.
 116                                         if (path == "") {       parameters["classify"] = inputDir + it->second;         }
 117                                 }
 118
 119                                 it = parameters.find("taxonomy");
 120                                 //user has given a template file
 121                                 if(it != parameters.end()){
 122                                         path = m->hasPath(it->second);
 123                                         //if the user has not given a path then, add inputdir. else leave path alone.
 124                                         if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
 125                                 }
 126
 127                                 it = parameters.find("pipeline");
 128                                 //user has given a template file
 129                                 if(it != parameters.end()){
 130                                         path = m->hasPath(it->second);
 131                                         //if the user has not given a path then, add inputdir. else leave path alone.
 132                                         if (path == "") {       parameters["pipeline"] = inputDir + it->second;         }
 133                                 }
 134                         }
 135
 136                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
 137
 138                         pipeFilename = validParameter.validFile(parameters, "pipeline", true);
 139                         if (pipeFilename == "not found") { pipeFilename = "";  }
 140                         else if (pipeFilename == "not open") { pipeFilename = ""; abort = true; }
 141
 142                         string temp = validParameter.validFile(parameters, "processors", false);        if (temp == "not found"){       temp = "1";                             }
 143                         convert(temp, processors);
 144
 145                         if (pipeFilename != "") {
 146                                 abort = readUsersPipeline();
 147                         }else{
 148                                 sffFile = validParameter.validFile(parameters, "sff", true);
 149                                 if (sffFile == "not found") { m->mothurOut("sff is a required parameter for the pipeline command."); m->mothurOutEndLine(); abort = true;  }
 150                                 else if (sffFile == "not open") { sffFile = ""; abort = true; }
 151
 152                                 oligosFile = validParameter.validFile(parameters, "oligos", true);
 153                                 if (oligosFile == "not found") { m->mothurOut("oligos is a required parameter for the pipeline command."); m->mothurOutEndLine(); abort = true;  }
 154                                 else if (oligosFile == "not open") { oligosFile = ""; abort = true; }
 155
 156                                 alignFile = validParameter.validFile(parameters, "align", true);
 157                                 if (alignFile == "not found") { m->mothurOut("align is a required parameter for the pipeline command. Please provide the template to align with."); m->mothurOutEndLine(); abort = true;  }
 158                                 else if (alignFile == "not open") { alignFile = ""; abort = true; }
 159
 160                                 chimeraFile = validParameter.validFile(parameters, "chimera", true);
 161                                 if (chimeraFile == "not found") { m->mothurOut("chimera is a required parameter for the pipeline command. Please provide the template to check for chimeras with."); m->mothurOutEndLine(); abort = true;  }
 162                                 else if (chimeraFile == "not open") { chimeraFile = ""; abort = true; }
 163
 164                                 classifyFile = validParameter.validFile(parameters, "classify", true);
 165                                 if (classifyFile == "not found") { m->mothurOut("classify is a required parameter for the pipeline command. Please provide the template to use with the classifier."); m->mothurOutEndLine(); abort = true;  }
 166                                 else if (classifyFile == "not open") { classifyFile = ""; abort = true; }
 167
 168                                 taxonomyFile = validParameter.validFile(parameters, "taxonomy", true);
 169                                 if (taxonomyFile == "not found") { m->mothurOut("taxonomy is a required parameter for the pipeline command."); m->mothurOutEndLine(); abort = true;  }
 170                                 else if (taxonomyFile == "not open") { taxonomyFile = ""; abort = true; }
 171                         }
 172                 }
 173
 174         }
 175         catch(exception& e) {
 176                 m->errorOut(e, "PipelineCommand", "PipelineCommand");
 177                 exit(1);
 178         }
 179 }
 180 //**********************************************************************************************************************
 181
 182 void PipelineCommand::help(){
 183         try {
 184                  m->mothurOut("The pipeline command is designed to guide you through your analysis using mothur.\n");
 185                  m->mothurOut("The pipeline command parameters are pipeline, sff, oligos, align, chimera, classify, taxonomy and processors.\n");
 186                  m->mothurOut("The sff parameter allows you to enter your sff file. It is required, if not using pipeline parameter.\n");
 187                  m->mothurOut("The oligos parameter allows you to enter your oligos file. It is required, if not using pipeline parameter.\n");
 188                  m->mothurOut("The align parameter allows you to enter a template to use with the aligner. It is required, if not using pipeline parameter.\n");
 189                  m->mothurOut("The chimera parameter allows you to enter a template to use for chimera detection. It is required, if not using pipeline parameter.\n");
 190                  m->mothurOut("The classify parameter allows you to enter a template to use for classification. It is required, if not using pipeline parameter.\n");
 191                  m->mothurOut("The taxonomy parameter allows you to enter a taxonomy file for the classify template to use for classification. It is required, if not using pipeline parameter.\n");
 192                  m->mothurOut("The processors parameter allows you to specify the number of processors to use. The default is 1.\n");
 193                  m->mothurOut("The pipeline parameter allows you to enter your own pipeline file. This file should look like a mothur batchfile, but where you would be using a mothur generated file, you can use mothurmade instead.\n");
 194                  m->mothurOut("Example: trim.seqs(processors=8, allfiles=T, maxambig=0, maxhomop=8, flip=T, bdiffs=1, pdiffs=2, qwindowaverage=35, qwindowsize=50, fasta=may1.v13.fasta, oligos=may1.v13.oligos, qfile=may1.v13.qual)\n");
 195                  m->mothurOut("then, you could enter unique.seqs(fasta=mothurmade), and mothur would use the .trim.fasta file from the trim.seqs command. \n");
 196                  m->mothurOut("then you could enter align.seqs(candidate=mothurmade, template=silva.v13.align, processors=8). , and mothur would use the .trim.unique.fasta file from the unique.seqs command. \n");
 197                  m->mothurOut("If no pipeline file is given then mothur will use Pat's pipeline. \n\n");
 198                  m->mothurOut("Here is a list of the commands used in Pat's pipeline.\n");
 199                  m->mothurOut("All paralellized commands will use the processors you entered.\n");
 200                  m->mothurOut("The sffinfo command takes your sff file and extracts the fasta and quality files.\n");
 201                  m->mothurOut("The trim.seqs command uses your oligos file and the quality and fasta files generated by sffinfo.\n");
 202                  m->mothurOut("The trim.seqs command sets the following parameters: allfiles=T, maxambig=0, maxhomop=8, flip=T, bdiffs=1, pdiffs=2, qwindowaverage=35, qwindowsize=50.\n");
 203                  m->mothurOut("The unique.seqs command uses the trimmed fasta file and removes redundant sequences, don't worry the names file generated by unique.seqs will be used in the pipeline to make sure they are included.\n");
 204                  m->mothurOut("The align.seqs command aligns the unique sequences using the aligners default options. \n");
 205                  m->mothurOut("The screen.seqs command screens the sequences using optimize=end-minlength. \n");
 206                  m->mothurOut("The pipeline uses chimera.slayer to detect chimeras using the default options. \n");
 207                  m->mothurOut("The pipeline removes all sequences determined to be chimeric by chimera.slayer. \n");
 208                  m->mothurOut("The filter.seqs command filters the sequences using vertical=T, trump=. \n");
 209                  m->mothurOut("The unique.seqs command uses the filtered fasta file and name file to remove sequences that have become redundant after filtering.\n");
 210                  m->mothurOut("The pre.cluster command clusters sequences that have no more than 2 differences.\n");
 211                  m->mothurOut("The dist.seqs command is used to generate a column and phylip formatted distance matrix using cutoff=0.20 for column.\n");
 212                  m->mothurOut("The pipeline uses cluster with method=average, hard=T. \n");
 213                  m->mothurOut("The classify.seqs command is used to classify the sequences using the bayesian method with a cutoff of 80.\n");
 214                  m->mothurOut("The phylotype command is used to cluster the sequences based on their classification.\n");
 215                  m->mothurOut("The clearcut command is used to generate a tree using neighbor=T. \n");
 216                  m->mothurOut("The summary.single and summary.shared commands are run on the otu files from cluster and phylotype commands. \n");
 217                  m->mothurOut("The summary.shared command uses calc=sharednseqs-sharedsobs-sharedchao-sharedace-anderberg-jclass-jest-kulczynski-kulczynskicody-lennon-ochiai-sorclass-sorest-whittaker-braycurtis-jabund-morisitahorn-sorabund-thetan-thetayc. \n");
 218                  m->mothurOut("The summary.single command uses calc=nseqs-sobs-coverage-bergerparker-chao-ace-jack-bootstrap-boneh-efron-shen-solow-shannon-npshannon-invsimpson-qstat-simpsoneven-shannoneven-heip-smithwilson. \n");
 219                  m->mothurOut("The classify.otu command is used to get the concensus taxonomy for otu files from cluster and phylotype commands. \n");
 220                  m->mothurOut("The phylo.diversity command run on the tree generated by clearcut with rarefy=T, iters=100. \n");
 221                  m->mothurOut("The unifrac commands are also run on the tree generated by clearcut with random=F, distance=T. \n");
 222                  m->mothurOut("\n\n");
 223         }
 224         catch(exception& e) {
 225                 m->errorOut(e, "PipelineCommand", "help");
 226                 exit(1);
 227         }
 228 }
 229
 230 //**********************************************************************************************************************
 231
 232 PipelineCommand::~PipelineCommand(){}
 233
 234 //**********************************************************************************************************************
 235
 236 int PipelineCommand::execute(){
 237         try {
 238                 if (abort == true) { return 0; }
 239
 240                 int start = time(NULL);
 241
 242                 if (pipeFilename == "") {
 243                         createPatsPipeline();
 244
 245                         //run Pats pipeline
 246                         for (int i = 0; i < commands.size(); i++) {
 247                                 m->mothurOutEndLine(); m->mothurOut("mothur > " + commands[i]); m->mothurOutEndLine();
 248
 249                                 if (m->control_pressed) { return 0; }
 250
 251                                 CommandOptionParser parser(commands[i]);
 252                                 string commandName = parser.getCommandString();
 253                                 string options = parser.getOptionString();
 254
 255                                 #ifdef USE_MPI
 256                                         int pid;
 257                                         MPI_Comm_rank(MPI_COMM_WORLD, &pid);
 258
 259                                         if ((cFactory->MPIEnabled(commandName)) || (pid == 0)) {
 260                                 #endif
 261
 262                                 //executes valid command
 263                                 Command* command = cFactory->getCommand(commandName, options, "pipe");
 264                                 command->execute();
 265
 266                                 //add output files to list
 267                                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 268                                 map<string, vector<string> >::iterator itMade;
 269                                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 270                                         vector<string> temp = itMade->second;
 271                                         for (int j = 0; j < temp.size(); j++) { outputNames.push_back(temp[j]); }
 272                                 }
 273
 274                                 #ifdef USE_MPI
 275                                         }
 276                                 #endif
 277                         }
 278
 279                 }else {  runUsersPipeline(); }
 280
 281                 if (m->control_pressed) { return 0; }
 282
 283                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run the pipeline analysis."); m->mothurOutEndLine(); m->mothurOutEndLine();
 284
 285                 m->mothurOutEndLine();
 286                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 287                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 288                 m->mothurOutEndLine();
 289
 290                 return 0;
 291         }
 292         catch(exception& e) {
 293                 m->errorOut(e, "PipelineCommand", "execute");
 294                 exit(1);
 295         }
 296 }
 297 //**********************************************************************************************************************
 298
 299 bool PipelineCommand::readUsersPipeline(){
 300         try {
 301
 302                 ifstream in;
 303                 m->openInputFile(pipeFilename, in);
 304
 305                 string nextCommand = "";
 306
 307                 map<string, vector<string> > mothurMadeFiles;
 308
 309                 while(!in.eof()) {
 310                         nextCommand = m->getline(in); m->gobble(in);
 311
 312                         if (nextCommand[0] != '#') {
 313                                 bool error = false;
 314
 315                                 string commandName, options;
 316                                 error = parseCommand(nextCommand, commandName, options);
 317
 318                                 if (error) { in.close(); return error; }
 319                                 if (commandName == "pipeline.pds") { m->mothurOut("Cannot run the pipeline.pds command from inside the pipeline.pds command."); m->mothurOutEndLine(); in.close(); return true; }
 320
 321                                 error = checkForValidAndRequiredParameters(commandName, options, mothurMadeFiles);
 322
 323                                 if (error) { in.close(); return error; }
 324                         }
 325                 }
 326
 327                 in.close();
 328
 329                 return false;
 330         }
 331         catch(exception& e) {
 332                 m->errorOut(e, "PipelineCommand", "readUsersPipeline");
 333                 exit(1);
 334         }
 335 }
 336 //**********************************************************************************************************************
 337
 338 bool PipelineCommand::parseCommand(string nextCommand, string& name, string& options){
 339         try {
 340                 CommandOptionParser parser(nextCommand);
 341                 name = parser.getCommandString();
 342                 options = parser.getOptionString();
 343
 344                 if (name == "") { return true; } //name == "" if () are not right
 345
 346                 return false;
 347         }
 348         catch(exception& e) {
 349                 m->errorOut(e, "PipelineCommand", "parseCommand");
 350                 exit(1);
 351         }
 352 }
 353 //**********************************************************************************************************************
 354
 355 bool PipelineCommand::checkForValidAndRequiredParameters(string name, string options, map<string, vector<string> >& mothurMadeFiles){
 356         try {
 357
 358                 if (name == "system") { return false; }
 359
 360                 if (name == "system") { return false; }
 361
 362                 //get shell of the command so we can check to make sure its valid without running it
 363                 Command* command = cFactory->getCommand(name);
 364
 365                 //check to make sure all parameters are valid for command
 366                 vector<string> validParameters = command->getValidParameters();
 367
 368                 OptionParser parser(options);
 369                 map<string, string> parameters = parser.getParameters();
 370
 371                 ValidParameters validParameter;
 372                 map<string, string>::iterator it;
 373                 map<string, vector<string> >::iterator itMade;
 374
 375                 for (it = parameters.begin(); it != parameters.end(); it++) {
 376
 377                         if (validParameter.isValidParameter(it->first, validParameters, it->second) != true) {  return true;  } // not valid
 378                         if (it->second == "mothurmade") {
 379                                 itMade = mothurMadeFiles.find(it->first);
 380
 381                                 if (itMade == mothurMadeFiles.end()) {
 382                                         if ((name == "align.seqs") && (it->first == "candidate")) {} //do nothing about candidate
 383                                         else {
 384                                                 m->mothurOut("You have the " + it->first + " listed as a mothurmade file for the " + name + " command, but it seems mothur will not make that file in your current pipeline, please correct."); m->mothurOutEndLine();
 385                                                 return true;
 386                                         }
 387                                 }
 388                         }
 389                 }
 390
 391                 //is the command missing any required
 392                 vector<string> requiredParameters = command->getRequiredParameters();
 393
 394                 //check for or
 395                 bool hasOr = false;
 396                 int numFound = 0;
 397                 if (requiredParameters.size() > 2) {
 398                         if (requiredParameters[(requiredParameters.size()-1)] == "or") { hasOr = true; }
 399                 }
 400
 401                 for (int i = 0; i < requiredParameters.size(); i++) {
 402                         it = parameters.find(requiredParameters[i]);
 403
 404                         if (it != parameters.end()) { numFound++; }
 405                         else {
 406                                 if (!hasOr) { m->mothurOut(name + " requires the " + requiredParameters[i] + " parameter, please correct."); m->mothurOutEndLine(); }
 407                         }
 408                 }
 409
 410                 // if all are needed and not all are found
 411                 if ((!hasOr) && (numFound != requiredParameters.size())) { return true; }
 412                 //if one is needed and none are found
 413                 else if ((hasOr) && (numFound == 0)) { return true; }
 414
 415                 //update MothurMade
 416                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 417                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 418                         mothurMadeFiles[itMade->first] = itMade->second; //adds any new types
 419                 }
 420
 421                 return false;
 422         }
 423         catch(exception& e) {
 424                 m->errorOut(e, "PipelineCommand", "checkForValidAndRequiredParameters");
 425                 exit(1);
 426         }
 427 }
 428 //**********************************************************************************************************************
 429 int PipelineCommand::runUsersPipeline(){
 430         try {
 431                 ifstream in;
 432                 m->openInputFile(pipeFilename, in);
 433
 434                 string nextCommand = "";
 435
 436                 map<string, vector<string> > mothurMadeFiles;
 437
 438                 while(!in.eof()) {
 439                         nextCommand = m->getline(in);  m->gobble(in);
 440
 441                         if (nextCommand[0] != '#') {
 442                                 CommandOptionParser parser(nextCommand);
 443                                 string commandName = parser.getCommandString();
 444                                 string options = parser.getOptionString();
 445
 446                                 if ((options != "") && (commandName != "system")) {
 447                                         bool error = fillInMothurMade(options, mothurMadeFiles);
 448                                         if (error) { in.close(); return 0; }
 449                                 }
 450
 451                                 m->mothurOutEndLine(); m->mothurOut("mothur > " + commandName + "(" + options + ")"); m->mothurOutEndLine();
 452
 453                                 if (m->control_pressed) { return 0; }
 454
 455                                 #ifdef USE_MPI
 456                                         int pid;
 457                                         MPI_Comm_rank(MPI_COMM_WORLD, &pid);
 458
 459                                         if ((cFactory->MPIEnabled(commandName)) || (pid == 0)) {
 460                                 #endif
 461
 462                                 //executes valid command
 463                                 Command* command = cFactory->getCommand(commandName, options, "pipe");
 464                                 command->execute();
 465
 466                                 //add output files to list
 467                                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 468                                 map<string, vector<string> >::iterator itMade;
 469                                 map<string, vector<string> >::iterator it;
 470                                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 471
 472                                         vector<string> temp = itMade->second;
 473                                         for (int k = 0; k < temp.size(); k++) { outputNames.push_back(temp[k]); }  //
 474
 475                                         //update Mothur Made for each file
 476                                         it = mothurMadeFiles.find(itMade->first);
 477
 478                                         if (it == mothurMadeFiles.end()) { //new type
 479
 480                                                 mothurMadeFiles[itMade->first] = temp;
 481
 482                                         }else{ //update existing type
 483                                                 vector<string> oldFileNames = it->second;
 484                                                 //look at new files, see if an old version of the file exists, if so update, else just add.
 485                                                 //for example you may have abrecovery.fasta and amazon.fasta as old files and you created a new amazon.trim.fasta.
 486
 487                                                 for (int k = 0; k < temp.size(); k++) {
 488
 489                                                         //get base name
 490                                                         string root = m->getSimpleName(temp[k]);
 491                                                         string individual = "";
 492                                                         for(int i=0;i<root.length();i++){
 493                                                                 if(root[i] == '.'){
 494                                                                         root = individual;
 495                                                                         break;
 496                                                                 }else{
 497                                                                         individual += root[i];
 498                                                                 }
 499                                                         }
 500
 501                                                         //look for that base name in oldfiles
 502                                                         int spot = -1;
 503                                                         for (int l = 0; l < oldFileNames.size(); l++) {
 504                                                                 int pos = oldFileNames[l].find(root);
 505                                                                 if (pos != string::npos) {
 506                                                                         spot = l;
 507                                                                         break;
 508                                                                 }
 509                                                         }
 510
 511                                                         //if you found it update it, else add it
 512                                                         if (spot != -1) {
 513                                                                 mothurMadeFiles[it->first][spot] = temp[k];
 514                                                         }else{
 515                                                                 mothurMadeFiles[it->first].push_back(temp[k]);
 516                                                         }
 517                                                 }
 518                                         }
 519                                 }
 520
 521                                 #ifdef USE_MPI
 522                                         }
 523                                 #endif
 524                         }
 525                 }
 526
 527                 in.close();
 528         }
 529         catch(exception& e) {
 530                 m->errorOut(e, "PipelineCommand", "runUsersPipeline");
 531                 exit(1);
 532         }
 533 }
 534 //**********************************************************************************************************************
 535 bool PipelineCommand::fillInMothurMade(string& options, map<string, vector<string> >& mothurMadeFiles){
 536         try {
 537
 538                 OptionParser parser(options);
 539                 map<string, string> parameters = parser.getParameters();
 540                 map<string, string>::iterator it;
 541                 map<string, vector<string> >::iterator itMade;
 542
 543                 options = "";
 544
 545                 //fill in mothurmade filenames
 546                 for (it = parameters.begin(); it != parameters.end(); it++) {
 547                         string paraType = it->first;
 548                         string tempOption = it->second;
 549
 550                         if (tempOption == "mothurmade") {
 551
 552                                 if (it->first == "candidate") { paraType = "fasta"; }
 553
 554                                 itMade = mothurMadeFiles.find(paraType);
 555
 556                                 if (itMade == mothurMadeFiles.end()) {
 557                                         m->mothurOut("Looking for a mothurmade " + paraType + " file, but it seems mothur has not made that file type in your current pipeline, please correct."); m->mothurOutEndLine();
 558                                         return true;
 559                                 }else{
 560                                         vector<string> temp = itMade->second;
 561
 562                                         if (temp.size() > 1) {
 563                                                 //ask user which file to use
 564                                                 m->mothurOut("More than one file has been created for the " + paraType + " parameter. "); m->mothurOutEndLine();
 565                                                 for (int i = 0; i < temp.size(); i++) {
 566                                                         m->mothurOut(toString(i) + " - " + temp[i]); m->mothurOutEndLine();
 567                                                 }
 568
 569                                                 m->mothurOut("Please select the number of the file you would like to use: ");
 570                                                 int num = 0;
 571                                                 cin >> num;
 572                                                 m->mothurOutJustToLog(toString(num)); m->mothurOutEndLine();
 573
 574                                                 if ((num < 0) || (num > (temp.size()-1))) { m->mothurOut("Not a valid response, quitting."); m->mothurOutEndLine(); return true; }
 575                                                 else {
 576                                                         tempOption = temp[num];
 577                                                 }
 578
 579                                                 //clears buffer so next command doesn't have error
 580                                                 string s;
 581                                                 getline(cin, s);
 582
 583                                                 vector<string> newTemp;
 584                                                 for (int i = 0; i < temp.size(); i++) {
 585                                                         if (i == num) { newTemp.push_back(temp[i]); }
 586                                                         else {
 587                                                                 m->mothurOut("Would you like to remove " + temp[i] + " as an option for " + paraType + ", (y/n): "); m->mothurOutEndLine();
 588                                                                 string response;
 589                                                                 cin >> response;
 590                                                                 m->mothurOutJustToLog(response); m->mothurOutEndLine();
 591
 592                                                                 if (response == "n") {  newTemp.push_back(temp[i]); }
 593
 594                                                                 //clears buffer so next command doesn't have error
 595                                                                 string s;
 596                                                                 getline(cin, s);
 597                                                         }
 598                                                 }
 599
 600                                                 mothurMadeFiles[paraType] = newTemp;
 601
 602
 603                                         }else if (temp.size() == 0){
 604                                                 m->mothurOut("Sorry, we seem to think you created a " + paraType + " file, but it seems mothur doesn't have a filename."); m->mothurOutEndLine();
 605                                                 return true;
 606                                         }else{
 607                                                 tempOption = temp[0];
 608                                         }
 609                                 }
 610                         }
 611
 612                         options += it->first + "=" + tempOption + ", ";
 613                 }
 614
 615                 //rip off extra comma
 616                 options = options.substr(0, (options.length()-2));
 617
 618                 return false;
 619         }
 620         catch(exception& e) {
 621                 m->errorOut(e, "PipelineCommand", "fillInMothurMade");
 622                 exit(1);
 623         }
 624 }
 625
 626 //**********************************************************************************************************************
 627 void PipelineCommand::createPatsPipeline(){
 628         try {
 629
 630                 //sff.info command
 631                 string thisCommand = "sffinfo(sff=" + sffFile + ")";
 632                 commands.push_back(thisCommand);
 633
 634                 //trim.seqs command
 635                 string fastaFile = m->getRootName(m->getSimpleName(sffFile)) + "fasta";
 636                 string qualFile = m->getRootName(m->getSimpleName(sffFile)) + "qual";
 637                 thisCommand = "trim.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", allfiles=T, maxambig=0, maxhomop=8, flip=T, bdiffs=1, pdiffs=2, qwindowaverage=35, qwindowsize=50, oligos=" + oligosFile + ", qfile=" + qualFile + ")";
 638                 commands.push_back(thisCommand);
 639
 640                 //unique.seqs
 641                 string groupFile = m->getRootName(m->getSimpleName(fastaFile)) + "groups";
 642                 qualFile =  m->getRootName(m->getSimpleName(fastaFile)) + "trim.qual";
 643                 fastaFile =  m->getRootName(m->getSimpleName(fastaFile)) + "trim.fasta";
 644                 thisCommand = "unique.seqs(fasta=" + fastaFile + ")";
 645                 commands.push_back(thisCommand);
 646
 647                 //align.seqs
 648                 string nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "names";
 649                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "unique" + m->getExtension(fastaFile);
 650                 thisCommand = "align.seqs(processors=" + toString(processors) + ", candidate=" + fastaFile + ", template=" + alignFile + ")";
 651                 commands.push_back(thisCommand);
 652
 653                 //screen.seqs
 654                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "align";
 655                 thisCommand = "screen.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", name=" + nameFile + ", group=" + groupFile + ", optimize=end-minlength)";
 656                 commands.push_back(thisCommand);
 657
 658                 //chimera.slayer
 659                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "good" + m->getExtension(fastaFile);
 660                 nameFile = m->getRootName(m->getSimpleName(nameFile)) + "good" + m->getExtension(nameFile);
 661                 groupFile = m->getRootName(m->getSimpleName(groupFile)) + "good" + m->getExtension(groupFile);
 662                 thisCommand = "chimera.slayer(processors=" + toString(processors) + ", fasta=" + fastaFile + ", template=" + chimeraFile + ")";
 663                 commands.push_back(thisCommand);
 664
 665                 //remove.seqs
 666                 string accnosFile = m->getRootName(m->getSimpleName(fastaFile))  + "slayer.accnos";
 667                 thisCommand = "remove.seqs(fasta=" + fastaFile + ", name=" + nameFile + ", group=" + groupFile + ", accnos=" + accnosFile + ", dups=T)";
 668                 commands.push_back(thisCommand);
 669
 670                 //filter.seqs
 671                 nameFile = m->getRootName(m->getSimpleName(nameFile)) + "pick" + m->getExtension(nameFile);
 672                 groupFile = m->getRootName(m->getSimpleName(groupFile)) + "pick" + m->getExtension(groupFile);
 673                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "pick" + m->getExtension(fastaFile);
 674                 thisCommand = "filter.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", vertical=T, trump=.)";
 675                 commands.push_back(thisCommand);
 676
 677                 //unique.seqs
 678                 fastaFile =  m->getRootName(m->getSimpleName(fastaFile)) + "filter.fasta";
 679                 thisCommand = "unique.seqs(fasta=" + fastaFile + ", name=" + nameFile + ")";
 680                 commands.push_back(thisCommand);
 681
 682                 //pre.cluster
 683                 nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "names";
 684                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "unique" + m->getExtension(fastaFile);
 685                 thisCommand = "pre.cluster(fasta=" + fastaFile + ", name=" + nameFile + ", diffs=2)";
 686                 commands.push_back(thisCommand);
 687
 688                 //dist.seqs
 689                 nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "precluster.names";
 690                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "precluster" + m->getExtension(fastaFile);
 691                 thisCommand = "dist.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", cutoff=0.20)";
 692                 commands.push_back(thisCommand);
 693
 694                 //dist.seqs
 695                 string columnFile = m->getRootName(m->getSimpleName(fastaFile)) + "dist";
 696                 thisCommand = "dist.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", output=lt)";
 697                 commands.push_back(thisCommand);
 698
 699                 //read.dist
 700                 string phylipFile = m->getRootName(m->getSimpleName(fastaFile)) + "phylip.dist";
 701                 thisCommand = "read.dist(column=" + columnFile + ", name=" + nameFile + ")";
 702                 commands.push_back(thisCommand);
 703
 704                 //cluster
 705                 thisCommand = "cluster(method=average, hard=T)";
 706                 commands.push_back(thisCommand);
 707
 708                 string listFile = m->getRootName(m->getSimpleName(columnFile)) + "an.list";
 709                 string rabundFile = m->getRootName(m->getSimpleName(columnFile)) + "an.rabund";
 710
 711                 //degap.seqs
 712                 thisCommand = "degap.seqs(fasta=" + fastaFile + ")";
 713                 commands.push_back(thisCommand);
 714
 715                 //classify.seqs
 716                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "ng.fasta";
 717                 thisCommand = "classify.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", name=" + nameFile + ", template=" + classifyFile + ", taxonomy=" + taxonomyFile + ", cutoff=80)";
 718                 commands.push_back(thisCommand);
 719
 720                 string RippedTaxName = m->getRootName(m->getSimpleName(taxonomyFile));
 721                 RippedTaxName = m->getExtension(RippedTaxName.substr(0, RippedTaxName.length()-1));
 722                 if (RippedTaxName[0] == '.') { RippedTaxName = RippedTaxName.substr(1, RippedTaxName.length()); }
 723                 RippedTaxName +=  ".";
 724
 725                 string fastaTaxFile = m->getRootName(m->getSimpleName(fastaFile)) + RippedTaxName + "taxonomy";
 726                 string taxSummaryFile = m->getRootName(m->getSimpleName(fastaFile)) + RippedTaxName + "tax.summary";
 727
 728                 //phylotype
 729                 thisCommand = "phylotype(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ")";
 730                 commands.push_back(thisCommand);
 731
 732                 string phyloListFile = m->getRootName(m->getSimpleName(fastaTaxFile)) + "tx.list";
 733                 string phyloRabundFile = m->getRootName(m->getSimpleName(fastaTaxFile)) + "tx.rabund";
 734
 735                 //clearcut
 736                 thisCommand = "clearcut(phylip=" + phylipFile + ", neighbor=T)";
 737                 commands.push_back(thisCommand);
 738
 739                 string treeFile = m->getRootName(m->getSimpleName(phylipFile)) + "tre";
 740
 741                 //read.otu
 742                 thisCommand = "read.otu(list=" + listFile + ", group=" + groupFile + ", label=0.03)";
 743                 commands.push_back(thisCommand);
 744
 745                 string sharedFile = m->getRootName(m->getSimpleName(listFile)) + "shared";
 746
 747                 //read.otu
 748                 thisCommand = "read.otu(list=" + phyloListFile + ", group=" + groupFile + ", label=1)";
 749                 commands.push_back(thisCommand);
 750
 751                 string phyloSharedFile = m->getRootName(m->getSimpleName(phyloListFile)) + "shared";
 752
 753                 //read.otu
 754                 thisCommand = "read.otu(shared=" + sharedFile + ")";
 755                 commands.push_back(thisCommand);
 756
 757                 //summary.single
 758                 thisCommand = "summary.single(calc=nseqs-sobs-coverage-bergerparker-chao-ace-jack-bootstrap-boneh-efron-shen-solow-shannon-npshannon-invsimpson-qstat-simpsoneven-shannoneven-heip-smithwilson, size=5000)";
 759                 commands.push_back(thisCommand);
 760
 761                 //summary.shared
 762                 thisCommand = "summary.shared(calc=sharednseqs-sharedsobs-sharedchao-sharedace-anderberg-jclass-jest-kulczynski-kulczynskicody-lennon-ochiai-sorclass-sorest-whittaker-braycurtis-jabund-morisitahorn-sorabund-thetan-thetayc)";
 763                 commands.push_back(thisCommand);
 764
 765                 //read.otu
 766                 thisCommand = "read.otu(rabund=" + rabundFile + ", label=0.03)";
 767                 commands.push_back(thisCommand);
 768
 769                 //summary.single
 770                 thisCommand = "summary.single(calc=nseqs-sobs-coverage-bergerparker-chao-ace-jack-bootstrap-boneh-efron-shen-solow-shannon-npshannon-invsimpson-qstat-simpsoneven-shannoneven-heip-smithwilson, size=5000)";
 771                 commands.push_back(thisCommand);
 772
 773                 //read.otu
 774                 thisCommand = "read.otu(shared=" + phyloSharedFile + ")";
 775                 commands.push_back(thisCommand);
 776
 777                 //summary.single
 778                 thisCommand = "summary.single(calc=nseqs-sobs-coverage-bergerparker-chao-ace-jack-bootstrap-boneh-efron-shen-solow-shannon-npshannon-invsimpson-qstat-simpsoneven-shannoneven-heip-smithwilson, size=5000)";
 779                 commands.push_back(thisCommand);
 780
 781                 //summary.shared
 782                 thisCommand = "summary.shared(calc=sharednseqs-sharedsobs-sharedchao-sharedace-anderberg-jclass-jest-kulczynski-kulczynskicody-lennon-ochiai-sorclass-sorest-whittaker-braycurtis-jabund-morisitahorn-sorabund-thetan-thetayc)";
 783                 commands.push_back(thisCommand);
 784
 785                 //read.otu
 786                 thisCommand = "read.otu(rabund=" + phyloRabundFile + ", label=1)";
 787                 commands.push_back(thisCommand);
 788
 789                 //summary.single
 790                 thisCommand = "summary.single(calc=nseqs-sobs-coverage-bergerparker-chao-ace-jack-bootstrap-boneh-efron-shen-solow-shannon-npshannon-invsimpson-qstat-simpsoneven-shannoneven-heip-smithwilson, size=5000)";
 791                 commands.push_back(thisCommand);
 792
 793                 //classify.otu
 794                 thisCommand = "classify.otu(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ", list=" + listFile + ", cutoff=51, label=0.03)";
 795                 commands.push_back(thisCommand);
 796
 797                 //classify.otu
 798                 thisCommand = "classify.otu(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ", list=" + phyloListFile + ", cutoff=51, label=1)";
 799                 commands.push_back(thisCommand);
 800
 801                 //read.tree
 802                 thisCommand = "read.tree(tree=" + treeFile + ", name=" + nameFile + ", group=" + groupFile + ")";
 803                 commands.push_back(thisCommand);
 804
 805                 //phylo.diversity
 806                 thisCommand = "phylo.diversity(iters=100,rarefy=T)";
 807                 commands.push_back(thisCommand);
 808
 809                 //unifrac.weighted
 810                 thisCommand = "unifrac.weighted(random=false, distance=true, groups=all, processors=" + toString(processors) + ")";
 811                 commands.push_back(thisCommand);
 812
 813                 //unifrac.unweighted
 814                 thisCommand = "unifrac.unweighted(random=false, distance=true, processors=" + toString(processors) + ")";
 815                 commands.push_back(thisCommand);
 816
 817
 818         }
 819         catch(exception& e) {
 820                 m->errorOut(e, "PipelineCommand", "createPatsPipeline");
 821                 exit(1);
 822         }
 823 }
 824
 825 //**********************************************************************************************************************