]> git.donarmstrong.com Git - mothur.git/blobdiff - blastdb.cpp
sffinfo bug with flow grams right index when clipQualRight=0
[mothur.git] / blastdb.cpp
index 438d90c30961dfbabe5436483f3d33a23fc537ae..2eced7e14dc6fa42e0147747fc622bfa088250f7 100644 (file)
 
 /**************************************************************************************************/
 
-BlastDB::BlastDB(string tag, float gO, float gE, float m, float mM) : Database(), 
-gapOpen(gO), gapExtend(gE), match(m), misMatch(mM) {
-       
-       count = 0;
-
-       int randNumber = rand();
-       dbFileName = tag + toString(randNumber) + ".template.unaligned.fasta";
-       queryFileName = tag + toString(randNumber) + ".candidate.unaligned.fasta";
-       blastFileName = tag + toString(randNumber) + ".blast";
+BlastDB::BlastDB(string tag, float gO, float gE, float mm, float mM, string b, int tid) : Database(), 
+gapOpen(gO), gapExtend(gE), match(mm), misMatch(mM) {
+       try {
+               count = 0;
+               path = b;
+               threadID = tid;
 
+               int randNumber = rand();
+               //int randNumber = 12345;
+               string pid = "";
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               pid += getpid();        
+#else
+               pid += toString(threadID);      
+#endif
+               
+               
+               dbFileName = tag + pid + toString(randNumber) + ".template.unaligned.fasta";
+               queryFileName = tag + pid + toString(randNumber) + ".candidate.unaligned.fasta";
+               blastFileName = tag + pid + toString(randNumber) + ".blast";
+               
+               //make sure blast exists in the write place
+               if (path == "") {
+                       path = m->argv;
+                       string tempPath = path;
+                       for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
+                       path = path.substr(0, (tempPath.find_last_of('m')));
+                       
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+                       path += "blast/bin/";   
+#else
+                       path += "blast\\bin\\";
+#endif                 
+               }
+               
+               
+               string formatdbCommand;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               formatdbCommand = path + "formatdb";    //      format the database, -o option gives us the ability
+#else
+               formatdbCommand = path + "formatdb.exe";
+#endif
+               
+               //test to make sure formatdb exists
+               ifstream in;
+               formatdbCommand = m->getFullPathName(formatdbCommand);
+               int ableToOpen = m->openInputFile(formatdbCommand, in, "no error"); in.close();
+               if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + formatdbCommand + " file does not exist. mothur requires formatdb.exe."); m->mothurOutEndLine(); m->control_pressed = true; }
+               
+               string blastCommand;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               blastCommand = path + "blastall";       //      format the database, -o option gives us the ability
+#else
+               blastCommand = path + "blastall.exe";
+               //wrap entire string in ""
+               //blastCommand = "\"" + blastCommand + "\"";
+#endif
+               
+               //test to make sure formatdb exists
+               ifstream in2;
+               blastCommand = m->getFullPathName(blastCommand);
+               ableToOpen = m->openInputFile(blastCommand, in2, "no error"); in2.close();
+               if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + blastCommand + " file does not exist. mothur requires blastall.exe."); m->mothurOutEndLine(); m->control_pressed = true; }
+               
+               
+               string megablastCommand;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               megablastCommand = path + "megablast";  //      format the database, -o option gives us the ability
+#else
+               megablastCommand = path + "megablast.exe";
+#endif
+               
+               //test to make sure formatdb exists
+               ifstream in3;
+               megablastCommand = m->getFullPathName(megablastCommand);
+               ableToOpen = m->openInputFile(megablastCommand, in3, "no error"); in3.close();
+               if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " +  megablastCommand + " file does not exist. mothur requires megablast.exe."); m->mothurOutEndLine(); m->control_pressed = true; }
+               
+       }
+       catch(exception& e) {
+               m->errorOut(e, "BlastDB", "BlastDB");
+               exit(1);
+       }
 }
 /**************************************************************************************************/
 
-BlastDB::BlastDB() : Database() {
+BlastDB::BlastDB(string b, int tid) : Database() {
        try {
                count = 0;
-
+               
+               path = b;
+               threadID = tid;
+               
+               //make sure blast exists in the write place
+               if (path == "") {
+                       path = m->argv;
+                       string tempPath = path;
+                       for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
+                       path = path.substr(0, (tempPath.find_last_of('m')));
+                       
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+                       path += "blast/bin/";   
+#else
+                       path += "blast\\bin\\";
+#endif                 
+               }
+               
                int randNumber = rand();
-               dbFileName = toString(randNumber) + ".template.unaligned.fasta";
-               queryFileName = toString(randNumber) + ".candidate.unaligned.fasta";
-               blastFileName = toString(randNumber) + ".blast";
+               string pid = "";
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               pid += getpid();        
+#else
+               pid += toString(threadID);      
+#endif
+               
+               dbFileName = pid + toString(randNumber) + ".template.unaligned.fasta";
+               queryFileName = pid + toString(randNumber) + ".candidate.unaligned.fasta";
+               blastFileName = pid + toString(randNumber) + ".blast";
+               
+               string formatdbCommand;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               formatdbCommand = path + "formatdb";    //      format the database, -o option gives us the ability
+#else
+               formatdbCommand = path + "formatdb.exe";
+               //wrap entire string in ""
+               //formatdbCommand = "\"" + formatdbCommand + "\"";
+#endif
+               
+               //test to make sure formatdb exists
+               ifstream in;
+               formatdbCommand = m->getFullPathName(formatdbCommand);
+               int ableToOpen = m->openInputFile(formatdbCommand, in, "no error"); in.close();
+               if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " +  formatdbCommand + " file does not exist. mothur requires formatdb.exe."); m->mothurOutEndLine(); m->control_pressed = true; }
+               
+               string blastCommand;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               blastCommand = path + "blastall";       //      format the database, -o option gives us the ability
+#else
+               blastCommand = path + "blastall.exe";
+               //wrap entire string in ""
+               //blastCommand = "\"" + blastCommand + "\"";
+#endif
+               
+               //test to make sure formatdb exists
+               ifstream in2;
+               blastCommand = m->getFullPathName(blastCommand);
+               ableToOpen = m->openInputFile(blastCommand, in2, "no error"); in2.close();
+               if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + blastCommand + " file does not exist. mothur requires blastall.exe."); m->mothurOutEndLine(); m->control_pressed = true; }
+               
+               
+               string megablastCommand;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               megablastCommand = path + "megablast";  //      format the database, -o option gives us the ability
+#else
+               megablastCommand = path + "megablast.exe";
+               //wrap entire string in ""
+               //megablastCommand = "\"" + megablastCommand + "\"";
+#endif
+               
+               //test to make sure formatdb exists
+               ifstream in3;
+               megablastCommand = m->getFullPathName(megablastCommand);
+               ableToOpen = m->openInputFile(megablastCommand, in3, "no error"); in3.close();
+               if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + megablastCommand + " file does not exist. mothur requires megablast.exe."); m->mothurOutEndLine(); m->control_pressed = true; }
+               
+               
        }
        catch(exception& e) {
                m->errorOut(e, "BlastDB", "BlastDB");
@@ -46,14 +191,14 @@ BlastDB::BlastDB() : Database() {
 
 BlastDB::~BlastDB(){
        try{
-               remove(queryFileName.c_str());                          //      let's clean stuff up and remove the temp files
-               remove(dbFileName.c_str());                                     //      let's clean stuff up and remove the temp files
-               remove((dbFileName+".nsq").c_str());                                    //      let's clean stuff up and remove the temp files
-               remove((dbFileName+".nsi").c_str());                                    //      let's clean stuff up and remove the temp files
-               remove((dbFileName+".nsd").c_str());                                    //      let's clean stuff up and remove the temp files
-               remove((dbFileName+".nin").c_str());                                    //      let's clean stuff up and remove the temp files
-               remove((dbFileName+".nhr").c_str());                                    //      let's clean stuff up and remove the temp files
-               remove(blastFileName.c_str());                          //      let's clean stuff up and remove the temp files
+               m->mothurRemove(queryFileName);                         //      let's clean stuff up and remove the temp files
+               m->mothurRemove(dbFileName);                                    //      let's clean stuff up and remove the temp files
+               m->mothurRemove((dbFileName+".nsq"));                                   //      let's clean stuff up and remove the temp files
+               m->mothurRemove((dbFileName+".nsi"));                                   //      let's clean stuff up and remove the temp files
+               m->mothurRemove((dbFileName+".nsd"));                                   //      let's clean stuff up and remove the temp files
+               m->mothurRemove((dbFileName+".nin"));                                   //      let's clean stuff up and remove the temp files
+               m->mothurRemove((dbFileName+".nhr"));                                   //      let's clean stuff up and remove the temp files
+               m->mothurRemove(blastFileName.c_str());                         //      let's clean stuff up and remove the temp files
        }
        catch(exception& e) {
                m->errorOut(e, "BlastDB", "~BlastDB");
@@ -67,7 +212,10 @@ vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
                vector<int> topMatches;
                
                ofstream queryFile;
-               m->openOutputFile((queryFileName+seq->getName()), queryFile);
+               int randNumber = rand();
+               string pid = scrubName(seq->getName());
+               
+               m->openOutputFile((queryFileName+pid+toString(randNumber)), queryFile);
                queryFile << '>' << seq->getName() << endl;
                queryFile << seq->getUnaligned() << endl;
                queryFile.close();
@@ -77,12 +225,22 @@ vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
                //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
                //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
                
-               string blastCommand = path + "blast/bin/blastall -p blastn -d " + dbFileName + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);;
-               blastCommand += (" -i " + (queryFileName+seq->getName()) + " -o " + blastFileName+seq->getName());
+               string blastCommand;
+               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+               
+                       blastCommand = path + "blastall -p blastn -d " + dbFileName + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);
+                       blastCommand += (" -i " + (queryFileName+pid+toString(randNumber)) + " -o " + blastFileName+pid+toString(randNumber));
+               #else
+                       blastCommand =  "\"" + path + "blastall\" -p blastn -d " + "\"" + dbFileName + "\"" + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);
+                       blastCommand += (" -i " + (queryFileName+pid+toString(randNumber)) + " -o " + blastFileName+pid+toString(randNumber));
+                       //wrap entire string in ""
+                       blastCommand = "\"" + blastCommand + "\"";
+               #endif
+               
                system(blastCommand.c_str());
                
                ifstream m8FileHandle;
-               m->openInputFile(blastFileName+seq->getName(), m8FileHandle, "no error");
+               m->openInputFile(blastFileName+pid+toString(randNumber), m8FileHandle, "no error");
                
                string dummy;
                int templateAccession;
@@ -98,8 +256,8 @@ vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
                        topMatches.push_back(templateAccession);
                }
                m8FileHandle.close();
-               remove((queryFileName+seq->getName()).c_str());
-               remove((blastFileName+seq->getName()).c_str());
+               m->mothurRemove((queryFileName+pid+toString(randNumber)));
+               m->mothurRemove((blastFileName+pid+toString(randNumber)));
 
                return topMatches;
        }
@@ -118,23 +276,36 @@ vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n, int minPerID) {
                Scores.clear();
                
                ofstream queryFile;
-
-               m->openOutputFile((queryFileName+seq->getName()), queryFile);
+               int randNumber = rand();
+               string pid = scrubName(seq->getName());
+               
+               m->openOutputFile((queryFileName+pid+toString(randNumber)), queryFile);
                queryFile << '>' << seq->getName() << endl;
                queryFile << seq->getUnaligned() << endl;
                queryFile.close();
-                               
+//             cout << seq->getUnaligned() << endl;
                //      the goal here is to quickly survey the database to find the closest match.  To do this we are using the default
                //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
                //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
 //7000004128189528left 0       100             66      0       0       1       66      61      126     1e-31    131    
-               
-               string blastCommand = path + "blast/bin/megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
-               blastCommand += (" -i " + (queryFileName+seq->getName()) + " -o " + blastFileName+seq->getName());
+               string blastCommand;
+               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+                       blastCommand = path + "megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
+                       blastCommand += (" -i " + (queryFileName+pid+toString(randNumber)) + " -o " + blastFileName+pid+toString(randNumber));
+               #else
+               //blastCommand = path + "blast\\bin\\megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
+               //blastCommand += (" -i " + (queryFileName+toString(randNumber)) + " -o " + blastFileName+toString(randNumber));
+
+                       blastCommand =  "\"" + path + "megablast\" -e 1e-10 -d " + "\"" + dbFileName + "\"" + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
+                       blastCommand += (" -i " + (queryFileName+pid+toString(randNumber)) + " -o " + blastFileName+pid+toString(randNumber));
+                       //wrap entire string in ""
+                       blastCommand = "\"" + blastCommand + "\"";
+
+               #endif
                system(blastCommand.c_str());
 
                ifstream m8FileHandle;
-               m->openInputFile(blastFileName+seq->getName(), m8FileHandle, "no error");
+               m->openInputFile(blastFileName+pid+toString(randNumber), m8FileHandle, "no error");
        
                string dummy, eScore;
                int templateAccession;
@@ -142,7 +313,7 @@ vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n, int minPerID) {
                
                while(!m8FileHandle.eof()){
                        m8FileHandle >> dummy >> templateAccession >> searchScore >> numBases >> mismatch >> gap >> startQuery >> endQuery >> startRef >> endRef >> eScore >> score;
-                       //cout << dummy << '\t' << templateAccession << '\t' << searchScore << '\t' << numBases << '\t' << mismatch << '\t' << gap << '\t' << startQuery << '\t' << endQuery << '\t' << startRef << '\t' << endRef << '\t' << eScore << '\t' << score << endl; 
+//                     cout << dummy << '\t' << templateAccession << '\t' << searchScore << '\t' << numBases << '\t' << mismatch << '\t' << gap << '\t' << startQuery << '\t' << endQuery << '\t' << startRef << '\t' << endRef << '\t' << eScore << '\t' << score << endl; 
                        
                        //get rest of junk in line
                        //while (!m8FileHandle.eof())   {       char c = m8FileHandle.get(); if (c == 10 || c == 13){   break;  }else{ cout << c; }     } //
@@ -155,8 +326,8 @@ vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n, int minPerID) {
 //cout << templateAccession << endl;
                }
                m8FileHandle.close();
-               remove((queryFileName+seq->getName()).c_str());
-               remove((blastFileName+seq->getName()).c_str());
+               m->mothurRemove((queryFileName+pid+toString(randNumber)));
+               m->mothurRemove((blastFileName+pid+toString(randNumber)));
 //cout << "\n" ;               
                return topMatches;
        }
@@ -189,13 +360,19 @@ void BlastDB::generateDB() {
        try {
        
                //m->mothurOut("Generating the temporary BLAST database...\t"); cout.flush();
+                       
+               string formatdbCommand;
                
-               path = m->argv;
-               string tempPath = path;
-               for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
-               path = path.substr(0, (tempPath.find_last_of('m')));
-       
-               string formatdbCommand = path + "blast/bin/formatdb -p F -o T -i " + dbFileName;        //      format the database, -o option gives us the ability
+               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+                       formatdbCommand = path + "formatdb -p F -o T -i " + dbFileName; //      format the database, -o option gives us the ability
+               #else
+                       //formatdbCommand = path + "blast\\bin\\formatdb -p F -o T -i " + dbFileName;   //      format the database, -o option gives us the ability
+
+                       formatdbCommand = "\"" + path + "formatdb\" -p F -o T -i " + "\"" +  dbFileName + "\"";
+                       //wrap entire string in ""
+                       formatdbCommand = "\"" + formatdbCommand + "\"";
+               #endif
+               //cout << formatdbCommand << endl;
                system(formatdbCommand.c_str());                                                                //      to get the right sequence names, i think. -p F
                                                                                                                                        //      option tells formatdb that seqs are DNA, not prot
                //m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); cout.flush();
@@ -206,6 +383,24 @@ void BlastDB::generateDB() {
        }
 }
 /**************************************************************************************************/
+string BlastDB::scrubName(string seqName) {
+       try {
+               
+               string cleanName = "";
+               
+               for (int i = 0; i < seqName.length(); i++) {
+                       if (isalnum(seqName[i])) { cleanName += seqName[i]; }
+                       else { cleanName += "_";  }
+               }
+               
+               return cleanName;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "BlastDB", "scrubName");
+               exit(1);
+       }
+}
+/**************************************************************************************************/
 
 /**************************************************************************************************/