]> git.donarmstrong.com Git - mothur.git/blobdiff - classify.cpp
added save parameter to align.seqs, chimera commands, classify.seqs, and seq.error...
[mothur.git] / classify.cpp
index 925dedacbf0edee93af4d89b4cf6386169b6cfc0..7770999e779917996a916e87836a03c3a2a8f26b 100644 (file)
 #include "suffixdb.hpp"
 #include "blastdb.hpp"
 #include "distancedb.hpp"
+#include "referencedb.h"
 
 /**************************************************************************************************/
 void Classify::generateDatabaseAndNames(string tfile, string tempFile, string method, int kmerSize, float gapOpen, float gapExtend, float match, float misMatch)  {            
        try {   
-               taxFile = tfile;
-               readTaxonomy(taxFile);  
+               ReferenceDB* rdb = ReferenceDB::getInstance();
                
-               templateFile = tempFile;        
+               if (tfile == "saved") { tfile = rdb->getSavedTaxonomy(); }
                
-               int start = time(NULL);
+               taxFile = tfile;
+               readTaxonomy(taxFile);  
                int numSeqs = 0;
                
-               m->mothurOut("Generating search database...    "); cout.flush();
-#ifdef USE_MPI 
-                       int pid, processors;
-                       vector<unsigned long int> positions;
-                       int tag = 2001;
-               
-                       MPI_Status status; 
-                       MPI_File inMPI;
-                       MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-                       MPI_Comm_size(MPI_COMM_WORLD, &processors);
-
-                       //char* inFileName = new char[tempFile.length()];
-                       //memcpy(inFileName, tempFile.c_str(), tempFile.length());
+               if (tempFile == "saved") {
+                       int start = time(NULL);
+                       m->mothurOutEndLine();  m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory.");        m->mothurOutEndLine();
+                       
+                       numSeqs = rdb->referenceSeqs.size();
+                       templateFile = rdb->getSavedReference();
+                       tempFile = rdb->getSavedReference();
+                       
+                       bool needToGenerate = true;
+                       string kmerDBName;
+                       if(method == "kmer")                    {       
+                               database = new KmerDB(tempFile, kmerSize);                      
+                               
+                               kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
+                               ifstream kmerFileTest(kmerDBName.c_str());
+                               if(kmerFileTest){       
+                                       bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion());
+                                       if (GoodFile) {  needToGenerate = false;        }
+                               }
+                       }
+                       else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
+                       else if(method == "blast")              {       database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch);   }
+                       else if(method == "distance")   {       database = new DistanceDB();    }
+                       else {
+                               m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
+                               m->mothurOutEndLine();
+                               database = new KmerDB(tempFile, 8);
+                       }
+                       
+                       if (needToGenerate) {
+                               for (int k = 0; k < rdb->referenceSeqs.size(); k++) {
+                                       Sequence temp(rdb->referenceSeqs[k].getName(), rdb->referenceSeqs[k].getAligned());
+                                       names.push_back(temp.getName());
+                                       database->addSequence(temp);    
+                               }
+                               database->generateDB();
+                       }else if ((method == "kmer") && (!needToGenerate)) {    
+                               ifstream kmerFileTest(kmerDBName.c_str());
+                               database->readKmerDB(kmerFileTest);     
+                               
+                               for (int k = 0; k < rdb->referenceSeqs.size(); k++) {
+                                       names.push_back(rdb->referenceSeqs[k].getName());
+                               }                       
+                       }       
+                       
+                       database->setNumSeqs(numSeqs);
+                       
+                       //sanity check
+                       bool okay = phyloTree->ErrorCheck(names);
+                       
+                       if (!okay) { m->control_pressed = true; }
                        
-                       char inFileName[1024];
-                       strcpy(inFileName, tempFile.c_str());
+                       m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences and generate the search databases.");m->mothurOutEndLine();  
+                       
+               }else {
+                       
+                       templateFile = tempFile;        
+                       
+                       int start = time(NULL);
+                       
+                       m->mothurOut("Generating search database...    "); cout.flush();
+       #ifdef USE_MPI  
+                               int pid, processors;
+                               vector<unsigned long int> positions;
+                               int tag = 2001;
+                       
+                               MPI_Status status; 
+                               MPI_File inMPI;
+                               MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
+                               MPI_Comm_size(MPI_COMM_WORLD, &processors);
+
+                               //char* inFileName = new char[tempFile.length()];
+                               //memcpy(inFileName, tempFile.c_str(), tempFile.length());
+                               
+                               char inFileName[1024];
+                               strcpy(inFileName, tempFile.c_str());
 
-                       MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
-                       //delete inFileName;
+                               MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
+                               //delete inFileName;
 
-                       if (pid == 0) { //only one process needs to scan file
-                               positions = m->setFilePosFasta(tempFile, numSeqs); //fills MPIPos, returns numSeqs
+                               if (pid == 0) { //only one process needs to scan file
+                                       positions = m->setFilePosFasta(tempFile, numSeqs); //fills MPIPos, returns numSeqs
 
-                               //send file positions to all processes
-                               for(int i = 1; i < processors; i++) { 
-                                       MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
-                                       MPI_Send(&positions[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
+                                       //send file positions to all processes
+                                       for(int i = 1; i < processors; i++) { 
+                                               MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
+                                               MPI_Send(&positions[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
+                                       }
+                               }else{
+                                       MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
+                                       positions.resize(numSeqs+1);
+                                       MPI_Recv(&positions[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
                                }
-                       }else{
-                               MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
-                               positions.resize(numSeqs+1);
-                               MPI_Recv(&positions[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
-                       }
+                               
+                               //create database
+                               if(method == "kmer")                    {       database = new KmerDB(tempFile, kmerSize);                      }
+                               else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
+                               else if(method == "blast")              {       database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch);   }
+                               else if(method == "distance")   {       database = new DistanceDB();    }
+                               else {
+                                       m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); m->mothurOutEndLine();
+                                       database = new KmerDB(tempFile, 8);
+                               }
+
+                               //read file 
+                               for(int i=0;i<numSeqs;i++){
+                                       //read next sequence
+                                       int length = positions[i+1] - positions[i];
+                                       char* buf4 = new char[length];
+                                       MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
+                                       
+                                       string tempBuf = buf4;
+                                       if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
+                                       delete buf4;
+                                       istringstream iss (tempBuf,istringstream::in);
+                                       
+                                       Sequence temp(iss);  
+                                       if (temp.getName() != "") {
+                                               if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
+                                               names.push_back(temp.getName());
+                                               database->addSequence(temp);    
+                                       }
+                               }
+                               
+                               database->generateDB();
+                               MPI_File_close(&inMPI);
+                               MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
+               #else
                        
-                       //create database
-                       if(method == "kmer")                    {       database = new KmerDB(tempFile, kmerSize);                      }
+                       //need to know number of template seqs for suffixdb
+                       if (method == "suffix") {
+                               ifstream inFASTA;
+                               m->openInputFile(tempFile, inFASTA);
+                               m->getNumSeqs(inFASTA, numSeqs);
+                               inFASTA.close();
+                       }
+
+                       bool needToGenerate = true;
+                       string kmerDBName;
+                       if(method == "kmer")                    {       
+                               database = new KmerDB(tempFile, kmerSize);                      
+                               
+                               kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
+                               ifstream kmerFileTest(kmerDBName.c_str());
+                               if(kmerFileTest){       
+                                       bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion());
+                                       if (GoodFile) {  needToGenerate = false;        }
+                               }
+                       }
                        else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
                        else if(method == "blast")              {       database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch);   }
                        else if(method == "distance")   {       database = new DistanceDB();    }
                        else {
-                               m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); m->mothurOutEndLine();
+                               m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
+                               m->mothurOutEndLine();
                                database = new KmerDB(tempFile, 8);
                        }
-
-                       //read file 
-                       for(int i=0;i<numSeqs;i++){
-                               //read next sequence
-                               int length = positions[i+1] - positions[i];
-                               char* buf4 = new char[length];
-                               MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
-                               
-                               string tempBuf = buf4;
-                               if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
-                               delete buf4;
-                               istringstream iss (tempBuf,istringstream::in);
+                       
+                       if (needToGenerate) {
+                               ifstream fastaFile;
+                               m->openInputFile(tempFile, fastaFile);
                                
-                               Sequence temp(iss);  
-                               if (temp.getName() != "") {
+                               while (!fastaFile.eof()) {
+                                       Sequence temp(fastaFile);
+                                       m->gobble(fastaFile);
+                                       
+                                       if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
+                                       
                                        names.push_back(temp.getName());
+                                                               
                                        database->addSequence(temp);    
                                }
-                       }
-                       
-                       database->generateDB();
-                       MPI_File_close(&inMPI);
-                       MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
-       #else
-               
-               //need to know number of template seqs for suffixdb
-               if (method == "suffix") {
-                       ifstream inFASTA;
-                       m->openInputFile(tempFile, inFASTA);
-                       m->getNumSeqs(inFASTA, numSeqs);
-                       inFASTA.close();
-               }
+                               fastaFile.close();
 
-               bool needToGenerate = true;
-               string kmerDBName;
-               if(method == "kmer")                    {       
-                       database = new KmerDB(tempFile, kmerSize);                      
+                               database->generateDB();
+                               
+                       }else if ((method == "kmer") && (!needToGenerate)) {    
+                               ifstream kmerFileTest(kmerDBName.c_str());
+                               database->readKmerDB(kmerFileTest);     
                        
-                       kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
-                       ifstream kmerFileTest(kmerDBName.c_str());
-                       if(kmerFileTest){       
-                               bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion());
-                               if (GoodFile) {  needToGenerate = false;        }
+                               ifstream fastaFile;
+                               m->openInputFile(tempFile, fastaFile);
+                               
+                               while (!fastaFile.eof()) {
+                                       Sequence temp(fastaFile);
+                                       m->gobble(fastaFile);
+                                       
+                                       if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
+                                       names.push_back(temp.getName());
+                               }
+                               fastaFile.close();
                        }
-               }
-               else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
-               else if(method == "blast")              {       database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch);   }
-               else if(method == "distance")   {       database = new DistanceDB();    }
-               else {
-                       m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
-                       m->mothurOutEndLine();
-                       database = new KmerDB(tempFile, 8);
-               }
+       #endif  
                
-               if (needToGenerate) {
-                       ifstream fastaFile;
-                       m->openInputFile(tempFile, fastaFile);
+                       database->setNumSeqs(names.size());
                        
-                       while (!fastaFile.eof()) {
-                               Sequence temp(fastaFile);
-                               m->gobble(fastaFile);
+                       //sanity check
+                       bool okay = phyloTree->ErrorCheck(names);
                        
-                               names.push_back(temp.getName());
-                                                       
-                               database->addSequence(temp);    
-                       }
-                       fastaFile.close();
-
-                       database->generateDB();
+                       if (!okay) { m->control_pressed = true; }
                        
-               }else if ((method == "kmer") && (!needToGenerate)) {    
-                       ifstream kmerFileTest(kmerDBName.c_str());
-                       database->readKmerDB(kmerFileTest);     
-               
-                       ifstream fastaFile;
-                       m->openInputFile(tempFile, fastaFile);
-                       
-                       while (!fastaFile.eof()) {
-                               Sequence temp(fastaFile);
-                               m->gobble(fastaFile);
-
-                               names.push_back(temp.getName());
-                       }
-                       fastaFile.close();
+                       m->mothurOut("DONE."); m->mothurOutEndLine();
+                       m->mothurOut("It took " + toString(time(NULL) - start) + " seconds generate search database. "); m->mothurOutEndLine();
                }
-#endif 
-       
-               database->setNumSeqs(names.size());
-               
-               //sanity check
-               bool okay = phyloTree->ErrorCheck(names);
-               
-               if (!okay) { m->control_pressed = true; }
-               
-               m->mothurOut("DONE."); m->mothurOutEndLine();
-               m->mothurOut("It took " + toString(time(NULL) - start) + " seconds generate search database. "); m->mothurOutEndLine();
 
        }
        catch(exception& e) {