Biopython學(xué)習(xí)筆記（三）Blast

使用BLAST通辰徊ィ可以分成2個(gè)步。這兩步都可以用上Biopython践付。第一步秦士，提交你的查詢 序列，運(yùn)行BLAST永高，并得到輸出數(shù)據(jù)隧土。第二步提针，用Python解析BLAST的輸出，并作進(jìn)一步 分析曹傀。

通過Internet運(yùn)行BLAST

使用 Bio.Blast.NCBIWWW 模塊的里qblast() 來調(diào)用在線版本的BLAST辐脖。有三個(gè)參數(shù)：
（１）用來搜索blast程序。目前只支持：blastn, blastp, blastx, tblast 和 tblastx.（用過NCBI里blast的同學(xué)一定不會(huì)對(duì)這些程序陌生的）
（２）指定要搜索的數(shù)據(jù)庫
（３）你要查詢序列的字符串皆愉。這個(gè)字符串可以是序列的本身 嗜价，后者fasta格式的的文件，或者是一個(gè)類似GI的id幕庐。

qblast 函數(shù)可以返回多種格式的BLAST結(jié)果久锥。你可以用format_type 來指定 "HTML", "Text", "ASN.1", 或者"XML"格式。默認(rèn)是"XML"异剥。

舉個(gè)例子瑟由，如果你有一個(gè)fasta文件（你可以隨便下載一個(gè)你感興趣的基因序列，關(guān)于如何得到基因的fasta文件請(qǐng)參照：https://zhuanlan.zhihu.com/p/31618808）冤寿，你想用BLASTN數(shù)據(jù)庫進(jìn)行比對(duì)歹苦，你可以這樣：

>>> from Bio.Blast import NCBIWWW
>>> from Bio import SeqIO
#這里我下載的是人類的snai1基因序列
>>> record = SeqIO.read("human_snai1.fasta", format="fasta")
>>> result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)

上面最重要的代碼就是最后一行的比對(duì)代碼。括號(hào)里是前面提到的最重要的三個(gè)參數(shù)督怜，你也可以加其他的過濾條件殴瘦。都有哪些可以加的過濾條件，可以通過看幫助文檔：

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

另外号杠，括號(hào)里中間的搜索數(shù)據(jù)庫蚪腋，也是非常重要的。而有哪些數(shù)據(jù)庫可以給你搜索究流，可以看一下NCBI里的blast辣吃，我截了張圖：

上面我用的數(shù)據(jù)庫就是stardard database(nr)里的动遭，其中nr里有14個(gè)小的分類芬探，可以用上圖下拉菜單里尋找。旁邊的問號(hào)可以告訴你選擇的數(shù)據(jù)庫都是什么厘惦。選好了數(shù)據(jù)庫偷仿，就可以替換上面的代碼里的第二個(gè)參數(shù)了。

上面的方法只提供序列宵蕉，意味著BLAST會(huì)自動(dòng)分配給你一個(gè)ID酝静。或者你還可以用 SeqRecord 對(duì)象的format方法來包裝一個(gè)fasta字符串羡玛，這個(gè)對(duì)象會(huì)包含fasta文件中已有的ID：

>>> from Bio.Blast import NCBIWWW
>>> from Bio import SeqIO
>>> record = SeqIO.read("human_snai1.fasta", format="fasta")
>>> result_handle = NCBIWWW.qblast("blastn", "nt", record.format("fasta"))

然后電腦會(huì)運(yùn)行一小會(huì)兒别智，比如我這個(gè)破電腦就差不多得兩三分鐘的樣子。運(yùn)行后它不會(huì)彈出來任何結(jié)果稼稿，這個(gè)運(yùn)行的結(jié)果是需要你手動(dòng)保存的：

>>> with open("my_blast.xml", "w") as out_handle:
...     out_handle.write(result_handle.read())
... 
4184652
>>> result_handle.close()
>>> result_handle = open("my_blast.xml")
#這些做好后薄榛，結(jié)果已經(jīng)存儲(chǔ)在 `my_blast.xml` 文件中讳窟，并且原先的handle中的數(shù)據(jù) 已經(jīng)被全部提取出來了(所以我們把它關(guān)閉了)。但是敞恋，BLAST解析器的 `parse` 函數(shù) 采用一個(gè)文件句柄類的對(duì)象丽啡，所以我們只需打開已經(jīng)保存的文件作為輸入。

在本地運(yùn)行

在本地運(yùn)行BLAST有2個(gè)主要優(yōu)點(diǎn)：
（１）本地運(yùn)行BLAST可能比通過internet運(yùn)行更快硬猫；
（２）本地運(yùn)行可以讓你用自己的數(shù)據(jù)庫來對(duì)序列進(jìn)行blast补箍。
但是本地運(yùn)行也有些缺點(diǎn) ：
（１）本地運(yùn)行BLAST需要你安裝相關(guān)命令行工具。
（２）本地運(yùn)行BLAST需要安裝一個(gè)很大的BLAST的數(shù)據(jù)庫（并且需要保持?jǐn)?shù)據(jù)更新）啸蜜。
這里我就不展開了坑雅，對(duì)于我來說還不太需要。有興趣的朋友可以具體看一下這一部分盔性。

怎么看BLAST結(jié)果的輸出

剛才在上面做了一個(gè)blast試了試霞丧，結(jié)果輸出一個(gè)xml文件，打開以后是這樣的冕香，還沒仔細(xì)看我就已經(jīng)暈菜了：

需要注意的是：上面你保存完你的xml文件蛹尝，一定要有最后一行的代碼，才可以分析你的blast結(jié)果：

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

如果你的比對(duì)只有一個(gè)悉尾，用下面的代碼進(jìn)行解析：

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

如果你的比對(duì)多個(gè)突那，則用：

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

但是當(dāng)你的比對(duì)成百上千條的時(shí)候怎么辦呢？用循環(huán)吧：

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

因?yàn)槲抑槐葘?duì)了一條序列构眯，所以我用的上面第一種方法進(jìn)行解析愕难。
那么解析以后怎么能看到數(shù)據(jù)庫里哪些序列和我的序列比對(duì)上了呢？看下面：

BLAST 記錄類

教材里給的代碼是：

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

但是我用了這個(gè)代碼運(yùn)行后惫霸，頓時(shí)后悔了猫缭。彈出了不知道多少條的比對(duì)結(jié)果。因?yàn)檫@個(gè)比對(duì)不僅僅是把高匹配的序列列出來壹店，就連不怎么相似的序列也列出來了猜丹。而且，是所有物種里和你的序列相似的基因硅卢！自己可以感受一下射窒，我就不貼運(yùn)行之后的結(jié)果圖了。

那么你也許會(huì)說将塑，怎么看那些匹配度特別高的呢脉顿？你可以通過改變 E_VALUE_THRESH這個(gè)值來篩選。一般如果匹配度很高点寥，基本上這個(gè)值會(huì)是：××e-**這種形式的艾疟。也就是說0.000000幾（省略好多個(gè)０）。比如我第二次嘗試，試了９e-100蔽莱。這次終于可以看到最前面幾條的記錄了：

****Alignment****
sequence: gi|1548994295|gb|CP034499.1| Eukaryotic synthetic construct chromosome 20 >gi|1549098600|gb|CP034524.1| Eukaryotic synthetic construct chromosome 20
length: 68480253
e value: 0.0
ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGCGACCCCAGTGCCTCGA...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||...
ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGCGACCCCAGTGCCTCGA...
****Alignment****
sequence: gi|9650757|emb|AL121712.27| Human DNA sequence from clone RP4-710H13 on chromosome 20, complete sequence
length: 79779
e value: 0.0
ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGCGACCCCAGTGCCTCGA...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||...
ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGCGACCCCAGTGCCTCGA...
****Alignment****
sequence: gi|5821735|gb|AF155233.1|AF155233 Homo sapiens snail zinc finger protein (SNAI1) gene, complete cds
length: 6258
e value: 0.0
ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGCGACCCCAGTGCCTCGA...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||...
ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGCGACCCCAGTGCCTCGA...
****Alignment****
sequence: gi|5725247|emb|AJ245659.1| Homo sapiens partial SNAI1 gene, exon 3
length: 1060
e value: 0.0
CCAGCCGTTGTCCCACGGCTCACTCGGCCTTTCTGGCGTTCTCTCCCCAGGCGAGAAGCCCTTCTCCTGTCCCCA...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||...
CCAGCCGTTGTCCCACGGCTCACTCGGCCTTTCTGGCGTTCTCTCCCCAGGCGAGAAGCCCTTCTCCTGTCCCCA...

這里我專門查了一下fasta格式的具體說明误褪，如果熟悉這部分的可以直接略過了。參考：https://blog.csdn.net/u011262253/article/details/51164756碾褂。

BLAST和其他序列搜索工具

隨著序列 數(shù)量的增加（匹配數(shù)也會(huì)隨之增加）兽间，將會(huì)有成百上千的可能匹配，解析這些結(jié)果 無疑變得越來越困難正塌。理所當(dāng)然嘀略，人工解析搜索結(jié)果就非常的困難了Biopython里有個(gè)模塊叫Bio.SearchIO。Bio.SearchIO 模塊使從搜索結(jié)果中提取信息變得簡單乓诽，并且可以處理不同工具的不同標(biāo)準(zhǔn)和規(guī)則帜羊。

（一）SearchIO對(duì)象模塊

SearchIO模塊里包含很多個(gè)對(duì)象，這些對(duì)象是嵌套分級(jí)的鸠天。怎么理解呢讼育？知道俄羅斯套娃么？最外面的是一個(gè)大盒子稠集，盒子的名字叫searchIO奶段，里面有一個(gè)小一點(diǎn)的盒子，叫QueryResult剥纷。這個(gè)小一點(diǎn)的盒子里又有一個(gè)更小的盒子痹籍，叫Hit。Hit盒子里是HSP晦鞋。不過一個(gè)Hit盒子里可以有好多個(gè)HSP蹲缠。HSP盒子里的是HSPFragment盒子。（所以有四個(gè)對(duì)象）
Bio.SearchIO有四個(gè)方法悠垛，分別是：read, parse, index, 和 index_db线定。read 用于單query對(duì)輸出文件進(jìn)行搜索并且返回一個(gè) QueryResult 對(duì)象。parse 用于多query對(duì)輸出文件進(jìn)行搜索并且返回一個(gè)可以產(chǎn)出 QueryResult 對(duì)象的輸出器确买。

下面來具體的看一下每一個(gè)對(duì)象：

（１）QueryResult

QueryResult斤讥，代表單query搜索，每個(gè) QueryResult 中有0個(gè)或多個(gè) Hit 對(duì)象∧赐铮現(xiàn)在看看上面我們得到的BLAST文件是什么樣的：

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read("my_blast2.xml", "blast-xml")
>>> print(blast_qresult)
Program: blastn (2.10.0+) #程序的名稱和版本
  Query: gi|568815578|ref|NC_000020.11|:49982980-49988886 (5907)
         Homo sapiens chromosome 20, GRCh38.p13 Primary Assembly #query的ID周偎，描述和序列長度
 Target: nt #搜索的目標(biāo)數(shù)據(jù)庫
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description          #hit結(jié)果的快速預(yù)覽抹剩。
         ----  -----  ----------------------------------------------------------
            0    976  gi|1548994295|gb|CP034499.1|  Eukaryotic synthetic cons...
            1     63  gi|9650757|emb|AL121712.27|  Human DNA sequence from cl...
            2      1  gi|5821735|gb|AF155233.1|AF155233  Homo sapiens snail z...
            3      1  gi|5725247|emb|AJ245659.1|  Homo sapiens partial SNAI1 ...
            4      3  gi|1519243938|ref|NM_005985.4|  Homo sapiens snail fami...
            5      3  gi|15277716|gb|BC012910.1|  Homo sapiens snail homolog ...
            6      3  gi|1753017152|ref|XM_004062349.3|  PREDICTED: Gorilla g...
            7      3  gi|1367236023|ref|XM_016938127.2|  PREDICTED: Pan trogl...
            8      3  gi|1367236020|ref|XM_016938125.2|  PREDICTED: Pan trogl...
            9      3  gi|1367236018|ref|XM_009437412.3|  PREDICTED: Pan trogl...
           10      3  gi|675701575|ref|XM_003809255.2|  PREDICTED: Pan panisc...
           11      3  gi|4325321|gb|AF125377.1|AF125377  Homo sapiens zinc fi...
           12      3  gi|1800046440|ref|XM_032142166.1|  PREDICTED: Hylobates...
           13      3  gi|1743154268|ref|XM_003252924.4|  PREDICTED: Nomascus ...
           14      3  gi|1351479674|ref|XM_002830412.3|  PREDICTED: Pongo abe...
           15      3  gi|1777284280|ref|XM_003904624.3|  PREDICTED: Papio anu...
           16      3  gi|1381483496|ref|XM_011722884.2|  PREDICTED: Macaca ne...
           17      3  gi|1059174534|ref|XM_017888131.1|  PREDICTED: Rhinopith...
           18      3  gi|795592937|ref|XM_012059057.1|  PREDICTED: Cercocebus...
           19      3  gi|1751209046|ref|XM_010383846.2|  PREDICTED: Rhinopith...
           20      3  gi|1622840906|ref|XM_001097698.4|  PREDICTED: Macaca mu...
           21      3  gi|1411087738|ref|XM_025399385.1|  PREDICTED: Theropith...
           22      3  gi|795350384|ref|XM_011928821.1|  PREDICTED: Colobus an...
           23      3  gi|635020274|ref|XM_008014341.1|  PREDICTED: Chlorocebu...
           24      3  gi|544466412|ref|XM_005569331.1|  PREDICTED: Macaca fas...
           25      3  gi|795204465|ref|XM_011993195.1|  PREDICTED: Mandrillus...
           26      3  gi|1788687212|ref|XM_023184155.2|  PREDICTED: Piliocolo...
           27      3  gi|817304701|ref|XM_012467438.1|  PREDICTED: Aotus nanc...
           28      3  gi|1044359970|ref|XM_017529925.1|  PREDICTED: Cebus cap...
           29      3  gi|1804439919|ref|XM_032256019.1|  PREDICTED: Sapajus a...
           ~~~
           47      1  gi|929047344|gb|KT583904.1|  Homo sapiens clone HsUT006...
           48      3  gi|1126231900|ref|XM_002912981.3|  PREDICTED: Ailuropod...
           49     99  gi|1353793171|gb|CP027081.1|  Bos mutus isolate yakQH1 ...

NOTE:注意撑帖， Bio.SearchIO 截?cái)嗔吮砀瘢伙@示0-29澳眷，然后是最后３個(gè)胡嘿。

那么QueryResult 到底是什么？QueryResult是混合了列表和字典的特性钳踊。就是我上面說的一個(gè)容器衷敌，一個(gè)盒子勿侯。QueryResult對(duì)象是可迭代的（iterable）。每一次迭代返回一個(gè)Hit對(duì)象：

>>> for hit in blast_qresult:
... hit
... 
Hit(id='gi|1548994295|gb|CP034499.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 976 hsps)
Hit(id='gi|9650757|emb|AL121712.27|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 63 hsps)
Hit(id='gi|5821735|gb|AF155233.1|AF155233', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 1 hsps)
Hit(id='gi|5725247|emb|AJ245659.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 1 hsps)
Hit(id='gi|1519243938|ref|NM_005985.4|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 3 hsps)
Hit(id='gi|15277716|gb|BC012910.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 3 hsps)
Hit(id='gi|1753017152|ref|XM_004062349.3|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 3 hsps)
Hit(id='gi|1367236023|ref|XM_016938127.2|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 3 hsps)
......#(此處省略n行)

要得到 QueryResult 對(duì)象有多少hit：

>>> len(blast_qresult)
50

可以用切片來獲得 QueryResult 對(duì)象的hit：（只看部分hit）

#查看最高的hit缴罗。索引為０助琐，這是之前提到的，python里的索引從０開始面氓。
>>> blast_qresult[0]  
Hit(id='gi|1548994295|gb|CP034499.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 976 hsps)
#查看最后一個(gè)hit兵钮，索引為－１。
>>> blast_qresult[-1]
Hit(id='gi|1353793171|gb|CP027081.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 99 hsps)

想要一次看幾個(gè)hit怎么辦舌界？

>>> blast_slice = blast_qresult[:3]
>>> print(blast_slice)
Program: blastn (2.10.0+)
  Query: gi|568815578|ref|NC_000020.11|:49982980-49988886 (5907)
         Homo sapiens chromosome 20, GRCh38.p13 Primary Assembly
 Target: nt
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description
         ----  -----  ----------------------------------------------------------
            0    976  gi|1548994295|gb|CP034499.1|  Eukaryotic synthetic cons...
            1     63  gi|9650757|emb|AL121712.27|  Human DNA sequence from cl...
            2      1  gi|5821735|gb|AF155233.1|AF155233  Homo sapiens snail z...

如果你知道一個(gè)特定的hit ID存在于一個(gè)搜索結(jié)果中時(shí)掘譬，特別有用：

>>> blast_qresult["gi|**********|ref|***********|"]

你可以用hits方法獲得完整的Hit對(duì)象，也可以用hit_keys方法獲得完整的Hit ID：

>>> blast_qresult.hits
[Hit(id='gi|1548994295|gb|CP034499.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 976 hsps), Hit(id='gi|9650757|emb|AL121712.27|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 63 hsps), Hit(id='gi|5821735|gb|AF155233.1|AF155233', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 1 hsps), Hit(id='gi|5725247|emb|AJ245659.1|', ......

>>> blast_qresult.hit_keys
['gi|1548994295|gb|CP034499.1|', 'gi|9650757|emb|AL121712.27|', 'gi|5821735|gb|AF155233.1|AF155233', 'gi|5725247|emb|AJ245659.1|', 'gi|1519243938|ref|NM_005985.4|', 'gi|15277716|gb|BC012910.1|', 'gi|1753017152|ref|XM_004062349.3|', ......

想確定一個(gè)特定的hit是否存在于查詢結(jié)果中該怎么做?

>>> "gi|262205317|ref|NR_030195.1|" in blast_qresult
False

有時(shí)候你會(huì)想知道某一個(gè)hit的排名（請(qǐng)注意這個(gè)排名是python的索引數(shù)字）：

>>> blast_qresult.index("gi|**********|ref|***********|")
22

如果原本的hit排序不合你意呻拌，可以用 QueryResult 對(duì)象的sort方法葱轩。這個(gè)方法基于每個(gè)hit 的完整序列長度。對(duì)于這個(gè)排序藐握，設(shè)置in_place參數(shù)為False 靴拱，這樣排序之后會(huì)返回一個(gè)新的QueryResult 對(duì)象，而原來的對(duì)象是未排序的猾普。同樣可以設(shè)置 reverse 參數(shù)為True以遞減排序：

>>> for hit in blast_qresult[:5]:
...     print("%s %i" % (hit.id, hit.seq_len))
... 
gi|1548994295|gb|CP034499.1| 68480253
gi|9650757|emb|AL121712.27| 79779
gi|5821735|gb|AF155233.1|AF155233 6258
gi|5725247|emb|AJ245659.1| 1060
gi|1519243938|ref|NM_005985.4| 1705
>>> 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|1353793171|gb|CP027081.1| 75983851
gi|1548994295|gb|CP034499.1| 68480253
gi|1051792985|ref|NG_051361.1| 234376
gi|4753226|gb|AC006385.3| 173508
gi|9650757|emb|AL121712.27| 79779

對(duì)于QueryResult對(duì)象有兩種更方便的方法：filter 和 map 方法缭嫡。filter方法有兩種：hit_filter 和 hsp_filter，分別過濾 QueryResult 對(duì)象的Hit對(duì)象或者HSP對(duì)象抬闷。同樣對(duì)于 map 也有兩種方法：hit_map 和 hsp_map妇蛀。下面看一下具體的使用：

從hit_filter開始。用 hit_filter篩選出多于10個(gè)HSP的Hit 對(duì)象：

>>> filter_func = lambda hit: len(hit.hsps) > 10
>>> len(blast_qresult) #過濾前的hit數(shù)量
50
>>> filtered_qresult = blast_qresult.hit_filter(filter_func)
>>> len(filtered_qresult) #過濾后的hit數(shù)量
5

hsp_filter和hit_filter功能一樣笤成，只不過它過濾每一個(gè)hit中的HSP對(duì)象评架，而不是Hit 。

對(duì)于map方法炕泳，返回的是修改過的Hit或HSP對(duì)象（取決于你是否使用 hit_map 或 hsp_map 方法）纵诞。

>>> def map_func(hit):
...     hit.id = hit.id.split("|")[3] # 把'gi|＠＠＠＠＠＠|ref|×××××|' 格式改為ref后面的 '×××××'
...     return hit
... 
>>> mapped_qresult = blast_qresult.hit_map(map_func)
>>> for hit in mapped_qresult[:5]:
...     print(hit.id)
... 
CP034499.1
AL121712.27
AF155233.1
AJ245659.1
NM_005985.4

hsp_map和hit_map作用相似, 但是作用于HSP對(duì)象而不是Hit對(duì)象。

（２） Hit對(duì)象

Hit對(duì)象代表從單個(gè)數(shù)據(jù)庫獲得所有查詢結(jié)果培遵。在Bio.SearchIO模塊等級(jí)中是二級(jí)容器浙芙。被包含在 QueryResult 對(duì)象中。

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read("my_blast2.xml", "blast-xml")
>>> blast_hit = blast_qresult[3] #查看索引為３的 hit
>>> print(blast_hit)
Query: gi|568815578|ref|NC_000020.11|:49982980-49988886
       Homo sapiens chromosome 20, GRCh38.p13 Primary Assembly
  Hit: gi|5725247|emb|AJ245659.1| (1060)
       Homo sapiens partial SNAI1 gene, exon 3
 HSPs: ----  --------  ---------  ------  ---------------  ---------------------
          #   E-value  Bit score    Span      Query range              Hit range
       ----  --------  ---------  ------  ---------------  ---------------------
          0         0    1912.86    1060      [4842:5902]               [0:1060]

Hit 對(duì)象也是可迭代的籽腕，每次迭代返回一個(gè)HSP對(duì)象：

>>> for hsp in blast_hit:
...     hsp
... 
HSP(hit_id='gi|5725247|emb|AJ245659.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 1 fragments)

你可以對(duì)Hit對(duì)象調(diào)用 len 方法查看它含有多少個(gè)HSP對(duì)象：

>>> len(blast_hit)
1

你可以對(duì)Hit對(duì)象作切片取得單個(gè)或多個(gè)HSP對(duì)象嗡呼，和QueryResult一樣，如果切取多個(gè) HSP皇耗，會(huì)返回包含被切片 HSP的一個(gè)新Hit對(duì)象南窗。因?yàn)槲疫@里的hit只有一個(gè)HSP對(duì)象，所以len里才會(huì)顯示０：

>>> blast_hit[0]
>>> sliced_hit = blast_hit[4:9]
>>> len(sliced_hit)
0
>>> print(sliced_hit)
Query: None
       Homo sapiens chromosome 20, GRCh38.p13 Primary Assembly
  Hit: None (1060)
       Homo sapiens partial SNAI1 gene, exon 3
 HSPs: ?
（３）HSP對(duì)象

HSP (高分片段)代表hit序列中的一個(gè)區(qū)域，該區(qū)域包含對(duì)于查詢序列有意義的比對(duì)万伤。它包含了你的查詢序列和一個(gè)數(shù)據(jù)庫條目之間精確的匹配窒悔。由于匹配取決于序列搜索工具的算法， HSP含有大部分統(tǒng)計(jì)信息敌买，這些統(tǒng)計(jì)是由搜索工具計(jì)算得到的简珠。這使得不同搜索工具的HSP對(duì)象之間的差異和你在QueryResult以及Hit對(duì)象看到的差異尤其明顯。

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast2.xml', 'blast-xml')
>>> blast_hsp = blast_qresult[0][0] #查看第一個(gè) hit虹钮，第一個(gè)HSP
>>> print(blast_hsp)
      Query: gi|568815578|ref|NC_000020.11|:49982980-49988886 Homo sapiens ch...
        Hit: gi|1548994295|gb|CP034499.1| Eukaryotic synthetic construct chro...
Query range: [0:5907] (1)
  Hit range: [48539516:48545423] (1)
Quick stats: evalue 0; bitscore 10653.80
  Fragments: 1 (5907 columns)
     Query - ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGC~~~AAAAA
             |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||~~~|||||
       Hit - ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGC~~~AAAAA

#查看比對(duì)的區(qū)間
>>> blast_hsp.query_range
(0, 5907)

#查看e值北救，因?yàn)槲冶葘?duì)的序列本來就是一個(gè)基因序列，所以hsp的e值是０芜抒，是完全比對(duì)上的
>>> blast_hsp.evalue
0.0

當(dāng)然珍策，你還可以查看的屬性有：

#這段hsp在hit上的起始位置在哪里
>>> blast_hsp.hit_start
48539516
#這段hsp的跨度（長度、多少個(gè)堿基）
>>> blast_hsp.query_span
5907
#這段hsp比對(duì)上的跨度
>>> blast_hsp.aln_span 
5907

你會(huì)發(fā)現(xiàn)宅倒，HSP里的query屬性和hit屬性不只是規(guī)律字符串：

>>> blast_hsp.query
SeqRecord(seq=Seq('ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAA...AAA', DNAAlphabet()), id='gi|568815578|ref|NC_000020.11|:49982980-49988886', name='aligned query sequence', description='Homo sapiens chromosome 20, GRCh38.p13 Primary Assembly', dbxrefs=[])
>>> blast_hsp.hit
SeqRecord(seq=Seq('ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAA...AAA', DNAAlphabet()), id='gi|1548994295|gb|CP034499.1|', name='aligned hit sequence', description='Eukaryotic synthetic construct chromosome 20', dbxrefs=[])
（４） HSP片段

HSPFragment代表query和hit之間單個(gè)連續(xù)匹配攘宙。其實(shí)我們應(yīng)該把它當(dāng)作搜索結(jié)果的核心，因?yàn)樗鼪Q定你的搜索是否有結(jié)果拐迁。在多數(shù)情況下蹭劈，它們僅僅包含直接與序列有關(guān)的屬性：正負(fù)鏈， 閱讀框线召，字母表铺韧，位置坐標(biāo)，序列本身以及它們的ID和描述缓淹。

>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast2.xml', 'blast-xml')
>>> blast_frag = blast_qresult[0][0][0] #第一個(gè)hit里的第一個(gè)hsp哈打，第一個(gè)hsp里的第一個(gè)fragment
>>> print(blast_frag)
      Query: gi|568815578|ref|NC_000020.11|:49982980-49988886 Homo sapiens ch...
        Hit: gi|1548994295|gb|CP034499.1| Eukaryotic synthetic construct chro...
Query range: [0:5907] (1)
  Hit range: [48539516:48545423] (1)
  Fragments: 1 (5907 columns)
     Query - ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGC~~~AAAAA
             |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||~~~|||||
       Hit - ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAATCGGC~~~AAAAA

#查詢這個(gè)frag比對(duì)的起始位置
>>> blast_frag.query_start 
0
#這個(gè)frag在hit上哪一條鏈
>>> blast_frag.hit_strand 
1
#所在的hit序列，是一個(gè)SeqRecord對(duì)象
>>> blast_frag.hit 
SeqRecord(seq=Seq('ATTGCGCCGCGGCACGGCCTAGCGAGTGGTTCTTCTGCGCTACTGCTGCGCGAA...AAA', DNAAlphabet()), id='gi|1548994295|gb|CP034499.1|', name='aligned hit sequence', description='Eukaryotic synthetic construct chromosome 20', dbxrefs=[])
Bio.SearchIO的慣例

（１）所有序列坐標(biāo)遵循Python 的坐標(biāo)風(fēng)格： 從0開始
（２）對(duì)于鏈值讯壶，只有四個(gè)可選值：1 (正鏈)料仗，-1 (負(fù)鏈)，0 (蛋白序列)伏蚊，和 None (無鏈)立轧。對(duì)于閱讀框， 可選值是從-3～3的整型以及 None 躏吊。

寫入和轉(zhuǎn)換搜索輸出文件

最后要說的是：你可以把xml文件轉(zhuǎn)換成其他文件氛改。目前只支持如下文件格式："blast-tab', 'blast-xml', 'blat-psl', 'hmmer3-tab', 'hmmscan3-domtab', 'hmmsearch3-domtab', 'phmmer3-domtab"

>>> from Bio import SearchIO
>>> qresults = SearchIO.read('my_blast.xml', 'blast-xml') 
#tab是表格文件
>>> SearchIO.write(qresults, 'results.tab', 'blast-tab')
(1, 50, 4050, 4050)
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請(qǐng)通過簡信或評(píng)論聯(lián)系作者比伏。
  • 人面猴
    序言：七十年代末胜卤，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子凳怨，更是在濱河造成了極大的恐慌瑰艘，老刑警劉巖，帶你破解...
    沈念sama閱讀 206,126評(píng)論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件肤舞，死亡現(xiàn)場離奇詭異紫新，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)李剖，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評(píng)論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門芒率，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人篙顺，你說我怎么就攤上這事偶芍。” “怎么了德玫？”我有些...
    開封第一講書人閱讀 152,445評(píng)論 0贊 341
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵匪蟀，是天一觀的道長。 經(jīng)常有香客問我宰僧，道長材彪，這世上最難降的妖魔是什么？ 我笑而不...
    開封第一講書人閱讀 55,185評(píng)論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任琴儿，我火速辦了婚禮段化，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘造成。我一直安慰自己显熏，他們只是感情好，可當(dāng)我...
    茶點(diǎn)故事閱讀 64,178評(píng)論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布晒屎。 她就那樣靜靜地躺著喘蟆，像睡著了一般。 火紅的嫁衣襯著肌膚如雪鼓鲁。 梳的紋絲不亂的頭發(fā)上履肃，一...
    開封第一講書人閱讀 48,970評(píng)論 1贊 284
  • 城市分裂傳說
    那天，我揣著相機(jī)與錄音坐桩，去河邊找鬼尺棋。 笑死，一個(gè)胖子當(dāng)著我的面吹牛绵跷，可吹牛的內(nèi)容都是我干的膘螟。 我是一名探鬼主播，決...
    沈念sama閱讀 38,276評(píng)論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼碾局，長吁一口氣：“原來是場噩夢啊……” “哼荆残！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起净当，我...
    開封第一講書人閱讀 36,927評(píng)論 0贊 259
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對(duì)情侶失蹤内斯，失蹤者是張志新（化名）和其女友劉穎蕴潦，沒想到半個(gè)月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體俘闯，經(jīng)...
    沈念sama閱讀 43,400評(píng)論 1贊 300
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡潭苞，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,883評(píng)論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了真朗。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片此疹。...
    茶點(diǎn)故事閱讀 37,997評(píng)論 1贊 333
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡，死狀恐怖遮婶，靈堂內(nèi)的尸體忽然破棺而出蝗碎，到底是詐尸還是另有隱情，我是刑警寧澤旗扑，帶...
    沈念sama閱讀 33,646評(píng)論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布蹦骑，位于F島的核電站，受9級(jí)特大地震影響臀防，放射性物質(zhì)發(fā)生泄漏脊串。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,213評(píng)論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一清钥、第九天 我趴在偏房一處隱蔽的房頂上張望琼锋。 院中可真熱鬧，春花似錦祟昭、人聲如沸缕坎。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,204評(píng)論 0贊 19
  • 一樁弒父案篡悟，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽谜叹。三九已至，卻和暖如春搬葬，著一層夾襖步出監(jiān)牢的瞬間荷腊，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,423評(píng)論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工急凰， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留女仰，地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 45,423評(píng)論 2贊 352
  • 代替公主和親
    正文 我出身青樓抡锈，卻偏偏與公主長得像疾忍，于是被迫代替她去往敵國和親。 傳聞我的和親對(duì)象是個(gè)殘疾皇子床三，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,722評(píng)論 2贊 345

推薦閱讀更多精彩內(nèi)容

  • 歡迎關(guān)注：oddxix 本章主要講序列分析與聯(lián)配 序列分析是生物信息學(xué)最主要的研究內(nèi)容之一一罩，它可以分為兩個(gè)主要部分...
    oddxix閱讀 2,208評(píng)論 0贊 10
  • Blast，全稱Basic Local Alignment Search Tool撇簿，即"基于局部比對(duì)算法的搜索工具...
    曉僉閱讀 14,865評(píng)論 1贊 26
  • 用過網(wǎng)頁版本 BLAST 的童鞋都會(huì)發(fā)現(xiàn)，提交的序列比對(duì)往往在幾分鐘汉嗽，甚至幾十秒就可以得到比對(duì)的結(jié)果欲逃；而通過調(diào)用 ...
    BioIT愛好者閱讀 1,936評(píng)論 0贊 0
  • 雖然高通量測序分析最常用的操作是將fastq比對(duì)到參考基因組得到BAM文件，但偶爾我們也需要提取BAM文件中特定區(qū)...
    xuzhougeng閱讀 25,882評(píng)論 4贊 31
  • DISC國慶班結(jié)束已經(jīng)5天了，至今才輸出一篇關(guān)于其主題文章竭望，不免有點(diǎn)鄙視自己的懶惰邪码。同時(shí)，也有些興奮咬清，至少闭专，...
    劉紅leona閱讀 1,275評(píng)論 1贊 4