Biopython学习笔记（三）Blast

使用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)
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
禁止转载，如需转载请通过简信或评论联系作者。
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 204,732评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 87,496评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 151,264评论 0赞 338
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,807评论 1赞 277
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,806评论 5赞 368
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,675评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,029评论 3赞 399
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,683评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 41,704评论 1赞 299
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,666评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,773评论 1赞 332
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,413评论 4赞 321
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,016评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,978评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,204评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,083评论 2赞 350
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,503评论 2赞 343

推荐阅读更多精彩内容

  • 【读书笔记】生信札记 第三章
    欢迎关注：oddxix 本章主要讲序列分析与联配 序列分析是生物信息学最主要的研究内容之一，它可以分为两个主要部分...
    oddxix阅读 2,205评论 0赞 10
  • Blast
    Blast，全称Basic Local Alignment Search Tool，即"基于局部比对算法的搜索工具...
    晓佥阅读 14,863评论 1赞 26
  • 为什么 Biopython 的在线 BLAST 这么慢？
    用过网页版本 BLAST 的童鞋都会发现，提交的序列比对往往在几分钟，甚至几十秒就可以得到比对的结果；而通过调用 ...
    BioIT爱好者阅读 1,934评论 0赞 0
  • 如何从BAM文件中提取fastq
    虽然高通量测序分析最常用的操作是将fastq比对到参考基因组得到BAM文件，但偶尔我们也需要提取BAM文件中特定区...
    xuzhougeng阅读 25,858评论 4赞 31
  • 什么是DISC？如何良好的应用？
    DISC国庆班结束已经5天了，至今才输出一篇关于其主题文章，不免有点鄙视自己的懒惰。同时，也有些兴奋，至少，...
    刘红leona阅读 1,274评论 1赞 4