第7章 BLAST

嗨,每个人都喜欢BLAST,对吧?我是指,通过BLAST把你的序列和世界上已知的序列 比较是多么简单方便啊。不过,这章当然不是讲Blast有多么酷,因为我们都已经 知道了。这章是来解决使用Blast的一些麻烦地方——处理大量的BLAST比对结果数据 通常是困难的,还有怎么自动运行BLAST序列比对。

幸运的是,Biopython社区的人早就了解了这些难处。所以,他们已经发展了很多 工具来简化BlAST使用和结果处理。这章会具体讲解怎么用这些工具来做些有用的 事情。

使用BLAST通常可以分成2个步。这两步都可以用上Biopython。第一步,提交你的查询 序列,运行BLAST,并得到输出数据。第二步,用Python解析BLAST的输出,并作进一步 分析。

你第一次接触并运行BLAST也许就是通过NCBI的web服务。事实上,你可以通过多种方式 (这些方式可以分成几类)来使用BLAST。这些方式最重要的区别在于你是在你的自己 电脑上运行一个本地BLAST,还是在远程服务器(另外一台电脑,通常是NCBI的服务器)上运行。 我们将在一个Python脚本里调用NCBI在线BLAST服务来开始这章的内容。

注意: 接下来的 第 8 章介绍的 Bio.SearchIO 是一个 Biopython实验性质的模块。我们准备最终用它来替换原来的 Bio.Blast 模块。 因为它提供了一个更为通用的序列搜索相关的框架。不过,在这个模块的稳定版本发布之前,在实际工作中的代码里,请继续用 Bio.Blast 模块来处理NCBI BLAST。

7.1 通过Internet运行BLAST

我们使用 Bio.Blast.NCBIWWW 模块的函数 qblast() 来调用在线版本的BLAST。 这个函数有3个必需的参数:

  • 第一个参数是用来搜索的blast程序,这是小写的字符串。对这个参数的选项和描述可以在 http://www.ncbi.nlm.nih.gov/BLAST/blast_program.shtml. 查到。目前 qblast 只支持 blastn, blastp, blastx, tblast 和 tblastx.
  • 第二个参数指定了将要搜索的数据库。同样地,这个参数的选项也可以在 http://www.ncbi.nlm.nih.gov/BLAST/blast_databases.shtml. 查到
  • 第三个参数是一个包含你要查询序列的字符串。这个字符串可以是序列的本身 (fasta格式的),或者是一个类似GI的id。

qblast 函数还可以接受许多其他选项和参数。这些参数基本上类似于你在BLAST网站页面 能够设置的参数。在这里我们只强调其中的一些:

  • qblast 函数可以返回多种格式的BLAST结果。你可以通过可选参数 format_type 指定格式关键字为:"HTML", "Text", "ASN.1", 或 "XML" 。默认 格式是 "XML" ,这是解析器期望的格式, 7.3 节对其有详细的描述。
  • 参数 expect 指定期望值,即阀值 e-value。

更多可选的BLAST参数,请参照NCBI的文档,或者是Biopython内置的文档。

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

请注意,NCBI BLAST 网站上的默认参数和QBLAST的默认参数不完全相同。如果你得到了 不同的结果,你就需要检查下参数设置 (比如,e-value阈值和gap值).

举个例子,如果你有条核酸序列,想使用BLAST对核酸数据库(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会自动分配给你一个ID。你可能更喜欢用 SeqRecord 对象的format方法来包装一个fasta字符串,因为这个对象会包含fasta文件中已有的ID

>>> 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"))

如果你的序列在一个非FASTA格式的文件中并且你用 Bio.SeqIO (看第 5 章) 把序列取出来了,那么这个方法更有用。

不论你给 qblast() 函数提供了什么参数,都应该返回一个handle object的结果( 默认是XML格式)。下一步就是将这个XML输出解析为代表BLAST搜索结果的Python 对象( 7.3 )。 不过,也许你想先把这个XML输出保存一个本地文件副本。当调试从BLAST结果提取信息的代码的时候,我发现这样做 尤其有用。(因为重新运行在线BLAST搜索很慢并且会浪费NCBI服务器的运行时间)。

这里我们需要注意下:因为用 result_handle.read() 来读取BLAST结果只能用一次 - 再次调用 result_handle.read() 会返回一个空的字符串.

>>> save_file = open("my_blast.xml", "w")
>>> save_file.write(result_handle.read())
>>> save_file.close()
>>> result_handle.close()

这些做好后,结果已经存储在 my_blast.xml 文件中了并且原先的handle中的数据 已经被全部提取出来了(所以我们把它关闭了)。但是,BLAST解析器的 parse 函数(描述见 7.3 ) 采用一个文件句柄类的对象,所以我们只需打开已经保存的文件作为输入。

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

既然现在已经把BLAST的结果又一次读回handle,我们可以分析下这些结果。所以我们正好可以去读 关于结果解析的章节(看下面 7.3 )。你现在也许想跳过去看吧 ...

7.2 本地运行BLAST

7.2.1 介绍

在本地运行BLAST(跟通过internet运行比,见 7.1 ) 至少有2个主要优点:

  • 本地运行BLAST可能比通过internet运行更快;
  • 本地运行可以让你建立自己的数据库来对序列进行搜索。

处理有版权的或者没有发表的序列数据也许是本地运行BLAST的另一个原因。你也许 不能泄露这些序列数据,所以没法提交给NCBI来BLAST。

不幸的是,本地运行也有些缺点 - 安装所有的东东并成功运行需要花些力气:

  • 本地运行BLAST需要你安装相关命令行工具。
  • 本地运行BLAST需要安装一个很大的BLAST的数据库(并且需要保持数据更新).

更令人困惑的是,至少有4种不同的BLAST安装程序包,并且还有其他的一些工具能 产生类似的BLAST 输出文件,比如BLAT。

7.2.2 单机版的NCBI老版本BLAST

NCBI “老版本” BLAST 包括命令行工具 blastallblastpgprpsblast 。 这是NCBI发布它的替代品BLAST+ 前使用最为广泛的单机版BLAST工具。

Bio.Blast.Applications 模块有个对老版本NCBI BLAST 工具像 blastallblastpgprpsblast 的封装, 并且在 Bio.Blast.NCBIStandalone 还有个辅助函数。 这些东东现在都被认为是过时的,并且当用户们迁移到BLAST+程序套件后,这些都会被弃用, 最终从Biopython删除。

为了减少你的困惑,我们在这个指南中不会提到怎么从Biopython调用这些老版本的工具。 如果你有兴趣,可以看下在Biopython 1.52中包含的基本指南。(看下 biopython-1.52.tar.gz 或者 biopython-1.52.zip 中Doc目录下的指南的PDF文件 或者 HTML 文件)。

7.2.3 单机版 NCBI BLAST+

NCBI “新版本”的 BLAST+ 在2009年发布。它替代了原来老版本的BLAST程序包。Bio.Blast.Applications 模块 包装了这些新工具像 blastn , blastp , blastx , tblastn , tblastx (这些以前都是由 blastall 处理)。而 rpsblastrpstblastn (替代了 原来的 rpsblast )。我们这里不包括对 makeblastdb 的包装,它在BLAST+中用于从FASTA文件 建立一个本地BLAST数据库,还有其在老版本BLAST中的等效工具 formatdb

这节将简要地介绍怎样在Python中使用这些工具。如果你已经阅读了并试过 6.4 节的序列联配(alignment)工具,下面介绍 的方法应该是很简单直接的。首先,我们构建一个命令行字符串(就像你使用单机版 BLAST的时候,在终端打入命令行一样)。然后,我们在Python中运行这个命令。

举个例子,你有个FASTA格式的核酸序列文件,你想用它通过BLASTX(翻译)来搜索 非冗余(NR)蛋白质数据库。如果你(或者你的系统管理员)下载并安装好了这个数据库, 那么你只要运行:

blastx -query opuntia.fasta -db nr -out opuntia.xml -evalue 0.001 -outfmt 5

这样就完成了运行BLASTX查找非冗余蛋白质数据库,用0.001的e值并产生XML格式的 输出结果文件(这样我们可以继续下一步解析)。在我的电脑上运行这条命令花了大约6分钟 - 这就是为什么我们需要保存输出到文件。这样我们就可以在需要时重复任何基于这个输出的分析。

在Biopython中,我们可以用NCBI BLASTX包装模块 Bio.Blast.Applications 来构建 命令行字符串并运行它:

>>> from Bio.Blast.Applications import NcbiblastxCommandline
>>> help(NcbiblastxCommandline)
...
>>> blastx_cline = NcbiblastxCommandline(query="opuntia.fasta", db="nr", evalue=0.001,
...                                      outfmt=5, out="opuntia.xml")
>>> blastx_cline
NcbiblastxCommandline(cmd='blastx', out='opuntia.xml', outfmt=5, query='opuntia.fasta',
db='nr', evalue=0.001)
>>> print blastx_cline
blastx -out opuntia.xml -outfmt 5 -query opuntia.fasta -db nr -evalue 0.001
>>> stdout, stderr = blastx_cline()

在这个例子中,终端里应该没有任何从BLASTX的输出,所以stdout和stderr是空的。 你可能想要检查下输出文件 opuntia.xml 是否已经创建。

如果你回想下这个指南的中的早先的例子,opuntia.fasta 包含7条序列, 所以BLAST XML 格式的结果输出文件应该包括多个结果。因此,我们在 下面的 7.3 节将用 Bio.Blast.NCBIXML.parse() 来 解析这个结果文件。

7.2.4 WU-BLAST 和 AB-BLAST

你也许会碰到 Washington University BLAST (WU-BLAST), 和它的后继版本`Advanced Biocomputing BLAST <http://blast.advbiocomp.com>`__ (AB-BLAST, 在2009年发布,免费但是没有开源)。这些程序包包括了命令工具行 wu-blastallab-blastall

Biopython 目前还没有提供调用这些工具的包装程序,但是应该可以解析它们 与NCBI兼容的输出结果。

7.3 解析BLAST 输出

就像上面提过的那样,BLAST能生成多种格式的输出,比如 XML, HTML 和纯文本格式。 以前,Biopython有针对HTML 和纯文本格式输出文件的解析器,因为当时只有这两种 格式的输出结果文件。不幸的是,这两种方式的BLAST 输出结果一直在变动,而每次 变动就会导致解析器失效。 所以,我们删除了针对HTML格式的解析器,不过纯文本格式 的解析还可以用(见 7.5 )。使用这个解析器 有一定的风险,它可能能工作也可能无效,依赖于你正在使用哪个BLAST版本。

跟上BLAST输出文件格式的改变很难,特别是当用户使用不同版本的BLAST的时候。 我们推荐使用XML格式的解析器。因为最近版本的BLAST能生成这种格式的输出结果。 XML格式的输出不仅比HTML 和纯文本格式的更稳定,而且解析起来更加容易自动化, 从而提高整个Biopython整体的稳定性。

你可以通过好几个途径来获得XML格式的BLAST输出文件。对解析器来说,不管你是 怎么生成输出的,只要是输出的格式是XML就行。

  • 你可以通过Biopython来运行因特网上的BLAST,就像 7.1 节描述的那样。
  • 你可以通过Biopython来运行本地的BLAST,就像 7.2 节描述的那样。
  • 你可以在通过浏览器在NCBI网站上进行BLAST搜索,然后保存结果文件。你需要选择输出 结果文件是XML格式的,并保存最终的结果网页(你知道,就是包含所有有趣结果的那个网页) 到文件。
  • 你也可以直接运行本地电脑上的BlAST,不通过Biopython,保存输出结果到文件。 同样的你也需要选择输出文件格式为XML。

关键点就是你不必用Biopython脚本来获取数据才能解析它。通过以上任何一种方式 获取了结果输出,你然后需要获得文件句柄来处理它。在Python中,一个文件句柄就是一种 用于描述到任何信息源的输入的良好通用的方式,以便于这些信息能够使用 read()readline() 函数(见章节 22.1 )来获取。

如果你一直跟着上几节用来和BLAST交互的代码的话,你已经有了个 result_handle ,一个用来得到BLAST的结果文件句柄。 比如通过GI号来进行一个在线BLAST搜索:

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

如果你通过其他方式运行了BLAST,并且XML格式的BLAST结果输出文件是 my_blast.xml , 那么你只需要打开文件来读:

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

好的,现在我们已经有了个文件句柄,可以解析输出结果了。解析结果的代码 很短。如果你想要一条BLAST输出结果(就是说,你只用了一条序列去搜索):

>>> from Bio.Blast import NCBIXML
>>> blast_record = NCBIXML.read(result_handle)

或者, 你有许多搜索结果(就是说,你用了多条序列去BLAST搜索)

>>> from Bio.Blast import NCBIXML
>>> blast_records = NCBIXML.parse(result_handle)

就像 Bio.SeqIOBio.AlignIO (参见第 5 和第 6 章), 我们有一对输入函数, readparse 。 当你只有一个输出结果的时候用 read 。当你有许多 输出结果的时候,可以用 parse 这个迭代器。 但是,我们调用函数获得结果 不是 SeqRecord 或者 MultipleSeqAlignment 对象,我们得到BLAST记录对象。

为了能处理BLAST结果文件很大有很多结果这种情况, NCBIXML.parse() 返回一个迭代器。简单来说,一个迭代器可以让你一个接着一个地获得BLAST 的搜索结果。

>>> from Bio.Blast import NCBIXML
>>> blast_records = NCBIXML.parse(result_handle)
>>> blast_record = blast_records.next()
# ... do something with blast_record
>>> blast_record = blast_records.next()
# ... do something with blast_record
>>> blast_record = blast_records.next()
# ... do something with blast_record
>>> blast_record = blast_records.next()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
StopIteration
# No further records

或者,你也可以使用 for - 循环

>>> for blast_record in blast_records:
...     # Do something with blast_record

注意对每个BLAST搜索结果只能迭代一次。通常,对于每个BLAST记录,你可能会保存你 感兴趣的信息。如果你想保存所有返回的BLAST记录,你可以把迭代 转换成列表。

>>> blast_records = list(blast_records)

现在,你可以像通常的做法通过索引从这个列表中获得每一条BLAST结果。 如果你的BLAST输出 结果文件很大,那么当把它们全部放入一个列表时,你也许会遇到内存不够的情况。

一般来说,你会一次运行一个BLAST搜索。然后,你只需提取第一条BLAST 搜索记录到 blast_records :

>>> from Bio.Blast import NCBIXML
>>> blast_records = NCBIXML.parse(result_handle)
>>> blast_record = blast_records.next()

or more elegantly:

或者更加优雅地:

>>> from Bio.Blast import NCBIXML
>>> blast_record = NCBIXML.read(result_handle)

我猜你现在在想BLAST搜索记录中到底有什么。

7.4 BLAST 记录类

一个BLAST搜索结果记录包括了所有你想要从中提取出来的信息。现在,我们将 用一个例子说明你怎么从BLAST搜索结果提取出一些信息。但是,如果你想从BLAST 搜索结果获得的信息没有在这里提到,你可以详细阅读BLAST搜索记录类, 并且可以参考下源代码 或者 是自动生成的文档 - 文档字符串里面包含了许多 关于各部分源代码是什么的很有用的信息。

继续我们的例子,让我们打印出所有大于某一特定阈值的BLAST命中结果的一些汇总信息。 代码如下:

>>> E_VALUE_THRESH = 0.04

>>> for alignment in blast_record.alignments:
...     for hsp in alignment.hsps:
...         if hsp.expect < E_VALUE_THRESH:
...             print '****Alignment****'
...             print 'sequence:', alignment.title
...             print 'length:', alignment.length
...             print 'e value:', hsp.expect
...             print hsp.query[0:75] + '...'
...             print hsp.match[0:75] + '...'
...             print hsp.sbjct[0:75] + '...'

上面代码会打印出如下图的总结报告:

****Alignment****
sequence: >gb|AF283004.1|AF283004 Arabidopsis thaliana cold acclimation protein WCOR413-like protein
alpha form mRNA, complete cds
length: 783
e value: 0.034
tacttgttgatattggatcgaacaaactggagaaccaacatgctcacgtcacttttagtcccttacatattcctc...
||||||||| | ||||||||||| || ||||  || || |||||||| |||||| |  | |||||||| ||| ||...
tacttgttggtgttggatcgaaccaattggaagacgaatatgctcacatcacttctcattccttacatcttcttc...

基本上,一旦你解析了BLAST搜索结果文件,你可以提取任何你需要的信息。 当然,这取决于你想要获得什么信息。但是希望这里的例子能够帮助你开始工作。

在用Biopython提取BLAST搜索结果信息的时候,重要的是你需要考虑到信息存储在什么 (Biopython)对象中。在Biopython中,解析器返回 Record 对象,这个对象 可以是 Blast 类型的,也可以是 PSIBlast 类型的,具体哪个取决你 解析什么。这些对象的定义都可以在 Bio.Blast.Record 找到 并且很完整。

下面是 我尝试画的 BlastPSIBlast 记录类的UML图。如果你对UML图很熟悉,不妨 看看下面的UML图是否有错误或者可以改进的地方,如果有,请联系我。 BLAST类图在这里 7.4

image1

PSIBlast 记录类是类似的,但是支持用在迭代器中的rounds方法。PSIBlast类图在这里 7.4

image2

7.5 废弃的BLAST 解析器

老版本的Biopython 有针对纯文本和HTML格式输出结果的解析器。但是经过几年 我们发现维护这些解析器很困难。基本上,任何BLAST输出的任何小改变都会导致 这些解析器失效。所以我们推荐你解析XML格式的BLAST输出结果,就像在 7.3 描述的那样。

取决于你使用Biopython的版本,纯文本格式的解析器也许有效也许失效。 用这个解析器的所带来的风险由你自己承担。

7.5.1 解析纯文本格式的BLAST输出

纯文本格式的解析器在 Bio.Blast.NCBIStandalone

和xml解析器类似, 我们也需要一个能够传给解析器的文件句柄。这个文件句柄必须 实现了 readline() 方法 。通常要获得这样文件句柄,既可以用Biopython提供的 blastallblastpgp 函数来调用本地的BLAST,或者从命令行运行本地的 BLAST, 并且如下处理:

>>> result_handle = open("my_file_of_blast_output.txt")

好了,既然现在得到了个文件句柄(我们称它是 result_handle ), 我们已经做好了解析它的准备。按下面的代码来解析:

>>> from Bio.Blast import NCBIStandalone
>>> blast_parser = NCBIStandalone.BlastParser()
>>> blast_record = blast_parser.parse(result_handle)

这样就能把BALST的搜索结果报告解析到Blast记录类中(取决你于你解析的对象, 解析结果可能返回一条 Blast 或者 PSIBlast记录)。这样你就可以从中提取 信息了。在我们的例子里,我们来打印出大于某个阈值的所有比对的一个总结 信息。

>>> E_VALUE_THRESH = 0.04
>>> for alignment in blast_record.alignments:
...     for hsp in alignment.hsps:
...         if hsp.expect < E_VALUE_THRESH:
...             print '****Alignment****'
...             print 'sequence:', alignment.title
...             print 'length:', alignment.length
...             print 'e value:', hsp.expect
...             print hsp.query[0:75] + '...'
...             print hsp.match[0:75] + '...'
...             print hsp.sbjct[0:75] + '...'

如果你已经读过 7.3 节关于解析XML格式的部分, 你将会发现上面的代码和那个章节的是一样的。一旦你把输出文件解析到记录类中, 你就能处理信息,不管你原来的BLAST输出格式是什么。很赞吧。

好,解析一条记录是不错,那么如果我有一个包含许多记录的BLAST文件 - 我该怎么处理它们呢?好吧,不要害怕,答案就在下个章节中。

7.5.2 解析包含多次BLAST结果的纯文本BLAST文件

我们可以用BLAST迭代器解析多次结果。为了得到一个迭代器,我们首先需要创建一个解析器,来 解析BLAST的搜索结果报告为Blast记录对象。

>>> from Bio.Blast import NCBIStandalone
>>> blast_parser = NCBIStandalone.BlastParser()

然后,我们假定我们有一个连接到一大堆blast记录的文件句柄,我们把这个文件句柄 叫做 result_handle 。 怎么得到一个文件句柄在上面blast解析章节有详细 描述。

好了,我们现在有了一个解析器和一个文件句柄,我们可以用以下命令来创建 一个迭代器。

>>> blast_iterator = NCBIStandalone.Iterator(result_handle, blast_parser)

第二个参数,解析器,是可选的。如果我们没有提供一个解析器,那么迭代器将会 一次返回一个原始的BLAST搜索结果。

现在我们已经有了个迭代器,就可以开始通过 next() 方法来获取BLAST 记录(由我们的解析器产生)。

>>> blast_record = blast_iterator.next()

每次调用next都会返回一条我们能处理的新记录。现在我们可以遍历所有记录,并打印一 个我们最爱、漂亮的、简洁的BLAST记录报告。

>>> for blast_record in blast_iterator:
...     E_VALUE_THRESH = 0.04
...     for alignment in blast_record.alignments:
...         for hsp in alignment.hsps:
...             if hsp.expect < E_VALUE_THRESH:
...                 print '****Alignment****'
...                 print 'sequence:', alignment.title
...                 print 'length:', alignment.length
...                 print 'e value:', hsp.expect
...                 if len(hsp.query) > 75:
...                     dots = '...'
...                 else:
...                     dots = ''
...                 print hsp.query[0:75] + dots
...                 print hsp.match[0:75] + dots
...                 print hsp.sbjct[0:75] + dots

迭代器允许你处理很多blast记录而不出现内存不足的问题。因为,它使一次处理 一个记录。我曾经用大处理过一个非常巨大的文件,没有出过任何问题。

7.5.3 在巨大的BLAST纯文本文件中发现不对的记录

当我开始解析一个巨大的blast 文件,有时候会碰到一个郁闷的问题就是解析器以一个 ValueError异常终止了。这是个严肃的问题。因为你无法分辨导致ValueError异常的是 解析器的问题还是BLAST的问题。更加糟糕是,你不知道在哪一行解析器失效了。所以, 你不能忽略这个错误。不然,可能会忽视一个重要的数据。

我们以前必须写一些小脚本来解决这个问题。不过,现在 Bio.Blast 模块包含了 BlastErrorParser ,可以更加简单地来解决这个问题。 BlastErrorParser 和常规的 BlastParser 类似,但是它加了特别一层来捕获由解析器产生的ValueErrors 异常,并尝试来诊断这些错误。

让我们来看看怎样用这个解析器 - 首先我们定义我们准备解析的文件和报告错误情况的 输出文件。

>>> import os
>>> blast_file = os.path.join(os.getcwd(), "blast_out", "big_blast.out")
>>> error_file = os.path.join(os.getcwd(), "blast_out", "big_blast.problems")

现在我们想要一个 BlastErrorParser

>>> from Bio.Blast import NCBIStandalone
>>> error_handle = open(error_file, "w")
>>> blast_error_parser = NCBIStandalone.BlastErrorParser(error_handle)

注意,解析器有个关于文件句柄的可选参数。如果传递了这个参数,那么解析器就会 把产生ValueError异常的记录写到这个文件句柄中。不然的话,这些错误记录就不会 被记录下来。

现在,我们可以像用常规的blast解析器一样地用 BlastErrorParser 。 特别的是,我们也许想要一个一次读入一个记录的迭代器并用 BlastErrorParser 来解析它。

>>> result_handle = open(blast_file)
>>> iterator = NCBIStandalone.Iterator(result_handle, blast_error_parser)

我们可以一次读一个记录,并且我们现在可以捕获并处理那些因为Blast引起的、 不是解析器本身导致的错误。

>>> try:
...     next_record = iterator.next()
... except NCBIStandalone.LowQualityBlastError, info:
...     print "LowQualityBlastError detected in id %s" % info[1]

.next() 方法通常被 for 循环间接地调用。现在, BlastErrorParser 能够捕获如下的错误:

  • ValueError - 这就是和常规BlastParser产生的一样的错误。这个错误产生 是因为解析器不能解析某个文件。通常是因为解析器有bug, 或者是 因为你使用解析器的版本和你BLAST命令的版本不一致。
  • LowQualityBlastError - 当Blast一条低质量的序列时(比如,一条 只有1个核苷酸的短序列),似乎Blast会终止并屏蔽掉整个序列,所有就没有什么可以 解析了。 这种情况下,Blast就会产生一个不完整的报告导致解析器出现ValueError 错误。 LowQualityBlastError 错误在这种情况下产生。这个错误返回如下 信息:
    • item[0] – The error message
    • item[0] - 错误消息
    • item[1] – The id of the input record that caused the error. This is really useful if you want to record all of the records that are causing problems.
    • item[1] - 导致错误产生的输入记录id。如果你想记录所有导致问题 记录的时候很有用。

就像上面提到的那样,BlastErrorParser 将会把有问题的记录写到指定的``error_handle``。 然后,你可以排查这些有问题记录。你可以针对某条记录来调试解析器,或者找到 你运行blast中的问题。无论哪种方式,这些都是有用的经验。

希望 BlastErrorParser 能帮你更简单的调试和处理一些数据巨大的Blast 文件。

7.6 处理PSI-BLAST

你可以通过 Bio.Blast.Applications 模块中的包装函数来运行单机版本的PSI-BLAST (老版本的NCBI命令工具行 blastpgp 或者它的替代程序 psiblast )。

在写这篇指南的时候,没有迹象表明NCBI将会支持通过internet来进行PSI-BLAST 搜索。

请注意 Bio.Blast.NCBIXML 解析器能读入并解析当前版本PSI-BLAST的、XML格式的 输出,但是像哪条序列在每个迭代循环中是新的还是复用的信息在XML格式输出中是没有的。 如果,你需要这些信息你应该用纯文本输出和 Bio.Blast.NCBIStandalone 模块的 PSIBlastParser

7.7 处理 RPS-BLAST

你可以通过 Bio.Blast.Applications 模块中的包装函数来运行单机版本的RPS-BLAST (或者老版本的NCBI命令工具行 rpsblast 或者同样名字的替代程序 )。

在写这篇指南的时候,没有迹象表明NCBI将会支持通过internet来进行RPS-BLAST 搜索

你可以通过 Bio.Blast.NCBIXML 这个解析器来读入并解析当前版本的RPS-BLAST的 XML格式的输出。