]> git.donarmstrong.com Git - mothur.git/blobdiff - filterseqscommand.cpp
added citation function to commands
[mothur.git] / filterseqscommand.cpp
index 4ee674856223d3ef6ec2e9605437ae915f6e3ac7..167dd3ce835f071e70c187ba9507fb2931a2e548 100644 (file)
 #include "filterseqscommand.h"
 #include "sequence.hpp"
 
-/**************************************************************************************/
 
+//**********************************************************************************************************************
+vector<string> FilterSeqsCommand::setParameters(){     
+       try {
+               CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pfasta);
+               CommandParameter phard("hard", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(phard);
+               CommandParameter ptrump("trump", "String", "", "*", "", "", "",false,false); parameters.push_back(ptrump);
+               CommandParameter psoft("soft", "Number", "", "0", "", "", "",false,false); parameters.push_back(psoft);
+               CommandParameter pvertical("vertical", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pvertical);
+               CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
+               CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
+               CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+               
+               vector<string> myArray;
+               for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
+               return myArray;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "FilterSeqsCommand", "setParameters");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+string FilterSeqsCommand::getHelpString(){     
+       try {
+               string helpString = "";
+               helpString += "The filter.seqs command reads a file containing sequences and creates a .filter and .filter.fasta file.\n";
+               helpString += "The filter.seqs command parameters are fasta, trump, soft, hard, processors and vertical. \n";
+               helpString += "The fasta parameter is required, unless you have a valid current fasta file. You may enter several fasta files to build the filter from and filter, by separating their names with -'s.\n";
+               helpString += "For example: fasta=abrecovery.fasta-amazon.fasta \n";
+               helpString += "The trump option will remove a column if the trump character is found at that position in any sequence of the alignment. Default=*, meaning no trump. \n";
+               helpString += "A soft mask removes any column where the dominant base (i.e. A, T, G, C, or U) does not occur in at least a designated percentage of sequences. Default=0.\n";
+               helpString += "The hard parameter allows you to enter a file containing the filter you want to use.\n";
+               helpString += "The vertical parameter removes columns where all sequences contain a gap character. The default is T.\n";
+               helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
+               helpString += "The filter.seqs command should be in the following format: \n";
+               helpString += "filter.seqs(fasta=yourFastaFile, trump=yourTrump) \n";
+               helpString += "Example filter.seqs(fasta=abrecovery.fasta, trump=.).\n";
+               helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
+               return helpString;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "FilterSeqsCommand", "getHelpString");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+FilterSeqsCommand::FilterSeqsCommand(){        
+       try {
+               abort = true; calledHelp = true; 
+               setParameters();
+               vector<string> tempOutNames;
+               outputTypes["fasta"] = tempOutNames;
+               outputTypes["filter"] = tempOutNames;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "FilterSeqsCommand", "FilterSeqsCommand");
+               exit(1);
+       }
+}
+/**************************************************************************************/
 FilterSeqsCommand::FilterSeqsCommand(string option)  {
        try {
-               abort = false;
+               abort = false; calledHelp = false;   
                filterFileName = "";
                
                //allow user to run help
-               if(option == "help") { help(); abort = true; }
+               if(option == "help") { help(); abort = true; calledHelp = true; }
+               else if(option == "citation") { citation(); abort = true; calledHelp = true;}
                
                else {
-                       //valid paramters for this command
-                       string Array[] =  {"fasta", "trump", "soft", "hard", "vertical", "outputdir","inputdir", "processors"};
-                       vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+                       vector<string> myArray = setParameters();
                        
                        OptionParser parser(option);
                        map<string,string> parameters = parser.getParameters();
@@ -36,6 +94,11 @@ FilterSeqsCommand::FilterSeqsCommand(string option)  {
                                if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
                        }
                        
+                       //initialize outputTypes
+                       vector<string> tempOutNames;
+                       outputTypes["fasta"] = tempOutNames;
+                       outputTypes["filter"] = tempOutNames;
+               
                        //if the user changes the input directory command factory will send this info to us in the output parameter 
                        string inputDir = validParameter.validFile(parameters, "inputdir", false);              
                        if (inputDir == "not found"){   inputDir = "";          }
@@ -60,42 +123,76 @@ FilterSeqsCommand::FilterSeqsCommand(string option)  {
                        
                        //check for required parameters
                        fasta = validParameter.validFile(parameters, "fasta", false);
-                       if (fasta == "not found") { m->mothurOut("fasta is a required parameter for the filter.seqs command."); m->mothurOutEndLine(); abort = true;  }
+                       if (fasta == "not found") {                             
+                               fasta = m->getFastaFile(); 
+                               if (fasta != "") { fastafileNames.push_back(fasta);  m->mothurOut("Using " + fasta + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
+                               else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
+                       }
                        else { 
                                m->splitAtDash(fasta, fastafileNames);
                                
                                //go through files and make sure they are good, if not, then disregard them
                                for (int i = 0; i < fastafileNames.size(); i++) {
-                                       if (inputDir != "") {
-                                               string path = m->hasPath(fastafileNames[i]);
-                                               //if the user has not given a path then, add inputdir. else leave path alone.
-                                               if (path == "") {       fastafileNames[i] = inputDir + fastafileNames[i];               }
-                                       }
-
-                                       ifstream in;
-                                       int ableToOpen = m->openInputFile(fastafileNames[i], in, "noerror");
-                               
-                                       //if you can't open it, try default location
-                                       if (ableToOpen == 1) {
-                                               if (m->getDefaultPath() != "") { //default path is set
-                                                       string tryPath = m->getDefaultPath() + m->getSimpleName(fastafileNames[i]);
-                                                       m->mothurOut("Unable to open " + fastafileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
-                                                       ableToOpen = m->openInputFile(tryPath, in, "noerror");
-                                                       fastafileNames[i] = tryPath;
+                                       
+                                       bool ignore = false;
+                                       if (fastafileNames[i] == "current") { 
+                                               fastafileNames[i] = m->getFastaFile(); 
+                                               if (fastafileNames[i] != "") {  m->mothurOut("Using " + fastafileNames[i] + " as input file for the fasta parameter where you had given current."); m->mothurOutEndLine(); }
+                                               else {  
+                                                       m->mothurOut("You have no current fastafile, ignoring current."); m->mothurOutEndLine(); ignore=true; 
+                                                       //erase from file list
+                                                       fastafileNames.erase(fastafileNames.begin()+i);
+                                                       i--;
                                                }
                                        }
-                                       in.close();
                                        
-                                       if (ableToOpen == 1) { 
-                                               m->mothurOut("Unable to open " + fastafileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
-                                               //erase from file list
-                                               fastafileNames.erase(fastafileNames.begin()+i);
-                                               i--;
-                                       }else{  
-                                               string simpleName = m->getSimpleName(fastafileNames[i]);
-                                               filterFileName += simpleName.substr(0, simpleName.find_first_of('.'));
+                                       if (!ignore) {
+                                               if (inputDir != "") {
+                                                       string path = m->hasPath(fastafileNames[i]);
+                                                       //if the user has not given a path then, add inputdir. else leave path alone.
+                                                       if (path == "") {       fastafileNames[i] = inputDir + fastafileNames[i];               }
+                                               }
+
+                                               ifstream in;
+                                               int ableToOpen = m->openInputFile(fastafileNames[i], in, "noerror");
+                                       
+                                               //if you can't open it, try default location
+                                               if (ableToOpen == 1) {
+                                                       if (m->getDefaultPath() != "") { //default path is set
+                                                               string tryPath = m->getDefaultPath() + m->getSimpleName(fastafileNames[i]);
+                                                               m->mothurOut("Unable to open " + fastafileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
+                                                               ifstream in2;
+                                                               ableToOpen = m->openInputFile(tryPath, in2, "noerror");
+                                                               in2.close();
+                                                               fastafileNames[i] = tryPath;
+                                                       }
+                                               }
+                                               
+                                               //if you can't open it, try default location
+                                               if (ableToOpen == 1) {
+                                                       if (m->getOutputDir() != "") { //default path is set
+                                                               string tryPath = m->getOutputDir() + m->getSimpleName(fastafileNames[i]);
+                                                               m->mothurOut("Unable to open " + fastafileNames[i] + ". Trying output directory " + tryPath); m->mothurOutEndLine();
+                                                               ifstream in2;
+                                                               ableToOpen = m->openInputFile(tryPath, in2, "noerror");
+                                                               in2.close();
+                                                               fastafileNames[i] = tryPath;
+                                                       }
+                                               }
+                                               
+                                               in.close();
+                                               
+                                               if (ableToOpen == 1) { 
+                                                       m->mothurOut("Unable to open " + fastafileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();
+                                                       //erase from file list
+                                                       fastafileNames.erase(fastafileNames.begin()+i);
+                                                       i--;
+                                               }else{  
+                                                       string simpleName = m->getSimpleName(fastafileNames[i]);
+                                                       filterFileName += simpleName.substr(0, simpleName.find_first_of('.'));
+                                               }
+                                               in.close();
                                        }
-                                       in.close();
                                }
                                
                                //make sure there is at least one valid file left
@@ -122,7 +219,8 @@ FilterSeqsCommand::FilterSeqsCommand(string option)  {
                        temp = validParameter.validFile(parameters, "soft", false);                             if (temp == "not found") { soft = 0; }
                        else {  soft = (float)atoi(temp.c_str()) / 100.0;  }
                        
-                       temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
+                       temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
+                       m->setProcessors(temp);
                        convert(temp, processors); 
                        
                        vertical = validParameter.validFile(parameters, "vertical", false);             
@@ -140,38 +238,12 @@ FilterSeqsCommand::FilterSeqsCommand(string option)  {
                exit(1);
        }
 }
-
-//**********************************************************************************************************************
-
-void FilterSeqsCommand::help(){
-       try {
-                               
-               m->mothurOut("The filter.seqs command reads a file containing sequences and creates a .filter and .filter.fasta file.\n");
-               m->mothurOut("The filter.seqs command parameters are fasta, trump, soft, hard and vertical. \n");
-               m->mothurOut("The fasta parameter is required. You may enter several fasta files to build the filter from and filter, by separating their names with -'s.\n");
-               m->mothurOut("For example: fasta=abrecovery.fasta-amazon.fasta \n");
-               m->mothurOut("The trump parameter .... The default is ...\n");
-               m->mothurOut("The soft parameter .... The default is ....\n");
-               m->mothurOut("The hard parameter allows you to enter a file containing the filter you want to use.\n");
-               m->mothurOut("The vertical parameter removes columns where all sequences contain a gap character. The default is T.\n");
-               m->mothurOut("The filter.seqs command should be in the following format: \n");
-               m->mothurOut("filter.seqs(fasta=yourFastaFile, trump=yourTrump) \n");
-               m->mothurOut("Example filter.seqs(fasta=abrecovery.fasta, trump=.).\n");
-               m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
-               
-       }
-       catch(exception& e) {
-               m->errorOut(e, "FilterSeqsCommand", "help");
-               exit(1);
-       }
-}
-
 /**************************************************************************************/
 
 int FilterSeqsCommand::execute() {     
        try {
        
-               if (abort == true) { return 0; }
+               if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
                ifstream inFASTA;
                m->openInputFile(fastafileNames[0], inFASTA);
@@ -187,7 +259,7 @@ int FilterSeqsCommand::execute() {
                
                m->mothurOutEndLine();  m->mothurOutEndLine();
                
-               if (m->control_pressed) { return 0; }
+               if (m->control_pressed) { outputTypes.clear(); return 0; }
                
                #ifdef USE_MPI
                        int pid;
@@ -198,11 +270,15 @@ int FilterSeqsCommand::execute() {
                
                ofstream outFilter;
                
-               string filterFile = outputDir + filterFileName + ".filter";
+               //prevent giantic file name
+               string filterFile;
+               if (fastafileNames.size() > 3) { filterFile = outputDir + "merge.filter"; }
+               else {  filterFile = outputDir + filterFileName + ".filter";  }
+               
                m->openOutputFile(filterFile, outFilter);
                outFilter << filter << endl;
                outFilter.close();
-               outputNames.push_back(filterFile);
+               outputNames.push_back(filterFile); outputTypes["filter"].push_back(filterFile);
                
                #ifdef USE_MPI
                        }
@@ -221,7 +297,7 @@ int FilterSeqsCommand::execute() {
                        if(filter[i] == '1'){   filteredLength++;       }
                }
                
-               if (m->control_pressed) {  for(int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }  return 0; }
+               if (m->control_pressed) {  outputTypes.clear(); for(int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }  return 0; }
 
                
                m->mothurOutEndLine();
@@ -230,6 +306,12 @@ int FilterSeqsCommand::execute() {
                m->mothurOut("Length of the original alignment: " + toString(alignmentLength)); m->mothurOutEndLine();
                m->mothurOut("Number of sequences used to construct filter: " + toString(numSeqs)); m->mothurOutEndLine();
                
+               //set fasta file as new current fastafile
+               string current = "";
+               itTypes = outputTypes.find("fasta");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
+               }
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();
@@ -256,7 +338,7 @@ int FilterSeqsCommand::filterSequences() {
                                
                                string filteredFasta = outputDir + m->getRootName(m->getSimpleName(fastafileNames[s])) + "filter.fasta";
 #ifdef USE_MPI 
-                               int pid, start, end, numSeqsPerProcessor, num; 
+                               int pid, numSeqsPerProcessor, num; 
                                int tag = 2001;
                                vector<unsigned long int>MPIPos;
                                                
@@ -265,7 +347,6 @@ int FilterSeqsCommand::filterSequences() {
                                MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
                                
                                MPI_File outMPI;
-                               MPI_File tempMPI;
                                MPI_File inMPI;
                                int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY; 
                                int inMode=MPI_MODE_RDONLY; 
@@ -305,8 +386,8 @@ int FilterSeqsCommand::filterSequences() {
                                        
                                        //wait on chidren
                                        for(int i = 1; i < processors; i++) { 
-                                               char buf[4];
-                                               MPI_Recv(buf, 4, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); 
+                                               char buf[5];
+                                               MPI_Recv(buf, 5, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status); 
                                        }
                                        
                                }else { //you are a child process
@@ -326,11 +407,11 @@ int FilterSeqsCommand::filterSequences() {
                                        
                                        if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPI);  return 0;  }
                                        
-                                       char buf[4];
+                                       char buf[5];
                                        strcpy(buf, "done"); 
                                        
                                        //tell parent you are done.
-                                       MPI_Send(buf, 4, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
+                                       MPI_Send(buf, 5, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
                                }
                                
                                MPI_File_close(&outMPI);
@@ -362,13 +443,13 @@ int FilterSeqsCommand::filterSequences() {
                                
                                if (m->control_pressed) {  return 1; }
                #else
-                               numFastaSeqs = driverRunFilter(filter, filteredFasta, fastafileNames[s], lines[0]);
+                               int numFastaSeqs = driverRunFilter(filter, filteredFasta, fastafileNames[s], lines[0]);
                                numSeqs += numFastaSeqs;
 
                                if (m->control_pressed) {  return 1; }
                #endif
 #endif
-                       outputNames.push_back(filteredFasta);
+                       outputNames.push_back(filteredFasta); outputTypes["fasta"].push_back(filteredFasta);
                }
 
                return 0;
@@ -485,8 +566,12 @@ int FilterSeqsCommand::driverRunFilter(string F, string outputFilename, string i
                                count++;
                        }
                        
-                       unsigned long int pos = in.tellg();
-                       if ((pos == -1) || (pos >= filePos->end)) { break; }
+                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+                               unsigned long int pos = in.tellg();
+                               if ((pos == -1) || (pos >= filePos->end)) { break; }
+                       #else
+                               if (in.eof()) { break; }
+                       #endif
                        
                        //report progress
                        if((count) % 100 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine();           }
@@ -533,7 +618,11 @@ int FilterSeqsCommand::createProcessesRunFilter(string F, string filename) {
                                out.close();
                                
                                exit(0);
-                       }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
+                       }else { 
+                               m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
+                               for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+                               exit(0);
+                       }
                }
                
                //force parent to wait until all the processes are done
@@ -662,7 +751,7 @@ string FilterSeqsCommand::createFilter() {
                                
                                if (m->control_pressed) {  return filterString; }
                #else
-                               numFastaSeqs = driverCreateFilter(F, fastafileNames[s], lines[0]);
+                               int numFastaSeqs = driverCreateFilter(F, fastafileNames[s], lines[0]);
                                numSeqs += numFastaSeqs;
                                if (m->control_pressed) {  return filterString; }
                #endif
@@ -786,8 +875,12 @@ int FilterSeqsCommand::driverCreateFilter(Filters& F, string filename, linePair*
                                        count++;
                        }
                        
-                       unsigned long int pos = in.tellg();
-                       if ((pos == -1) || (pos >= filePos->end)) { break; }
+                       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+                               unsigned long int pos = in.tellg();
+                               if ((pos == -1) || (pos >= filePos->end)) { break; }
+                       #else
+                               if (in.eof()) { break; }
+                       #endif
                        
                        //report progress
                        if((count) % 100 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine();           }
@@ -855,7 +948,7 @@ int FilterSeqsCommand::MPICreateFilter(int start, int num, Filters& F, MPI_File&
 int FilterSeqsCommand::createProcessesCreateFilter(Filters& F, string filename) {
        try {
 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
-               int process = 0;
+               int process = 1;
                int num = 0;
                processIDS.clear();
                
@@ -893,16 +986,23 @@ int FilterSeqsCommand::createProcessesCreateFilter(Filters& F, string filename)
                                out.close();
                                
                                exit(0);
-                       }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
+                       }else { 
+                               m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); 
+                               for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+                               exit(0);
+                       }
                }
                
+               //parent do your part
+               num = driverCreateFilter(F, filename, lines[0]);
+               
                //force parent to wait until all the processes are done
-               for (int i=0;i<processors;i++) { 
+               for (int i=0;i<(processors-1);i++) { 
                        int temp = processIDS[i];
                        wait(&temp);
                }
                
-               //parent reads in and combine Filter info
+               //parent reads in and combines Filter info
                for (int i = 0; i < processIDS.size(); i++) {
                        string tempFilename = filename + toString(processIDS[i]) + "filterValues.temp";
                        ifstream in;