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(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         dbFileName = toString(randNumber) + ".template.unaligned.fasta";
  24         queryFileName = toString(randNumber) + ".candidate.unaligned.fasta";
  25         blastFileName = toString(randNumber) + ".blast";
  26
  27 }
  28 /**************************************************************************************************/
  29
  30 BlastDB::BlastDB() : Database() {
  31         try {
  32                 count = 0;
  33
  34                 int randNumber = rand();
  35                 dbFileName = toString(randNumber) + ".template.unaligned.fasta";
  36                 queryFileName = toString(randNumber) + ".candidate.unaligned.fasta";
  37                 blastFileName = toString(randNumber) + ".blast";
  38         }
  39         catch(exception& e) {
  40                 m->errorOut(e, "BlastDB", "BlastDB");
  41                 exit(1);
  42         }
  43 }
  44
  45 /**************************************************************************************************/
  46
  47 BlastDB::~BlastDB(){
  48         try{
  49                 remove(queryFileName.c_str());                          //      let's clean stuff up and remove the temp files
  50                 remove(dbFileName.c_str());                                     //      let's clean stuff up and remove the temp files
  51                 remove((dbFileName+".nsq").c_str());                                    //      let's clean stuff up and remove the temp files
  52                 remove((dbFileName+".nsi").c_str());                                    //      let's clean stuff up and remove the temp files
  53                 remove((dbFileName+".nsd").c_str());                                    //      let's clean stuff up and remove the temp files
  54                 remove((dbFileName+".nin").c_str());                                    //      let's clean stuff up and remove the temp files
  55                 remove((dbFileName+".nhr").c_str());                                    //      let's clean stuff up and remove the temp files
  56                 remove(blastFileName.c_str());                          //      let's clean stuff up and remove the temp files
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "BlastDB", "~BlastDB");
  60                 exit(1);
  61         }
  62 }
  63 /**************************************************************************************************/
  64 //assumes you have added all the template sequences using the addSequence function and run generateDB.
  65 vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
  66         try{
  67                 vector<int> topMatches;
  68
  69                 ofstream queryFile;
  70                 m->openOutputFile((queryFileName+seq->getName()), queryFile);
  71                 queryFile << '>' << seq->getName() << endl;
  72                 queryFile << seq->getUnaligned() << endl;
  73                 queryFile.close();
  74
  75
  76                 //      the goal here is to quickly survey the database to find the closest match.  To do this we are using the default
  77                 //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
  78                 //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
  79
  80                 string blastCommand = path + "blast/bin/blastall -p blastn -d " + dbFileName + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);;
  81                 blastCommand += (" -i " + (queryFileName+seq->getName()) + " -o " + blastFileName+seq->getName());
  82                 system(blastCommand.c_str());
  83
  84                 ifstream m8FileHandle;
  85                 m->openInputFile(blastFileName+seq->getName(), m8FileHandle, "no error");
  86
  87                 string dummy;
  88                 int templateAccession;
  89                 m->gobble(m8FileHandle);
  90
  91                 while(!m8FileHandle.eof()){
  92                         m8FileHandle >> dummy >> templateAccession >> searchScore;
  93
  94                         //get rest of junk in line
  95                         while (!m8FileHandle.eof())     {       char c = m8FileHandle.get(); if (c == 10 || c == 13){   break;  }       }
  96
  97                         m->gobble(m8FileHandle);
  98                         topMatches.push_back(templateAccession);
  99                 }
 100                 m8FileHandle.close();
 101                 remove((queryFileName+seq->getName()).c_str());
 102                 remove((blastFileName+seq->getName()).c_str());
 103
 104                 return topMatches;
 105         }
 106         catch(exception& e) {
 107                 m->errorOut(e, "BlastDB", "findClosestSequences");
 108                 exit(1);
 109         }
 110
 111 }
 112 /**************************************************************************************************/
 113 //assumes you have added all the template sequences using the addSequence function and run generateDB.
 114 vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n) {
 115         try{
 116                 vector<int> topMatches;
 117
 118                 ofstream queryFile;
 119                 m->openOutputFile((queryFileName+seq->getName()), queryFile);
 120                 queryFile << '>' << seq->getName() << endl;
 121                 queryFile << seq->getUnaligned() << endl;
 122                 queryFile.close();
 123
 124                 //      the goal here is to quickly survey the database to find the closest match.  To do this we are using the default
 125                 //      wordsize used in megablast.  I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
 126                 //      long.  With this setting, it seems comparable in speed to the suffix tree approach.
 127
 128                 string blastCommand = path + "blast/bin/megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
 129                 blastCommand += (" -i " + (queryFileName+seq->getName()) + " -o " + blastFileName+seq->getName());
 130                 system(blastCommand.c_str());
 131
 132                 ifstream m8FileHandle;
 133                 m->openInputFile(blastFileName+seq->getName(), m8FileHandle, "no error");
 134
 135                 string dummy;
 136                 int templateAccession;
 137                 m->gobble(m8FileHandle);
 138
 139                 while(!m8FileHandle.eof()){
 140                         m8FileHandle >> dummy >> templateAccession >> searchScore;
 141                         //cout << templateAccession << '\t' << searchScore << endl;
 142
 143                         //get rest of junk in line
 144                         while (!m8FileHandle.eof())     {       char c = m8FileHandle.get(); if (c == 10 || c == 13){   break;  }       }
 145
 146                         m->gobble(m8FileHandle);
 147                         topMatches.push_back(templateAccession);
 148 //cout << templateAccession << endl;
 149                 }
 150                 m8FileHandle.close();
 151                 remove((queryFileName+seq->getName()).c_str());
 152                 remove((blastFileName+seq->getName()).c_str());
 153 //cout << "\n" ;
 154                 return topMatches;
 155         }
 156         catch(exception& e) {
 157                 m->errorOut(e, "BlastDB", "findClosest");
 158                 exit(1);
 159         }
 160 }
 161 /**************************************************************************************************/
 162 void BlastDB::addSequence(Sequence seq) {
 163         try {
 164
 165                 ofstream unalignedFastaFile;
 166                 m->openOutputFileAppend(dbFileName, unalignedFastaFile);
 167
 168                 //      generating a fasta file with unaligned template
 169                 unalignedFastaFile << '>' << count << endl;                                     //      sequences, which will be input to formatdb
 170                 unalignedFastaFile << seq.getUnaligned() << endl;
 171                 unalignedFastaFile.close();
 172
 173                 count++;
 174         }
 175         catch(exception& e) {
 176                 m->errorOut(e, "BlastDB", "addSequence");
 177                 exit(1);
 178         }
 179 }
 180 /**************************************************************************************************/
 181 void BlastDB::generateDB() {
 182         try {
 183
 184                 //m->mothurOut("Generating the temporary BLAST database...\t"); cout.flush();
 185
 186                 path = m->argv;
 187                 path = path.substr(0, (path.find_last_of('m')));
 188
 189                 string formatdbCommand = path + "blast/bin/formatdb -p F -o T -i " + dbFileName;        //      format the database, -o option gives us the ability
 190                 system(formatdbCommand.c_str());                                                                //      to get the right sequence names, i think. -p F
 191                                                                                                                                         //      option tells formatdb that seqs are DNA, not prot
 192                 //m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); cout.flush();
 193         }
 194         catch(exception& e) {
 195                 m->errorOut(e, "BlastDB", "generateDB");
 196                 exit(1);
 197         }
 198 }
 199 /**************************************************************************************************/
 200
 201 /**************************************************************************************************/
 202