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.\n");
 187                  m->mothurOut("The oligos parameter allows you to enter your oligos file. It is required.\n");
 188                  m->mothurOut("The align parameter allows you to enter a template to use with the aligner. It is required.\n");
 189                  m->mothurOut("The chimera parameter allows you to enter a template to use for chimera detection. It is required.\n");
 190                  m->mothurOut("The classify parameter allows you to enter a template to use for classification. It is required.\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.\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("First column contains the command name, and the second column contains the parameter options or 'defaults', meaning use defaults. You may leave out file options.\n");
 195                  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");
 196                  m->mothurOut("then, you could enter unique.seqs(fasta=mothurmade), and mothur would use the .trim.fasta file from the trim.seqs command. \n");
 197                  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");
 198                  m->mothurOut("If no pipeline file is given then mothur will use Pat's pipeline. \n\n");
 199                  m->mothurOut("Here is a list of the commands used in Pat's pipeline.\n");
 200                  m->mothurOut("All paralellized commands will use the processors you entered.\n");
 201                  m->mothurOut("The sffinfo command takes your sff file and extracts the fasta and quality files.\n");
 202                  m->mothurOut("The trim.seqs command uses your oligos file and the quality and fasta files generated by sffinfo.\n");
 203                  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");
 204                  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");
 205                  m->mothurOut("The align.seqs command aligns the unique sequences using the aligners default options. \n");
 206                  m->mothurOut("The screen.seqs command screens the sequences using optimize=end-minlength. \n");
 207                  m->mothurOut("The pipeline uses chimera.slayer to detect chimeras using the default options. \n");
 208                  m->mothurOut("The pipeline removes all sequences determined to be chimeric by chimera.slayer. \n");
 209                  m->mothurOut("The filter.seqs command filters the sequences using vertical=T, trump=. \n");
 210                  m->mothurOut("The unique.seqs command uses the filtered fasta file and name file to remove sequences that have become redundant after filtering.\n");
 211                  m->mothurOut("The pre.cluster command clusters sequences that have no more than 2 differences.\n");
 212                  m->mothurOut("The dist.seqs command is used to generate a column and phylip formatted distance matrix using cutoff=0.20 for column.\n");
 213                  m->mothurOut("The pipeline uses cluster with method=average, hard=T. \n");
 214                  m->mothurOut("The classify.seqs command is used to classify the sequences using the bayesian method with a cutoff of 80.\n");
 215                  m->mothurOut("The phylotype command is used to cluster the sequences based on their classification.\n");
 216                  m->mothurOut("The clearcut command is used to generate a tree using neighbor=T. \n");
 217                  m->mothurOut("The summary.single and summary.shared commands are run on the otu files from cluster and phylotype commands. \n");
 218                  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");
 219                  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");
 220                  m->mothurOut("The classify.otu command is used to get the concensus taxonomy for otu files from cluster and phylotype commands. \n");
 221                  m->mothurOut("The phylo.diversity command run on the tree generated by clearcut with rarefy=T, iters=100. \n");
 222                  m->mothurOut("The unifrac commands are also run on the tree generated by clearcut with random=F, distance=T. \n");
 223                  m->mothurOut("\n\n");
 224         }
 225         catch(exception& e) {
 226                 m->errorOut(e, "PipelineCommand", "help");
 227                 exit(1);
 228         }
 229 }
 230
 231 //**********************************************************************************************************************
 232
 233 PipelineCommand::~PipelineCommand(){}
 234
 235 //**********************************************************************************************************************
 236
 237 int PipelineCommand::execute(){
 238         try {
 239                 if (abort == true) { return 0; }
 240
 241                 int start = time(NULL);
 242
 243                 if (pipeFilename == "") {
 244                         createPatsPipeline();
 245
 246                         //run Pats pipeline
 247                         for (int i = 0; i < commands.size(); i++) {
 248                                 m->mothurOutEndLine(); m->mothurOut("mothur > " + commands[i]); m->mothurOutEndLine();
 249
 250                                 if (m->control_pressed) { return 0; }
 251
 252                                 CommandOptionParser parser(commands[i]);
 253                                 string commandName = parser.getCommandString();
 254                                 string options = parser.getOptionString();
 255
 256                                 #ifdef USE_MPI
 257                                         int pid;
 258                                         MPI_Comm_rank(MPI_COMM_WORLD, &pid);
 259
 260                                         if ((cFactory->MPIEnabled(commandName)) || (pid == 0)) {
 261                                 #endif
 262
 263                                 //executes valid command
 264                                 Command* command = cFactory->getCommand(commandName, options, "pipe");
 265                                 command->execute();
 266
 267                                 //add output files to list
 268                                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 269                                 map<string, vector<string> >::iterator itMade;
 270                                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 271                                         vector<string> temp = itMade->second;
 272                                         for (int j = 0; j < temp.size(); j++) { outputNames.push_back(temp[j]); }
 273                                 }
 274
 275                                 #ifdef USE_MPI
 276                                         }
 277                                 #endif
 278                         }
 279
 280                 }else {  runUsersPipeline(); }
 281
 282                 if (m->control_pressed) { return 0; }
 283
 284                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run the pipeline analysis."); m->mothurOutEndLine(); m->mothurOutEndLine();
 285
 286                 m->mothurOutEndLine();
 287                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 288                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 289                 m->mothurOutEndLine();
 290
 291                 return 0;
 292         }
 293         catch(exception& e) {
 294                 m->errorOut(e, "PipelineCommand", "execute");
 295                 exit(1);
 296         }
 297 }
 298 //**********************************************************************************************************************
 299
 300 bool PipelineCommand::readUsersPipeline(){
 301         try {
 302
 303                 ifstream in;
 304                 m->openInputFile(pipeFilename, in);
 305
 306                 string nextCommand = "";
 307
 308                 map<string, vector<string> > mothurMadeFiles;
 309
 310                 while(!in.eof()) {
 311                         nextCommand = m->getline(in); m->gobble(in);
 312
 313                         if (nextCommand[0] != '#') {
 314                                 bool error = false;
 315
 316                                 string commandName, options;
 317                                 error = parseCommand(nextCommand, commandName, options);
 318
 319                                 if (error) { in.close(); return error; }
 320                                 if (commandName == "pipeline.pds") { m->mothurOut("Cannot run the pipeline.pds command from inside the pipeline.pds command."); m->mothurOutEndLine(); in.close(); return true; }
 321
 322                                 error = checkForValidAndRequiredParameters(commandName, options, mothurMadeFiles);
 323
 324                                 if (error) { in.close(); return error; }
 325                         }
 326                 }
 327
 328                 in.close();
 329
 330                 return false;
 331         }
 332         catch(exception& e) {
 333                 m->errorOut(e, "PipelineCommand", "readUsersPipeline");
 334                 exit(1);
 335         }
 336 }
 337 //**********************************************************************************************************************
 338
 339 bool PipelineCommand::parseCommand(string nextCommand, string& name, string& options){
 340         try {
 341                 CommandOptionParser parser(nextCommand);
 342                 name = parser.getCommandString();
 343                 options = parser.getOptionString();
 344
 345                 if (name == "") { return true; } //name == "" if () are not right
 346
 347                 return false;
 348         }
 349         catch(exception& e) {
 350                 m->errorOut(e, "PipelineCommand", "parseCommand");
 351                 exit(1);
 352         }
 353 }
 354 //**********************************************************************************************************************
 355
 356 bool PipelineCommand::checkForValidAndRequiredParameters(string name, string options, map<string, vector<string> >& mothurMadeFiles){
 357         try {
 358
 359                 //get shell of the command so we can check to make sure its valid without running it
 360                 Command* command = cFactory->getCommand(name);
 361
 362                 //check to make sure all parameters are valid for command
 363                 vector<string> validParameters = command->getValidParameters();
 364
 365                 OptionParser parser(options);
 366                 map<string, string> parameters = parser.getParameters();
 367
 368                 ValidParameters validParameter;
 369                 map<string, string>::iterator it;
 370                 map<string, vector<string> >::iterator itMade;
 371
 372                 for (it = parameters.begin(); it != parameters.end(); it++) {
 373                         if (validParameter.isValidParameter(it->first, validParameters, it->second) != true) {  return true;  } // not valid
 374                         if (it->second == "mothurmade") {
 375                                 itMade = mothurMadeFiles.find(it->first);
 376
 377                                 if (itMade == mothurMadeFiles.end()) {
 378                                         if ((name == "align.seqs") && (it->first == "candidate")) {} //do nothing about candidate
 379                                         else {
 380                                                 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();
 381                                                 return true;
 382                                         }
 383                                 }
 384                         }
 385                 }
 386
 387                 //is the command missing any required
 388                 vector<string> requiredParameters = command->getRequiredParameters();
 389
 390                 //check for or
 391                 bool hasOr = false;
 392                 int numFound = 0;
 393                 if (requiredParameters.size() > 2) {
 394                         if (requiredParameters[(requiredParameters.size()-1)] == "or") { hasOr = true; }
 395                 }
 396
 397                 for (int i = 0; i < requiredParameters.size(); i++) {
 398                         it = parameters.find(requiredParameters[i]);
 399
 400                         if (it != parameters.end()) { numFound++; }
 401                         else {
 402                                 if (!hasOr) { m->mothurOut(name + " requires the " + requiredParameters[i] + " parameter, please correct."); m->mothurOutEndLine(); }
 403                         }
 404                 }
 405
 406                 // if all are needed and not all are found
 407                 if ((!hasOr) && (numFound != requiredParameters.size())) { return true; }
 408                 //if one is needed and none are found
 409                 else if ((hasOr) && (numFound == 0)) { return true; }
 410
 411                 //update MothurMade
 412                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 413                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 414                         mothurMadeFiles[itMade->first] = itMade->second; //adds any new types
 415                 }
 416
 417                 return false;
 418         }
 419         catch(exception& e) {
 420                 m->errorOut(e, "PipelineCommand", "checkForValidAndRequiredParameters");
 421                 exit(1);
 422         }
 423 }
 424 //**********************************************************************************************************************
 425 int PipelineCommand::runUsersPipeline(){
 426         try {
 427                 ifstream in;
 428                 m->openInputFile(pipeFilename, in);
 429
 430                 string nextCommand = "";
 431
 432                 map<string, vector<string> > mothurMadeFiles;
 433
 434                 while(!in.eof()) {
 435                         nextCommand = m->getline(in);  m->gobble(in);
 436
 437                         if (nextCommand[0] != '#') {
 438                                 CommandOptionParser parser(nextCommand);
 439                                 string commandName = parser.getCommandString();
 440                                 string options = parser.getOptionString();
 441
 442                                 if (options != "") {
 443                                         bool error = fillInMothurMade(options, mothurMadeFiles);
 444                                         if (error) { in.close(); return 0; }
 445                                 }
 446
 447                                 m->mothurOutEndLine(); m->mothurOut("mothur > " + commandName + "(" + options + ")"); m->mothurOutEndLine();
 448
 449                                 if (m->control_pressed) { return 0; }
 450
 451                                 #ifdef USE_MPI
 452                                         int pid;
 453                                         MPI_Comm_rank(MPI_COMM_WORLD, &pid);
 454
 455                                         if ((cFactory->MPIEnabled(commandName)) || (pid == 0)) {
 456                                 #endif
 457
 458                                 //executes valid command
 459                                 Command* command = cFactory->getCommand(commandName, options, "pipe");
 460                                 command->execute();
 461
 462                                 //add output files to list
 463                                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 464                                 map<string, vector<string> >::iterator itMade;
 465                                 map<string, vector<string> >::iterator it;
 466                                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 467
 468                                         vector<string> temp = itMade->second;
 469                                         for (int k = 0; k < temp.size(); k++) { outputNames.push_back(temp[k]); }  //
 470
 471                                         //update Mothur Made for each file
 472                                         it = mothurMadeFiles.find(itMade->first);
 473
 474                                         if (it == mothurMadeFiles.end()) { //new type
 475
 476                                                 mothurMadeFiles[itMade->first] = temp;
 477
 478                                         }else{ //update existing type
 479                                                 vector<string> oldFileNames = it->second;
 480                                                 //look at new files, see if an old version of the file exists, if so update, else just add.
 481                                                 //for example you may have abrecovery.fasta and amazon.fasta as old files and you created a new amazon.trim.fasta.
 482
 483                                                 for (int k = 0; k < temp.size(); k++) {
 484
 485                                                         //get base name
 486                                                         string root = m->getSimpleName(temp[k]);
 487                                                         string individual = "";
 488                                                         for(int i=0;i<root.length();i++){
 489                                                                 if(root[i] == '.'){
 490                                                                         root = individual;
 491                                                                         break;
 492                                                                 }else{
 493                                                                         individual += root[i];
 494                                                                 }
 495                                                         }
 496
 497                                                         //look for that base name in oldfiles
 498                                                         int spot = -1;
 499                                                         for (int l = 0; l < oldFileNames.size(); l++) {
 500                                                                 int pos = oldFileNames[l].find(root);
 501                                                                 if (pos != string::npos) {
 502                                                                         spot = l;
 503                                                                         break;
 504                                                                 }
 505                                                         }
 506
 507                                                         //if you found it update it, else add it
 508                                                         if (spot != -1) {
 509                                                                 mothurMadeFiles[it->first][spot] = temp[k];
 510                                                         }else{
 511                                                                 mothurMadeFiles[it->first].push_back(temp[k]);
 512                                                         }
 513                                                 }
 514                                         }
 515                                 }
 516
 517                                 #ifdef USE_MPI
 518                                         }
 519                                 #endif
 520                         }
 521                 }
 522
 523                 in.close();
 524         }
 525         catch(exception& e) {
 526                 m->errorOut(e, "PipelineCommand", "runUsersPipeline");
 527                 exit(1);
 528         }
 529 }
 530 //**********************************************************************************************************************
 531 bool PipelineCommand::fillInMothurMade(string& options, map<string, vector<string> > mothurMadeFiles){
 532         try {
 533                 OptionParser parser(options);
 534                 map<string, string> parameters = parser.getParameters();
 535                 map<string, string>::iterator it;
 536                 map<string, vector<string> >::iterator itMade;
 537
 538                 options = "";
 539
 540                 //fill in mothurmade filenames
 541                 for (it = parameters.begin(); it != parameters.end(); it++) {
 542                         string paraType = it->first;
 543
 544                         if (it->second == "mothurmade") {
 545
 546                                 if (it->first == "candidate") { paraType = "fasta"; }
 547
 548                                 itMade = mothurMadeFiles.find(paraType);
 549
 550                                 if (itMade == mothurMadeFiles.end()) {
 551                                         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();
 552                                         return true;
 553                                 }else{
 554                                         vector<string> temp = itMade->second;
 555
 556                                         if (temp.size() > 1) {
 557                                                 //ask user which file to use
 558                                                 m->mothurOut("More than one file has been created for the " + paraType + " parameter. "); m->mothurOutEndLine();
 559                                                 for (int i = 0; i < temp.size(); i++) {
 560                                                         m->mothurOut(toString(i) + " - " + temp[i]); m->mothurOutEndLine();
 561                                                 }
 562                                                 m->mothurOut("Please select the number of the file you would like to use: ");
 563                                                 int num = 0;
 564                                                 cin >> num;
 565                                                 m->mothurOutJustToLog(toString(num)); m->mothurOutEndLine();
 566
 567                                                 if ((num < 0) || (num > (temp.size()-1))) { m->mothurOut("Not a valid response, quitting."); m->mothurOutEndLine(); return true; }
 568                                                 else {
 569                                                         parameters[paraType] = temp[num];
 570                                                 }
 571
 572                                                 //clears buffer so next command doesn't have error
 573                                                 string s;
 574                                                 getline(cin, s);
 575
 576                                         }else if (temp.size() == 0){
 577                                                 m->mothurOut("Sorry, we seem to think you created a " + paraType + " file, but it seems mothur doesn't have a filename."); m->mothurOutEndLine();
 578                                                 return true;
 579                                         }else{
 580                                                 parameters[paraType] = temp[0];
 581                                         }
 582                                 }
 583                         }
 584
 585                         options += it->first + "=" + parameters[paraType] + ", ";
 586                 }
 587
 588                 //rip off extra comma
 589                 options = options.substr(0, (options.length()-2));
 590
 591                 return false;
 592         }
 593         catch(exception& e) {
 594                 m->errorOut(e, "PipelineCommand", "fillInMothurMade");
 595                 exit(1);
 596         }
 597 }
 598
 599 //**********************************************************************************************************************
 600 void PipelineCommand::createPatsPipeline(){
 601         try {
 602
 603                 //sff.info command
 604                 string thisCommand = "sffinfo(sff=" + sffFile + ")";
 605                 commands.push_back(thisCommand);
 606
 607                 //trim.seqs command
 608                 string fastaFile = m->getRootName(m->getSimpleName(sffFile)) + "fasta";
 609                 string qualFile = m->getRootName(m->getSimpleName(sffFile)) + "qual";
 610                 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 + ")";
 611                 commands.push_back(thisCommand);
 612
 613                 //unique.seqs
 614                 string groupFile = m->getRootName(m->getSimpleName(fastaFile)) + "groups";
 615                 qualFile =  m->getRootName(m->getSimpleName(fastaFile)) + "trim.qual";
 616                 fastaFile =  m->getRootName(m->getSimpleName(fastaFile)) + "trim.fasta";
 617                 thisCommand = "unique.seqs(fasta=" + fastaFile + ")";
 618                 commands.push_back(thisCommand);
 619
 620                 //align.seqs
 621                 string nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "names";
 622                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "unique" + m->getExtension(fastaFile);
 623                 thisCommand = "align.seqs(processors=" + toString(processors) + ", candidate=" + fastaFile + ", template=" + alignFile + ")";
 624                 commands.push_back(thisCommand);
 625
 626                 //screen.seqs
 627                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "align";
 628                 thisCommand = "screen.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", name=" + nameFile + ", group=" + groupFile + ", optimize=end-minlength)";
 629                 commands.push_back(thisCommand);
 630
 631                 //chimera.slayer
 632                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "good" + m->getExtension(fastaFile);
 633                 nameFile = m->getRootName(m->getSimpleName(nameFile)) + "good" + m->getExtension(nameFile);
 634                 groupFile = m->getRootName(m->getSimpleName(groupFile)) + "good" + m->getExtension(groupFile);
 635                 thisCommand = "chimera.slayer(processors=" + toString(processors) + ", fasta=" + fastaFile + ", template=" + chimeraFile + ")";
 636                 commands.push_back(thisCommand);
 637
 638                 //remove.seqs
 639                 string accnosFile = m->getRootName(m->getSimpleName(fastaFile))  + "slayer.accnos";
 640                 thisCommand = "remove.seqs(fasta=" + fastaFile + ", name=" + nameFile + ", group=" + groupFile + ", accnos=" + accnosFile + ", dups=T)";
 641                 commands.push_back(thisCommand);
 642
 643                 //filter.seqs
 644                 nameFile = m->getRootName(m->getSimpleName(nameFile)) + "pick" + m->getExtension(nameFile);
 645                 groupFile = m->getRootName(m->getSimpleName(groupFile)) + "pick" + m->getExtension(groupFile);
 646                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "pick" + m->getExtension(fastaFile);
 647                 thisCommand = "filter.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", vertical=T, trump=.)";
 648                 commands.push_back(thisCommand);
 649
 650                 //unique.seqs
 651                 fastaFile =  m->getRootName(m->getSimpleName(fastaFile)) + "filter.fasta";
 652                 thisCommand = "unique.seqs(fasta=" + fastaFile + ", name=" + nameFile + ")";
 653                 commands.push_back(thisCommand);
 654
 655                 //pre.cluster
 656                 nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "names";
 657                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "unique" + m->getExtension(fastaFile);
 658                 thisCommand = "pre.cluster(fasta=" + fastaFile + ", name=" + nameFile + ", diffs=2)";
 659                 commands.push_back(thisCommand);
 660
 661                 //dist.seqs
 662                 nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "precluster.names";
 663                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "precluster" + m->getExtension(fastaFile);
 664                 thisCommand = "dist.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", cutoff=0.20)";
 665                 commands.push_back(thisCommand);
 666
 667                 //dist.seqs
 668                 string columnFile = m->getRootName(m->getSimpleName(fastaFile)) + "dist";
 669                 thisCommand = "dist.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", output=lt)";
 670                 commands.push_back(thisCommand);
 671
 672                 //read.dist
 673                 string phylipFile = m->getRootName(m->getSimpleName(fastaFile)) + "phylip.dist";
 674                 thisCommand = "read.dist(column=" + columnFile + ", name=" + nameFile + ")";
 675                 commands.push_back(thisCommand);
 676
 677                 //cluster
 678                 thisCommand = "cluster(method=average, hard=T)";
 679                 commands.push_back(thisCommand);
 680
 681                 string listFile = m->getRootName(m->getSimpleName(columnFile)) + "an.list";
 682                 string rabundFile = m->getRootName(m->getSimpleName(columnFile)) + "an.rabund";
 683
 684                 //degap.seqs
 685                 thisCommand = "degap.seqs(fasta=" + fastaFile + ")";
 686                 commands.push_back(thisCommand);
 687
 688                 //classify.seqs
 689                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "ng.fasta";
 690                 thisCommand = "classify.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", name=" + nameFile + ", template=" + classifyFile + ", taxonomy=" + taxonomyFile + ", cutoff=80)";
 691                 commands.push_back(thisCommand);
 692
 693                 string RippedTaxName = m->getRootName(m->getSimpleName(taxonomyFile));
 694                 RippedTaxName = m->getExtension(RippedTaxName.substr(0, RippedTaxName.length()-1));
 695                 if (RippedTaxName[0] == '.') { RippedTaxName = RippedTaxName.substr(1, RippedTaxName.length()); }
 696                 RippedTaxName +=  ".";
 697
 698                 string fastaTaxFile = m->getRootName(m->getSimpleName(fastaFile)) + RippedTaxName + "taxonomy";
 699                 string taxSummaryFile = m->getRootName(m->getSimpleName(fastaFile)) + RippedTaxName + "tax.summary";
 700
 701                 //phylotype
 702                 thisCommand = "phylotype(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ")";
 703                 commands.push_back(thisCommand);
 704
 705                 string phyloListFile = m->getRootName(m->getSimpleName(fastaTaxFile)) + "tx.list";
 706                 string phyloRabundFile = m->getRootName(m->getSimpleName(fastaTaxFile)) + "tx.rabund";
 707
 708                 //clearcut
 709                 thisCommand = "clearcut(phylip=" + phylipFile + ", neighbor=T)";
 710                 commands.push_back(thisCommand);
 711
 712                 string treeFile = m->getRootName(m->getSimpleName(phylipFile)) + "tre";
 713
 714                 //read.otu
 715                 thisCommand = "read.otu(list=" + listFile + ", group=" + groupFile + ", label=0.03)";
 716                 commands.push_back(thisCommand);
 717
 718                 string sharedFile = m->getRootName(m->getSimpleName(listFile)) + "shared";
 719
 720                 //read.otu
 721                 thisCommand = "read.otu(list=" + phyloListFile + ", group=" + groupFile + ", label=1)";
 722                 commands.push_back(thisCommand);
 723
 724                 string phyloSharedFile = m->getRootName(m->getSimpleName(phyloListFile)) + "shared";
 725
 726                 //read.otu
 727                 thisCommand = "read.otu(shared=" + sharedFile + ")";
 728                 commands.push_back(thisCommand);
 729
 730                 //summary.single
 731                 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)";
 732                 commands.push_back(thisCommand);
 733
 734                 //summary.shared
 735                 thisCommand = "summary.shared(calc=sharednseqs-sharedsobs-sharedchao-sharedace-anderberg-jclass-jest-kulczynski-kulczynskicody-lennon-ochiai-sorclass-sorest-whittaker-braycurtis-jabund-morisitahorn-sorabund-thetan-thetayc)";
 736                 commands.push_back(thisCommand);
 737
 738                 //read.otu
 739                 thisCommand = "read.otu(rabund=" + rabundFile + ", label=0.03)";
 740                 commands.push_back(thisCommand);
 741
 742                 //summary.single
 743                 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)";
 744                 commands.push_back(thisCommand);
 745
 746                 //read.otu
 747                 thisCommand = "read.otu(shared=" + phyloSharedFile + ")";
 748                 commands.push_back(thisCommand);
 749
 750                 //summary.single
 751                 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)";
 752                 commands.push_back(thisCommand);
 753
 754                 //summary.shared
 755                 thisCommand = "summary.shared(calc=sharednseqs-sharedsobs-sharedchao-sharedace-anderberg-jclass-jest-kulczynski-kulczynskicody-lennon-ochiai-sorclass-sorest-whittaker-braycurtis-jabund-morisitahorn-sorabund-thetan-thetayc)";
 756                 commands.push_back(thisCommand);
 757
 758                 //read.otu
 759                 thisCommand = "read.otu(rabund=" + phyloRabundFile + ", label=1)";
 760                 commands.push_back(thisCommand);
 761
 762                 //summary.single
 763                 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)";
 764                 commands.push_back(thisCommand);
 765
 766                 //classify.otu
 767                 thisCommand = "classify.otu(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ", list=" + listFile + ", cutoff=51, label=0.03)";
 768                 commands.push_back(thisCommand);
 769
 770                 //classify.otu
 771                 thisCommand = "classify.otu(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ", list=" + phyloListFile + ", cutoff=51, label=1)";
 772                 commands.push_back(thisCommand);
 773
 774                 //read.tree
 775                 thisCommand = "read.tree(tree=" + treeFile + ", name=" + nameFile + ", group=" + groupFile + ")";
 776                 commands.push_back(thisCommand);
 777
 778                 //phylo.diversity
 779                 thisCommand = "phylo.diversity(iters=100,rarefy=T)";
 780                 commands.push_back(thisCommand);
 781
 782                 //unifrac.weighted
 783                 thisCommand = "unifrac.weighted(random=false, distance=true, groups=all, processors=" + toString(processors) + ")";
 784                 commands.push_back(thisCommand);
 785
 786                 //unifrac.unweighted
 787                 thisCommand = "unifrac.unweighted(random=false, distance=true, processors=" + toString(processors) + ")";
 788                 commands.push_back(thisCommand);
 789
 790
 791         }
 792         catch(exception& e) {
 793                 m->errorOut(e, "PipelineCommand", "createPatsPipeline");
 794                 exit(1);
 795         }
 796 }
 797
 798 //**********************************************************************************************************************