用过网页版本 BLAST 的童鞋都会发现,提交的序列比对往往在几分钟,甚至几十秒就可以得到比对的结果;而通过调用 API 却要花费几十分钟或者更长的时间!这到底是为什么呢?
首先,我们来看一下提供了基于 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
(biopython==1.7.4)仅适用于 blastn,blastp,blastx,tblast 和 tblastx。qblast
函数还接受许多其他选项参数,这些参数基本上类似于我们可以在 BLAST 网页上设置的不同参数。我们在这里只重点介绍其中一些:
url_base
是设置用于在 Internet 上运行 BLAST 的基本 URL。默认情况下,它连接到 NCBI(即 url_base='https://blast.ncbi.nlm.nih.gov/Blast.cgi'),但是可以使用它连接到云端运行的 NCBI BLAST 实例。更多详细信息请参阅 qblast 功能的文档。format_type
关键字进行选择:“HTML”,“Text”,"ASN.1” 或 "XML"。默认值为 “XML”,因为这是解析器期望的格式。expect
用于设置期望值或 e-value 阈值。有关可选的 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 的实现前,我们先来看一下 NCBI BLAST 对于 API 使用的一些说明:
对于 API 的使用准则:
我们再来看一下 NCBIWWW 在源码层面的处理:
可以看到 NCBIWWW 从 20 秒的延迟开始,然后开始每隔一分钟执行一次 request 轮询,直至任务完成或者任务出现异常。
所以,总的来说,NCBI BLAST API 的使用准则,加上 NCBI BLAST 对用户请求的任务队列处理,甚至 NCBI BLAST 服务器共享资源的限制,以及总用户请求数,这些都可能成为 NCBIWWW.qblast()
异常耗时的原因,这其中还不算个人服务器的网络影响。综上种种原因,如果考虑使用 NCBIWWW.qblast()
执行频繁的序列在线批处理,或许不是一个好的解决方案。
最后,基于 Python 的 NCBI BLAST 在线批处理,如果你有更好的方法,欢迎留言交流。