我有一个执行 BLAST 查询的脚本 (bl2seq)
脚本是这样工作的：
- 获取序列a，序列b
- 将序列 a 写入文件
- 将序列 b 写入文件 b
- 运行命令'bl2seq -i filea -j fileb -n blastn'
- 从 STDOUT 获取输出，解析
- 重复 2000 万次
bl2seq 程序不支持管道。 有没有办法做到这一点并避免写入/读取硬盘？
我正在使用 Python BTW。
【问题讨论】：
balstall 的 1500 万次查询 - 也是来自 Blast2 的命令）并通过 Google'ing 到达这里。请考虑在biostars.org重新提问
标签： python perl unix shell bioinformatics
我使用 R 脚本调用 blast2：
....
system("mkfifo seq1")
system("mkfifo seq2")
system("echo sequence1 > seq1"), wait = FALSE)
system("echo sequence2 > seq2"), wait = FALSE)
system("blast2 -p blastp -i seq1 -j seq2 -m 8", intern = TRUE)
....
这比从硬盘写入和读取要慢 2 倍（！）！
【讨论】：
哇。我想通了。
答案是使用python的子进程模块和管道！
编辑：忘了提到我正在使用 支持管道的 blast2。
（这是类的一部分）
def _query(self):
from subprocess import Popen, PIPE, STDOUT
pipe = Popen([BLAST,
'-p', 'blastn',
'-d', self.database,
'-m', '8'],
stdin=PIPE,
stdout=PIPE)
pipe.stdin.write('%s\n' % self.sequence)
print pipe.communicate()[0]
其中 self.database 是包含数据库文件名的字符串，即 'nt.fa' self.sequence 是一个包含查询序列的字符串
这会将输出打印到屏幕上，但您可以轻松地对其进行解析。没有慢速磁盘 I/O。没有缓慢的 XML 解析。我要为此写一个模块，放到github上。
此外，我还没有走到这一步，但我认为您可以执行多个查询，这样就不需要为每个查询读取爆炸数据库并将其加载到 RAM 中。
【讨论】：
你怎么知道 bl2seq 不支持管道？顺便说一句，管道是操作系统功能，而不是程序。如果您的 bl2seq 程序输出某些内容，无论是 STDOUT 还是文件，您都应该能够解析输出。检查 bl2seq 的帮助文件以获取输出到文件的选项，例如 -o 选项。然后就可以解析文件了。
此外，由于您使用的是 Python，因此您可以使用 BioPython 模块。
【讨论】：
根据您运行的操作系统，您可能可以使用bash's process substitution 之类的内容。我不确定您如何在 Python 中进行设置，但您基本上使用的是命名管道（或命名文件描述符）。如果bl2seq 尝试在文件中查找，这将不起作用，但如果它只是按顺序读取它们应该会起作用。
【讨论】：
这是来自BioPerl 的bl2seq 程序吗？如果是这样，看起来您无法对其进行管道处理。但是，您可以使用 Bio::Tools::Run::AnalysisFactory::Pise 编写自己的 hack，这是推荐的方法。不过，您必须在 Perl 中进行操作。
如果这是另一个bl2seq，则忽略该消息。无论如何，您可能应该提供更多细节。
【讨论】：