用过网页版本 BLAST 的童鞋都会发现,提交的序列比对往往在几分钟,甚至几十秒就可以得到比对的结果;而通过调用 API 却要花费几十分钟或者更长的时间!这到底是为什么呢?

ncbi_blast.png


NCBIWWW 基本用法


首先,我们来看一下提供了基于 API 在线比对的 Biopython 模块。


Biopython 中的 BLAST 提供了 over the Internet 和 locally 两种选择:Bio.Blast.NCBIWWW 主要是基于 NCBI BLAST API 用于在线比对;Bio.Blast.Applications 模块则是调用本地安装好的 BLAST 程序以及数据库执行比对。在这里我们来重点看一下 Bio.Blast.NCBIWWW


Bio.Blast.NCBIWWW 模块中主要是通过 qblast() 函数来调用 BLAST 的在线版本。它具有三个非可选参数:


qblast 函数还接受许多其他选项参数,这些参数基本上类似于我们可以在 BLAST 网页上设置的不同参数。我们在这里只重点介绍其中一些:


有关可选的 BLAST 参数的更多信息,请参考 NCBI 自己的文档或 Biopython 内置的文档:

>>> from Bio.Blast import NCBIWWW
>>> help(NCBIWWW.qblast)


请注意,NCBI BLAST 网站上的默认设置与 qblast 上的默认设置不太相同。如果获得不同的结果,则需要检查参数(例如,e-value 值和 gap 值)。


例如,如果您要使用 BLASTN 在核苷酸数据库(nt)中搜索核苷酸序列,并且知道查询序列的 GI 号,则可以使用:

>>> from Bio.Blast import NCBIWWW
>>> result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")


另外,如果我们的查询序列已经存在于 FASTA 格式的文件中,则只需打开文件并以字符串形式读取此记录,然后将其用作查询参数:

>>> from Bio.Blast import NCBIWWW
>>> fasta_string = open("m_cold.fasta").read()
>>> result_handle = NCBIWWW.qblast("blastn", "nt", fasta_string)


我们还可以将 FASTA 文件作为 SeqRecord 对象进行读取,然后仅提供序列本身进行比对:

>>> from Bio.Blast import NCBIWWW
>>> from Bio import SeqIO
>>> record = SeqIO.read("m_cold.fasta", format="fasta")
>>> result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)


仅提供序列意味着 BLAST 将自动为您的序列分配一个标识符。您可能更喜欢使用 SeqRecord 对象的 format 方法来制作 FASTA 字符串(其中将包含现有标识符):

>>> from Bio.Blast import NCBIWWW
>>> from Bio import SeqIO
>>> record = SeqIO.read("m_cold.fasta", format="fasta")
>>> result_handle = NCBIWWW.qblast("blastn", "nt", record.format("fasta"))


无论给 qblast() 函数提供什么参数,都应在 handle 对象(默认为 XML 格式)中返回结果。下一步是将 XML 输出解析为表示搜索结果的 Python 对象,但是您可能想先保存输出文件的本地副本。在调试从 BLAST 结果中提取信息的代码时,我发现这特别有用(因为重新运行在线搜索速度很慢,并且浪费了 NCBI 计算机时间)。


我们需要小心一点,因为我们只能使用 result_handle.read() 读取一次 BLAST 输出——再次调用 result_handle.read() 会返回一个空字符串。

>>> with open("my_blast.xml", "w") as out_handle:
...     out_handle.write(result_handle.read())
...
>>> result_handle.close()


完成上面的操作后,结果将保存在文件 my_blast.xml 中,并且原始句柄已提取了所有数据(因此我们将其关闭了)。 但是,BLAST 解析器的解析功能采用了类似于文件句柄的对象,因此我们可以打开保存的文件进行输入:

>>> result_handle = open("my_blast.xml")


现在我们已经将 BLAST 结果重新放回了句柄中,下一步,如果我们准备对它们进行处理,我们可以参考 Biopython 中 Parsing BLAST output 部分的内容,这里不再说明。



NCBIWWW 实现


在了解 NCBIWWW 的实现前,我们先来看一下 NCBI BLAST 对于 API 使用的一些说明:

https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=DeveloperInfo

biopython.png


对于 API 的使用准则:


我们再来看一下 NCBIWWW 在源码层面的处理:

image.png


可以看到 NCBIWWW 从 20秒的延迟开始,然后开始每隔一分钟执行一次 request 轮询,直至任务完成或者任务出现异常。


所以,总的来说,NCBI BLAST API 的使用准则,加上 NCBI BLAST 对用户请求的任务队列处理,甚至 NCBI BLAST 服务器共享资源的限制,以及总用户请求数,这些都可能成为 NCBIWWW.qblast() 异常耗时的原因,这其中还不算个人服务器的网络影响。综上种种原因,如果考虑使用 NCBIWWW.qblast() 执行频繁的序列在线批处理,或许不是一个好的解决方案。


最后,基于 Python 的 NCBI BLAST 在线批处理,如果你有更好的方法,欢迎留言交流。