使用BLAST通常可以分成2个步。这两步都可以用上Biopython。第一步，提交你的查询 序列，运行BLAST，并得到输出数据。第二步，用Python解析BLAST的输出，并作进一步 分析。
通过Internet运行BLAST
使用 Bio.Blast.NCBIWWW 模块的里qblast() 来调用在线版本的BLAST。有三个参数：
（１）用来搜索blast程序。目前只支持：blastn, blastp, blastx, tblast 和 tblastx.（用过NCBI里blast的同学一定不会对这些程序陌生的）
（２）指定要搜索的数据库
（３）你要查询序列的字符串。这个字符串可以是序列的本身 ，后者fasta格式的的文件，或者是一个类似GI的id。
qblast 函数可以返回多种格式的BLAST结果。你可以用format_type 来指定 "HTML", "Text", "ASN.1", 或者"XML"格式。默认是"XML"。
举个例子，如果你有一个fasta文件（你可以随便下载一个你感兴趣的基因序列，关于如何得到基因的fasta文件请参照：https://zhuanlan.zhihu.com/p/31618808），你想用BLASTN数据库进行比对，你可以这样：
>>> 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)
上面最重要的代码就是最后一行的比对代码。括号里是前面提到的最重要的三个参数，你也可以加其他的过滤条件。都有哪些可以加的过滤条件，可以通过看帮助文档：
>>> from Bio.Blast import NCBIWWW
>>> help(NCBIWWW.qblast)
另外，括号里中间的搜索数据库，也是非常重要的。而有哪些数据库可以给你搜索，可以看一下NCBI里的blast，我截了张图：
上面我用的数据库就是stardard database(nr)里的，其中nr里有14个小的分类，可以用上图下拉菜单里寻找。旁边的问号可以告诉你选择的数据库都是什么。选好了数据库，就可以替换上面的代码里的第二个参数了。
上面的方法只提供序列，意味着BLAST会自动分配给你一个ID。或者你还可以用 SeqRecord 对象的format方法来包装一个fasta字符串，这个对象会包含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"))
然后电脑会运行一小会儿，比如我这个破电脑就差不多得两三分钟的样子。运行后它不会弹出来任何结果，这个运行的结果是需要你手动保存的：
>>> 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")
#这些做好后，结果已经存储在 `my_blast.xml` 文件中，并且原先的handle中的数据 已经被全部提取出来了(所以我们把它关闭了)。但是，BLAST解析器的 `parse` 函数 采用一个文件句柄类的对象，所以我们只需打开已经保存的文件作为输入。
在本地运行
在本地运行BLAST有2个主要优点：
（１）本地运行BLAST可能比通过internet运行更快；
（２）本地运行可以让你用自己的数据库来对序列进行blast。
但是本地运行也有些缺点 ：
（１）本地运行BLAST需要你安装相关命令行工具。
（２）本地运行BLAST需要安装一个很大的BLAST的数据库（并且需要保持数据更新）。
这里我就不展开了，对于我来说还不太需要。有兴趣的朋友可以具体看一下这一部分。
怎么看BLAST结果的输出
刚才在上面做了一个blast试了试，结果输出一个xml文件，打开以后是这样的，还没仔细看我就已经晕菜了：
需要注意的是：上面你保存完你的xml文件，一定要有最后一行的代码，才可以分析你的blast结果：
>>> result_handle = open("my_blast.xml")
如果你的比对只有一个，用下面的代码进行解析：
>>> from Bio.Blast import NCBIXML
>>> blast_record = NCBIXML.read(result_handle)
如果你的比对多个，则用：
>>> from Bio.Blast import NCBIXML
>>> blast_records = NCBIXML.parse(result_handle)
但是当你的比对成百上千条的时候怎么办呢？用循环吧：
>>> for blast_record in blast_records:
... # Do something with blast_record
因为我只比对了一条序列，所以我用的上面第一种方法进行解析。
那么解析以后怎么能看到数据库里哪些序列和我的序列比对上了呢？看下面：
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] + "...")
但是我用了这个代码运行后，顿时后悔了。弹出了不知道多少条的比对结果。因为这个比对不仅仅是把高匹配的序列列出来，就连不怎么相似的序列也列出来了。而且，是所有物种里和你的序列相似的基因！自己可以感受一下，我就不贴运行之后的结果图了。
那么你也许会说，怎么看那些匹配度特别高的呢？你可以通过改变 E_VALUE_THRESH这个值来筛选。一般如果匹配度很高，基本上这个值会是：××e-**这种形式的。也就是说0.000000几（省略好多个０）。比如我第二次尝试，试了９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和其他序列搜索工具
随着序列 数量的增加（匹配数也会随之增加），将会有成百上千的可能匹配，解析这些结果 无疑变得越来越困难。理所当然，人工解析搜索结果就非常的困难了Biopython里有个模块叫Bio.SearchIO。Bio.SearchIO 模块使从搜索结果中提取信息变得简单，并且可以处理不同工具的不同标准和规则。
（一）SearchIO对象模块
SearchIO模块里包含很多个对象，这些对象是嵌套分级的。怎么理解呢？知道俄罗斯套娃么？最外面的是一个大盒子，盒子的名字叫searchIO，里面有一个小一点的盒子，叫QueryResult。这个小一点的盒子里又有一个更小的盒子，叫Hit。Hit盒子里是HSP。不过一个Hit盒子里可以有好多个HSP。HSP盒子里的是HSPFragment盒子。（所以有四个对象）
Bio.SearchIO有四个方法，分别是：read, parse, index, 和 index_db。read 用于单query对输出文件进行搜索并且返回一个 QueryResult 对象。parse 用于多query对输出文件进行搜索并且返回一个可以产出 QueryResult 对象的输出器。
下面来具体的看一下每一个对象：
（１）QueryResult
QueryResult，代表单query搜索，每个 QueryResult 中有0个或多个 Hit 对象。现在看看上面我们得到的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 #搜索的目标数据库
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description #hit结果的快速预览。
---- ----- ----------------------------------------------------------
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 截断了表格，只显示0-29，然后是最后３个。
那么QueryResult 到底是什么？QueryResult是混合了列表和字典的特性。就是我上面说的一个容器，一个盒子。QueryResult对象是可迭代的（iterable）。每一次迭代返回一个Hit对象：
>>> 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 对象有多少hit：
>>> len(blast_qresult)
50
可以用切片来获得 QueryResult 对象的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)
#查看最后一个hit，索引为－１。
>>> blast_qresult[-1]
Hit(id='gi|1353793171|gb|CP027081.1|', query_id='gi|568815578|ref|NC_000020.11|:49982980-49988886', 99 hsps)
想要一次看几个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...
如果你知道一个特定的hit ID存在于一个搜索结果中时，特别有用：
>>> blast_qresult["gi|**********|ref|***********|"]
你可以用hits方法获得完整的Hit对象，也可以用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|', ......
想确定一个特定的hit是否存在于查询结果中该怎么做?
>>> "gi|262205317|ref|NR_030195.1|" in blast_qresult
False
有时候你会想知道某一个hit的排名（请注意这个排名是python的索引数字）：
>>> blast_qresult.index("gi|**********|ref|***********|")
22
如果原本的hit排序不合你意，可以用 QueryResult 对象的sort方法。这个方法基于每个hit 的完整序列长度。对于这个排序，设置in_place参数为False ，这样排序之后会返回一个新的QueryResult 对象，而原来的对象是未排序的。同样可以设置 reverse 参数为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
对于QueryResult对象有两种更方便的方法：filter 和 map 方法。filter方法有两种：hit_filter 和 hsp_filter，分别过滤 QueryResult 对象的Hit对象或者HSP对象。同样对于 map 也有两种方法：hit_map 和 hsp_map。下面看一下具体的使用：
从hit_filter开始。用 hit_filter筛选出多于10个HSP的Hit 对象：
>>> filter_func = lambda hit: len(hit.hsps) > 10
>>> len(blast_qresult) #过滤前的hit数量
50
>>> filtered_qresult = blast_qresult.hit_filter(filter_func)
>>> len(filtered_qresult) #过滤后的hit数量
5
hsp_filter和hit_filter功能一样，只不过它过滤每一个hit中的HSP对象，而不是Hit 。
对于map方法，返回的是修改过的Hit或HSP对象（取决于你是否使用 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对象而不是Hit对象。
（２） Hit对象
Hit对象代表从单个数据库获得所有查询结果。在Bio.SearchIO模块等级中是二级容器。被包含在 QueryResult 对象中。
>>> 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 对象也是可迭代的，每次迭代返回一个HSP对象：
>>> 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)
你可以对Hit对象调用 len 方法查看它含有多少个HSP对象：
>>> len(blast_hit)
1
你可以对Hit对象作切片取得单个或多个HSP对象，和QueryResult一样，如果切取多个 HSP，会返回包含被切片 HSP的一个新Hit对象。因为我这里的hit只有一个HSP对象，所以len里才会显示０：
>>> 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对象
HSP (高分片段)代表hit序列中的一个区域，该区域包含对于查询序列有意义的比对。它包含了你的查询序列和一个数据库条目之间精确的匹配。由于匹配取决于序列搜索工具的算法， HSP含有大部分统计信息，这些统计是由搜索工具计算得到的。这使得不同搜索工具的HSP对象之间的差异和你在QueryResult以及Hit对象看到的差异尤其明显。
>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast2.xml', 'blast-xml')
>>> blast_hsp = blast_qresult[0][0] #查看第一个 hit，第一个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
#查看比对的区间
>>> blast_hsp.query_range
(0, 5907)
#查看e值，因为我比对的序列本来就是一个基因序列，所以hsp的e值是０，是完全比对上的
>>> blast_hsp.evalue
0.0
当然，你还可以查看的属性有：
#这段hsp在hit上的起始位置在哪里
>>> blast_hsp.hit_start
48539516
#这段hsp的跨度（长度、多少个碱基）
>>> blast_hsp.query_span
5907
#这段hsp比对上的跨度
>>> blast_hsp.aln_span
5907
你会发现，HSP里的query属性和hit属性不只是规律字符串：
>>> 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之间单个连续匹配。其实我们应该把它当作搜索结果的核心，因为它决定你的搜索是否有结果。在多数情况下，它们仅仅包含直接与序列有关的属性：正负链， 阅读框，字母表，位置坐标，序列本身以及它们的ID和描述。
>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast2.xml', 'blast-xml')
>>> blast_frag = blast_qresult[0][0][0] #第一个hit里的第一个hsp，第一个hsp里的第一个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
#查询这个frag比对的起始位置
>>> blast_frag.query_start
0
#这个frag在hit上哪一条链
>>> blast_frag.hit_strand
1
#所在的hit序列，是一个SeqRecord对象
>>> 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的惯例
（１）所有序列坐标遵循Python 的坐标风格： 从0开始
（２）对于链值，只有四个可选值：1 (正链)，-1 (负链)，0 (蛋白序列)，和 None (无链)。对于阅读框， 可选值是从-3～3的整型以及 None 。
写入和转换搜索输出文件
最后要说的是：你可以把xml文件转换成其他文件。目前只支持如下文件格式："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)
网友评论