第8章 BLAST和其他序列搜索工具(实验性质的代码)¶

WARNING: 这章教程介绍了Biopython中一个 实验的 模块。它正在被加入到 Biopython中，并且以一个预尾期的状态整理到教程当中，这样在我们发布稳定版的之前可以收到一系列的反馈和并作改进。那时有些细节可能会改变，并且用到当前 Bio.SearchIO 模块的脚本也需要更新。切记！为了与NCBI BLAST有关的代码可以稳定工作，请继续使用第 7 章介绍的Bio.Blast。

生物序列的鉴定是生物信息工作的主要部分。有几个工具像BLAST（可能是最流行的），FASTA ，HMMER还有许多其它的都有这个功能，每个工具都有独特的算法和途径。一般来说，这些工具都是用你的序列在数据库中搜索可能的匹配。随着序列数量的增加（匹配数也会随之增加），将会有成百上千的可能匹配，解析这些结果无疑变得越来越困难。理所当然，人工解析搜索结果变得不可能。而且你可能会同时用几种不同的搜索工具，每种工具都有独特的统计方法、规则和输出格式。可以想象，同时用多种工具搜索多条序列是多么恐怖的事。

我们对此非常了解，所以我们在Biopython构建了 Bio.SearchIO 亚模块。Bio.SearchIO 模块使从搜索结果中提取信息变得简单，并且可以处理不同工具的不同标准和规则。SearchIO 和BioPerl中模块名称一致。

在本章中，我们将学习 Bio.SearchIO 的主要功能，了解它可以做什么。我们将使用两个主要的搜索工具：BLAST和FASTA。它们只是用来阐明思路，你可以轻易地把工作流程应用到 Bio.SearchIO 支持的其他工具中。欢迎你使用我们将要用到的搜索结果文件。BLAST搜索结果文件可以在这里下载。 BLAT输出结果文件可以在这里下载。两个结果文件都是用下面这条序列搜索产生的：

>mystery_seq
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG

BLAST的XML结果是用 blastn 搜索NCBI的 refseq_rna 数据库得到的。对于 BLAT，数据库是2009年2月份的 hg19 人类基因组草图，输出格式是PSL。

我们从 Bio.SearchIO 的对象模型的介绍开始。这个模型代表你的搜索结果，因此它是 Bio.SearchIO 的核心。然后，我们会介绍 Bio.SearchIO 常用的主要功能。

现在一切就绪，让我们开始第一步：介绍核心对象模型。

8.1 SearchIO对象模型¶

尽管多数搜索工具的输出风格极为不同，但是它们蕴含的理念很相似：

输出文件可能包含一条或更多的搜索查询的结果。
在每次查询中，你会在给定的数据库中得到一个或更多的hit（命中）。
在每个hit中，你会得到一个或更多包含query（查询)序列和数据库序列实际比对的区域。
一些工具如BLAT和Exonerate可能会把这些区域分成几个比对片段（或在BLAT中称为区块，在Exonerate中称为可能外显子）。这并不是很常见，像BLAST和 HMMER就不这么做。

知道这些共性之后，我们决定把它们作为创造 Bio.SearchIO 对象模型的基础。对象模型是Python对象组成的嵌套分级系统，每个对象都代表一个上面列出的概念。这些对象是：

QueryResult，代表单个搜索查询。
Hit，代表单个的数据库hit。Hit 对象包含在 QueryResult 中，每个 QueryResult 中有0个或多个 Hit 对象。
HSP (high-scoring pair（高分片段）)，代表query和hit序列中有意义比对的区域。HSP 对象包含在 Hit 对象中，而且每个 Hit 有一个或更多的 HSP 对象。
HSPFragment，代表query和hit序列中单个的邻近比对。 HSPFragment 对象包含在 HSP 对象中。多数的搜索工具如BLAST和HMMER把 HSP 和 HSPFragment 合并，因为一个 HSP 只含有一个 HSPFragment。但是像BLAT和Exonerate会产生含有多个 HSPFragment 的 HSP 。似乎有些困惑？不要紧，稍后我们将详细介绍这两个对象。

这四个对象是当你用 Bio.SearchIO 会碰到的。 Bio.SearchIO 四个主要方法： read ，parse，index 或 index_db 中任意一个都可以产生这四个对象。这些方法的会在后面的部分详细介绍。这部分只会用到 read 和 parse ，这两个方法和 Bio.SeqIO 以及 Bio.AlignIO 中的 read 和 parse 方法功能相似：

read 用于单query对输出文件进行搜索并且返回一个 QueryResult 对象。
parse 用于多query对输出文件进行搜索并且返回一个可以yield QueryResult 对象的generator。

完成这些之后，让我们开始学习每个 Bio.SearchIO 对象，从 QueryResult 开始。

8.1.1 QueryResult¶

QueryResult，代表单query搜索，每个 QueryResult 中有0个或多个 Hit 对象。我们来看看BLAST文件是什么样的：

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> print blast_qresult
Program: blastn (2.2.27+)
  Query: 42291 (61)
         mystery_seq
 Target: refseq_rna
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description
         ----  -----  ----------------------------------------------------------
            0      1  gi|262205317|ref|NR_030195.1|  Homo sapiens microRNA 52...
            1      1  gi|301171311|ref|NR_035856.1|  Pan troglodytes microRNA...
            2      1  gi|270133242|ref|NR_032573.1|  Macaca mulatta microRNA ...
            3      2  gi|301171322|ref|NR_035857.1|  Pan troglodytes microRNA...
            4      1  gi|301171267|ref|NR_035851.1|  Pan troglodytes microRNA...
            5      2  gi|262205330|ref|NR_030198.1|  Homo sapiens microRNA 52...
            6      1  gi|262205302|ref|NR_030191.1|  Homo sapiens microRNA 51...
            7      1  gi|301171259|ref|NR_035850.1|  Pan troglodytes microRNA...
            8      1  gi|262205451|ref|NR_030222.1|  Homo sapiens microRNA 51...
            9      2  gi|301171447|ref|NR_035871.1|  Pan troglodytes microRNA...
           10      1  gi|301171276|ref|NR_035852.1|  Pan troglodytes microRNA...
           11      1  gi|262205290|ref|NR_030188.1|  Homo sapiens microRNA 51...
           12      1  gi|301171354|ref|NR_035860.1|  Pan troglodytes microRNA...
           13      1  gi|262205281|ref|NR_030186.1|  Homo sapiens microRNA 52...
           14      2  gi|262205298|ref|NR_030190.1|  Homo sapiens microRNA 52...
           15      1  gi|301171394|ref|NR_035865.1|  Pan troglodytes microRNA...
           16      1  gi|262205429|ref|NR_030218.1|  Homo sapiens microRNA 51...
           17      1  gi|262205423|ref|NR_030217.1|  Homo sapiens microRNA 52...
           18      1  gi|301171401|ref|NR_035866.1|  Pan troglodytes microRNA...
           19      1  gi|270133247|ref|NR_032574.1|  Macaca mulatta microRNA ...
           20      1  gi|262205309|ref|NR_030193.1|  Homo sapiens microRNA 52...
           21      2  gi|270132717|ref|NR_032716.1|  Macaca mulatta microRNA ...
           22      2  gi|301171437|ref|NR_035870.1|  Pan troglodytes microRNA...
           23      2  gi|270133306|ref|NR_032587.1|  Macaca mulatta microRNA ...
           24      2  gi|301171428|ref|NR_035869.1|  Pan troglodytes microRNA...
           25      1  gi|301171211|ref|NR_035845.1|  Pan troglodytes microRNA...
           26      2  gi|301171153|ref|NR_035838.1|  Pan troglodytes microRNA...
           27      2  gi|301171146|ref|NR_035837.1|  Pan troglodytes microRNA...
           28      2  gi|270133254|ref|NR_032575.1|  Macaca mulatta microRNA ...
           29      2  gi|262205445|ref|NR_030221.1|  Homo sapiens microRNA 51...
           ~~~
           97      1  gi|356517317|ref|XM_003527287.1|  PREDICTED: Glycine ma...
           98      1  gi|297814701|ref|XM_002875188.1|  Arabidopsis lyrata su...
           99      1  gi|397513516|ref|XM_003827011.1|  PREDICTED: Pan panisc...

虽然我们才接触对象模型的皮毛，但是你已经可以看到一些有用的信息了。通过调用``QueryResult`` 对象的 print 方法，你可以看到：

程序的名称和版本 (blastn version 2.2.27+)
query的ID，描述和序列长度(ID是42291，描述是 ‘mystery_seq’，长度是61)
搜索的目标数据库 (refseq_rna)
hit结果的快速预览。对于我们的查询序列，有100个可能的hit（表格中表示的是 0-99）对于每个hit，我们可以看到它包含的高分比对片段（HSP)，ID和一个片段描述。注意， Bio.SearchIO 截断了表格，只显示0-29，然后是97-99。

现在让我们用同样的步骤来检查BLAT的结果：

>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> print blat_qresult
Program: blat (<unknown version>)
  Query: mystery_seq (61)
         <unknown description>
 Target: <unknown target>
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description
         ----  -----  ----------------------------------------------------------
            0     17  chr19  <unknown description>

马上可以看到一些不同点。有些是由于BLAT使用PSL格式储存它的信息，稍后会看到。其余是由于BLAST和BLAT搜索的程序和数据库之间明显的差异造成的：

程序名称和版本。 Bio.SearchIO 知道程序是BLAST，但是在输出文件中没有信息显示程序版本，所以默认是 ‘<unknown version>’。
query的ID，描述和序列的长度。注意，这些细节和BLAST的细节只有细小的差别， ID是 ‘mystery_seq’ 而不是42991，缺少描述，但是序列长度仍是61。这实际上是文件格式本身导致的差异。BLAST有时创建自己的query ID并且用你的原始ID作为序列描述。
目标数据库是未知的，因为BLAT输出文件没提到相关信息。
最后，hit列表完全不同。这里，我们的查询序列只命中到 ‘chr19’ 数据库条目，但是我们可以看到它含有17个HSP区域。这并不让人诧异，考虑到我们使用的是不同的程序，并且这些程序都有自己的数据库。

所有通过调用 print 方法看到的信息都可以单独地用Python的对象属性获得（又叫点标记法）。同样还可以用相同的方法获得其他格式特有的属性。

>>> print "%s %s" % (blast_qresult.program, blast_qresult.version)
blastn 2.2.27+
>>> print "%s %s" % (blat_qresult.program, blat_qresult.version)
blat <unknown version>
>>> blast_qresult.param_evalue_threshold    # blast-xml specific
10.0

想获得一个可访问属性的完整列表，可以查询每个格式特有的文档。这些是 BLAST BLAT.

已经知道了在 QueryResult 对象上可以调用 print 方法，让我们研究的更深一些。 QueryResult 到底是什么？就Python对象来说， QueryResult 混合了列表和字典的特性。换句话说，也就是一个包含了列表和字典功能的容器对象。

和列表以及字典一样， QueryResult 对象是可迭代的。每次迭代返回一个hit 对象：

>>> for hit in blast_qresult:
...     hit
Hit(id='gi|262205317|ref|NR_030195.1|', query_id='42291', 1 hsps)
Hit(id='gi|301171311|ref|NR_035856.1|', query_id='42291', 1 hsps)
Hit(id='gi|270133242|ref|NR_032573.1|', query_id='42291', 1 hsps)
Hit(id='gi|301171322|ref|NR_035857.1|', query_id='42291', 2 hsps)
Hit(id='gi|301171267|ref|NR_035851.1|', query_id='42291', 1 hsps)
...

要得到 QueryResult 对象有多少hit，可以简单调用Python的 len 方法：

>>> len(blast_qresult)
100
>>> len(blat_qresult)
1

同列表类似，你可以用切片来获得 QueryResult 对象的hit：

>>> blast_qresult[0]        # retrieves the top hit
Hit(id='gi|262205317|ref|NR_030195.1|', query_id='42291', 1 hsps)
>>> blast_qresult[-1]       # retrieves the last hit
Hit(id='gi|397513516|ref|XM_003827011.1|', query_id='42291', 1 hsps)

要得到多个hit，你同样可以对 QueryResult 对象作切片。这种情况下，返回一个包含被切hit的新 QueryResult 对象：

>>> blast_slice = blast_qresult[:3]     # slices the first three hits
>>> print blast_slice
Program: blastn (2.2.27+)
  Query: 42291 (61)
         mystery_seq
 Target: refseq_rna
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description
         ----  -----  ----------------------------------------------------------
            0      1  gi|262205317|ref|NR_030195.1|  Homo sapiens microRNA 52...
            1      1  gi|301171311|ref|NR_035856.1|  Pan troglodytes microRNA...
            2      1  gi|270133242|ref|NR_032573.1|  Macaca mulatta microRNA ...

同字典类似，可以通过ID获取hit。如果你知道一个特定的hit ID存在于一个搜索结果中时，特别有用：

>>> blast_qresult['gi|262205317|ref|NR_030195.1|']
Hit(id='gi|262205317|ref|NR_030195.1|', query_id='42291', 1 hsps)

你可以用 hits 方法获得完整的 Hit 对象，也可以用 hit_keys 方法获得完整的 Hit ID：

>>> blast_qresult.hits
[...]       # list of all hits
>>> blast_qresult.hit_keys
[...]       # list of all hit IDs

如果你想确定一个特定的hit是否存在于查询结果中该怎么做呢？可以用 in 关键字作一个简单的成员检验：

>>> 'gi|262205317|ref|NR_030195.1|' in blast_qresult
True
>>> 'gi|262205317|ref|NR_030194.1|' in blast_qresult
False

有时候，只知道一个hit是否存在是不够的，你可能也会想知道hit的排名。 index 方法可以帮助你：

>>> blast_qresult.index('gi|301171437|ref|NR_035870.1|')
22

记住，我们用的是Python风格的索引，是从0开始。这代表hit的排名是23而不是22。

同样，注意你看的hit排名是基于原始搜索输出文件的本来顺序。不同的搜索工具可能会基于不同的标准排列hit。

如果原本的hit排序不合你意，可以用 QueryResult 对象的 sort 方法。它和Python的 list.sort 方法很相似，只是有个是否创建一个新的排序后的 QueryResult 对象的选项。

这里有个用 QueryResult.sort 方法对hit排序的例子，这个方法基于每个hit 的完整序列长度。对于这个特殊的排序，我们设置 in_place 参数等于 False ，这样排序方法会返回一个新的 QueryResult 对象，而原来的对象是未排序的。我们同样可以设置 reverse 参数等于 `` True `` 以递减排序。

>>> for hit in blast_qresult[:5]:   # id and sequence length of the first five hits
...     print hit.id, hit.seq_len
...
gi|262205317|ref|NR_030195.1| 61
gi|301171311|ref|NR_035856.1| 60
gi|270133242|ref|NR_032573.1| 85
gi|301171322|ref|NR_035857.1| 86
gi|301171267|ref|NR_035851.1| 80

>>> sort_key = lambda hit: hit.seq_len
>>> sorted_qresult = blast_qresult.sort(key=sort_key, reverse=True, in_place=False)
>>> for hit in sorted_qresult[:5]:
...     print hit.id, hit.seq_len
...
gi|397513516|ref|XM_003827011.1| 6002
gi|390332045|ref|XM_776818.2| 4082
gi|390332043|ref|XM_003723358.1| 4079
gi|356517317|ref|XM_003527287.1| 3251
gi|356543101|ref|XM_003539954.1| 2936

有 in_place 参数的好处是可以保留原本的顺序，后面会用到。注意这不是 QueryResult.sort 的默认行为，需要我们明确地设置 in_place 为 True 。

现在，你已经知道使用 QueryResult 对象。但是，在我们学习 Bio.SearchIO 模块下个对象前，先了解下可以让 QueryResult 对象更易使用的两个方法： filter 和 map 方法。

如果你对Python的列表推导式、generator表达式或内建的 filter 和 map 很熟悉，就知道（不知道就是看看吧!)它们在处理list-like的对象时有多有用。你可以用这些内建的方法来操作 QueryResult 对象，但是这只对正常list有效，并且可操作性也会受到限制。

这就是为什么 QueryResult 对象提供自己特有的 filter 和 map 方法。对于 filter 有相似的 hit_filter 和 hsp_filter 方法，从名称就可以看出，这些方法过滤 QueryResult 对象的 Hit 对象或者 HSP 对象。同样的，对于 map ， QueryResult 对象同样提供相似的 hit_map 和 hsp_map 方法。这些方法分别应用于 QueryResult 对象 hit 或者 HSP 对象。

让我们来看看这些方法的功能，从 hit_filter 开始。这个方法接受一个回调函数，这个函数检验给定的 Hit 是否符合你设定的条件。换句话说，这个方法必须接受一个单独 Hit 对象作为参数并且返回 True 或 False 。

这里有个用 hit_filter 筛选出只有一个HSP的 Hit 对象的例子：

>>> filter_func = lambda hit: len(hit.hsps) > 1     # the callback function
>>> len(blast_qresult)      # no. of hits before filtering
100
>>> filtered_qresult = blast_qresult.hit_filter(filter_func)
>>> len(filtered_qresult)   # no. of hits after filtering
37
>>> for hit in filtered_qresult[:5]:    # quick check for the hit lengths
...     print hit.id, len(hit.hsps)
gi|301171322|ref|NR_035857.1| 2
gi|262205330|ref|NR_030198.1| 2
gi|301171447|ref|NR_035871.1| 2
gi|262205298|ref|NR_030190.1| 2
gi|270132717|ref|NR_032716.1| 2

hsp_filter 和 hit_filter 功能相同，只是它过滤每个hit中的 HSP 对象，而不是 Hit 。

对于 map 方法，同样接受一个回调函数作为参数。但是回调函数返回修改过的 Hit 或 HSP 对象（取决于你是否使用 hit_map 或 hsp_map 方法），而不是返回 True 或 False。

来看一个用 hit_map 方法来重命名hit ID的例子：

>>> def map_func(hit):
...     hit.id = hit.id.split('|')[3]   # renames 'gi|301171322|ref|NR_035857.1|' to 'NR_035857.1'
...     return hit
...
>>> mapped_qresult = blast_qresult.hit_map(map_func)
>>> for hit in mapped_qresult[:5]:
...     print hit.id
NR_030195.1
NR_035856.1
NR_032573.1
NR_035857.1
NR_035851.1

同样的， hsp_map 和 hit_map 作用相似, 但是作用于 HSP 对象而不是 Hit 对象。

8.1.2 Hit¶

Hit 对象代表从单个数据库获得所有查询结果。在 Bio.SearchIO 对象等级中是二级容器。它们被包含在 QueryResult 对象中，同时它们又包含 HSP 对象。

看看它们是什么样的，从我们的BLAST搜索开始：

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> blast_hit = blast_qresult[3]    # fourth hit from the query result

>>> print blast_hit
Query: 42291
       mystery_seq
  Hit: gi|301171322|ref|NR_035857.1| (86)
       Pan troglodytes microRNA mir-520c (MIR520C), microRNA
 HSPs: ----  --------  ---------  ------  ---------------  ---------------------
          #   E-value  Bit score    Span      Query range              Hit range
       ----  --------  ---------  ------  ---------------  ---------------------
          0   8.9e-20     100.47      60           [1:61]                [13:73]
          1   3.3e-06      55.39      60           [0:60]                [13:73]

可以看到我们获得了必要的信息：

query的ID和描述信息。一个hit总是和一个query绑定，所以我们同样希望记录原始 query。这些值可以通过 query_id 和 query_description 属性从hit 中获取。
我们同样得到了hit的ID、描述和序列全长。它们可以分别通过 id， description，和 seq_len 获取。
最后，有一个hit含有的HSP的简短信息表。在每行中，HSP重要信息被列出来：HSP索引，e值，得分，长度（包括gap），query坐标和hit坐标。

现在，和BLAT结果作对比。记住，在BLAT搜索结果中，我们发现有一个含有17HSP的hit。

>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> blat_hit = blat_qresult[0]      # the only hit
>>> print blat_hit
Query: mystery_seq
       <unknown description>
  Hit: chr19 (59128983)
       <unknown description>
 HSPs: ----  --------  ---------  ------  ---------------  ---------------------
          #   E-value  Bit score    Span      Query range              Hit range
       ----  --------  ---------  ------  ---------------  ---------------------
          0         ?          ?       ?           [0:61]    [54204480:54204541]
          1         ?          ?       ?           [0:61]    [54233104:54264463]
          2         ?          ?       ?           [0:61]    [54254477:54260071]
          3         ?          ?       ?           [1:61]    [54210720:54210780]
          4         ?          ?       ?           [0:60]    [54198476:54198536]
          5         ?          ?       ?           [0:61]    [54265610:54265671]
          6         ?          ?       ?           [0:61]    [54238143:54240175]
          7         ?          ?       ?           [0:60]    [54189735:54189795]
          8         ?          ?       ?           [0:61]    [54185425:54185486]
          9         ?          ?       ?           [0:60]    [54197657:54197717]
         10         ?          ?       ?           [0:61]    [54255662:54255723]
         11         ?          ?       ?           [0:61]    [54201651:54201712]
         12         ?          ?       ?           [8:60]    [54206009:54206061]
         13         ?          ?       ?          [10:61]    [54178987:54179038]
         14         ?          ?       ?           [8:61]    [54212018:54212071]
         15         ?          ?       ?           [8:51]    [54234278:54234321]
         16         ?          ?       ?           [8:61]    [54238143:54238196]

我们得到了和前面看到的BLAST hit详细程度相似的结果。但是有些不同点需要解释：

e-value和bit score列的值。因为BLAT HSP没有e-values和bit scores，默认显示‘?’.
span列是怎么回事呢？span值本来是显示完整的比对长度，包含所有的残基和 gap。但是PSL格式目前还不支持这些信息并且 Bio.SearchIO 也不打算去猜它到底是多少，所有我们得到了和e-value以及bit score列相同的 ‘?’。

就Python对象来说， Hit 和列表行为最相似，但是额外含有 HSP 。如果你对列表熟悉，在使用 Hit 对象是不会遇到困难。

和列表一样， Hit 对象是可迭代的，并且每次迭代返回一个 HSP 对象：

>>> for hsp in blast_hit:
...     hsp
HSP(hit_id='gi|301171322|ref|NR_035857.1|', query_id='42291', 1 fragments)
HSP(hit_id='gi|301171322|ref|NR_035857.1|', query_id='42291', 1 fragments)

你可以对 Hit 对象调用 len 方法查看它含有多少个 HSP 对象：

>>> len(blast_hit)
2
>>> len(blat_hit)
17

你可以对 Hit 对象作切片取得单个或多个 HSP 对象，和 QueryResult 一样，如果切取多个 HSP ，会返回包含被切片 HSP 的一个新 Hit 对象。

>>> blat_hit[0]                 # retrieve single items
HSP(hit_id='chr19', query_id='mystery_seq', 1 fragments)
>>> sliced_hit = blat_hit[4:9]  # retrieve multiple items
>>> len(sliced_hit)
5
>>> print sliced_hit
Query: mystery_seq
       <unknown description>
  Hit: chr19 (59128983)
       <unknown description>
 HSPs: ----  --------  ---------  ------  ---------------  ---------------------
          #   E-value  Bit score    Span      Query range              Hit range
       ----  --------  ---------  ------  ---------------  ---------------------
          0         ?          ?       ?           [0:60]    [54198476:54198536]
          1         ?          ?       ?           [0:61]    [54265610:54265671]
          2         ?          ?       ?           [0:61]    [54238143:54240175]
          3         ?          ?       ?           [0:60]    [54189735:54189795]
          4         ?          ?       ?           [0:61]    [54185425:54185486]

你同样可以对一个 Hit 里的 HSP 排序，和你在 QueryResult 对象中看到的方法一样。

最后，同样可以对 Hit 对象使用 filter 和 map 方法。和 QueryResult 不同， Hit 对象只有一种 filter (Hit.filter) 和一种 map (Hit.map)。

8.1.3 HSP¶

HSP (高分片段)代表hit序列中的一个区域，该区域包含对于查询序列有意义的比对。它包含了你的查询序列和一个数据库条目之间精确的匹配。由于匹配取决于序列搜索工具的算法， HSP 含有大部分统计信息，这些统计是由搜索工具计算得到的。这使得不同搜索工具的 HSP 对象之间的差异和你在 QueryResult 以及 Hit 对象看到的差异更加明显。

我们来看看BLAST和BLAT搜索的例子。先看BLAST HSP：

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> blast_hsp = blast_qresult[0][0]    # first hit, first hsp

>>> print blast_hsp
      Query: 42291 mystery_seq
        Hit: gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 520b (MIR520...
Query range: [0:61] (1)
  Hit range: [0:61] (1)
Quick stats: evalue 4.9e-23; bitscore 111.29
  Fragments: 1 (61 columns)
     Query - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
             |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
       Hit - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG

和 QueryResult 以及 Hit 类似，调用 HSP 的 print 方法, 显示细节：

有query和hit的ID以及描述。我们需要这些来识别我们的 HSP 。
我们同样得到了query和hit序列的匹配范围。这里用的的切片标志着范围的表示是使用Python的索引风格（从0开始，半开区间）。圆括号里的数字表示正负链。这里，两条序列都是正链。
还有一些简短统计：e-value和bitscore。
还有一些HSP片段的信息。现在可以忽略，稍后会解释。
最后，还有query和hit的比对。

这些信息可以用点标记从它们本身获得，和 Hit 以及 QueryResult 相同：

>>> blast_hsp.query_range
(0, 61)

>>> blast_hsp.evalue
4.91307e-23

它们并不是仅有的属性， HSP 对象有一系列的属性，使得获得它们的具体信息更加容易。下面是一些例子：

>>> blast_hsp.hit_start         # start coordinate of the hit sequence
0
>>> blast_hsp.query_span        # how many residues in the query sequence
61
>>> blast_hsp.aln_span          # how long the alignment is
61

查看 HSP 文档获取完整的属性列表。

不仅如此，每个搜索工具通常会对它的 HSP 对象作统计学或其他细节计算。例如，一个 XML BLAST搜索同样输出gap以及相同的残基数量。这些属性可以像这样被获取：

>>> blast_hsp.gap_num       # number of gaps
0
>>> blast_hsp.ident_num     # number of identical residues
61

这些细节是格式特异的；它们可能不会出现在其他的格式中。要知道哪些细节在给定的序列搜索工具中是存在的，你应该查看那种格式的在 Bio.SearchIO 中的文档。或者可以用 .__dict__.keys() 获得快速列表：

>>> blast_hsp.__dict__.keys()
['bitscore', 'evalue', 'ident_num', 'gap_num', 'bitscore_raw', 'pos_num', '_items']

最后，你可能已经注意到了，我们HSP的 query 和 hit 属性不只是规律字符串：

>>> blast_hsp.query
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG', DNAAlphabet()), id='42291', name='aligned query sequence', description='mystery_seq', dbxrefs=[])
>>> blast_hsp.hit
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG', DNAAlphabet()), id='gi|262205317|ref|NR_030195.1|', name='aligned hit sequence', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])

它们是你已经在第 4 章看到过的 SeqRecord 对象！意味着你可以对 SeqRecord 对象做的各种有趣的事同样适用于 HSP.query 和 HSP.hit 对象。

现在 HSP 对象有个 alignment 属性（一个 MultipleSeqAlignment 对象）应该不会让你感到惊讶：

>>> print blast_hsp.aln
DNAAlphabet() alignment with 2 rows and 61 columns
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAG...GGG 42291
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAG...GGG gi|262205317|ref|NR_030195.1|

探索完BLAST HSP对象，让我们看看来自BLAT结果的不一样的HSP。我们将对它调用 print 方法：

>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> blat_hsp = blat_qresult[0][0]       # first hit, first hsp
>>> print blat_hsp
      Query: mystery_seq <unknown description>
        Hit: chr19 <unknown description>
Query range: [0:61] (1)
  Hit range: [54204480:54204541] (1)
Quick stats: evalue ?; bitscore ?
  Fragments: 1 (? columns)

一些输出你应该已经猜到了。我们得到了查询序列、hit ID、描述以及序列坐标。 evalue和bitscore的值是 ‘?’ ，因为BLAT HSP并没有这些属性。但是最大的不同是你看不到任何的序列比对展示。如果你看的更仔细，PSL格式本身并没有任何的 hit和query序列，所以 Bio.SearchIO 不会创建任何序列或者比对对象。如果你尝试获取 HSP.query ，HSP.hit ，或者 HSP.aln 属性会怎么样呢？你会得到这些属性的默认值 None ：

>>> blat_hsp.hit is None
True
>>> blat_hsp.query is None
True
>>> blat_hsp.aln is None
True

这并不影响其他的属性。例如，你仍然可以获取query和hit比对的长度。尽管不显示任何的属性，但是PSL格式还是有这些信息的，所以 Bio.SearchIO 可以抽提出这些信息。

>>> blat_hsp.query_span     # length of query match
61
>>> blat_hsp.hit_span       # length of hit match
61

其他格式特异的属性同样被展示出来：

>>> blat_hsp.score          # PSL score
61
>>> blat_hsp.mismatch_num   # the mismatch column
0

到目前为止，一切还不错？当你看到BLAT结果中不同的HSP时，事情变得更有趣了。你可能会回想起在BLAT搜索中，有时我们把结果分成 ‘blocks’ 。这些区块是必需比对片段，可能会有些内含子在它们之间。

让我们看看 Bio.SearchIO 怎么处理包含多个区块的BLAT HSP：

>>> blat_hsp2 = blat_qresult[0][1]      # first hit, second hsp
>>> print blat_hsp2
      Query: mystery_seq <unknown description>
        Hit: chr19 <unknown description>
Query range: [0:61] (1)
  Hit range: [54233104:54264463] (1)
Quick stats: evalue ?; bitscore ?
  Fragments: ---  --------------  ----------------------  ----------------------
               #            Span             Query range               Hit range
             ---  --------------  ----------------------  ----------------------
               0               ?                  [0:18]     [54233104:54233122]
               1               ?                 [18:61]     [54264420:54264463]

怎么回事？我们仍然得到了一些必要的信息：ID，描述信息，坐标和快速统计，和你前面看到的一样。但是片段信息完全不同。我们得到了有两行数据的表格，而不是显示 ‘Fragment: 1’。

这就是 Bio.SearchIO 处理含有多片段HSP的方式。和前面提到的一样，一个 HSP比对可能会被内含子分成多个片段。内含子不是query-hit匹配的一部分，所以它们不能被当成query或hit序列的一部分。但是，它们确实影响我们处理序列坐标，所以我们不能忽视。

看看上面的HSP的hit坐标。在 Hit range 区域，我们看到坐标是 [54233104:54264463]。但是看看表格中的行，我们发现不是坐标跨度的所有区域都能匹配我们的query。特殊的是，间断区域从 54233122 到 54264420 。

你可能会问，为什么query坐标好像是邻近的?这是很好的。在这个例子中，query是连续的（无间断区域），但是hit却不是。

所有的这些属性都是可以直接从HSP获取的，通过这样的方式：

>>> blat_hsp2.hit_range         # hit start and end coordinates of the entire HSP
(54233104, 54264463)
>>> blat_hsp2.hit_range_all     # hit start and end coordinates of each fragment
[(54233104, 54233122), (54264420, 54264463)]
>>> blat_hsp2.hit_span          # hit span of the entire HSP
31359
>>> blat_hsp2.hit_span_all      # hit span of each fragment
[18, 43]
>>> blat_hsp2.hit_inter_ranges  # start and end coordinates of intervening regions in the hit sequence
[(54233122, 54264420)]
>>> blat_hsp2.hit_inter_spans   # span of intervening regions in the hit sequence
[31298]

这些属性中大多数都不能简单地从PSL文件获得，但是当你分析PSL文件时， Bio.SearchIO 会动态地帮你计算。它需要的只是每个片段的开始和结束坐标。

query， hit，和 aln 属性又是什么情况？如果HSP含有多个片段，你就不能使用这些属性，因为它们只取回单个 SeqRecord 或 MultipleSeqAlignment 对象。但是，你可以用相应的 *_all 方法： query_all， hit_all，和 aln_all。这些属性会返回包含每个HSP 片段的 SeqRecord 或 MultipleSeqAlignment 对象的列表。还有其他相同功能的属性，也就是只对只有一个片段的HSP有效。查看 HSP 文档获得完整的列表。

最后，想要检查是否是多片段HSP，你可以用 is_fragmented 属性：

>>> blat_hsp2.is_fragmented     # BLAT HSP with 2 fragments
True
>>> blat_hsp.is_fragmented      # BLAT HSP from earlier, with one fragment
False

在进入下部分之前，你只需要了解我们可以对 HSP 对象使用切片，和 QueryResult 或 Hit 对象一样。当你使用切片的时候，会返回一个 HSPFragment 对象。

8.1.4 HSP片段¶

HSPFragment 代表query和hit之间单个连续匹配。应该把它当作对象模型和搜索结果的核心，因为它决定你的搜索是否有结果。

在多数情况下，你不必直接处理 HSPFragment 对象，因为没有那么多搜索工具断裂它们的HSP。当你确实需要处理它们时，需要记住的是 HSPFragment 对象要被写地尽量压缩。在多数情况下，它们仅仅包含直接与序列有关的属性：正负链，阅读框，字母表，位置坐标，序列本身以及它们的ID和描述。

当你对 HSPFragment 对象调用 print 方法时，这些属性可以非常简单地显示出来。这里有个从我们BLAST搜索得到的例子：

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> blast_frag = blast_qresult[0][0][0]    # first hit, first hsp, first fragment
>>> print blast_frag
      Query: 42291 mystery_seq
        Hit: gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 520b (MIR520...
Query range: [0:61] (1)
  Hit range: [0:61] (1)
  Fragments: 1 (61 columns)
     Query - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
             |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
       Hit - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG

在这个水平上，BLAT和BLAST片段看起来非常相似，除了没有出现的query和hit序列：

>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> blat_frag = blat_qresult[0][0][0]    # first hit, first hsp, first fragment
>>> print blat_frag
      Query: mystery_seq <unknown description>
        Hit: chr19 <unknown description>
Query range: [0:61] (1)
  Hit range: [54204480:54204541] (1)
  Fragments: 1 (? columns)

在所有情况下，这些属性都可以通过我们最爱的点标记访问。一些例子：

>>> blast_frag.query_start      # query start coordinate
0
>>> blast_frag.hit_strand       # hit sequence strand
1
>>> blast_frag.hit              # hit sequence, as a SeqRecord object
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG', DNAAlphabet()), id='gi|262205317|ref|NR_030195.1|', name='aligned hit sequence', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])

8.2 一个关于标准和惯例的注意事项¶

在我们进入到主要功能前，你需要知道 Bio.SearchIO 使用的一些标准。如果你已经接触过多序列搜索工具，你可能必须面对每个程序处理事情方式不同的问题，如序列位置坐标。这可能不是一个令人高兴的经历，因为这些搜索工具通常有它们自己的标准。例如，一种工具可能使用“从1开始”(one-based)的坐标，而其他工具使用“从0开始”(zero-based)的坐标。或者，一种程序在处理负链时，可能会反转开始和结束坐标，而其他程序确不会。简而言之，会产生一些必须要处理的混乱。

我们意识到这种问题，并且打算在 Bio.SearchIO 中解决。毕竟， Bio.SearchIO 的目标之一就是创建一个通用简单的接口来处理多种不同的搜索输出文件。意味着要制定一个超越你所见的对象模型的标准。

现在，你可能抱怨，”不要又来一个标准“。好吧，最后我们必须选择一个标准，这是必须的。并且，我们并不是创造一个全新的事物；只是采用一个我们觉得对Python 使用者最好的标准（这是Biopython，毕竟）。

在使用 Bio.SearchIO 时你可以认为有个三个隐含的标准：

第一个适用于序列坐标。在 Bio.SearchIO 模块中，所有序列坐标遵循Python 的坐标风格：从0开始，半开区间。例如，在一个BLAST XML输出文件中，HSP的起始和结束坐标是10和28，它们在 Bio.SearchIO 中将变成9和28。起始坐标变成9因为Python 中索引是从0开始，而结束坐标仍然是28因为Python索引删除了区间中最后一个项目。
第二个是关于序列坐标顺序。在 Bio.SearchIO 中，开始坐标总是小于或等于结束坐标。但是这不是在所有的序列搜索工具中都始终适用。因为当序列为负链时，起始坐标会更大一些。
最后一个标准是关于链和阅读框的值。对于链值，只有四个可选值： 1 (正链)， -1 (负链)， 0 (蛋白序列)，和 None (无链)。对于阅读框，可选值是从 -3 至 3 的整型以及 None 。

注意，这些标准只是存在于 Bio.SearchIO 对象中。如果你把 Bio.SearchIO 对象写入一种输出格式， Bio.SearchIO 会使用该格式的标准来输出。它并不强加它的标准到你的输出文件。

8.3 读取搜索输出文件¶

有两个方法，你可以用来读取搜索输出文件到 Bio.SearchIO 对象： read 和 parse。它们和其他亚模块如 Bio.SeqIO 或 Bio.AlignIO 中的 read 和 parse 方法在本质上是相似的。你都需要提供搜索输出文件名和文件格式名，都是Python字符串类型。你可以查阅文档来获得 Bio.SearchIO 可以识别的格式清单。

Bio.SearchIO.read 用于读取单query的搜索输出文件并且返回一个 QueryResult 对象。你在前面的例子中已经看到过 read 的使用了。你没看到的是， read 同样接受额外的关键字参数，取决于文件的格式。

这里有一些例子。在第一个例子中，我们和前面一样用 read 读BLAST表格输出文件。在第二个例子中，我们用一个关键字来修饰，所以它分析带有注释的BLAST 表格变量。

>>> from Bio import SearchIO
>>> qresult = SearchIO.read('tab_2226_tblastn_003.txt', 'blast-tab')
>>> qresult
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)
>>> qresult2 = SearchIO.read('tab_2226_tblastn_007.txt', 'blast-tab', comments=True)
>>> qresult2
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)

这些关键字在不同的文件格式中是不一样的。查看格式文档，看看它是否有关键字参数来控制它的分析器行为。

对于 Bio.SearchIO.parse，是用来读取含有任意数量query的搜索输出文件。这个方法返回一个generator对象，在每次迭代中yield一个 QueryResult 对象。和 Bio.SearchIO.read 一样，它同样接受格式特异的关键字参数：

>>> from Bio import SearchIO
>>> qresults = SearchIO.parse('tab_2226_tblastn_001.txt', 'blast-tab')
>>> for qresult in qresults:
...     print qresult.id
gi|16080617|ref|NP_391444.1|
gi|11464971:4-101
>>> qresults2 = SearchIO.parse('tab_2226_tblastn_005.txt', 'blast-tab', comments=True)
>>> for qresult in qresults2:
...     print qresult.id
random_s00
gi|16080617|ref|NP_391444.1|
gi|11464971:4-101

8.4 用索引处理含有大量搜索输出的文件¶

有时，你得到了一个包含成百上千个query的搜索输出文件要分析，你当然可以使用 Bio.SearchIO.parse 来处理，但是如果你仅仅需要访问少数query的话，效率是及其低下的。这是因为 parse 会分析所有的query，直到找到你感兴趣。

在这种情况下，理想的选择是用 Bio.SearchIO.index 或 Bio.SearchIO.index_db 来索引文件。如果名字听起来很熟悉，是因为你之前已经见过了，在章节 5.4.2 。这些方法和 Bio.SeqIO 中相应的方法行为很相似，只是多了些格式特异的关键字参数。

这里有一些例子。你可以只用文件名和格式名来 index

>>> from Bio import SearchIO
>>> idx = SearchIO.index('tab_2226_tblastn_001.txt', 'blast-tab')
>>> sorted(idx.keys())
['gi|11464971:4-101', 'gi|16080617|ref|NP_391444.1|']
>>> idx['gi|16080617|ref|NP_391444.1|']
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)

或者依旧使用格式特异的关键字参数：

>>> idx = SearchIO.index('tab_2226_tblastn_005.txt', 'blast-tab', comments=True)
>>> sorted(idx.keys())
['gi|11464971:4-101', 'gi|16080617|ref|NP_391444.1|', 'random_s00']
>>> idx['gi|16080617|ref|NP_391444.1|']
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)

或者使用 key_function 参数，和 Bio.SeqIO 中一样：

>>> key_function = lambda id: id.upper()    # capitalizes the keys
>>> idx = SearchIO.index('tab_2226_tblastn_001.txt', 'blast-tab', key_function=key_function)
>>> sorted(idx.keys())
['GI|11464971:4-101', 'GI|16080617|REF|NP_391444.1|']
>>> idx['GI|16080617|REF|NP_391444.1|']
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)

Bio.SearchIO.index_db 和 index 作用差不多，不同的只是它把query 偏移量写入一个SQLite数据库文件中。

8.5 写入和转换搜索输出文件¶

有时候，读取一个搜索输出文件，作些调整并写到一个新的文件是很有用的。 Bio.SearchIO 提供了一个 write 方法，让你可以准确地完成这种工作。它需要的参数是：一个可迭代返回 QueryResult 的对象，输出文件名，输出文件格式和一些可选的格式特异的关键字参数。它返回一个4项目的元组，分别代表被写入的 QueryResult， Hit， HSP，和 HSPFragment 对象的数量。

>>> from Bio import SearchIO
>>> qresults = SearchIO.parse('mirna.xml', 'blast-xml')     # read XML file
>>> SearchIO.write(qresults, 'results.tab', 'blast-tab')    # write to tabular file
(3, 239, 277, 277)

你应该注意，不同的文件格式需要 QueryResult， Hit， HSP 和 HSPFragment 对象的不同属性。如果这些属性不存在，那么将不能写入。也就是，你想写入的格式可能有时也会失效。举个例子，如果你读取一个BLASTXML文件，你就不能将结果写入PSL文件，因为PSL文件需要一些属性，而这些属性BLAST却不能提供（如重复匹配的数量）。如果你确实想写到PSL，可以手工设置这些属性。

和 read， parse， index 和 index_db 相似， write 同样接受格式特异的关键字参数。查阅文档获得 Bio.SearchIO 可写格式和这些格式的参数的完整清单。

最后， Bio.SearchIO 同样提供一个 convert 方法，可以理解为 Bio.SearchIO.parse 和 Bio.SearchIO.write 的简单替代方法。使用 convert 方法的例子如下：

>>> from Bio import SearchIO
>>> SearchIO.convert('mirna.xml', 'blast-xml', 'results.tab', 'blast-tab')
(3, 239, 277, 277)

因为 convert 使用 write 方法，所以只有所有需要的属性都存在时，格式转换才能正常工作。这里由于BLAST XML文件提供BLAST 表格文件所需的所有默认值，格式转换才能正常完成。但是，其他格式转换就可能不会正常工作，因为你需要先手工指定所需的属性。