]> git.donarmstrong.com Git - mothur.git/blobdiff - classifyseqscommand.cpp
some bug fixes
[mothur.git] / classifyseqscommand.cpp
index 38378491c072c1ed33cc257ec954de6a81f3bb7d..0635fdeb4cd9a72de87d786d6a0aa30056b99a7d 100644 (file)
@@ -32,7 +32,7 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                        OptionParser parser(option);
                        map<string, string> parameters = parser.getParameters(); 
                        
-                       ValidParameters validParameter;
+                       ValidParameters validParameter("classify.seqs");
                        map<string, string>::iterator it;
                        
                        //check to make sure all parameters are valid for command
@@ -82,10 +82,7 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                        }
                        else if (templateFileName == "not open") { abort = true; }      
                        
-                       groupfile = validParameter.validFile(parameters, "group", true);
-                       if (groupfile == "not open") { abort = true; }  
-                       else if (groupfile == "not found") { groupfile = ""; }
-                       
+                                               
                        fastaFileName = validParameter.validFile(parameters, "fasta", false);
                        if (fastaFileName == "not found") { m->mothurOut("fasta is a required parameter for the classify.seqs command."); m->mothurOutEndLine(); abort = true;  }
                        else { 
@@ -110,7 +107,17 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                                        #endif
                                        
                                        ifstream in;
-                                       ableToOpen = openInputFile(fastaFileNames[i], in);
+                                       ableToOpen = 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() + getSimpleName(fastaFileNames[i]);
+                                                       m->mothurOut("Unable to open " + fastaFileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
+                                                       ableToOpen = openInputFile(tryPath, in, "noerror");
+                                                       fastaFileNames[i] = tryPath;
+                                               }
+                                       }
                                        in.close();
                                        
                                        #ifdef USE_MPI  
@@ -125,7 +132,7 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                                        #endif
                                        
                                        if (ableToOpen == 1) { 
-                                               m->mothurOut(fastaFileNames[i] + " will be disregarded."); m->mothurOutEndLine(); 
+                                               m->mothurOut("Unable to open " + fastaFileNames[i] + ". It will be disregarded."); m->mothurOutEndLine(); 
                                                //erase from file list
                                                fastaFileNames.erase(fastaFileNames.begin()+i);
                                                i--;
@@ -171,7 +178,17 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                                        #endif
 
                                        ifstream in;
-                                       ableToOpen = openInputFile(namefileNames[i], in);
+                                       ableToOpen = openInputFile(namefileNames[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() + getSimpleName(namefileNames[i]);
+                                                       m->mothurOut("Unable to open " + namefileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
+                                                       ableToOpen = openInputFile(tryPath, in, "noerror");
+                                                       namefileNames[i] = tryPath;
+                                               }
+                                       }
                                        in.close();
                                        
                                        #ifdef USE_MPI  
@@ -184,8 +201,14 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                                                }
                                                
                                        #endif
-                                       if (ableToOpen == 1) {  m->mothurOut("Unable to match name file with fasta file."); m->mothurOutEndLine(); abort = true;        }
                                        
+                                       if (ableToOpen == 1) { 
+                                               m->mothurOut("Unable to open " + namefileNames[i] + ". It will be disregarded."); m->mothurOutEndLine();  abort = true;
+                                               //erase from file list
+                                               namefileNames.erase(namefileNames.begin()+i);
+                                               i--;
+                                       }
+
                                }
                        }
 
@@ -193,6 +216,68 @@ ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
                                if (namefileNames.size() != fastaFileNames.size()) { abort = true; m->mothurOut("If you provide a name file, you must have one for each fasta file."); m->mothurOutEndLine(); }
                        }
                        
+                       groupfile = validParameter.validFile(parameters, "group", false);
+                       if (groupfile == "not found") { groupfile = "";  }
+                       else { 
+                               splitAtDash(groupfile, groupfileNames);
+                               
+                               //go through files and make sure they are good, if not, then disregard them
+                               for (int i = 0; i < groupfileNames.size(); i++) {
+                                       if (inputDir != "") {
+                                               string path = hasPath(groupfileNames[i]);
+                                               //if the user has not given a path then, add inputdir. else leave path alone.
+                                               if (path == "") {       groupfileNames[i] = inputDir + groupfileNames[i];               }
+                                       }
+                                       int ableToOpen;
+                                       
+                                       #ifdef USE_MPI  
+                                               int pid;
+                                               MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
+                                               MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
+                               
+                                               if (pid == 0) {
+                                       #endif
+
+                                       ifstream in;
+                                       ableToOpen = openInputFile(groupfileNames[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() + getSimpleName(groupfileNames[i]);
+                                                       m->mothurOut("Unable to open " + groupfileNames[i] + ". Trying default " + tryPath); m->mothurOutEndLine();
+                                                       ableToOpen = openInputFile(tryPath, in, "noerror");
+                                                       groupfileNames[i] = tryPath;
+                                               }
+                                       }
+                                       in.close();
+                                       
+                                       #ifdef USE_MPI  
+                                                       for (int j = 1; j < processors; j++) {
+                                                               MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD); 
+                                                       }
+                                               }else{
+                                                       MPI_Status status;
+                                                       MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+                                               }
+                                               
+                                       #endif
+                                       
+                                       if (ableToOpen == 1) { 
+                                               m->mothurOut("Unable to open " + groupfileNames[i] + ". It will be disregarded."); m->mothurOutEndLine(); groupfileNames[i] = "";
+                                               //erase from file list
+                                               groupfileNames.erase(groupfileNames.begin()+i);
+                                               i--;
+                                       }
+                               }
+                       }
+
+                       if (groupfile != "") {
+                               if (groupfileNames.size() != fastaFileNames.size()) { abort = true; m->mothurOut("If you provide a group file, you must have one for each fasta file."); m->mothurOutEndLine(); }
+                       }else {
+                               for (int i = 0; i < fastaFileNames.size(); i++) {  groupfileNames.push_back("");  }
+                       }
+                       
                        //check for optional parameter and set defaults
                        // ...at some point should added some additional type checking...
                        string temp;
@@ -276,7 +361,7 @@ void ClassifySeqsCommand::help(){
                m->mothurOut("The gapextend parameter allows you to specify the penalty for extending a gap in an alignment.  The default is -1.0.\n");
                m->mothurOut("The numwanted parameter allows you to specify the number of sequence matches you want with the knn method.  The default is 10.\n");
                m->mothurOut("The cutoff parameter allows you to specify a bootstrap confidence threshold for your taxonomy.  The default is 0.\n");
-               m->mothurOut("The probs parameter shut off the bootstrapping results for the bayesian method. The default is true, meaning you want the bootstrapping to be run.\n");
+               m->mothurOut("The probs parameter shuts off the bootstrapping results for the bayesian method. The default is true, meaning you want the bootstrapping to be shown.\n");
                m->mothurOut("The iters parameter allows you to specify how many iterations to do when calculating the bootstrap confidence score for your taxonomy with the bayesian method.  The default is 100.\n");
                m->mothurOut("The classify.seqs command should be in the following format: \n");
                m->mothurOut("classify.seqs(template=yourTemplateFile, fasta=yourFastaFile, method=yourClassificationMethod, search=yourSearchmethod, ksize=yourKmerSize, taxonomy=yourTaxonomyFile, processors=yourProcessors) \n");
@@ -329,7 +414,7 @@ int ClassifySeqsCommand::execute(){
 #ifdef USE_MPI 
                                int pid, end, numSeqsPerProcessor; 
                                int tag = 2001;
-                               vector<long> MPIPos;
+                               vector<unsigned long int> MPIPos;
                                
                                MPI_Status status; 
                                MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
@@ -369,21 +454,22 @@ int ClassifySeqsCommand::execute(){
                                //delete inFileName;
 
                                if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPINewTax);   MPI_File_close(&outMPITempTax);  delete classify; return 0;  }
-
-                               if(namefile != "") {  MPIReadNamesFile(namefileNames[s]);  }
                                
                                if (pid == 0) { //you are the root process 
                                        
                                        MPIPos = setFilePosFasta(fastaFileNames[s], numFastaSeqs); //fills MPIPos, returns numSeqs
                                        
                                        //send file positions to all processes
-                                       MPI_Bcast(&numFastaSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD);  //send numSeqs
-                                       MPI_Bcast(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos   
+                                       for(int i = 1; i < processors; i++) { 
+                                               MPI_Send(&numFastaSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
+                                               MPI_Send(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
+                                       }
                                        
                                        //figure out how many sequences you have to align
                                        numSeqsPerProcessor = numFastaSeqs / processors;
-                                       if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
                                        int startIndex =  pid * numSeqsPerProcessor;
+                                       if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
+                                       
                                
                                        //align your part
                                        driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPINewTax, outMPITempTax, MPIPos);
@@ -395,14 +481,15 @@ int ClassifySeqsCommand::execute(){
                                                MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
                                        }
                                }else{ //you are a child process
-                                       MPI_Bcast(&numFastaSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs
+                                       MPI_Recv(&numFastaSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
                                        MPIPos.resize(numFastaSeqs+1);
-                                       MPI_Bcast(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions
+                                       MPI_Recv(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
                                        
                                        //figure out how many sequences you have to align
                                        numSeqsPerProcessor = numFastaSeqs / processors;
-                                       if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
                                        int startIndex =  pid * numSeqsPerProcessor;
+                                       if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
+                                       
                                        
                                        //align your part
                                        driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPINewTax, outMPITempTax, MPIPos);
@@ -417,28 +504,14 @@ int ClassifySeqsCommand::execute(){
                                MPI_File_close(&inMPI);
                                MPI_File_close(&outMPINewTax);
                                MPI_File_close(&outMPITempTax);
+                               MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
                                
 #else
-                       //read namefile
-                       if(namefile != "") {
-                               nameMap.clear(); //remove old names
-                               
-                               ifstream inNames;
-                               openInputFile(namefileNames[s], inNames);
-                               
-                               string firstCol, secondCol;
-                               while(!inNames.eof()) {
-                                       inNames >> firstCol >> secondCol; gobble(inNames);
-                                       nameMap[firstCol] = getNumNames(secondCol);  //ex. seq1 seq1,seq3,seq5 -> seq1 = 3.
-                               }
-                               inNames.close();
-                       }
-
-       #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
                        if(processors == 1){
                                ifstream inFASTA;
                                openInputFile(fastaFileNames[s], inFASTA);
-                               numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+                               getNumSeqs(inFASTA, numFastaSeqs);
                                inFASTA.close();
                                
                                lines.push_back(new linePair(0, numFastaSeqs));
@@ -446,7 +519,7 @@ int ClassifySeqsCommand::execute(){
                                driver(lines[0], newTaxonomyFile, tempTaxonomyFile, fastaFileNames[s]);
                        }
                        else{
-                               vector<int> positions;
+                               vector<unsigned long int> positions;
                                processIDS.resize(0);
                                
                                ifstream inFASTA;
@@ -456,7 +529,7 @@ int ClassifySeqsCommand::execute(){
                                while(!inFASTA.eof()){
                                        input = getline(inFASTA);
                                        if (input.length() != 0) {
-                                               if(input[0] == '>'){    int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1);       }
+                                               if(input[0] == '>'){    unsigned long int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1); }
                                        }
                                }
                                inFASTA.close();
@@ -464,9 +537,9 @@ int ClassifySeqsCommand::execute(){
                                numFastaSeqs = positions.size();
                                
                                int numSeqsPerProcessor = numFastaSeqs / processors;
-                               
+       
                                for (int i = 0; i < processors; i++) {
-                                       int startPos = positions[ i * numSeqsPerProcessor ];
+                                       unsigned long int startPos = positions[ i * numSeqsPerProcessor ];
                                        if(i == processors - 1){
                                                numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
                                        }
@@ -488,7 +561,7 @@ int ClassifySeqsCommand::execute(){
        #else
                        ifstream inFASTA;
                        openInputFile(fastaFileNames[s], inFASTA);
-                       numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
+                       getNumSeqs(inFASTA, numFastaSeqs);
                        inFASTA.close();
                        
                        lines.push_back(new linePair(0, numFastaSeqs));
@@ -497,18 +570,48 @@ int ClassifySeqsCommand::execute(){
        #endif  
 #endif
 
+               m->mothurOutEndLine();
+               m->mothurOut("It took " + toString(time(NULL) - start) + " secs to classify " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
+               start = time(NULL);
+
+
                #ifdef USE_MPI  
                        if (pid == 0) {  //this part does not need to be paralellized
+                       
+                               if(namefile != "") { m->mothurOut("Reading " + namefileNames[s] + "..."); cout.flush();  MPIReadNamesFile(namefileNames[s]);  m->mothurOut("  Done."); m->mothurOutEndLine(); }
+               #else
+                       //read namefile
+                       if(namefile != "") {
+                       
+                           m->mothurOut("Reading " + namefileNames[s] + "..."); cout.flush();
+                               
+                               nameMap.clear(); //remove old names
+                               
+                               ifstream inNames;
+                               openInputFile(namefileNames[s], inNames);
+                               
+                               string firstCol, secondCol;
+                               while(!inNames.eof()) {
+                                       inNames >> firstCol >> secondCol; gobble(inNames);
+                                       
+                                       vector<string> temp;
+                                       splitAtComma(secondCol, temp);
+                       
+                                       nameMap[firstCol] = temp;  
+                               }
+                               inNames.close();
+                               
+                               m->mothurOut("  Done."); m->mothurOutEndLine();
+                       }
                #endif
 
-                       m->mothurOutEndLine();
-                       m->mothurOut("It took " + toString(time(NULL) - start) + " secs to classify " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
-                       start = time(NULL);
+                       string group = "";
+                       if (groupfile != "") {  group = groupfileNames[s]; }
                        
-                       PhyloSummary taxaSum(taxonomyFileName, groupfile);
+                       PhyloSummary taxaSum(taxonomyFileName, group);
                        
                        if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } delete classify; return 0; }
-                       
+               
                        if (namefile == "") {  taxaSum.summarize(tempTaxonomyFile);  }
                        else {
                                ifstream in;
@@ -516,6 +619,7 @@ int ClassifySeqsCommand::execute(){
                                
                                //read in users taxonomy file and add sequences to tree
                                string name, taxon;
+                               
                                while(!in.eof()){
                                        in >> name >> taxon; gobble(in);
                                        
@@ -524,9 +628,11 @@ int ClassifySeqsCommand::execute(){
                                        if (itNames == nameMap.end()) { 
                                                m->mothurOut(name + " is not in your name file please correct."); m->mothurOutEndLine(); exit(1);
                                        }else{
-                                               for (int i = 0; i < itNames->second; i++) { 
-                                                       taxaSum.addSeqToTree(name+toString(i), taxon);  //add it as many times as there are identical seqs
+                                               for (int i = 0; i < itNames->second.size(); i++) { 
+                                                       taxaSum.addSeqToTree(itNames->second[i], taxon);  //add it as many times as there are identical seqs
                                                }
+                                               itNames->second.clear();
+                                               nameMap.erase(itNames->first);
                                        }
                                }
                                in.close();
@@ -570,7 +676,7 @@ int ClassifySeqsCommand::execute(){
                        rename(unclass.c_str(), newTaxonomyFile.c_str());
                        
                        m->mothurOutEndLine();
-                       m->mothurOut("It took " + toString(time(NULL) - start) + " secs to create the summary file for  " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
+                       m->mothurOut("It took " + toString(time(NULL) - start) + " secs to create the summary file for " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
                        
                        #ifdef USE_MPI  
                                }
@@ -580,7 +686,6 @@ int ClassifySeqsCommand::execute(){
                        m->mothurOut("Output File Names: "); m->mothurOutEndLine();
                        for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
                        m->mothurOutEndLine();
-                       
                }
                
                delete classify;
@@ -697,9 +802,9 @@ int ClassifySeqsCommand::driver(linePair* line, string taxFName, string tempTFNa
 
                for(int i=0;i<line->numSeqs;i++){
                        if (m->control_pressed) { return 0; }
-                       
-                       Sequence* candidateSeq = new Sequence(inFASTA);
-                       
+               
+                       Sequence* candidateSeq = new Sequence(inFASTA); gobble(inFASTA);
+               
                        if (candidateSeq->getName() != "") {
                                taxonomy = classify->getTaxonomy(candidateSeq);
                                
@@ -736,7 +841,7 @@ int ClassifySeqsCommand::driver(linePair* line, string taxFName, string tempTFNa
 }
 //**********************************************************************************************************************
 #ifdef USE_MPI
-int ClassifySeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& newFile, MPI_File& tempFile, vector<long>& MPIPos){
+int ClassifySeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& newFile, MPI_File& tempFile, vector<unsigned long int>& MPIPos){
        try {
                MPI_Status statusNew; 
                MPI_Status statusTemp; 
@@ -838,7 +943,11 @@ int ClassifySeqsCommand::MPIReadNamesFile(string nameFilename){
                string firstCol, secondCol;
                while(!iss.eof()) {
                        iss >> firstCol >> secondCol; gobble(iss);
-                       nameMap[firstCol] = getNumNames(secondCol);  //ex. seq1 seq1,seq3,seq5 -> seq1 = 3.
+                       
+                       vector<string> temp;
+                       splitAtComma(secondCol, temp);
+                       
+                       nameMap[firstCol] = temp;  
                }
        
                MPI_File_close(&inMPI);