blastdb.cpp

   1 /*
   2  *  blastdb.cpp
   3  *
   4  *
   5  *  Created by Pat Schloss on 12/22/08.
   6  *  Copyright 2008 Patrick D. Schloss. All rights reserved.
   7  *
   8  */
   9
  10
  11 #include "database.hpp"
  12 #include "sequence.hpp"
  13 #include "blastdb.hpp"
  14
  15 /**************************************************************************************************/
  16
  17 BlastDB::BlastDB(string tag, float gO, float gE, float m, float mM) : Database(),
  18 gapOpen(gO), gapExtend(gE), match(m), misMatch(mM) {
  19
  20         count = 0;
  21
  22         int randNumber = rand();
  23         //int randNumber = 12345;
  24         dbFileName = tag + toString(randNumber) + ".template.unaligned.fasta";
  25         queryFileName = tag + toString(randNumber) + ".candidate.unaligned.fasta";
  26         blastFileName = tag + toString(randNumber) + ".blast";
  27
  28 }
  29 /**************************************************************************************************/
  30
  31 BlastDB::BlastDB() : Database() {
  32         try {
  33                 count = 0;
  34
  35                 int randNumber = rand();
  36                 //int randNumber = 12345;
  37                 dbFileName = toString(randNumber) + ".template.unaligned.fasta";
  38                 queryFileName = toString(randNumber) + ".candidate.unaligned.fasta";
  39                 blastFileName = toString(randNumber) + ".blast";
  40         }
  41         catch(exception& e) {
  42                 m->errorOut(e, "BlastDB", "BlastDB");
  43                 exit(1);
  44         }
  45 }
  46
  47 /**************************************************************************************************/
  48
  49 BlastDB::~BlastDB(){
  50         try{
  51                 remove(queryFileName.c_str());                          //      let's clean stuff up and remove the temp files
  52                 remove(dbFileName.c_str());                                     //      let's clean stuff up and remove the temp files
  53                 remove((dbFileName+".nsq").c_str());                                    //      let's clean stuff up and remove the temp files
  54                 remove((dbFileName+".nsi").c_str());                                    //      let's clean stuff up and remove the temp files
  55                 remove((dbFileName+".nsd").c_str());                                    //      let's clean stuff up and remove the temp files
  56                 remove((dbFileName+".nin").c_str());                                    //      let's clean stuff up and remove the temp files
  57                 remove((dbFileName+".nhr").c_str());                                    //      let's clean stuff up and remove the temp files
  58                 remove(blastFileName.c_str());                          //      let's clean stuff up and remove the temp files
  59         }
  60         catch(exception& e) {
  61                 m->errorOut(e, "BlastDB", "~BlastDB");
  62                 exit(1);
  63         }
  64 }
  65 /**************************************************************************************************/
  66 //assumes you have added all the template sequences using the addSequence function and run generateDB.
  67 vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
  68         try{
  69                 vector<int> topMatches;
  70
  71                 ofstream queryFile;
  72                 int randNumber = rand();
  73                 m->openOutputFile((queryFileName+toString(randNumber)), queryFile);
  74                 queryFile << '>' << seq->getName() << endl;
  75                 queryFile << seq->getUnaligned() << endl;
  76                 queryFile.close();
  77
  78
  79                 //      the goal here is to quickly survey the database to find the closest match.  To do this we are using the default
  80                 //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
  81                 //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
  82
  83                 string blastCommand;
  84                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
  85
  86                         blastCommand = path + "blast/bin/blastall -p blastn -d " + dbFileName + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);
  87                         blastCommand += (" -i " + (queryFileName+toString(randNumber)) + " -o " + blastFileName+toString(randNumber));
  88                 #else
  89                         blastCommand =  "\"" + path + "blast\\bin\\blastall\" -p blastn -d " + "\"" + dbFileName + "\"" + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);
  90                         blastCommand += (" -i " + (queryFileName+toString(randNumber)) + " -o " + blastFileName+toString(randNumber));
  91                         //wrap entire string in ""
  92                         blastCommand = "\"" + blastCommand + "\"";
  93                 #endif
  94
  95                 system(blastCommand.c_str());
  96
  97                 ifstream m8FileHandle;
  98                 m->openInputFile(blastFileName+toString(randNumber), m8FileHandle, "no error");
  99
 100                 string dummy;
 101                 int templateAccession;
 102                 m->gobble(m8FileHandle);
 103
 104                 while(!m8FileHandle.eof()){
 105                         m8FileHandle >> dummy >> templateAccession >> searchScore;
 106
 107                         //get rest of junk in line
 108                         while (!m8FileHandle.eof())     {       char c = m8FileHandle.get(); if (c == 10 || c == 13){   break;  }       }
 109
 110                         m->gobble(m8FileHandle);
 111                         topMatches.push_back(templateAccession);
 112                 }
 113                 m8FileHandle.close();
 114                 remove((queryFileName+toString(randNumber)).c_str());
 115                 remove((blastFileName+toString(randNumber)).c_str());
 116
 117                 return topMatches;
 118         }
 119         catch(exception& e) {
 120                 m->errorOut(e, "BlastDB", "findClosestSequences");
 121                 exit(1);
 122         }
 123
 124 }
 125 /**************************************************************************************************/
 126 //assumes you have added all the template sequences using the addSequence function and run generateDB.
 127 vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n, int minPerID) {
 128         try{
 129                 vector<int> topMatches;
 130                 float numBases, mismatch, gap, startQuery, endQuery, startRef, endRef, score;
 131                 Scores.clear();
 132
 133                 ofstream queryFile;
 134                 int randNumber = rand();
 135                 //int randNumber = 12345;
 136                 m->openOutputFile((queryFileName+toString(randNumber)), queryFile);
 137                 queryFile << '>' << seq->getName() << endl;
 138                 queryFile << seq->getUnaligned() << endl;
 139                 queryFile.close();
 140 //              cout << seq->getUnaligned() << endl;
 141                 //      the goal here is to quickly survey the database to find the closest match.  To do this we are using the default
 142                 //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
 143                 //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
 144 //7000004128189528left  0       100             66      0       0       1       66      61      126     1e-31    131
 145                 string blastCommand;
 146                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 147                         blastCommand = path + "blast/bin/megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
 148                         blastCommand += (" -i " + (queryFileName+toString(randNumber)) + " -o " + blastFileName+toString(randNumber));
 149                 #else
 150                 //blastCommand = path + "blast\\bin\\megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
 151                 //blastCommand += (" -i " + (queryFileName+toString(randNumber)) + " -o " + blastFileName+toString(randNumber));
 152
 153                         blastCommand =  "\"" + path + "blast\\bin\\megablast\" -e 1e-10 -d " + "\"" + dbFileName + "\"" + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
 154                         blastCommand += (" -i " + (queryFileName+toString(randNumber)) + " -o " + blastFileName+toString(randNumber));
 155                         //wrap entire string in ""
 156                         blastCommand = "\"" + blastCommand + "\"";
 157
 158                 #endif
 159                 //cout << blastCommand << endl;
 160                 system(blastCommand.c_str());
 161
 162                 ifstream m8FileHandle;
 163                 m->openInputFile(blastFileName+toString(randNumber), m8FileHandle, "no error");
 164
 165                 string dummy, eScore;
 166                 int templateAccession;
 167                 m->gobble(m8FileHandle);
 168
 169                 while(!m8FileHandle.eof()){
 170                         m8FileHandle >> dummy >> templateAccession >> searchScore >> numBases >> mismatch >> gap >> startQuery >> endQuery >> startRef >> endRef >> eScore >> score;
 171 //                      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;
 172
 173                         //get rest of junk in line
 174                         //while (!m8FileHandle.eof())   {       char c = m8FileHandle.get(); if (c == 10 || c == 13){   break;  }else{ cout << c; }     } //
 175                                 //cout << endl;
 176                         m->gobble(m8FileHandle);
 177                         if (searchScore >= minPerID) {
 178                                 topMatches.push_back(templateAccession);
 179                                 Scores.push_back(searchScore);
 180                         }
 181 //cout << templateAccession << endl;
 182                 }
 183                 m8FileHandle.close();
 184                 remove((queryFileName+toString(randNumber)).c_str());
 185                 remove((blastFileName+toString(randNumber)).c_str());
 186 //cout << "\n" ;
 187                 return topMatches;
 188         }
 189         catch(exception& e) {
 190                 m->errorOut(e, "BlastDB", "findClosestMegaBlast");
 191                 exit(1);
 192         }
 193 }
 194 /**************************************************************************************************/
 195 void BlastDB::addSequence(Sequence seq) {
 196         try {
 197
 198                 ofstream unalignedFastaFile;
 199                 m->openOutputFileAppend(dbFileName, unalignedFastaFile);
 200
 201                 //      generating a fasta file with unaligned template
 202                 unalignedFastaFile << '>' << count << endl;                                     //      sequences, which will be input to formatdb
 203                 unalignedFastaFile << seq.getUnaligned() << endl;
 204                 unalignedFastaFile.close();
 205
 206                 count++;
 207         }
 208         catch(exception& e) {
 209                 m->errorOut(e, "BlastDB", "addSequence");
 210                 exit(1);
 211         }
 212 }
 213 /**************************************************************************************************/
 214 void BlastDB::generateDB() {
 215         try {
 216
 217                 //m->mothurOut("Generating the temporary BLAST database...\t"); cout.flush();
 218
 219                 path = m->argv;
 220                 string tempPath = path;
 221                 for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
 222                 path = path.substr(0, (tempPath.find_last_of('m')));
 223
 224                 string formatdbCommand;
 225
 226                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 227                         formatdbCommand = path + "blast/bin/formatdb -p F -o T -i " + dbFileName;       //      format the database, -o option gives us the ability
 228                 #else
 229                         //formatdbCommand = path + "blast\\bin\\formatdb -p F -o T -i " + dbFileName;   //      format the database, -o option gives us the ability
 230
 231                         formatdbCommand = "\"" + path + "blast\\bin\\formatdb\" -p F -o T -i " + "\"" +  dbFileName + "\"";
 232                         //wrap entire string in ""
 233                         formatdbCommand = "\"" + formatdbCommand + "\"";
 234                 #endif
 235                 //cout << formatdbCommand << endl;
 236                 system(formatdbCommand.c_str());                                                                //      to get the right sequence names, i think. -p F
 237                                                                                                                                         //      option tells formatdb that seqs are DNA, not prot
 238                 //m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); cout.flush();
 239         }
 240         catch(exception& e) {
 241                 m->errorOut(e, "BlastDB", "generateDB");
 242                 exit(1);
 243         }
 244 }
 245 /**************************************************************************************************/
 246
 247 /**************************************************************************************************/
 248