5 * Created by Pat Schloss on 12/22/08.
6 * Copyright 2008 Patrick D. Schloss. All rights reserved.
13 #include "database.hpp"
14 #include "sequence.hpp"
15 #include "blastdb.hpp"
17 /**************************************************************************************************/
19 BlastDB::BlastDB(string fastaFileName, float gO, float gE, float m, float mM) : Database(fastaFileName),
20 gapOpen(gO), gapExtend(gE), match(m), misMatch(mM) {
22 cout << "Generating the temporary BLAST database...\t"; cout.flush();
24 int randNumber = rand();
25 dbFileName = toString(randNumber) + ".template.unaligned.fasta";
26 queryFileName = toString(randNumber) + ".candidate.unaligned.fasta";
27 blastFileName = toString(randNumber) + ".blast";
30 ofstream unalignedFastaFile;
31 openOutputFile(dbFileName, unalignedFastaFile);
33 for(int i=0;i<numSeqs;i++){ // generating a fasta file with unaligned template
34 unalignedFastaFile << '>' << i << endl; // sequences, which will be input to formatdb
35 unalignedFastaFile << templateSequences[i]->getUnaligned() << endl;
37 unalignedFastaFile.close();
39 string formatdbCommand = "~/Pipeline/src/cpp/production/blast/bin/formatdb -p F -o T -i " + dbFileName; // format the database, -o option gives us the ability
40 system(formatdbCommand.c_str()); // to get the right sequence names, i think. -p F
41 // option tells formatdb that seqs are DNA, not prot
42 cout << "DONE." << endl << endl; cout.flush();
43 emptySequence = new Sequence();
44 emptySequence->setName("no_match");
45 emptySequence->setUnaligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
46 emptySequence->setAligned("XXXXXXXXXXXXXXXXXXXXXXXXXXXXX");
51 /**************************************************************************************************/
54 remove(queryFileName.c_str()); // let's clean stuff up and remove the temp files
55 remove(dbFileName.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
59 /**************************************************************************************************/
61 Sequence* BlastDB::findClosestSequence(Sequence* candidate){
64 openOutputFile(queryFileName, queryFile);
65 queryFile << '>' << candidate->getName() << endl;
66 queryFile << candidate->getUnaligned() << endl;
70 // the goal here is to quickly survey the database to find the closest match. To do this we are using the default
71 // wordsize used in megablast. I'm sure we're sacrificing accuracy for speed, but anyother way would take way too
72 // long. With this setting, it seems comparable in speed to the suffix tree approach.
74 string blastCommand = "~/Pipeline/src/cpp/production/blast/bin/blastall -p blastn -d " + dbFileName + " -b 1 -m 8 -W 28";
75 blastCommand += (" -i " + queryFileName + " -o " + blastFileName);
76 system(blastCommand.c_str());
78 ifstream m8FileHandle;
79 openInputFile(blastFileName, m8FileHandle);
82 int templateAccession;
84 if(!m8FileHandle.eof()){
85 m8FileHandle >> dummy >> templateAccession >> searchScore;
87 return templateSequences[templateAccession];
95 /**************************************************************************************************/