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 結果是一個 XML 檔案，使用 blastn 針對 NCBI refseq_rna 資料庫產生。對於 BLAT，序列資料庫是 2009 年 2 月的 hg19 人類基因組草案，而輸出格式為 PSL。
我們將從介紹 Bio.SearchIO 物件模型開始。該模型是您的搜尋結果的表示，因此它是 Bio.SearchIO 本身的核心。之後，我們將檢查您可能經常使用的 Bio.SearchIO 中的主要函數。
現在我們都準備好了，讓我們開始第一步：介紹核心物件模型。
SearchIO 物件模型
儘管許多序列搜尋工具的輸出樣式差異很大，但事實證明它們的基本概念是相似的
輸出檔案可能包含一個或多個搜尋查詢的結果。
在每個搜尋查詢中，您會看到來自給定搜尋資料庫的一個或多個命中。
在每個資料庫命中中，您會看到一個或多個區域，其中包含您的查詢序列與資料庫序列之間的實際序列比對。
像 BLAT 或 Exonerate 這樣的一些程式可能會進一步將這些區域分割成多個比對片段（或 BLAT 中的區塊，以及 exonerate 中可能的 exons）。這並不是您總是會看到的，因為像 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 有多少個項目 (hits)，您可以直接調用 Python 的 len 方法。
>>> len(blast_qresult)
100
>>> len(blat_qresult)
1
如同 Python 列表，您可以使用切片表示法從 QueryResult 中擷取項目 (hits)。
>>> 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)
若要擷取多個 hits，您也可以使用切片表示法對 QueryResult 物件進行切片。在此情況下，切片將傳回一個新的 QueryResult 物件，其中僅包含切片的 hits。
>>> 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 字典，您也可以使用 hit 的 ID 來擷取 hits。如果您知道給定的 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 關鍵字執行簡單的 Python 成員資格測試。
>>> "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 的索引樣式，它是從零開始的。這表示我們上面的 hit 排名第 23，而不是第 22。
另外，請注意您在此看到的 hit 排名是基於原始搜尋輸出檔案中存在的原生 hit 順序。不同的搜尋工具可能會根據不同的標準對這些 hit 進行排序。
如果原生 hit 順序不符合您的偏好，您可以使用 QueryResult 物件的 sort 方法。它與 Python 的 list.sort 方法非常相似，但額外提供了一個選項，可建立新的已排序 QueryResult 物件，或不建立。
以下範例說明如何使用 QueryResult.sort 根據每個 hit 的完整序列長度對 hits 進行排序。對於此特定排序，我們會將 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 物件中的所有 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("%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 物件，而是對每個 hit 中的 HSP 物件執行篩選。
至於 map 方法，它們也接受回呼函式作為其引數。但是，回呼函式不是傳回 True 或 False，而是必須傳回修改過的 Hit 或 HSP 物件 (取決於您使用的是 hit_map 還是 hsp_map)。
讓我們看看一個範例，其中我們使用 hit_map 來重新命名 hit 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 和描述都存在。hit 始終與查詢相關聯，因此我們也想追蹤原始查詢。可以使用
query_id和query_description屬性從 hit 存取這些值。我們也有唯一的 hit ID、描述和完整序列長度。可以使用
id、description和seq_len分別存取它們。最後，還有一個表格，其中包含此 hit 所含 HSP 的快速資訊。在每一列中，我們列出了重要的 HSP 詳細資訊：HSP 索引、其 e 值、其位元分數、其跨距 (包含間隙的比對長度)、其查詢座標和其 hit 座標。
現在，讓我們將其與 BLAT 搜尋進行比較。請記住，在 BLAT 搜尋中，我們有一個包含 17 個 HSP 的 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 值和位元分數欄的值。由於 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]
您也可以使用與在 QueryResult 物件中看到的排序方法完全相同的參數，來排序 Hit 內的 HSP。
最後，您還可以在 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 值和位元分數。
有關於 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 物件！這表示您可以使用 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|
在探討完 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 檔案直接取得，但是當您剖析 PSL 檔案時，Bio.SearchIO 會即時為您計算它們。它只需要每個片段的開始和結束座標。
那麼 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 物件上使用切片 (slice) 符號，就像 QueryResult 或 Hit 物件一樣。當您使用這種符號時，您將會得到一個 HSPFragment 物件作為回傳值，它是物件模型的最後一個組件。
HSPFragment
HSPFragment 代表查詢序列和命中序列之間單一、連續的匹配。您可以將其視為物件模型和搜尋結果的核心，因為這些片段的存在與否決定了您的搜尋是否有結果。
在大多數情況下，您不必直接處理 HSPFragment 物件，因為沒有那麼多序列搜尋工具會將它們的 HSP 分割成片段。當您必須處理它們時，您應該記住 HSPFragment 物件的設計盡可能地精簡。在大多數情況下，它們只包含直接與序列相關的屬性：股 (strands)、讀碼框架 (reading frames)、分子類型、座標、序列本身，以及它們的 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 索引從零開始，而結束座標保持 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 表格檔案所需的所有預設值，因此它可以正常運作。然而，其他格式轉換不太可能運作，因為您需要先手動指定必要的屬性。