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