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 mm, float mM) : Database(),
  18 gapOpen(gO), gapExtend(gE), match(mm), misMatch(mM) {
  19         try {
  20                 count = 0;
  21
  22                 int randNumber = rand();
  23                 //int randNumber = 12345;
  24                 dbFileName = tag + toString(getpid()) + toString(randNumber) + ".template.unaligned.fasta";
  25                 queryFileName = tag + toString(getpid()) + toString(randNumber) + ".candidate.unaligned.fasta";
  26                 blastFileName = tag + toString(getpid()) + toString(randNumber) + ".blast";
  27
  28                 //make sure blast exists in the write place
  29                 path = m->argv;
  30                 string tempPath = path;
  31                 for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
  32                 path = path.substr(0, (tempPath.find_last_of('m')));
  33
  34                 string formatdbCommand;
  35 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
  36                 formatdbCommand = path + "blast/bin/formatdb";  //      format the database, -o option gives us the ability
  37 #else
  38                 formatdbCommand = path + "blast\\bin\\formatdb.exe";
  39                 //wrap entire string in ""
  40                 //formatdbCommand = "\"" + formatdbCommand + "\"";
  41 #endif
  42
  43                 //test to make sure formatdb exists
  44                 ifstream in;
  45                 formatdbCommand = m->getFullPathName(formatdbCommand);
  46                 int ableToOpen = m->openInputFile(formatdbCommand, in, "no error"); in.close();
  47                 if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + formatdbCommand + " file does not exist. mothur requires formatdb.exe to be in the ./blast/bin folder relative to the mothur.exe location."); m->mothurOutEndLine(); m->control_pressed = true; }
  48
  49                 string blastCommand;
  50 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
  51                 blastCommand = path + "blast/bin/blastall";     //      format the database, -o option gives us the ability
  52 #else
  53                 blastCommand = path + "blast\\bin\\blastall.exe";
  54                 //wrap entire string in ""
  55                 //blastCommand = "\"" + blastCommand + "\"";
  56 #endif
  57
  58                 //test to make sure formatdb exists
  59                 ifstream in2;
  60                 blastCommand = m->getFullPathName(blastCommand);
  61                 ableToOpen = m->openInputFile(blastCommand, in2, "no error"); in2.close();
  62                 if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + blastCommand + " file does not exist. mothur requires blastall.exe to be in the ./blast/bin folder relative to the mothur.exe location."); m->mothurOutEndLine(); m->control_pressed = true; }
  63
  64
  65                 string megablastCommand;
  66 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
  67                 megablastCommand = path + "blast/bin/megablast";        //      format the database, -o option gives us the ability
  68 #else
  69                 megablastCommand = path + "blast\\bin\\megablast.exe";
  70                 //wrap entire string in ""
  71                 //megablastCommand = "\"" + megablastCommand + "\"";
  72 #endif
  73
  74                 //test to make sure formatdb exists
  75                 ifstream in3;
  76                 megablastCommand = m->getFullPathName(megablastCommand);
  77                 ableToOpen = m->openInputFile(megablastCommand, in3, "no error"); in3.close();
  78                 if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " +  megablastCommand + " file does not exist. mothur requires megablast.exe to be in the ./blast/bin folder relative to the mothur.exe location."); m->mothurOutEndLine(); m->control_pressed = true; }
  79
  80         }
  81         catch(exception& e) {
  82                 m->errorOut(e, "BlastDB", "BlastDB");
  83                 exit(1);
  84         }
  85 }
  86 /**************************************************************************************************/
  87
  88 BlastDB::BlastDB() : Database() {
  89         try {
  90                 count = 0;
  91
  92                 int randNumber = rand();
  93                 //int randNumber = 12345;
  94                 dbFileName = toString(randNumber) + toString(getpid()) + ".template.unaligned.fasta";
  95                 queryFileName = toString(randNumber) + toString(getpid()) + ".candidate.unaligned.fasta";
  96                 blastFileName = toString(randNumber) + toString(getpid()) + ".blast";
  97
  98                 //make sure blast exists in the write place
  99                 path = m->argv;
 100                 string tempPath = path;
 101                 for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
 102                 path = path.substr(0, (tempPath.find_last_of('m')));
 103
 104                 string formatdbCommand;
 105 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 106                 formatdbCommand = path + "blast/bin/formatdb";  //      format the database, -o option gives us the ability
 107 #else
 108                 formatdbCommand = path + "blast\\bin\\formatdb.exe";
 109                 //wrap entire string in ""
 110                 //formatdbCommand = "\"" + formatdbCommand + "\"";
 111 #endif
 112
 113                 //test to make sure formatdb exists
 114                 ifstream in;
 115                 formatdbCommand = m->getFullPathName(formatdbCommand);
 116                 int ableToOpen = m->openInputFile(formatdbCommand, in, "no error"); in.close();
 117                 if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " +  formatdbCommand + " file does not exist. mothur requires formatdb.exe to be in the ./blast/bin folder relative to the mothur.exe location."); m->mothurOutEndLine(); m->control_pressed = true; }
 118
 119                 string blastCommand;
 120 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 121                 blastCommand = path + "blast/bin/blastall";     //      format the database, -o option gives us the ability
 122 #else
 123                 blastCommand = path + "blast\\bin\\blastall.exe";
 124                 //wrap entire string in ""
 125                 //blastCommand = "\"" + blastCommand + "\"";
 126 #endif
 127
 128                 //test to make sure formatdb exists
 129                 ifstream in2;
 130                 blastCommand = m->getFullPathName(blastCommand);
 131                 ableToOpen = m->openInputFile(blastCommand, in2, "no error"); in2.close();
 132                 if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + blastCommand + " file does not exist. mothur requires blastall.exe to be in the ./blast/bin folder relative to the mothur.exe location."); m->mothurOutEndLine(); m->control_pressed = true; }
 133
 134
 135                 string megablastCommand;
 136 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 137                 megablastCommand = path + "blast/bin/megablast";        //      format the database, -o option gives us the ability
 138 #else
 139                 megablastCommand = path + "blast\\bin\\megablast.exe";
 140                 //wrap entire string in ""
 141                 //megablastCommand = "\"" + megablastCommand + "\"";
 142 #endif
 143
 144                 //test to make sure formatdb exists
 145                 ifstream in3;
 146                 megablastCommand = m->getFullPathName(megablastCommand);
 147                 ableToOpen = m->openInputFile(megablastCommand, in3, "no error"); in3.close();
 148                 if(ableToOpen == 1) {   m->mothurOut("[ERROR]: " + megablastCommand + " file does not exist. mothur requires megablast.exe to be in the ./blast/bin folder relative to the mothur.exe location."); m->mothurOutEndLine(); m->control_pressed = true; }
 149
 150
 151         }
 152         catch(exception& e) {
 153                 m->errorOut(e, "BlastDB", "BlastDB");
 154                 exit(1);
 155         }
 156 }
 157
 158 /**************************************************************************************************/
 159
 160 BlastDB::~BlastDB(){
 161         try{
 162                 remove(queryFileName.c_str());                          //      let's clean stuff up and remove the temp files
 163                 remove(dbFileName.c_str());                                     //      let's clean stuff up and remove the temp files
 164                 remove((dbFileName+".nsq").c_str());                                    //      let's clean stuff up and remove the temp files
 165                 remove((dbFileName+".nsi").c_str());                                    //      let's clean stuff up and remove the temp files
 166                 remove((dbFileName+".nsd").c_str());                                    //      let's clean stuff up and remove the temp files
 167                 remove((dbFileName+".nin").c_str());                                    //      let's clean stuff up and remove the temp files
 168                 remove((dbFileName+".nhr").c_str());                                    //      let's clean stuff up and remove the temp files
 169                 remove(blastFileName.c_str());                          //      let's clean stuff up and remove the temp files
 170         }
 171         catch(exception& e) {
 172                 m->errorOut(e, "BlastDB", "~BlastDB");
 173                 exit(1);
 174         }
 175 }
 176 /**************************************************************************************************/
 177 //assumes you have added all the template sequences using the addSequence function and run generateDB.
 178 vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
 179         try{
 180                 vector<int> topMatches;
 181
 182                 ofstream queryFile;
 183                 int randNumber = rand();
 184                 m->openOutputFile((queryFileName+toString(getpid())+toString(randNumber)), queryFile);
 185                 queryFile << '>' << seq->getName() << endl;
 186                 queryFile << seq->getUnaligned() << endl;
 187                 queryFile.close();
 188
 189
 190                 //      the goal here is to quickly survey the database to find the closest match.  To do this we are using the default
 191                 //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
 192                 //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
 193
 194                 string blastCommand;
 195                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 196
 197                         blastCommand = path + "blast/bin/blastall -p blastn -d " + dbFileName + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);
 198                         blastCommand += (" -i " + (queryFileName+toString(getpid())+toString(randNumber)) + " -o " + blastFileName+toString(getpid())+toString(randNumber));
 199                 #else
 200                         blastCommand =  "\"" + path + "blast\\bin\\blastall\" -p blastn -d " + "\"" + dbFileName + "\"" + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);
 201                         blastCommand += (" -i " + (queryFileName+toString(getpid())+toString(randNumber)) + " -o " + blastFileName+toString(getpid())+toString(randNumber));
 202                         //wrap entire string in ""
 203                         blastCommand = "\"" + blastCommand + "\"";
 204                 #endif
 205
 206                 system(blastCommand.c_str());
 207
 208                 ifstream m8FileHandle;
 209                 m->openInputFile(blastFileName+toString(getpid())+toString(randNumber), m8FileHandle, "no error");
 210
 211                 string dummy;
 212                 int templateAccession;
 213                 m->gobble(m8FileHandle);
 214
 215                 while(!m8FileHandle.eof()){
 216                         m8FileHandle >> dummy >> templateAccession >> searchScore;
 217
 218                         //get rest of junk in line
 219                         while (!m8FileHandle.eof())     {       char c = m8FileHandle.get(); if (c == 10 || c == 13){   break;  }       }
 220
 221                         m->gobble(m8FileHandle);
 222                         topMatches.push_back(templateAccession);
 223                 }
 224                 m8FileHandle.close();
 225                 remove((queryFileName+toString(getpid())+toString(randNumber)).c_str());
 226                 remove((blastFileName+toString(getpid())+toString(randNumber)).c_str());
 227
 228                 return topMatches;
 229         }
 230         catch(exception& e) {
 231                 m->errorOut(e, "BlastDB", "findClosestSequences");
 232                 exit(1);
 233         }
 234
 235 }
 236 /**************************************************************************************************/
 237 //assumes you have added all the template sequences using the addSequence function and run generateDB.
 238 vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n, int minPerID) {
 239         try{
 240                 vector<int> topMatches;
 241                 float numBases, mismatch, gap, startQuery, endQuery, startRef, endRef, score;
 242                 Scores.clear();
 243
 244                 ofstream queryFile;
 245                 int randNumber = rand();
 246                 //int randNumber = 12345;
 247                 m->openOutputFile((queryFileName+toString(getpid())+toString(randNumber)), queryFile);
 248                 queryFile << '>' << seq->getName() << endl;
 249                 queryFile << seq->getUnaligned() << endl;
 250                 queryFile.close();
 251 //              cout << seq->getUnaligned() << endl;
 252                 //      the goal here is to quickly survey the database to find the closest match.  To do this we are using the default
 253                 //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
 254                 //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
 255 //7000004128189528left  0       100             66      0       0       1       66      61      126     1e-31    131
 256                 string blastCommand;
 257                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 258                         blastCommand = path + "blast/bin/megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
 259                         blastCommand += (" -i " + (queryFileName+toString(getpid())+toString(randNumber)) + " -o " + blastFileName+toString(getpid())+toString(randNumber));
 260                 #else
 261                 //blastCommand = path + "blast\\bin\\megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
 262                 //blastCommand += (" -i " + (queryFileName+toString(randNumber)) + " -o " + blastFileName+toString(randNumber));
 263
 264                         blastCommand =  "\"" + path + "blast\\bin\\megablast\" -e 1e-10 -d " + "\"" + dbFileName + "\"" + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
 265                         blastCommand += (" -i " + (queryFileName+toString(getpid())+toString(randNumber)) + " -o " + blastFileName+toString(getpid())+toString(randNumber));
 266                         //wrap entire string in ""
 267                         blastCommand = "\"" + blastCommand + "\"";
 268
 269                 #endif
 270                 system(blastCommand.c_str());
 271
 272                 ifstream m8FileHandle;
 273                 m->openInputFile(blastFileName+toString(getpid())+toString(randNumber), m8FileHandle, "no error");
 274
 275                 string dummy, eScore;
 276                 int templateAccession;
 277                 m->gobble(m8FileHandle);
 278
 279                 while(!m8FileHandle.eof()){
 280                         m8FileHandle >> dummy >> templateAccession >> searchScore >> numBases >> mismatch >> gap >> startQuery >> endQuery >> startRef >> endRef >> eScore >> score;
 281 //                      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;
 282
 283                         //get rest of junk in line
 284                         //while (!m8FileHandle.eof())   {       char c = m8FileHandle.get(); if (c == 10 || c == 13){   break;  }else{ cout << c; }     } //
 285                                 //cout << endl;
 286                         m->gobble(m8FileHandle);
 287                         if (searchScore >= minPerID) {
 288                                 topMatches.push_back(templateAccession);
 289                                 Scores.push_back(searchScore);
 290                         }
 291 //cout << templateAccession << endl;
 292                 }
 293                 m8FileHandle.close();
 294                 remove((queryFileName+toString(getpid())+toString(randNumber)).c_str());
 295                 remove((blastFileName+toString(getpid())+toString(randNumber)).c_str());
 296 //cout << "\n" ;
 297                 return topMatches;
 298         }
 299         catch(exception& e) {
 300                 m->errorOut(e, "BlastDB", "findClosestMegaBlast");
 301                 exit(1);
 302         }
 303 }
 304 /**************************************************************************************************/
 305 void BlastDB::addSequence(Sequence seq) {
 306         try {
 307
 308                 ofstream unalignedFastaFile;
 309                 m->openOutputFileAppend(dbFileName, unalignedFastaFile);
 310
 311                 //      generating a fasta file with unaligned template
 312                 unalignedFastaFile << '>' << count << endl;                                     //      sequences, which will be input to formatdb
 313                 unalignedFastaFile << seq.getUnaligned() << endl;
 314                 unalignedFastaFile.close();
 315
 316                 count++;
 317         }
 318         catch(exception& e) {
 319                 m->errorOut(e, "BlastDB", "addSequence");
 320                 exit(1);
 321         }
 322 }
 323 /**************************************************************************************************/
 324 void BlastDB::generateDB() {
 325         try {
 326
 327                 //m->mothurOut("Generating the temporary BLAST database...\t"); cout.flush();
 328
 329                 path = m->argv;
 330                 string tempPath = path;
 331                 for (int i = 0; i < path.length(); i++) { tempPath[i] = tolower(path[i]); }
 332                 path = path.substr(0, (tempPath.find_last_of('m')));
 333
 334                 string formatdbCommand;
 335
 336                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 337                         formatdbCommand = path + "blast/bin/formatdb -p F -o T -i " + dbFileName;       //      format the database, -o option gives us the ability
 338                 #else
 339                         //formatdbCommand = path + "blast\\bin\\formatdb -p F -o T -i " + dbFileName;   //      format the database, -o option gives us the ability
 340
 341                         formatdbCommand = "\"" + path + "blast\\bin\\formatdb\" -p F -o T -i " + "\"" +  dbFileName + "\"";
 342                         //wrap entire string in ""
 343                         formatdbCommand = "\"" + formatdbCommand + "\"";
 344                 #endif
 345                 //cout << formatdbCommand << endl;
 346                 system(formatdbCommand.c_str());                                                                //      to get the right sequence names, i think. -p F
 347                                                                                                                                         //      option tells formatdb that seqs are DNA, not prot
 348                 //m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); cout.flush();
 349         }
 350         catch(exception& e) {
 351                 m->errorOut(e, "BlastDB", "generateDB");
 352                 exit(1);
 353         }
 354 }
 355 /**************************************************************************************************/
 356
 357 /**************************************************************************************************/
 358