基因组作图软件GenomeDiagram.pdf

《基因组作图软件GenomeDiagram.pdf》由会员分享，可在线阅读，更多相关《基因组作图软件GenomeDiagram.pdf（35页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第第 17 章章 Graphics 模块中的基模块中的基因组可视化包因组可视化包GenomeDiagramBio.Graphics模块基于 Python 的第三方扩展库 ReportLab，ReportLab 主要生成PDF 文件，同时也能生成 EPS（Encapsulated Postscript）文件和 SVG文件。ReportLa 可以导出矢量图，如果安装依赖关系（Dependencies），比如 PIL(PythonImaging Library)，ReportLab 也可以导出 JPEG,PNG,GIF,BMP和 PICT 格式的位图(Bitmap image)。17.1 基因组可视

2、化包基因组可视化包GenomeDiagram17.1.1 GenomeDiagram简介简介Bio.Graphics.GenomeDiagram包被整合到 Biopython 1.50版之前，就已经是 Biopython的独立模块。GenomeDiagram包首次出现在 2006 年 Pritchard 等人在Bioinformatics 杂志的一篇文章 :ref:2 ，文中展示了一些图像示例，更多图像示例请查看 GenomeDiagram手册 http:/biopython.org/DIST/docs/GenomeDiagram/userguide.pdf。正如“GenomeDiagram”

3、名称所指，它主要用于可视化全基因组(特别是原核生物基因组)，即可绘制线型图也可绘制环形图，Toth 等人在 2006 年发表的文章 :ref:3 中图 2 就是一个示例。Van der Auwera 等人在 2009 年发表的文章:ref:4 中图 1 和图 2 也进一步说明，GenomeDiagram适用于噬菌体、质粒和线粒体等微小基因组的可视化。如果存储基因组信息的是从 GenBank文件中下载的SeqRecord话，它会包含许多SeqFeature，那么用这个模块处理就很简单（详见 第:ref:4 章和第:ref:5 章）。17.1.2 图形，轨迹图形，轨迹,特征集和特征特征集和特征Ge

4、nomeDiagram使用一组嵌套的对象，图层中沿着水平轴或圆圈的图形对象（diagram object）表示一个序列（sequence）或序列区域（sequence region）。一个图形可以包含多个轨迹（track），呈现为横向排列或者环形放射图。这些轨迹的长度通常相等，代表相同的序列区域。可用一个轨迹表示基因的位置，另一个轨迹表示调节区域，第三个轨迹表示 GC 含量。可将最常用轨迹的特征打包为一个特征集（feature-sets）。CDS的特征可以用一个特征集，而tRNA 的特征可以用另外一个特征集。这不是强制性的要求，你可以在 diagram 中用同样的特征集。如果diagram 中

5、用不同的特征集，修改一个特征会很容易，比如把所有 tRNA 的特征都变为红色，你只需选择 tRNA 的特征就行。新建图形主要有两种方式。第一种是自上而下的方法（Top-Down），首先新建diagram 对象，然后用diagram 的方法来添加 track(s)，最后用track的方法添加特征。第二种是自下而上的方法（Bottom-Up），首先单独新建对象，然后再将其进行组合。17.1.3 自上而下的实例自上而下的实例我们用一个从 GenBank文件中读取出来的SeqRecord来绘制全基因组（详见第:ref:5 章）。这里用鼠疫杆菌 Yersinia pestis biovarMicrotu

6、s的 pPCP1质粒，元数据文件 NC_005816.gb在 Biopython中 GenBank的 tests目录下，NC_005816.gb也可下载fromfrom reportlab.lib importimport colorsfromfrom reportlab.lib.units importimport cmfromfrom Bio.Graphics importimport GenomeDiagramfromfrom Bio importimport SeqIOrecord=SeqIO.read(NC_005816.gb,genbank)这里用自上而下的方法，导入目标序列后，新

7、建一个 diagram，然后新建一个 track，最后新建一个特征集（feature set）：gd_diagram=GenomeDiagram.Diagram(Yersinia pestis biovar Microtus plasmidpPCP1)gd_track_for_features=gd_diagram.new_track(1,name=Annotated Features)gd_feature_set=gd_track_for_features.new_set()接下来的部分最有趣，提取SeqRecord中每个基因的SeqFeature对象，就会为 diagram生成一个相应的特

8、征（feature），将其颜色设置为蓝色，分别用深蓝和浅蓝表示。forfor feature inin record.features:ifif feature.type!=!=gene:#Exclude this featurecontinuecontinueifif len(gd_feature_set)%2=0:color=colors.blueelseelse:color=colors.lightblue gd_feature_set.add_feature(feature,color=color,label=True)创建导出文件需要两步，首先是draw方法，它用ReportLab

9、对象生成全部图形。然后是write方法，将图形存储到格式文件。注意：输出文件格式不止一种。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.write(plasmid_linear.eps,EPS)gd_diagram.write(plasmid_linear.svg,SVG)如果安装了依赖关系（Dependencies），也可以生成位图（Bitm

10、ap image），代码如下：gd_diagram.write(plasmid_linear.png,PNG)注意，我们将代码中的fragments变量设置为“4”，基因组就会被分为“4”个片段。如果想要环形图，可以试试以下的代码：gd_diagram.draw(format=circular,circular=True,pagesize=(20*cm,20*cm),start=0,end=len(record),circle_core=0.7)gd_diagram.write(plasmid_circular.pdf,PDF)示例图不是非常精彩，但这仅仅是精彩的开始。17.1.4 自下而上的

11、实例自下而上的实例现在，用“自下而上”的方法来创建相同的图形。首先新建不同的对象（可以是任何顺序），然后将其组合。fromfrom reportlab.lib importimport colorsfromfrom reportlab.lib.units importimport cmfromfrom Bio.Graphics importimport GenomeDiagramfromfrom Bio importimport SeqIOrecord=SeqIO.read(NC_005816.gb,genbank)#Create the feature set and its feature

12、 objects,gd_feature_set=GenomeDiagram.FeatureSet()forfor feature inin record.features:ifif feature.type!=!=gene:#Exclude this featurecontinuecontinueifif len(gd_feature_set)%2=0:color=colors.blueelseelse:color=colors.lightblue gd_feature_set.add_feature(feature,color=color,label=True)#(this for loop

13、 is the same as in the previous example)#Create a track,and a diagramgd_track_for_features=GenomeDiagram.Track(name=Annotated Features)gd_diagram=GenomeDiagram.Diagram(Yersinia pestis biovar Microtus plasmidpPCP1)#Now have to glue the bits together.gd_track_for_features.add_set(gd_feature_set)gd_dia

14、gram.add_track(gd_track_for_features,1)同样，利用draw和write方法来创建线形图或者环形图，结果应该完全相同（“draw”和“write”部分的代码见 17.1.3）。17.1.5 简单的简单的 Feature以上示例中，创建 diagram 使用的SeqRecord的SeqFeature对象（详见:ref:4.3章节）。如果你不需要SeqFeature对象，只将目标 feature定位在坐标轴，仅需要创建 minimalSeqFeature对象，方法很简单，代码如下：fromfrom Bio.SeqFeature importimport SeqF

15、eature,FeatureLocationmy_seq_feature=SeqFeature(FeatureLocation(50,100),strand=+=+1)对于序列来说，+1代表正向，-1代表反向，None代表两者都有，下面举个简单的示例：fromfrom Bio.SeqFeature importimport SeqFeature,FeatureLocationfromfrom Bio.Graphics importimport GenomeDiagramfromfrom reportlab.lib.units importimport cmgdd=GenomeDiagram.D

16、iagram(Test Diagram)gdt_features=gdd.new_track(1,greytrack=False)gds_features=gdt_features.new_set()#Add three features to show the strand options,feature=SeqFeature(FeatureLocation(25,125),strand=+=+1)gds_features.add_feature(feature,name=Forward,label=True)feature=SeqFeature(FeatureLocation(150,25

17、0),strand=None)gds_features.add_feature(feature,name=Strandless,label=True)feature=SeqFeature(FeatureLocation(275,375),strand=1)gds_features.add_feature(feature,name=Reverse,label=True)gdd.draw(format=linear,pagesize=(15*cm,4*cm),fragments=1,start=0,end=400)gdd.write(GD_labels_default.pdf,pdf)图形示例结果

18、请见下一节图中的第一个图，缺省的 feature为浅绿色。注意，这里用name参数作为 feature的“说明文本”（caption text）。下文将会讲述更多细节。17.1.6 Feature说明说明下面代码中，feature作为SeqFeature的对象添加到 diagram。gd_feature_set.add_feature(feature,color=color,label=True)前面的示例用SeqFeature的注释为 feature做了恰当的文字说明。SeqFeature对象的限定符（qualifiers dictionary）缺省值是：gene,label,name,lo

19、cus_tag,和product。简单地说，你可以定义一个名称：gd_feature_set.add_feature(feature,color=color,label=True,name=My Gene)每个 feature标签的说明文本可以设置字体、位置和方向。说明文本默认的位置在图形符号（sigil）的左边，可选择在中间或者右边，线形图中文本的默认方向是45旋转。#Large font,parallel with the trackgd_feature_set.add_feature(feature,label=True,color=green,label_size=25,label_a

20、ngle=0)#Very small font,perpendicular to the track(towards it)gd_feature_set.add_feature(feature,label=True,color=purple,label_position=end,label_size=4,label_angle=90)#Small font,perpendicular to the track(away from it)gd_feature_set.add_feature(feature,label=True,color=blue,label_position=middle,l

21、abel_size=6,label_angle=90)用前面示例的代码将这三个片段组合之后应该可以得到如下的结果：除此之外，还可以设置“label_color”来调节标签的颜色（第:ref:17.1.9 节也将用到这一步），这里没有进行演示。示例中默认的字体很小，这是比较明智的，因为通常我们会把许多 Feature同时展示，而不像这里只展示了几个比较大的 feature。17.1.7 表示表示 Feature 的图形符号的图形符号以上示例中 Feature的图形符号（sigil）默认是一个方框（plain box），GenomeDiagram第一版中只有这一选项，后来 GenomeDiagra

22、m被整合到Biopython1.50时，新增了箭头状的图形符号（sigil）。#Default uses a BOX sigilgd_feature_set.add_feature(feature)#You can make this explicit:gd_feature_set.add_feature(feature,sigil=BOX)#Or opt for an arrow:gd_feature_set.add_feature(feature,sigil=ARROW)Biopython 1.61又新增 3 个图形形状（sigil）。#Box with corners cut off(m

23、aking it an octagon)gd_feature_set.add_feature(feature,sigil=OCTO)#Box with jagged edges(useful for showing breaks in contains)gd_feature_set.add_feature(feature,sigil=JAGGY)#Arrow which spans the axis with strand used only for directiongd_feature_set.add_feature(feature,sigil=BIGARROW)下面就是这些新增的图形形状

24、（sigil），多数的图形形状都在边界框（bounding box）内部，在坐标轴的上/下位置代表序列（Strand）方向的正/反向，或者上下跨越坐标轴，高度是其他图形形状的两倍。“BIGARROW”有所不同，它总是跨越坐标轴，方向由 feature的序列决定。17.1.8 箭头形状箭头形状上一部分我们简单引出了箭头形状。还有两个选项可以对箭头形状进行设置：首先根据边界框的高度比例来设置箭杆宽度。#Full height shafts,giving pointed boxes:gd_feature_set.add_feature(feature,sigil=ARROW,color=brown,

25、arrowshaft_height=1.0)#Or,thin shafts:gd_feature_set.add_feature(feature,sigil=ARROW,color=teal,arrowshaft_height=0.2)#Or,very thin shafts:gd_feature_set.add_feature(feature,sigil=ARROW,color=darkgreen,arrowshaft_height=0.1)结果见下图：其次，根据边界框的高度比例设置箭头长度（默认为 0.5或 50%）：#Short arrow heads:gd_feature_set.ad

26、d_feature(feature,sigil=ARROW,color=blue,arrowhead_length=0.25)#Or,longer arrow heads:gd_feature_set.add_feature(feature,sigil=ARROW,color=orange,arrowhead_length=1)#Or,very very long arrow heads(i.e.all head,no shaft,so triangles):gd_feature_set.add_feature(feature,sigil=ARROW,color=red,arrowhead_l

27、ength=10000)结果见下图：Biopython1.61新增BIGARROW箭头形状，它经常跨越坐标轴，箭头指向”左边“代表”反向“，指向”右边“代表”正向“。#A large arrow straddling the axis:gd_feature_set.add_feature(feature,sigil=BIGARROW)上述ARROW形状中的箭杆和箭头设置选项都适用于BIGARROW。17.1.9 完美示例完美示例回到”自上而下的示例 Section:ref:17.1.3 中鼠疫杆菌 Yersinia pestis biovar Microtus的 pPCP1质粒，现在使用”图形

28、符号“的高级选项。箭头表示基因，窄框穿越箭头表示限制性内切酶的切割位点。fromfrom reportlab.lib importimport colorsfromfrom reportlab.lib.units importimport cmfromfrom Bio.Graphics importimport GenomeDiagramfromfrom Bio importimport SeqIOfromfrom Bio.SeqFeature importimport SeqFeature,FeatureLocationrecord=SeqIO.read(NC_005816.gb,genba

29、nk)gd_diagram=GenomeDiagram.Diagram(record.id)gd_track_for_features=gd_diagram.new_track(1,name=Annotated Features)gd_feature_set=gd_track_for_features.new_set()forfor feature inin record.features:ifif feature.type!=!=gene:#Exclude this featurecontinuecontinueifif len(gd_feature_set)%2=0:color=color

30、s.blueelseelse:color=colors.lightblue gd_feature_set.add_feature(feature,sigil=ARROW,color=color,label=True,label_size=14,label_angle=0)#I want to include some strandless features,so for an example#will use EcoRI recognition sites etc.forfor site,name,color inin(GAATTC,EcoRI,colors.green),(CCCGGG,Sm

31、aI,colors.orange),(AAGCTT,HindIII,colors.red),(GGATCC,BamHI,colors.purple):index=0whilewhile True:index=record.seq.find(site,start=index)ifif index=1:breakbreak feature=SeqFeature(FeatureLocation(index,index+len(site)gd_feature_set.add_feature(feature,color=color,name=name,label=True,label_size=10,l

32、abel_color=color)index+=+=len(site)gd_diagram.draw(format=linear,pagesize=A4,fragments=4,start=0,end=len(record)gd_diagram.write(plasmid_linear_nice.pdf,PDF)gd_diagram.write(plasmid_linear_nice.eps,EPS)gd_diagram.write(plasmid_linear_nice.svg,SVG)gd_diagram.draw(format=circular,circular=True,pagesiz

33、e=(20*cm,20*cm),start=0,end=len(record),circle_core=0.5)gd_diagram.write(plasmid_circular_nice.pdf,PDF)gd_diagram.write(plasmid_circular_nice.eps,EPS)gd_diagram.write(plasmid_circular_nice.svg,SVG)输出结果见下图：17.1.10 多重轨迹前面实例中都是单独的 track，我们可以创建多个 track，比如，一个 track展示基因，另一个 track展示重复序列。Proux 等人 2002 年报道的文

34、章 5 中图 6 是一个很好的范例，下面我们将三个噬菌体基因组依次进行展示。首先需要三个噬菌体的GenBank文件。NC_002703 Lactococcus phage Tuc2009,全基因组大小(38347 bp)AF323668 Bacteriophage bIL285,全基因组大小(35538 bp)NC_003212 Listeria innocua Clip11262,我们将仅关注前噬菌体 5 的全基因组(长度大体相同).这三个文件可以从 Entrez下载，详情请查阅:ref:9.6。从三个噬菌体基因组文件中分离（slice）提取相关 Features信息（请查阅:ref:4.6

35、），保证前两个噬菌体的反向互补链与其起始点对齐，再次保存 Feature(详情请查阅:ref:4.8 )。fromfrom Bio importimport SeqIOA_rec=SeqIO.read(NC_002703.gbk,gb)B_rec=SeqIO.read(AF323668.gbk,gb)C_rec=SeqIO.read(NC_003212.gbk,gb)2587879:2625807.reverse_complement(name=True)图像中用不同颜色表示基因功能的差异。这需要编辑 GenBank文件中每一个feature的颜色参数就像用 Sangers Artemis e

36、ditor 处理 才能被GenomeDiagram识别。但是，这里只需要硬编码（hard code）三个颜色列表。上述 GenBank文件中的注释信息与 Proux 所用的文件信息并不完全相同，他们还添加了一些未注释的基因。fromfrom reportlab.lib.colors importimport red,grey,orange,green,brown,blue,lightblue,purpleA_colors=red*5+grey*7+orange*2+grey*2+orange+grey*11+green*4+grey+green*2+grey,green+brown*5+blu

37、e*4+lightblue*5+grey,lightblue+purple*2+greyB_colors=red*6+grey*8+orange*2+grey+orange+grey*21+green*5+grey+brown*4+blue*3+lightblue*3+grey*5+purple*2C_colors=grey*30+green*5+brown*4+blue*2+grey,blue+lightblue*2+grey*5接下来是“draw”方法，给 diagram 添加 3 个 track。我们在示例中设置不同的开始/结束值来体现它们之间长度不等（Biopython 1.59及更高

38、级的版本）。fromfrom Bio.Graphics importimport GenomeDiagramname=Proux Fig 6gd_diagram=GenomeDiagram.Diagram(name)max_len=0forfor record,gene_colors inin zip(A_rec,B_rec,C_rec,A_colors,B_colors,C_colors):max_len=max(max_len,len(record)gd_track_for_features=gd_diagram.new_track(1,name=record.name,greytrack

39、=True,start=0,end=len(record)gd_feature_set=gd_track_for_features.new_set()i=0forfor feature inin record.features:ifif feature.type!=!=gene:#Exclude this featurecontinuecontinue gd_feature_set.add_feature(feature,sigil=ARROW,color=gene_colorsi,label=True,name=str(i+1),label_position=start,label_size

40、=6,label_angle=0)i+=+=1gd_diagram.draw(format=linear,pagesize=A4,fragments=1,start=0,end=max_len)gd_diagram.write(name+.pdf,PDF)gd_diagram.write(name+.eps,EPS)gd_diagram.write(name+.svg,SVG)结果如图所示：在示例图中底部的噬菌体没有红色或橙色的基因标记。另外，三个噬菌体可视化图的长度不同，这是因为它们的比例相同，长度却不同。另外有一点不同，不同噬菌体的同源蛋白质之间用有颜色的 links 相连，下一部分将解决

41、这个问题。17.1.11 不同不同 Track 之间的之间的 Cross-LinksBiopython 1.59新增绘制不同 track之间 Cross-Links的功能，这个功能可用于将要展示的简单线形图中，也可用于将线形图分割为短片段（fragments）和环形图。我们接着模仿 Proux 等人 :ref:5 的图像，我们需要一个包含基因之间的“cross links”、“得分”或“颜色”的列表。实际应用中，可以从 BLAST文件自动提取这些信息，这里是手动输入的。噬菌体的名称同样表示为 A，B 和 C。这里将要展示的是A 与 B 之间的 links，噬菌体 A 和 B 基因的相似百分比存

42、储在元组中。#Tuc2009(NC_002703)vs bIL285(AF323668)A_vs_B=(99,Tuc2009_01,int),(33,Tuc2009_03,orf4),(94,Tuc2009_05,orf6),(100,Tuc2009_06,orf7),(97,Tuc2009_07,orf8),(98,Tuc2009_08,orf9),(98,Tuc2009_09,orf10),(100,Tuc2009_10,orf12),(100,Tuc2009_11,orf13),(94,Tuc2009_12,orf14),(87,Tuc2009_13,orf15),(94,Tuc2009

43、_14,orf16),(94,Tuc2009_15,orf17),(88,Tuc2009_17,rusA),(91,Tuc2009_18,orf20),(93,Tuc2009_19,orf22),(71,Tuc2009_20,orf23),(51,Tuc2009_22,orf27),(97,Tuc2009_23,orf28),(88,Tuc2009_24,orf29),(26,Tuc2009_26,orf38),(19,Tuc2009_46,orf52),(77,Tuc2009_48,orf54),(91,Tuc2009_49,orf55),(95,Tuc2009_52,orf60),对噬菌体

44、 B 和 C 做同样的处理：#bIL285(AF323668)vs Listeria innocua prophage 5(in NC_003212)B_vs_C=(42,orf39,lin2581),(31,orf40,lin2580),(49,orf41,lin2579),#terL (54,orf42,lin2578),#portal (55,orf43,lin2577),#protease (33,orf44,lin2576),#mhp (51,orf46,lin2575),(33,orf47,lin2574),(40,orf48,lin2573),(25,orf49,lin2572)

45、,(50,orf50,lin2571),(48,orf51,lin2570),(24,orf52,lin2568),(30,orf53,lin2567),(28,orf54,lin2566),噬菌体 A 和 C 的标识符（Identifiers）是基因座标签（locus tags），噬菌体B 没有基因座标签，这里用基因名称来代替。以下的辅助函数可用基因座标签或基因名称来寻找 Feature。defdef get_featureget_feature(features,id,tags=locus_tag,gene):Search list of SeqFeature objects for an

46、 identifier under the giventags.forfor f inin features:forfor key inin tags:#tag may not be present in this featureforfor x inin f.qualifiers.get(key,):ifif x=id:returnreturn fraiseraise KeyErrorKeyError(id)现在将这些标识符对（identifier pairs）的列表转换为“SeqFeature”列表，因此来查找它们的坐标定位。现在将下列代码添加到上段代码中gd_diagram.draw(.

47、)这一行之前，将 cross-links 添加到图像中。示例中的脚本文件 Proux_et_al_2002_Figure_6.py 在 Biopython源程序文件夹的Doc/examples目录下。fromfrom Bio.Graphics.GenomeDiagram importimport CrossLinkfromfrom reportlab.lib importimport colors#Note it might have been clearer to assign the track numbers explicitly.forfor rec_X,tn_X,rec_Y,tn_Y

48、,X_vs_Y inin(A_rec,3,B_rec,2,A_vs_B),(B_rec,2,C_rec,1,B_vs_C):track_X=gd_diagram.trackstn_X track_Y=gd_diagram.trackstn_Yforfor score,id_X,id_Y inin X_vs_Y:feature_X=get_feature(rec_X.features,id_X)feature_Y=get_feature(rec_Y.features,id_Y)color=colors.linearlyInterpolatedColor(colors.white,colors.f

49、irebrick,0,100,score)link_xy=CrossLink(track_X,feature_X.location.start,feature_X.location.end),(track_Y,feature_Y.location.start,feature_Y.location.end),color,colors.lightgrey)gd_diagram.cross_track_links.append(link_xy)这段代码有几个要点，第一，GenomeDiagram对象有一个cross_track_links属性，这个属性只是CrossLink对象的一组数据。每个Cro

50、ssLink对象有两个 track-specific坐标，示例中用元组（tuples）来展现，可用GenomeDiagram.Feature对象来代替。可选择添加颜色和边框颜色，还可以说明这个 link 是否需要翻转，这个功能易于表现染色体异位。你也可以看我们是如何将 BLAST中特征百分比（Percentage Identity Score）转换为白-红的渐变色（白-0%，红-100%）。这个实例中没有 cross-links 的重叠，如果有 links 重叠可以用 ReportLab 库中的透明度（transparency）来解决，通过设置颜色的 alpha通道来使用。然而，若同时使用边框