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