5 * Created by Pat Schloss on 12/22/08.
6 * Copyright 2008 Patrick D. Schloss. All rights reserved.
11 #include "database.hpp"
12 #include "sequence.hpp"
13 #include "blastdb.hpp"
15 /**************************************************************************************************/
17 BlastDB::BlastDB(float gO, float gE, float m, float mM) : Database(),
18 gapOpen(gO), gapExtend(gE), match(m), misMatch(mM) {
20 globaldata = GlobalData::getInstance();
23 int randNumber = rand();
24 dbFileName = toString(randNumber) + ".template.unaligned.fasta";
25 queryFileName = toString(randNumber) + ".candidate.unaligned.fasta";
26 blastFileName = toString(randNumber) + ".blast";
30 /**************************************************************************************************/
33 remove(queryFileName.c_str()); // let's clean stuff up and remove the temp files
34 remove(dbFileName.c_str()); // let's clean stuff up and remove the temp files
35 remove(blastFileName.c_str()); // let's clean stuff up and remove the temp files
37 /**************************************************************************************************/
38 //assumes you have added all the template sequences using the addSequence function and run generateDB.
39 vector<int> BlastDB::findClosestSequences(Sequence* seq, int n) {
41 vector<int> topMatches;
44 openOutputFile(queryFileName, queryFile);
45 queryFile << '>' << seq->getName() << endl;
46 queryFile << seq->getUnaligned() << endl;
49 // the goal here is to quickly survey the database to find the closest match. To do this we are using the default
50 // wordsize used in megablast. I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
51 // long. With this setting, it seems comparable in speed to the suffix tree approach.
53 string blastCommand = path + "blast/bin/blastall -p blastn -d " + dbFileName + " -m 8 -W 28 -v " + toString(n) + " -b " + toString(n);;
54 blastCommand += (" -i " + queryFileName + " -o " + blastFileName);
55 system(blastCommand.c_str());
57 ifstream m8FileHandle;
58 openInputFile(blastFileName, m8FileHandle);
61 int templateAccession;
64 while(!m8FileHandle.eof()){
65 m8FileHandle >> dummy >> templateAccession >> searchScore;
67 //get rest of junk in line
68 while (!m8FileHandle.eof()) { char c = m8FileHandle.get(); if (c == 10 || c == 13){ break; } }
71 topMatches.push_back(templateAccession);
78 errorOut(e, "BlastDB", "findClosestSequences");
83 /**************************************************************************************************/
84 //assumes you have added all the template sequences using the addSequence function and run generateDB.
85 vector<int> BlastDB::findClosestMegaBlast(Sequence* seq, int n) {
87 vector<int> topMatches;
90 openOutputFile(queryFileName, queryFile);
91 queryFile << '>' << seq->getName() << endl;
92 queryFile << seq->getUnaligned() << endl;
95 // the goal here is to quickly survey the database to find the closest match. To do this we are using the default
96 // wordsize used in megablast. I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
97 // long. With this setting, it seems comparable in speed to the suffix tree approach.
99 string blastCommand = path + "blast/bin/megablast -e 1e-10 -d " + dbFileName + " -m 8 -b " + toString(n) + " -v " + toString(n); //-W 28 -p blastn
100 blastCommand += (" -i " + queryFileName + " -o " + blastFileName);
101 system(blastCommand.c_str());
103 ifstream m8FileHandle;
104 openInputFile(blastFileName, m8FileHandle);
107 int templateAccession;
108 gobble(m8FileHandle);
110 while(!m8FileHandle.eof()){
111 m8FileHandle >> dummy >> templateAccession >> searchScore;
113 //get rest of junk in line
114 while (!m8FileHandle.eof()) { char c = m8FileHandle.get(); if (c == 10 || c == 13){ break; } }
116 gobble(m8FileHandle);
117 topMatches.push_back(templateAccession);
119 m8FileHandle.close();
123 catch(exception& e) {
124 errorOut(e, "BlastDB", "findClosest");
128 /**************************************************************************************************/
129 void BlastDB::addSequence(Sequence seq) {
132 ofstream unalignedFastaFile;
133 openOutputFileAppend(dbFileName, unalignedFastaFile);
135 // generating a fasta file with unaligned template
136 unalignedFastaFile << '>' << count << endl; // sequences, which will be input to formatdb
137 unalignedFastaFile << seq.getUnaligned() << endl;
138 unalignedFastaFile.close();
142 catch(exception& e) {
143 errorOut(e, "BlastDB", "addSequence");
147 /**************************************************************************************************/
148 void BlastDB::generateDB() {
151 //mothurOut("Generating the temporary BLAST database...\t"); cout.flush();
153 path = globaldata->argv;
154 path = path.substr(0, (path.find_last_of('m')));
156 string formatdbCommand = path + "blast/bin/formatdb -p F -o T -i " + dbFileName; // format the database, -o option gives us the ability
157 system(formatdbCommand.c_str()); // to get the right sequence names, i think. -p F
158 // option tells formatdb that seqs are DNA, not prot
159 //mothurOut("DONE."); mothurOutEndLine(); mothurOutEndLine(); cout.flush();
161 catch(exception& e) {
162 errorOut(e, "BlastDB", "generateDB");
167 /**************************************************************************************************/