AMPS和其他序列搜索工具
生物序列识别是生物信息学不可或缺的一部分。有多种工具可用于此，每种工具都有自己的算法和方法，例如AMPS（可以说是最流行的）、FASTA、HMMER等。一般来说，这些工具通常使用您的序列来搜索潜在匹配的数据库。随着已知序列数量的不断增加（因此潜在匹配数量的不断增加），解释结果变得越来越困难，因为可能有数百甚至数千个潜在匹配。当然，手动解释这些搜索结果是不可能的。此外，您经常需要使用多个序列搜索工具，每个工具都有自己的统计数据、约定和输出格式。想象一下，当您需要使用多种搜索工具处理多个序列时，这将是多么令人畏惧。
我们自己也很清楚这一点，这就是为什么我们创造了 Bio.SearchIO Biopython的子模块。 Bio.SearchIO 允许您以方便的方式从搜索结果中提取信息，同时还可以处理不同搜索工具使用的不同标准和惯例。名称 SearchIO 是对BioPerl同名模块的致敬。
在本章中，我们将介绍的主要功能 Bio.SearchIO to show what it can do for you. We’ll use two popular search tools along the way: BLAST and BLAT. They are used merely for illustrative purposes, and you should be able to adapt the workflow to any other search tools supported by Bio.SearchIO in a breeze. You’re very welcome to follow along with the search output files we’ll be using. The BLAST output file can be downloaded here ，以及BLAT输出文件 here 或包含在Biopython源代码中 Doc/examples/ 文件夹.两个输出文件都是使用以下序列生成的：
>mystery_seq
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
BST结果是使用以下方式生成的ML文件 blastn 针对ncbi refseq_rna 数据库对于BLAT，序列数据库是2009年2月 hg19 人类基因组草案，输出格式为PSL。
我们将从介绍 Bio.SearchIO 对象模型。该模型是搜索结果的代表，因此它是 Bio.SearchIO 本身之后，我们将查看中的主要功能 Bio.SearchIO 您可能经常使用的。
现在我们都准备好了，让我们进入第一步：引入核心对象模型。
SearchIO对象模型
尽管许多序列搜索工具的输出风格差异很大，但事实证明它们的基本概念是相似的：
输出文件可能包含一个或多个搜索查询的结果。
在每个搜索查询中，您将看到来自给定搜索数据库的一个或多个命中。
在每个数据库命中中，您将看到一个或多个包含查询序列和数据库序列之间的实际序列比对的区域。
一些程序，如BLAT或Exonerate可能会进一步将这些区域分裂成几个比对片段（或BLAT中的片段，可能还有exonerate中的exonerate）。这并不是您经常看到的事情，因为像AMPS和HMMER这样的程序不会做到这一点。
意识到这种普遍性，我们决定使用它作为创建 Bio.SearchIO 对象模型。对象模型由Python对象的嵌套层次结构组成，每个对象代表上面概述的一个概念。这些对象是：
QueryResult，以表示单个搜索查询。Hit，以表示单个数据库命中。Hit对象包含在QueryResult和每个QueryResult有零或更多Hit对象HSP（高分对的缩写），代表查询和命中序列之间的重要比对区域。HSP对象包含在Hit对象和每个Hit具有一个或多HSP对象HSPFragment，以表示查询和命中序列之间的单个连续比对。HSPFragment对象包含在HSP对象大多数序列搜索工具，如BST和HMMER统一HSP和HSPFragment对象作为每个HSP只有一个HSPFragment.然而，有一些工具，如BLAT和Exonerate可以产生HSP包含多HSPFragment.如果现在看起来有点混乱，请不要担心，稍后我们将详细说明这两个对象。
这四个对象是您使用时将与之互动的对象 Bio.SearchIO .它们是使用主要的 Bio.SearchIO 方法： read , parse , index ，或者 index_db .这些方法的详细信息将在后面的部分中提供。对于本节，我们将仅使用读取和解析。这些函数的行为与它们的类似 Bio.SeqIO 和 Bio.AlignIO 同行：
read用于通过单个查询搜索输出文件并返回QueryResult对象parse用于搜索具有多个查询的输出文件，并返回生成器QueryResult对象
解决了问题，让我们开始探索每个 Bio.SearchIO 对象，以 QueryResult .
QueryResult
CredyResponse对象表示单个搜索查询并包含零个或多个Hit对象。让我们看看使用我们拥有的BST文件是什么样子：
>>> 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...
我们刚刚开始触及对象模型的表面，但您可以看到已经有一些有用的信息。通过调用 print 上 QueryResult 对象，您可以看到：
程序名称和版本（blackmail版本2.2.27+）
查询ID、描述及其序列长度（ID为42291，描述为“神秘_seq”，长度为61个核苷酸）
要搜索的目标数据库（refseq_rna）
最终点击的快速概述。对于我们的查询序列，有100个潜在命中（表中编号为0-99）。对于每个点击，我们还可以看到它包含多少个HSPs、它的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格式存储其详细信息的方式引起的。其余的是由于我们的BST和BLAT搜索之间的真实程序和目标数据库差异造成的：
程序名称和版本。
Bio.SearchIO知道该程序是BLAT，但在输出文件中没有有关程序版本的信息，因此默认为“<unknown version>”。查询ID、描述及其序列长度。请注意，这些细节与我们在AMPS中看到的略有不同。ID为“mysty_seq”而不是42991，没有已知的描述，但查询长度仍然是61。这实际上是文件格式本身引入的差异。AMPS有时会创建自己的查询ID并使用您的原始ID作为序列描述。
目标数据库未知，因为BLT输出文件中没有说明它。
最后，我们拥有的热门歌曲列表完全不同。在这里，我们看到我们的查询序列仅命中“chr19”数据库条目，但在其中我们看到17个热休克蛋白区域。然而，这并不奇怪，因为我们正在使用不同的程序，每个程序都有自己的目标数据库。
调用时您看到的所有详细信息 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 .
研究过使用 print 对 QueryResult 对象，让我们深入挖掘。究竟什么是 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)
...
要检查有多少项（点击）a 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)
您还可以获得完整列表 Hit 对象使用 hits 以及完整的列表 Hit ID使用 hit_keys :
>>> blast_qresult.hits
[...] # list of all hits
>>> blast_qresult.hit_keys
[...] # list of all hit IDs
如果您只是想检查查询结果中是否存在特定的命中，该怎么办？您可以使用进行简单的Python成员资格测试 in 关键字：
>>> "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位。
另外，请注意，您在这里看到的点击率排名基于原始搜索输出文件中存在的本地点击率顺序。不同的搜索工具可能会根据不同的标准对这些点击进行排序。
如果原生热门订购不适合您的口味，您可以使用 sort 方法 QueryResult object.它与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 方法.顾名思义，这些方法过滤 QueryResult 对象要么在其 Hit 对象或 HSP 对象同样，类似于 map , QueryResult 对象还提供 hit_map 和 hsp_map 方法.这些方法将给定的函数应用于一个或多个匹配中的所有命中或HSP。 QueryResult 分别是对象。
让我们看看这些方法的实际应用，首先 hit_filter .此方法接受回调函数，该函数检查给定的 Hit 对象是否通过您设置的条件。换句话说，该函数必须接受一个参数 Hit 对象并返回 True 或 False .
这是一个使用的例子 hit_filter 以滤除 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 对象表示来自单个数据库条目的所有查询结果。它们是中的二级容器 Bio.SearchIO 对象层次结构。你已经看到它们被包含在 QueryResult 物体，但它们本身包含 HSP 对象
让我们看看它们是什么样子，从我们的AMPS搜索开始：
>>> 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，分别。最后，有一个表格包含有关此热门内容包含的热休克蛋白的快速信息。在每一行中，我们都列出了重要的热休克详细信息：热休克指数、其e值、其位得分、其跨度（包括间隙的对齐长度）、其查询坐标和其命中坐标。
现在让我们将其与BLAT搜索进行比较。请记住，在BLAT搜索中，我们发现了一次与17个热休克蛋白的匹配。
>>> 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]
在这里，我们获得了与我们之前看到的井喷类似的细节水平。不过，也有一些差异值得解释：
e-value和bit score列值。由于BLAT HSP没有e值和位评分，因此显示默认为“？”。
跨度柱怎么样？跨度值旨在显示完整的对齐长度，其由所有残基和可能存在的任何间隙组成。PSL格式没有现成的此信息，
Bio.SearchIO没有尝试猜测它是什么，所以我们得到一个“？”类似于e值和bit score列。
就Python对象而言， Hit 行为几乎与Python列表相同，但包含 HSP 唯一的对象。如果您熟悉列表，那么使用 Hit object.
就像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)
您可以调用 len 上 Hit 看到有多少 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]
您还可以对 HSP 内 Hit ，使用与您在 QueryResult object.
Finally, there are also the filter and map methods you can use
on Hit objects. Unlike in the QueryResult object, Hit
objects only have one variant of filter (Hit.filter) and one
variant of map (Hit.map). Both of Hit.filter and Hit.map
work on the HSP objects a Hit has.
HSP
HSP （高分对）代表命中序列中包含与查询序列的显着比对的区域。它包含您的查询序列和数据库条目之间的实际匹配。由于该匹配是由序列搜索工具的算法确定的，因此 HSP 对象包含搜索工具计算的大部分统计数据。这也区分了 HSP 与您所看到的差异相比，来自不同搜索工具的对象更加明显 QueryResult 或 Hit 对象
让我们看看我们的AMPS和BLAT搜索中的一些示例。我们首先来看看AMPS：
>>> 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 ，援引 print 上 HSP 显示了其一般详细信息：
有查询和点击ID和描述。我们需要这些来识别我们的
HSP.我们还获得了查询和命中序列的匹配范围。我们在这里使用的切片符号表明该范围是使用Python的索引风格（从零开始，半开）显示的。括号内的数字表示链。在这种情况下，两个序列都有正链。
有一些快速统计数据：e值和bitscore。
有有关热休克蛋白片段的信息。暂时忽略这个;稍后会解释。
最后，我们有了查询和命中序列比对本身。
这些详细信息可以使用点符号自行访问，就像在 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
查看 HSP 下的文件 Bio.SearchIO 获取这些预定义属性的完整列表。
此外，每个序列搜索工具通常计算其自己的统计数据/详细信息 HSP 对象例如，HTML搜索还输出空位和相同残基的数量。这些属性可以像这样访问：
>>> 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']
最后，您可能已经注意到 query 和 hit 我们的SPP的属性不仅仅是常规字符串：
>>> 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 您之前在部分中看到的对象 顺序注释对象 !这意味着您可以做各种有趣的事情 SeqRecord 上的对象 HSP.query 和/或 HSP.hit .
现在你不应该感到惊讶 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|
在探索了AMPS之后，现在让我们来看看不同类型的热休克蛋白的BLAT结果中的热休克蛋白。像往常一样，我们将首先调用 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结果中存在的另一种热休克蛋白“变体”时，事情会变得更加有趣。您可能还记得，在BLAT搜索中，有时我们会将结果分成“块”。这些区块本质上是排列片段，它们之间可能有一些插入序列。
让我们来看看包含多个块的BLAT SPP，看看如何 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和描述、坐标和快速统计数据与您之前看到的相似。但碎片细节却各不相同。我们现在不再显示“Fragments：1”，而是拥有一个包含两个数据行的表。
就是这样 Bio.SearchIO 处理具有多个片段的热休克蛋白。如前所述，热休克蛋白比对可以通过插入序列分离成片段。中间序列不是查询命中匹配的一部分，因此它们不应被视为查询或命中序列的一部分。然而，它们确实会影响我们处理序列坐标的方式，所以我们不能忽视它们。
看看上面的热休克蛋白的命中坐标。在 Hit range: 字段，我们看到坐标是 [54233104:54264463] .但查看表行，我们发现该坐标跨越的整个区域并不匹配我们的查询。具体来说，中间区域跨越自 54233122 到 54264420 .
那么，为什么查询坐标似乎是连续的，您会问？这完全没问题。在这种情况下，这意味着查询匹配是连续的（没有插入区域），而命中匹配不是。
顺便说一句，所有这些属性都可以直接从SPP访问：
>>> 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 属性？如果SPP有多个片段，您将无法使用这些属性，因为它们只能获取单个片段 SeqRecord 或 MultipleSeqAlignment 对象但是，您可以使用他们的 *_all 同行： query_all , hit_all ，而且 aln_all .这些属性将返回一个列表，其中包含 SeqRecord 或 MultipleSeqAlignment 来自每个热休克蛋白片段的对象。还有其他属性的行为类似，即它们仅适用于具有一个片段的HSPs。查看 HSP 下的文件 Bio.SearchIO 查看完整列表。
最后，要检查是否有多个片段，可以使用 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 直接对象，因为没有那么多序列搜索工具片段化它们的热休克蛋白。当你确实必须与他们打交道时，你应该记住的是 HSPFragment 对象的编写是为了尽可能紧凑。在大多数情况下，它们只包含与序列直接相关的属性：链、阅读框、分子类型、坐标、序列本身以及它们的ID和描述。
当您调用 print 上 HSPFragment .这是一个例子，取自我们的AMPS搜索：
>>> 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 使用.如果您使用过多个序列搜索工具，那么您可能必须处理每个程序处理序列坐标等事情的许多不同方式。这可能不是一次愉快的体验，因为这些搜索工具通常有自己的标准。例如，一种工具可能使用以一为基的坐标，而另一种工具可能使用以零为基的坐标。或者，如果链为负，一个程序可能会颠倒开始和结束坐标，而其他程序则不会。简而言之，这些往往会造成不必要的混乱，必须处理。
我们自己意识到这个问题，并打算在 Bio.SearchIO .毕竟，目标之一 Bio.SearchIO 是创建一个通用的、易于使用的界面来处理各种搜索输出文件。这意味着创建超出您刚刚看到的对象模型的标准。
现在，您可能会抱怨：“不是另一个标准！".好吧，最终我们必须选择一个惯例或另一个惯例，所以这是必要的。另外，我们并不是在这里创建全新的东西;只是采用我们认为最适合Python程序员的标准（毕竟是Biopython）。
在工作时，您可以期待三个隐含标准 Bio.SearchIO :
第一个涉及序列坐标。在
Bio.SearchIO，所有序列坐标都遵循Python的坐标风格：从零开始、半开。例如，如果在BLASTML输出文件中，SPP的开始和结束坐标是10和28，那么它们将在中变成9和28Bio.SearchIO.开始坐标变为9，因为Python索引从零开始，而结束坐标保持28，因为Python切片在一定时间内省略了最后一项。第二个是序列坐标顺序。在
Bio.SearchIO，开始坐标始终小于或等于结束坐标。并非所有序列搜索工具都是如此，因为当序列链为负时，其中一些工具具有更大的开始坐标。最后一个是在链上并阅读框架值。对于股，只有四个有效选择：
1(plus股），-1（负链），0（蛋白质序列），以及None(no strand）。对于阅读框架，有效选择是从-3到3和None.
请注意，这些标准仅存在于 Bio.SearchIO 对象如果你写 Bio.SearchIO 对象转换为输出格式， Bio.SearchIO 将使用该格式的标准进行输出。它不会将其标准强加给您的输出文件。
读取搜索输出文件
您可以使用两个功能将搜索输出文件读取到 Bio.SearchIO 对象： read 和 parse .它们本质上类似于 read 和 parse 其他子模块中的功能，例如 Bio.SeqIO 或 Bio.AlignIO .在这两种情况下，您都需要提供搜索输出文件名和文件格式名称，两者都作为Python字符串。您可以检查文档中的格式名称列表 Bio.SearchIO 认识到。
Bio.SearchIO.read 用于仅使用一个查询读取搜索输出文件并返回 QueryResult object.你见过 read 在我们之前的例子中使用过。你没有看到的是 read 也可以接受附加的关键字参数，这取决于文件格式。
这里有一些例子。在第一个中，我们使用 read 就像以前一样读取AMPS表格输出文件。在第二个中，我们使用关键字参数进行修改，以便它解析包含注释的BST表格变体：
>>> 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 对象如果这些属性不存在，写作就行不通。换句话说，您不能总是写入您想要的输出格式。例如，如果您读取了一个RST的HTML文件，则无法将结果写入PSL文件，因为PSL文件需要未由RST计算的属性（例如重复匹配的数量）。不过，如果您真的想写入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 ，它仅限于具有所有所需属性的格式转换。在这里，AMPS HTML文件提供了BST表格文件所需的所有默认值，因此它工作得很好。然而，其他格式转换不太可能起作用，因为您需要首先手动分配所需的属性。