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                 if (name == "system") { return false; }
 360
 361                 //get shell of the command so we can check to make sure its valid without running it
 362                 Command* command = cFactory->getCommand(name);
 363
 364                 //check to make sure all parameters are valid for command
 365                 vector<string> validParameters = command->getValidParameters();
 366
 367                 OptionParser parser(options);
 368                 map<string, string> parameters = parser.getParameters();
 369
 370                 ValidParameters validParameter;
 371                 map<string, string>::iterator it;
 372                 map<string, vector<string> >::iterator itMade;
 373
 374                 for (it = parameters.begin(); it != parameters.end(); it++) {
 375
 376                         if (validParameter.isValidParameter(it->first, validParameters, it->second) != true) {  return true;  } // not valid
 377                         if (it->second == "mothurmade") {
 378                                 itMade = mothurMadeFiles.find(it->first);
 379
 380                                 if (itMade == mothurMadeFiles.end()) {
 381                                         if ((name == "align.seqs") && (it->first == "candidate")) {} //do nothing about candidate
 382                                         else {
 383                                                 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();
 384                                                 return true;
 385                                         }
 386                                 }
 387                         }
 388                 }
 389
 390                 //is the command missing any required
 391                 vector<string> requiredParameters = command->getRequiredParameters();
 392
 393                 //check for or
 394                 bool hasOr = false;
 395                 int numFound = 0;
 396                 if (requiredParameters.size() > 2) {
 397                         if (requiredParameters[(requiredParameters.size()-1)] == "or") { hasOr = true; }
 398                 }
 399
 400                 for (int i = 0; i < requiredParameters.size(); i++) {
 401                         it = parameters.find(requiredParameters[i]);
 402
 403                         if (it != parameters.end()) { numFound++; }
 404                         else {
 405                                 if (!hasOr) { m->mothurOut(name + " requires the " + requiredParameters[i] + " parameter, please correct."); m->mothurOutEndLine(); }
 406                         }
 407                 }
 408
 409                 // if all are needed and not all are found
 410                 if ((!hasOr) && (numFound != requiredParameters.size())) { return true; }
 411                 //if one is needed and none are found
 412                 else if ((hasOr) && (numFound == 0)) { return true; }
 413
 414                 //update MothurMade
 415                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 416                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 417                         mothurMadeFiles[itMade->first] = itMade->second; //adds any new types
 418                 }
 419
 420                 return false;
 421         }
 422         catch(exception& e) {
 423                 m->errorOut(e, "PipelineCommand", "checkForValidAndRequiredParameters");
 424                 exit(1);
 425         }
 426 }
 427 //**********************************************************************************************************************
 428 int PipelineCommand::runUsersPipeline(){
 429         try {
 430                 ifstream in;
 431                 m->openInputFile(pipeFilename, in);
 432
 433                 string nextCommand = "";
 434
 435                 map<string, vector<string> > mothurMadeFiles;
 436
 437                 while(!in.eof()) {
 438                         nextCommand = m->getline(in);  m->gobble(in);
 439
 440                         if (nextCommand[0] != '#') {
 441                                 CommandOptionParser parser(nextCommand);
 442                                 string commandName = parser.getCommandString();
 443                                 string options = parser.getOptionString();
 444
 445                                 if ((options != "") && (commandName != "system")) {
 446                                         bool error = fillInMothurMade(options, mothurMadeFiles);
 447                                         if (error) { in.close(); return 0; }
 448                                 }
 449
 450                                 m->mothurOutEndLine(); m->mothurOut("mothur > " + commandName + "(" + options + ")"); m->mothurOutEndLine();
 451
 452                                 if (m->control_pressed) { return 0; }
 453
 454                                 #ifdef USE_MPI
 455                                         int pid;
 456                                         MPI_Comm_rank(MPI_COMM_WORLD, &pid);
 457
 458                                         if ((cFactory->MPIEnabled(commandName)) || (pid == 0)) {
 459                                 #endif
 460
 461                                 //executes valid command
 462                                 Command* command = cFactory->getCommand(commandName, options, "pipe");
 463                                 command->execute();
 464
 465                                 //add output files to list
 466                                 map<string, vector<string> > thisCommandsFile = command->getOutputFiles();
 467                                 map<string, vector<string> >::iterator itMade;
 468                                 map<string, vector<string> >::iterator it;
 469                                 for (itMade = thisCommandsFile.begin(); itMade != thisCommandsFile.end(); itMade++) {
 470
 471                                         vector<string> temp = itMade->second;
 472                                         for (int k = 0; k < temp.size(); k++) { outputNames.push_back(temp[k]); }  //
 473
 474                                         //update Mothur Made for each file
 475                                         it = mothurMadeFiles.find(itMade->first);
 476
 477                                         if (it == mothurMadeFiles.end()) { //new type
 478
 479                                                 mothurMadeFiles[itMade->first] = temp;
 480
 481                                         }else{ //update existing type
 482                                                 vector<string> oldFileNames = it->second;
 483                                                 //look at new files, see if an old version of the file exists, if so update, else just add.
 484                                                 //for example you may have abrecovery.fasta and amazon.fasta as old files and you created a new amazon.trim.fasta.
 485
 486                                                 for (int k = 0; k < temp.size(); k++) {
 487
 488                                                         //get base name
 489                                                         string root = m->getSimpleName(temp[k]);
 490                                                         string individual = "";
 491                                                         for(int i=0;i<root.length();i++){
 492                                                                 if(root[i] == '.'){
 493                                                                         root = individual;
 494                                                                         break;
 495                                                                 }else{
 496                                                                         individual += root[i];
 497                                                                 }
 498                                                         }
 499
 500                                                         //look for that base name in oldfiles
 501                                                         int spot = -1;
 502                                                         for (int l = 0; l < oldFileNames.size(); l++) {
 503                                                                 int pos = oldFileNames[l].find(root);
 504                                                                 if (pos != string::npos) {
 505                                                                         spot = l;
 506                                                                         break;
 507                                                                 }
 508                                                         }
 509
 510                                                         //if you found it update it, else add it
 511                                                         if (spot != -1) {
 512                                                                 mothurMadeFiles[it->first][spot] = temp[k];
 513                                                         }else{
 514                                                                 mothurMadeFiles[it->first].push_back(temp[k]);
 515                                                         }
 516                                                 }
 517                                         }
 518                                 }
 519
 520                                 #ifdef USE_MPI
 521                                         }
 522                                 #endif
 523                         }
 524                 }
 525
 526                 in.close();
 527         }
 528         catch(exception& e) {
 529                 m->errorOut(e, "PipelineCommand", "runUsersPipeline");
 530                 exit(1);
 531         }
 532 }
 533 //**********************************************************************************************************************
 534 bool PipelineCommand::fillInMothurMade(string& options, map<string, vector<string> >& mothurMadeFiles){
 535         try {
 536
 537                 OptionParser parser(options);
 538                 map<string, string> parameters = parser.getParameters();
 539                 map<string, string>::iterator it;
 540                 map<string, vector<string> >::iterator itMade;
 541
 542                 options = "";
 543
 544                 //fill in mothurmade filenames
 545                 for (it = parameters.begin(); it != parameters.end(); it++) {
 546                         string paraType = it->first;
 547                         string tempOption = it->second;
 548
 549                         if (tempOption == "mothurmade") {
 550
 551                                 if (it->first == "candidate") { paraType = "fasta"; }
 552
 553                                 itMade = mothurMadeFiles.find(paraType);
 554
 555                                 if (itMade == mothurMadeFiles.end()) {
 556                                         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();
 557                                         return true;
 558                                 }else{
 559                                         vector<string> temp = itMade->second;
 560
 561                                         if (temp.size() > 1) {
 562                                                 //ask user which file to use
 563                                                 m->mothurOut("More than one file has been created for the " + paraType + " parameter. "); m->mothurOutEndLine();
 564                                                 for (int i = 0; i < temp.size(); i++) {
 565                                                         m->mothurOut(toString(i) + " - " + temp[i]); m->mothurOutEndLine();
 566                                                 }
 567
 568                                                 m->mothurOut("Please select the number of the file you would like to use: ");
 569                                                 int num = 0;
 570                                                 cin >> num;
 571                                                 m->mothurOutJustToLog(toString(num)); m->mothurOutEndLine();
 572
 573                                                 if ((num < 0) || (num > (temp.size()-1))) { m->mothurOut("Not a valid response, quitting."); m->mothurOutEndLine(); return true; }
 574                                                 else {
 575                                                         tempOption = temp[num];
 576                                                 }
 577
 578                                                 //clears buffer so next command doesn't have error
 579                                                 string s;
 580                                                 getline(cin, s);
 581
 582                                                 vector<string> newTemp;
 583                                                 for (int i = 0; i < temp.size(); i++) {
 584                                                         if (i == num) { newTemp.push_back(temp[i]); }
 585                                                         else {
 586                                                                 m->mothurOut("Would you like to remove " + temp[i] + " as an option for " + paraType + ", (y/n): "); m->mothurOutEndLine();
 587                                                                 string response;
 588                                                                 cin >> response;
 589                                                                 m->mothurOutJustToLog(response); m->mothurOutEndLine();
 590
 591                                                                 if (response == "n") {  newTemp.push_back(temp[i]); }
 592
 593                                                                 //clears buffer so next command doesn't have error
 594                                                                 string s;
 595                                                                 getline(cin, s);
 596                                                         }
 597                                                 }
 598
 599                                                 mothurMadeFiles[paraType] = newTemp;
 600
 601
 602                                         }else if (temp.size() == 0){
 603                                                 m->mothurOut("Sorry, we seem to think you created a " + paraType + " file, but it seems mothur doesn't have a filename."); m->mothurOutEndLine();
 604                                                 return true;
 605                                         }else{
 606                                                 tempOption = temp[0];
 607                                         }
 608                                 }
 609                         }
 610
 611                         options += it->first + "=" + tempOption + ", ";
 612                 }
 613
 614                 //rip off extra comma
 615                 options = options.substr(0, (options.length()-2));
 616
 617                 return false;
 618         }
 619         catch(exception& e) {
 620                 m->errorOut(e, "PipelineCommand", "fillInMothurMade");
 621                 exit(1);
 622         }
 623 }
 624
 625 //**********************************************************************************************************************
 626 void PipelineCommand::createPatsPipeline(){
 627         try {
 628
 629                 //sff.info command
 630                 string thisCommand = "sffinfo(sff=" + sffFile + ")";
 631                 commands.push_back(thisCommand);
 632
 633                 //trim.seqs command
 634                 string fastaFile = m->getRootName(m->getSimpleName(sffFile)) + "fasta";
 635                 string qualFile = m->getRootName(m->getSimpleName(sffFile)) + "qual";
 636                 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 + ")";
 637                 commands.push_back(thisCommand);
 638
 639                 //unique.seqs
 640                 string groupFile = m->getRootName(m->getSimpleName(fastaFile)) + "groups";
 641                 qualFile =  m->getRootName(m->getSimpleName(fastaFile)) + "trim.qual";
 642                 fastaFile =  m->getRootName(m->getSimpleName(fastaFile)) + "trim.fasta";
 643                 thisCommand = "unique.seqs(fasta=" + fastaFile + ")";
 644                 commands.push_back(thisCommand);
 645
 646                 //align.seqs
 647                 string nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "names";
 648                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "unique" + m->getExtension(fastaFile);
 649                 thisCommand = "align.seqs(processors=" + toString(processors) + ", candidate=" + fastaFile + ", template=" + alignFile + ")";
 650                 commands.push_back(thisCommand);
 651
 652                 //screen.seqs
 653                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "align";
 654                 thisCommand = "screen.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", name=" + nameFile + ", group=" + groupFile + ", optimize=end-minlength)";
 655                 commands.push_back(thisCommand);
 656
 657                 //chimera.slayer
 658                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "good" + m->getExtension(fastaFile);
 659                 nameFile = m->getRootName(m->getSimpleName(nameFile)) + "good" + m->getExtension(nameFile);
 660                 groupFile = m->getRootName(m->getSimpleName(groupFile)) + "good" + m->getExtension(groupFile);
 661                 thisCommand = "chimera.slayer(processors=" + toString(processors) + ", fasta=" + fastaFile + ", template=" + chimeraFile + ")";
 662                 commands.push_back(thisCommand);
 663
 664                 //remove.seqs
 665                 string accnosFile = m->getRootName(m->getSimpleName(fastaFile))  + "slayer.accnos";
 666                 thisCommand = "remove.seqs(fasta=" + fastaFile + ", name=" + nameFile + ", group=" + groupFile + ", accnos=" + accnosFile + ", dups=T)";
 667                 commands.push_back(thisCommand);
 668
 669                 //filter.seqs
 670                 nameFile = m->getRootName(m->getSimpleName(nameFile)) + "pick" + m->getExtension(nameFile);
 671                 groupFile = m->getRootName(m->getSimpleName(groupFile)) + "pick" + m->getExtension(groupFile);
 672                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "pick" + m->getExtension(fastaFile);
 673                 thisCommand = "filter.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", vertical=T, trump=.)";
 674                 commands.push_back(thisCommand);
 675
 676                 //unique.seqs
 677                 fastaFile =  m->getRootName(m->getSimpleName(fastaFile)) + "filter.fasta";
 678                 thisCommand = "unique.seqs(fasta=" + fastaFile + ", name=" + nameFile + ")";
 679                 commands.push_back(thisCommand);
 680
 681                 //pre.cluster
 682                 nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "names";
 683                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "unique" + m->getExtension(fastaFile);
 684                 thisCommand = "pre.cluster(fasta=" + fastaFile + ", name=" + nameFile + ", diffs=2)";
 685                 commands.push_back(thisCommand);
 686
 687                 //dist.seqs
 688                 nameFile = m->getRootName(m->getSimpleName(fastaFile)) + "precluster.names";
 689                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "precluster" + m->getExtension(fastaFile);
 690                 thisCommand = "dist.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", cutoff=0.20)";
 691                 commands.push_back(thisCommand);
 692
 693                 //dist.seqs
 694                 string columnFile = m->getRootName(m->getSimpleName(fastaFile)) + "dist";
 695                 thisCommand = "dist.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", output=lt)";
 696                 commands.push_back(thisCommand);
 697
 698                 //read.dist
 699                 string phylipFile = m->getRootName(m->getSimpleName(fastaFile)) + "phylip.dist";
 700                 thisCommand = "read.dist(column=" + columnFile + ", name=" + nameFile + ")";
 701                 commands.push_back(thisCommand);
 702
 703                 //cluster
 704                 thisCommand = "cluster(method=average, hard=T)";
 705                 commands.push_back(thisCommand);
 706
 707                 string listFile = m->getRootName(m->getSimpleName(columnFile)) + "an.list";
 708                 string rabundFile = m->getRootName(m->getSimpleName(columnFile)) + "an.rabund";
 709
 710                 //degap.seqs
 711                 thisCommand = "degap.seqs(fasta=" + fastaFile + ")";
 712                 commands.push_back(thisCommand);
 713
 714                 //classify.seqs
 715                 fastaFile = m->getRootName(m->getSimpleName(fastaFile)) + "ng.fasta";
 716                 thisCommand = "classify.seqs(processors=" + toString(processors) + ", fasta=" + fastaFile + ", name=" + nameFile + ", template=" + classifyFile + ", taxonomy=" + taxonomyFile + ", cutoff=80)";
 717                 commands.push_back(thisCommand);
 718
 719                 string RippedTaxName = m->getRootName(m->getSimpleName(taxonomyFile));
 720                 RippedTaxName = m->getExtension(RippedTaxName.substr(0, RippedTaxName.length()-1));
 721                 if (RippedTaxName[0] == '.') { RippedTaxName = RippedTaxName.substr(1, RippedTaxName.length()); }
 722                 RippedTaxName +=  ".";
 723
 724                 string fastaTaxFile = m->getRootName(m->getSimpleName(fastaFile)) + RippedTaxName + "taxonomy";
 725                 string taxSummaryFile = m->getRootName(m->getSimpleName(fastaFile)) + RippedTaxName + "tax.summary";
 726
 727                 //phylotype
 728                 thisCommand = "phylotype(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ")";
 729                 commands.push_back(thisCommand);
 730
 731                 string phyloListFile = m->getRootName(m->getSimpleName(fastaTaxFile)) + "tx.list";
 732                 string phyloRabundFile = m->getRootName(m->getSimpleName(fastaTaxFile)) + "tx.rabund";
 733
 734                 //clearcut
 735                 thisCommand = "clearcut(phylip=" + phylipFile + ", neighbor=T)";
 736                 commands.push_back(thisCommand);
 737
 738                 string treeFile = m->getRootName(m->getSimpleName(phylipFile)) + "tre";
 739
 740                 //read.otu
 741                 thisCommand = "read.otu(list=" + listFile + ", group=" + groupFile + ", label=0.03)";
 742                 commands.push_back(thisCommand);
 743
 744                 string sharedFile = m->getRootName(m->getSimpleName(listFile)) + "shared";
 745
 746                 //read.otu
 747                 thisCommand = "read.otu(list=" + phyloListFile + ", group=" + groupFile + ", label=1)";
 748                 commands.push_back(thisCommand);
 749
 750                 string phyloSharedFile = m->getRootName(m->getSimpleName(phyloListFile)) + "shared";
 751
 752                 //read.otu
 753                 thisCommand = "read.otu(shared=" + sharedFile + ")";
 754                 commands.push_back(thisCommand);
 755
 756                 //summary.single
 757                 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)";
 758                 commands.push_back(thisCommand);
 759
 760                 //summary.shared
 761                 thisCommand = "summary.shared(calc=sharednseqs-sharedsobs-sharedchao-sharedace-anderberg-jclass-jest-kulczynski-kulczynskicody-lennon-ochiai-sorclass-sorest-whittaker-braycurtis-jabund-morisitahorn-sorabund-thetan-thetayc)";
 762                 commands.push_back(thisCommand);
 763
 764                 //read.otu
 765                 thisCommand = "read.otu(rabund=" + rabundFile + ", label=0.03)";
 766                 commands.push_back(thisCommand);
 767
 768                 //summary.single
 769                 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)";
 770                 commands.push_back(thisCommand);
 771
 772                 //read.otu
 773                 thisCommand = "read.otu(shared=" + phyloSharedFile + ")";
 774                 commands.push_back(thisCommand);
 775
 776                 //summary.single
 777                 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)";
 778                 commands.push_back(thisCommand);
 779
 780                 //summary.shared
 781                 thisCommand = "summary.shared(calc=sharednseqs-sharedsobs-sharedchao-sharedace-anderberg-jclass-jest-kulczynski-kulczynskicody-lennon-ochiai-sorclass-sorest-whittaker-braycurtis-jabund-morisitahorn-sorabund-thetan-thetayc)";
 782                 commands.push_back(thisCommand);
 783
 784                 //read.otu
 785                 thisCommand = "read.otu(rabund=" + phyloRabundFile + ", label=1)";
 786                 commands.push_back(thisCommand);
 787
 788                 //summary.single
 789                 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)";
 790                 commands.push_back(thisCommand);
 791
 792                 //classify.otu
 793                 thisCommand = "classify.otu(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ", list=" + listFile + ", cutoff=51, label=0.03)";
 794                 commands.push_back(thisCommand);
 795
 796                 //classify.otu
 797                 thisCommand = "classify.otu(taxonomy=" + fastaTaxFile + ", name=" + nameFile + ", list=" + phyloListFile + ", cutoff=51, label=1)";
 798                 commands.push_back(thisCommand);
 799
 800                 //read.tree
 801                 thisCommand = "read.tree(tree=" + treeFile + ", name=" + nameFile + ", group=" + groupFile + ")";
 802                 commands.push_back(thisCommand);
 803
 804                 //phylo.diversity
 805                 thisCommand = "phylo.diversity(iters=100,rarefy=T)";
 806                 commands.push_back(thisCommand);
 807
 808                 //unifrac.weighted
 809                 thisCommand = "unifrac.weighted(random=false, distance=true, groups=all, processors=" + toString(processors) + ")";
 810                 commands.push_back(thisCommand);
 811
 812                 //unifrac.unweighted
 813                 thisCommand = "unifrac.unweighted(random=false, distance=true, processors=" + toString(processors) + ")";
 814                 commands.push_back(thisCommand);
 815
 816
 817         }
 818         catch(exception& e) {
 819                 m->errorOut(e, "PipelineCommand", "createPatsPipeline");
 820                 exit(1);
 821         }
 822 }
 823
 824 //**********************************************************************************************************************