BLAST 和其他序列搜索工具
生物序列识别是生物信息学中不可或缺的一部分。有几种工具可用于此，每种工具都有自己的算法和方法，例如 BLAST（可以说是最流行的），FASTA，HMMER 等等。一般来说，这些工具通常使用您的序列搜索潜在匹配的数据库。随着已知序列数量的不断增加（因此潜在匹配的数量也在不断增加），解释结果变得越来越困难，因为可能存在数百甚至数千个潜在匹配。当然，手动解释这些搜索结果是不可能的。此外，您通常需要使用多个序列搜索工具，每个工具都有自己的统计数据、约定和输出格式。想象一下，当您需要使用多个搜索工具处理多个序列时，这将是多么艰巨的任务。
我们自己也深知这一点，这就是我们在 Biopython 中创建 Bio.SearchIO 子模块的原因。Bio.SearchIO 允许您以便捷的方式从搜索结果中提取信息，同时处理不同搜索工具使用的不同标准和约定。名称 SearchIO 是对 BioPerl 中同名模块的致敬。
在本章中，我们将介绍 Bio.SearchIO 的主要功能，以展示它可以为您做什么。在此过程中，我们将使用两种流行的搜索工具：BLAST 和 BLAT。它们仅用于说明目的，您应该能够轻松地将工作流程适应 Bio.SearchIO 支持的任何其他搜索工具。欢迎您使用我们将要使用的搜索输出文件进行操作。BLAST 输出文件可以从 这里 下载，BLAT 输出文件可以从 这里 下载，或者包含在 Biopython 源代码中，位于 Doc/examples/ 文件夹下。这两个输出文件都是使用以下序列生成的
>mystery_seq
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
BLAST 结果是使用 blastn 对 NCBI refseq_rna 数据库进行搜索生成的 XML 文件。对于 BLAT，序列数据库是 2009 年 2 月的人类基因组草图 hg19，输出格式是 PSL。
我们将从介绍 Bio.SearchIO 对象模型开始。该模型是搜索结果的表示，因此它是 Bio.SearchIO 本身的核心。之后，我们将查看 Bio.SearchIO 中您可能经常使用的一些主要功能。
现在我们已经准备就绪，让我们进入第一步：介绍核心对象模型。
SearchIO 对象模型
尽管许多序列搜索工具的输出风格差异很大，但事实证明它们的底层概念是相似的
输出文件可能包含来自一个或多个搜索查询的结果。
在每个搜索查询中，您将看到来自给定搜索数据库的一个或多个匹配项。
在每个数据库匹配项中，您将看到一个或多个区域，其中包含查询序列和数据库序列之间的实际序列比对。
某些程序，如 BLAT 或 Exonerate，可能会将这些区域进一步细分为多个比对片段（或 BLAT 中的块，以及 exonerate 中可能的内含子）。这不是您始终会看到的情况，因为像 BLAST 和 HMMER 这样的程序不会这样做。
意识到这种普遍性，我们决定将其用作创建 Bio.SearchIO 对象模型的基础。该对象模型由嵌套层次结构的 Python 对象组成，每个对象都代表上面概述的一个概念。这些对象是
QueryResult，用于表示单个搜索查询。Hit，用于表示单个数据库匹配项。Hit对象包含在QueryResult中，并且每个QueryResult中有零个或多个Hit对象。HSP（高分对的缩写），用于表示查询和匹配序列之间重要比对的区域。HSP对象包含在Hit对象中，并且每个Hit有一个或多个HSP对象。HSPFragment，用于表示查询和匹配序列之间单个连续比对。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用于具有单个查询的搜索输出文件，并返回一个QueryResult对象parse用于具有多个查询的搜索输出文件，并返回一个生成器，生成QueryResult对象
有了这些，让我们开始探查每个 Bio.SearchIO 对象，从 QueryResult 开始。
QueryResult
QueryResult 对象表示单个搜索查询，并包含零个或多个 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 版本 2.2.27+）
查询 ID、描述及其序列长度（ID 为 42291，描述为 ‘mystery_seq’，长度为 61 个核苷酸）
要搜索的目标数据库（refseq_rna）
结果匹配项的快速概述。对于我们的查询序列，有 100 个潜在匹配项（在表格中编号为 0-99）。对于每个匹配项，我们还可以看到它包含多少个 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>
您会立即注意到有一些区别。其中一些是由 PSL 格式存储详细信息的方式造成的，您将看到。其余部分是由于我们的 BLAST 和 BLAT 搜索之间的真实程序和目标数据库差异造成的
程序名称和版本。
Bio.SearchIO知道该程序是 BLAT，但在输出文件中没有关于程序版本的任何信息，因此它默认为 ‘<unknown version>’。查询 ID、描述及其序列长度。请注意，这些详细信息与我们在 BLAST 中看到的略有不同。ID 是 ‘mystery_seq’ 而不是 42991，没有已知描述，但查询长度仍然是 61。这实际上是文件格式本身引入的差异。BLAST 有时会创建它自己的查询 ID，并将您的原始 ID 用作序列描述。
目标数据库未知，因为 BLAT 输出文件中没有说明。
最后，我们拥有的匹配项列表完全不同。这里，我们看到我们的查询序列只匹配 ‘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
有关可访问属性的完整列表，您可以查看每个特定于格式的文档。例如，Bio.SearchIO.BlastIO 和 Bio.SearchIO.BlatIO。
在研究了对 QueryResult 对象使用 print 之后，让我们深入探究。QueryResult 究竟是什么？在 Python 对象方面，QueryResult 是列表和字典的混合体。换句话说，它是一个包含列表和字典所有便利功能的容器对象。
与 Python 列表和字典类似，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 包含多少个项目（命中），您只需调用 Python 的 len 方法。
>>> len(blast_qresult)
100
>>> len(blat_qresult)
1
与 Python 列表类似，您可以使用切片符号从 QueryResult 中检索项目（命中）。
>>> 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)
要检索多个命中，您也可以使用切片符号来切片 QueryResult 对象。在这种情况下，切片将返回一个新的 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 ...
与 Python 字典类似，您还可以使用命中的 ID 来检索命中。如果您知道给定的命中 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
如果您只想检查特定命中是否出现在查询结果中，您可以使用 in 关键字进行简单的 Python 成员资格测试。
>>> "gi|262205317|ref|NR_030195.1|" in blast_qresult
True
>>> "gi|262205317|ref|NR_030194.1|" in blast_qresult
False
有时，仅知道命中是否存在还不够；您还想了解命中的排名。在这里，index 方法可以提供帮助。
>>> blast_qresult.index("gi|301171437|ref|NR_035870.1|")
22
请记住，我们在这里使用 Python 的索引样式，它是从零开始的。这意味着我们上面的命中排名在第 23 位，而不是第 22 位。
此外，请注意，您在这里看到的命中排名基于原始搜索输出文件中的原生命中排序。不同的搜索工具可能会根据不同的标准对这些命中进行排序。
如果原生命中排序不符合您的喜好，您可以使用 QueryResult 对象的 sort 方法。它与 Python 的 list.sort 方法非常相似，只是增加了是否创建新的排序 QueryResult 对象的选择。
以下是如何使用 QueryResult.sort 根据每个命中的完整序列长度对命中进行排序的示例。对于这种特定排序，我们将 in_place 标志设置为 False，以便排序将返回一个新的 QueryResult 对象，并将我们的初始对象保持为未排序状态。我们还将 reverse 标志设置为 True，以便我们以降序进行排序。
>>> for hit in blast_qresult[:5]: # id and sequence length of the first five hits
... print("%s %i" % (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("%s %i" % (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 的列表推导式、生成器表达式或内置的 filter 和 map 函数，您就会知道它们对于使用类似列表的对象有多么有用（如果您不熟悉，请查看它们！）。您可以使用这些内置方法来操作 QueryResult 对象，但您最终会得到普通的 Python 列表，并且会失去执行更多有趣操作的能力。
这就是 QueryResult 对象提供其自身的 filter 和 map 方法版本的原因。类似于 filter，存在 hit_filter 和 hsp_filter 方法。顾名思义，这些方法会根据其 Hit 对象或 HSP 对象对 QueryResult 对象进行过滤。类似地，类似于 map，QueryResult 对象还提供 hit_map 和 hsp_map 方法。这些方法分别将给定的函数应用于 QueryResult 对象中的所有命中或 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("%s %i" % (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 对象执行过滤。
至于 map 方法，它们也接受回调函数作为参数。但是，它们不是返回 True 或 False，而是回调函数必须返回修改后的 Hit 或 HSP 对象（取决于您使用的是 hit_map 还是 hsp_map）。
让我们看一个使用 hit_map 重命名命中 ID 的示例。
>>> def map_func(hit):
... # renames "gi|301171322|ref|NR_035857.1|" to "NR_035857.1"
... hit.id = hit.id.split("|")[3]
... 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 对象。
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]
您会看到我们在这里涵盖了基本内容。
查询 ID 和描述都存在。命中始终与查询相关联，因此我们也希望跟踪源查询。可以使用
query_id和query_description属性从命中中访问这些值。我们还有唯一的命中 ID、描述和完整序列长度。可以使用
id、description和seq_len分别访问它们。最后，还有一个表格包含有关此命中包含的 HSP 的快速信息。在每一行中，我们列出了重要的 HSP 细节：HSP 索引、其 e 值、其位得分、其跨度（包括间隙的比对长度）、其查询坐标和其命中坐标。
现在让我们将其与 BLAT 搜索进行对比。请记住，在 BLAT 搜索中，我们有一个包含 17 个 HSP 的命中。
>>> 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 命中类似的详细程度。但是，有一些差异值得解释。
e 值和位得分列值。由于 BLAT HSP 没有 e 值和位得分，因此显示默认值为“？”。
跨度列呢？跨度值用于显示完整的比对长度，包括所有残基和可能存在的任何间隙。PSL 格式没有此信息可以直接获取，并且
Bio.SearchIO不会尝试猜测它是什么，因此我们得到了与 e 值和位得分列类似的“？”。
从 Python 对象的角度来看，Hit 的行为几乎与 Python 列表相同，但只包含 HSP 对象。如果您熟悉列表，那么在使用 Hit 对象时应该不会遇到任何困难。
就像 Python 列表一样，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 还是多个 HSP 对象。与 QueryResult 一样，如果您切片多个 HSP，则会返回一个新的 Hit 对象，其中只包含切片后的 HSP 对象。
>>> 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)。Hit.filter 和 Hit.map 都作用于 Hit 所包含的 HSP 对象。
HSP
HSP（高得分对）表示命中序列中与查询序列存在显著比对的区域。它包含了您的查询序列和数据库条目之间的实际匹配。由于这种匹配是由序列搜索工具的算法决定的，因此 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 会显示它的常规详细信息。
有查询和命中 ID 和描述。我们需要这些来识别我们的
HSP。我们也得到了查询和命中序列的匹配范围。我们在这里使用的切片符号表明范围使用 Python 的索引风格显示（从零开始，半开）。括号内的数字表示链。在本例中，两个序列都具有正链。
一些快速统计数据可用：e 值和 bitscore。
有关于 HSP 片段的信息。现在忽略它；稍后会解释。
最后，我们有查询和命中序列比对本身。
这些详细信息可以使用点符号单独访问，就像在 QueryResult 和 Hit 中一样。
>>> 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
查看 Bio.SearchIO 下的 HSP 文档以获取这些预定义属性的完整列表。
此外，每个序列搜索工具通常会为其 HSP 对象计算自己的统计数据/详细信息。例如，XML BLAST 搜索还会输出间隙数和相同残基数。这些属性可以通过以下方式访问。
>>> 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'), id='42291', name='aligned query sequence', description='mystery_seq', dbxrefs=[])
>>> blast_hsp.hit
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG'), id='gi|262205317|ref|NR_030195.1|', name='aligned hit sequence', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])
它们是您在第 序列注释对象 节中看到的 SeqRecord 对象！这意味着您可以在 HSP.query 和/或 HSP.hit 上对 SeqRecord 对象执行各种有趣的操作。
您现在应该不会感到惊讶，HSP 对象具有一个 alignment 属性，它是一个 MultipleSeqAlignment 对象。
>>> print(blast_hsp.aln)
Alignment with 2 rows and 61 columns
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAG...GGG 42291
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAG...GGG gi|262205317|ref|NR_030195.1|
在探究了 BLAST HSP 之后，我们现在来看看 BLAT 结果中的 HSP，了解另一种 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)
您可能已经猜到了一些输出。我们有查询和命中 ID 和描述以及序列坐标。evalue 和 bitscore 的值为“？”，因为 BLAT HSP 没有这些属性。但这里最大的区别是，您没有看到任何序列比对显示。如果您仔细观察，PSL 格式本身没有任何命中或查询序列，因此 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
但这不会影响其他属性。例如，您仍然可以访问查询或命中比对的长度。尽管没有显示任何属性，但 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 搜索中，有时我们会将结果分成“块”。这些块本质上是比对片段，它们之间可能存在一些间隔序列。
让我们看看一个包含多个块的 BLAT HSP，看看 Bio.SearchIO 如何处理它。
>>> 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 和描述、坐标以及快速统计数据与您之前看到的类似。但片段详细信息完全不同。我们现在没有显示“片段：1”，而是一个包含两行数据的表格。
这就是 Bio.SearchIO 处理具有多个片段的 HSP 的方式。如前所述，HSP 比对可能被间隔序列分成片段。间隔序列不是查询-命中匹配的一部分，因此不应被视为查询或命中序列的一部分。但是，它们确实会影响我们处理序列坐标的方式，因此我们不能忽略它们。
看看上面 HSP 的命中坐标。在 Hit range: 字段中，我们看到坐标是 [54233104:54264463]。但查看表格行，我们看到，并非此坐标跨越的整个区域都与我们的查询匹配。具体来说，间隔区域从 54233122 跨越到 54264420。
那么，为什么查询坐标看起来是连续的呢？这完全没问题。在这种情况下，这意味着查询匹配是连续的（没有间隔区域），而命中匹配不是。
顺便说一下，所有这些属性都可以直接从 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 文件中获取，但 Bio.SearchIO 会在您解析 PSL 文件时动态地为您计算它们。它只需要每个片段的起始和结束坐标。
那么 query、hit 和 aln 属性呢？如果 HSP 包含多个片段，则无法使用这些属性，因为它们仅获取单个 SeqRecord 或 MultipleSeqAlignment 对象。但是，可以使用它们的 *_all 对应项：query_all、hit_all 和 aln_all。这些属性将返回一个列表，其中包含来自每个 HSP 片段的 SeqRecord 或 MultipleSeqAlignment 对象。还有一些其他属性的行为类似，即它们仅适用于包含一个片段的 HSP。查看 Bio.SearchIO 下的 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 对象，它是对象模型的最后一个组件。
HSPFragment
HSPFragment 代表查询序列和匹配序列之间的一个连续匹配。您可以将其视为对象模型和搜索结果的核心，因为这些片段的存在决定了您的搜索是否有结果。
在大多数情况下，您不必直接处理 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 片段非常相似，除了查询序列和匹配序列不在其中
>>> 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'), id='gi|262205317|ref|NR_030195.1|', name='aligned hit sequence', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])
关于标准和约定的说明
在我们继续进行主要功能之前，您应该了解 Bio.SearchIO 使用的标准。如果您使用过多个序列搜索工具，您可能不得不处理每个程序处理序列坐标的方式各不相同。这可能不是一种愉快的体验，因为这些搜索工具通常有自己的标准。例如，一个工具可能使用基于 1 的坐标，而另一个工具使用基于 0 的坐标。或者，一个程序可能会在链为负数的情况下反转起始坐标和结束坐标，而其他程序则不会。简而言之，这些通常会造成不必要的混乱，必须加以处理。
我们自己也意识到了这个问题，我们打算在 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 将使用该格式的标准进行输出。它不会将自己的标准强加到您的输出文件上。
读取搜索输出文件
您可以使用两个函数将搜索输出文件读入 Bio.SearchIO 对象：read 和 parse。它们本质上类似于其他子模块（如 Bio.SeqIO 或 Bio.AlignIO）中的 read 和 parse 函数。在这两种情况下，您都需要提供搜索输出文件名和文件格式名称，两者都作为 Python 字符串。您可以查看文档以了解 Bio.SearchIO 识别的格式名称列表。
Bio.SearchIO.read 用于读取只有一个查询的搜索输出文件，并返回一个 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，它用于读取包含任意数量查询的搜索输出文件。该函数返回一个生成器对象，它在每次迭代时生成一个 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
使用索引处理大型搜索输出文件
有时，您会获得一个包含数百或数千个查询的搜索输出文件，您需要对其进行解析。当然，您可以使用 Bio.SearchIO.parse 来处理此文件，但这在您只需要访问其中几个查询时效率极低。这是因为 parse 会在获取您感兴趣的查询之前解析它看到的所有查询。
在这种情况下，理想的选择是使用 Bio.SearchIO.index 或 Bio.SearchIO.index_db 对文件进行索引。如果这些名称听起来很熟悉，那是因为您之前在第 序列文件作为字典 - 索引文件 节中见过它们。这些函数的行为也类似于它们的 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.close()
或者也可以使用特定于格式的关键字参数
>>> 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)
>>> idx.close()
或者使用 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)
>>> idx.close()
Bio.SearchIO.index_db 的工作原理与 index 相似，只是它将查询偏移量写入 SQLite 数据库文件。
写入和转换搜索输出文件
有时需要能够操作输出文件中的搜索结果并将其写入新文件。 Bio.SearchIO 提供了一个 write 函数，它允许你执行此操作。它接受一个返回 QueryResult 对象的可迭代对象、要写入的输出文件名、要写入的格式名称以及可选的一些特定于格式的关键字参数。它返回一个包含四个元素的元组，表示写入的 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 对象的不同属性。如果这些属性不存在，则写入将无法进行。换句话说，你不能总是写入你想要的输出格式。例如，如果你读取一个 BLAST XML 文件，你将无法将结果写入 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 表格文件所需的所有默认值，因此它可以正常工作。但是，其他格式转换不太可能起作用，因为你需要首先手动分配必需的属性。