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