使用biopython可视化染色体和基因元件

欢迎关注”生信修炼手册”!

基因组结构元件的可视化有多种方式，比如IGV等基因组浏览器中以track为单位的展示形式，亦或以circos为代表的圈图形式，比如在细胞器基因组组装中，基因元件常用圈图形式展示，示例如下

在biopython中，通过BiolGraphics子模块可以对基因组结构进行可视化，支持线性和圈图两种可视化方式。其中，基因组结构信息存储在genebank格式的文件中，首先通过Bio.SeqIO读取结构信息，然后通过Bio.Graphics模块进行可视化。

以下列数据为例，先来看下可视化的用法

>https://www.ncbi.nlm.nih.gov/nuccore/NC_005816

首先是读取gb文件，代码如下

>>> from reportlab.lib import colors
>>> from reportlab.lib.units import cm
>>> from Bio.Graphics import GenomeDiagram
>>> from Bio import SeqIO
>>> record = SeqIO.read("sequence.gb", "genbank")

接下来提取gb文件中的feature信息，构建用于绘图的数据结构，代码如下

>>> gd_diagram = GenomeDiagram.Diagram("Yersinia pestis biovar Microtus plasmid pPCP1")
>>> gd_track_for_features = gd_diagram.new_track(1, name="Annotated Features")
>>> gd_feature_set = gd_track_for_features.new_set()
>>> for feature in record.features:
...     if feature.type != "gene":
...         continue
...     if len(gd_feature_set) % 2 == 0:
...         color = colors.blue
...     else:
...         color = colors.lightblue
...     gd_feature_set.add_feature(feature, color=color, label=True)

最后进行绘图即可，代码如下

>>> gd_diagram.draw(format="linear", orientation="landscape", pagesize='A4',
... fragments=4, start=0, end=len(record))
>>>
>>> gd_diagram.write("plasmid_linear.pdf", "PDF")

输出结果如下

对于圈图，只需要修改简单修改绘图的参数即可，代码如下

>>> gd_diagram.draw(format="circular", circular=True, pagesize=(20*cm,20*cm),
... start=0, end=len(record), circle_core=0.7)
>>> gd_diagr