最近在学习 Bioinformatics with python cookbook 这本书第六章 Phylogenetics 的内容，了解到python中与系统发育相关的两个模块 Dendropy和 ete3 (A Python framework for the analysis and visualization of trees)，浏览ete3的文档的时候发现了很多非常漂亮的图片，第一感觉是和R语言里的ggtree功能很相似，所以觉得还是有必要学习一下。以下内容记录自己重复ete3文档中漂亮图片的过程。（题外话：个人感觉python绘图系统的默认配色比R语言的配色漂亮一点）
- 第一步 安装
自己 windows 的电脑按住了Anaconda3,直接在DOS命令行下使用easy_install即可安装相应的python模块.(正常应该使用pip install安装也是可以的，但是自己尝试的时候遇到了报错，没有搞清楚是什么原因)
easy_install ete3
- 第一个简单的小例子
读入树文件，查看，然后保存为pdf文件
from ete3 import Tree
t = Tree("../../Desktop/Malus.output.fasta.treefile")
t.show()
运行完 t.show() 会跳出来一个ETE Tree Browser
25.PNG
有点像figtree
未完待续......
更新
将读入的树文件写入到新文件中
from ete3 import Tree
t = Tree("(A:1,(B:1,(E:1,D:1)Internal_1:0.5)Internal_2:0.5)Root;")
t.write() #输出到屏幕
t.write(outfile="new_tree.nex") #写入到文件中
文档的内容有些枯燥，还是先从重复美图开始吧
t.show()函数运行后会跳出来ETE Tree Browser窗口，将树显示到桌面上
t.render()函数可以将树输出到图片里，可以生成png,pdf,svg格式
一个简单的小例子
from ete3 import Tree, TreeStyle
t = Tree()
t.render("mytree.png",w=183,units="mm")
mytree.png
- 第二个简单的小例子
from ete3 import Tree
from ete3 import TreeStyle
t = Tree()
t.populate(10)
ts.show_leaf_name = True
ts.mode = "c"
ts.arc_start = -180
ts.arc_span = 180
t.show(tree_style=ts)
t.render("tree.png",tree_style=ts)
tree.png
- 3、第三个简单的小例子
from ete3 import Tree
t = Tree("((((a1,a2),a3), ((b1,b2),(b3,b4))), ((c1,c2),c3));")
t.render("46.png")
46.png
from ete3 import Tree
from ete3 import NodeStyle
t = Tree("((((a1,a2),a3), ((b1,b2),(b3,b4))), ((c1,c2),c3));")
n1 = t.get_common_ancestor("a1","a2","a3")
nst1 = NodeStyle()
nst1["bgcolor"] = "LightSteelBlue"
n1.set_style(nst1)
t.render("47.png")
47.png
from ete3 import Tree
from ete3 import NodeStyle
from ete3 import AttrFace
from ete3 import faces
from ete3 import TreeStyle
t = Tree("((((a1,a2),a3), ((b1,b2),(b3,b4))), ((c1,c2),c3));")
n1 = t.get_common_ancestor("a1","a2","a3")
nst1 = NodeStyle()
nst1["bgcolor"] = "LightSteelBlue"
n1.set_style(nst1)
n2 = t.get_common_ancestor("b1","b2","b3","b4")
nst2 = NodeStyle()
nst2["bgcolor"] = "DarkSeaGreen"
n2.set_style(nst2)
def lauout(node):
if node.is_leaf():
N = AttrFace("name",fsize=30)
faces.add_face_to_node(N,node,0,position="aligned")
ts = TreeStyle()
ts.layout_fn = layout
ts.show_leaf_name = False
ts.render(tree_style = ts,file_name = "48.png")
48.png
rom ete3 import Tree
from ete3 import NodeStyle
from ete3 import AttrFace
from ete3 import faces
from ete3 import TreeStyle
t = Tree("((((a1,a2),a3), ((b1,b2),(b3,b4))), ((c1,c2),c3));")
for n in t.traverse():
n.dist = 2
n1 = t.get_common_ancestor("a1","a2","a3")
nst1 = NodeStyle()
nst1["bgcolor"] = "LightSteelBlue"
n1.set_style(nst1)
n2 = t.get_common_ancestor("b1","b2","b3","b4")
nst2 = NodeStyle()
nst2["bgcolor"] = "Moccasin"
n2.set_style(nst2)
n2 = t.get_common_ancestor("c1","c2","c3")
nst3 = NodeStyle()
nst3["bgcolor"] = "DarkSeaGreen"
n2.set_style(nst3)
ts = TreeStyle()
ts.mode = "c"
t.render(tree_style=ts,file_name="49.png",w=1000,h=1000,dpi=300)
49.png
- 第4个小例子
from ete3 import Tree
from ete3 import TreeStyle
from ete3 import faces
from ete3 import AttrFace
from ete3 import PieChartFace
from ete3 import COLOR_SCHEMES
from random import sample
from random import randint
t = Tree("((((a1,a2),a3), ((b1,b2),(b3,b4))), ((c1,c2),c3));")
ts = TreeStyle()
def layout(node):
if node.is_leaf():
N = AttrFace("name",fsize=20)
faces.add_face_to_node(N,node,column=0,position="branch-right")
pieF = PieChartFace([10,20,60,10],colors=COLOR_SCHEMES[sample(schema_names,1)[0]],width=40,height=40)
faces.add_face_to_node(pieF,node,column=0,position="aligned")
else:
node.img_style["size"] = randint(3,6)
node.img_style["shape"] = "square"
node.img_style["fgcolor" ] = "green"
ts.layout_fn = layout
ts.show_leaf_name = False
ts.show_scale = False
t.render(tree_style=ts,file_name = "50.png",w=500,h=500)
50.png
- 第五个小例子
from ete3 import Tree
from ete3 import TreeStyle
from ete3 import faces
from ete3 import TextFace
from ete3 import AttrFace
from ete3 import CircleFace
from random import randint
t = Tree("((((a1,a2),a3), ((b1,b2),(b3,b4))), ((c1,c2),c3));")
def layout(node):
if node.is_leaf():
N = AttrFace("name",fsize=20)
faces.add_face_to_node(N,node,column=0,position="branch-right")
node.img_style["size"] = 0
else:
node.img_style['size'] = randint(5,8)
node.img_style["shape"] = "square"
node.img_style["fgcolor"] = "green"
bubble_face = CircleFace(randint(5,10),'steelblue','sphere')
bubble_face.opacity = 0.6
faces.add_face_to_node(bubble_face,node,column=0,position="float-behind")
faces.add_face_to_node(AttrFace("dist",fsize=7,fgcolor="red"),node,column=0,position="branch-top")
if node.up and not node.up.up:
node.img_style['vt_line_width'] = 3
node.img_style['hz_line_width'] = 4
ts = TreeStyle()
ts.lsyout _fn = layout
ts.show_leaf_name = False
ts.show_scale = False
ts.mode = 'c'
ts.arc_start = 270
ts.arc_span = 185
t.show(tree_style=ts)
t.render(tree_style=ts,w=800,file_name="51.png")
51.png
更新 Dendropy 模块的内容
比对格式之间的转化，比较常用的比如从fasta格式转换成newick格式，或者newick转换成nexus格式，自己之前遇到此类问题一直使用的是在线工具 http://sing.ei.uvigo.es/ALTER/ 。今天浏览dendropy文档时发现这个模块也可以实现格式转换，多了一种选择，简单记录。（具体都可以转换那些格式自己还不是很清楚，自己目前知道的是fasta,newick,nexus,phylip）使用到的示例文件
import dendropy
ccsA = dendropy.DnaCharacterMatrix.get(path = "ccsA_KaKs_pra.fas", schema = "fasta")
ccsA.write(path = "ccsA.phy",schema = "phylip")
ccsA.write(path = "ccsA.newick", schema = "newick")
ccsA.write(path = "ccsA.nexus", schema = "nexus")
使用mega利用上一步的比对文件建一棵树，导出为newick格式，然后利用dendropy模块转化为nexus格式（converting a single tree from Newick schema to nexus）
import dendropy
ccsA = dendropy.Tree.get(path = "ccsA.newick",schema = "newick")
ccsA.write(path="ccsA.nex",schema = "nexus")
查看树（viewing and displaying trees）
两种方式
- print_plot()可以查看拓扑结构
- as_string()可以查看文本形式的树
import dendropy
t = dendropy.Tree.get(path = "ccsA.newick",schema = "newick")
t.print_plot()
print(t.as_string(schema="newick"))
print(t.as_string(schema="nexus"))
自genbank数据库下载fasta格式的数据（这部分是重复Bioinformatics with python cookbook 这本书第六章 Phylogenetics 的内容第一步：下载诶博拉病毒的基因组数据，之前尝试了好多次一直没有看懂书中的代码，尝试原封不动的重复一直遇到错误，今天浏览dendropy的文档的过程中找到了一直遇到报错的原因：dendropy的部分代码已经更新，书中提到的部分代码已经不再使用）
先重复文档中的两个小例子
import dendropy
from dendropy.interop import genbank
gb_dna = genbank.GenBankDna(ids = ['EU105474','EU105475'])
#如果序列号之间是连续的，还可以换一种写法
gb_dna = genbank.GenBankDna(id_range=(74,75),prefix="EU1054")
for gb in gb_dna:
print(gb)
char_mat = gb_dna.generate_char_matrix()
#输出到屏幕
print(char_mat.as_string("fasta"))
#写到文件里
fw = open("dendropy_practice_1.fasta","w")
char_mat.write_to_stream(fw,'fasta')
fw.close()
接下来重复书中下载序列用到的的部分代码（书中的内容还涉及到了 yield 函数，自己还没有太搞懂这个函数的用法 ，可以参考 https://www.ibm.com/developerworks/cn/opensource/os-cn-python-yield/）
import dendropy
from dendropy.interop import genbank
def get_other_ebolavirus_sources():
yield 'BDBV', genbank.GenBankDna(id_range=(3,6),prefix='KC54539')
yield 'BDBV', genbank.GenBankDna(ids=['FJ217161'])
yield 'RESTV', genbank.GenBankDna(ids=['AB050936','JX477165','JX477166','FJ621583','FJ621584','FJ621585'])
yield 'SUDV', genbank.GenBankDna(ids=['KC242783','AY729654','EU338380','JN638998','FJ968794','KC589025','JN638998'])
yield 'SUDV', genbank.GenBankDna(id_range=(89,92),prefix='KC5453')
yield 'TAFV', genbank.GenBankDna(ids=['FJ217162'])
#原书中需要更新的代码
#这部分代码自己也不是太明白，反正目的是将序列的名字改成自己想要的格式
def gb_to_taxon(gb,taxon_namespace):
label = species + "_" + gb.accession
taxon = taxon_namespace.require_taxon(label=label)
return taxon
taxon_namespace = dendropy.TaxonNamespace()
other = open('other.fasta','w')
for species, recs in get_other_ebolavirus_sources():
char_mat = recs.generate_char_matrix(taxon_namespace = taxon_namespace,gb_to_taxon_fn = gb_to_taxon)
print(char_mat.as_string("fasta"))
char_mat.write_to_stream(other,'fasta')
other.close()
下载所有序列用到的完整代码(小插曲：第一次试运行遇到了报错，仔细检查才发现把序列号中的数字0错看成了字母O)
import dendropy
from dendropy.interop import genbank
def get_other_ebolavirus_sources():
yield 'BDBV', genbank.GenBankDna(id_range=(3,6),prefix='KC54539')
yield 'BDBV', genbank.GenBankDna(ids=['FJ217161'])
yield 'RESTV', genbank.GenBankDna(ids=['AB050936','JX477165','JX477166','FJ621583','FJ621584','FJ621585'])
yield 'SUDV', genbank.GenBankDna(ids=['KC242783','AY729654','EU338380','JN638998','FJ968794','KC589025','JN638998'])
yield 'SUDV', genbank.GenBankDna(id_range=(89,92),prefix='KC5453')
yield 'TAFV', genbank.GenBankDna(ids=['FJ217162'])
def get_ebov_2014_sources():
yield 'EBOV_2014', genbank.GenBankDna(id_range=(233036,233118),prefix="KM")
yield 'EBOV_2014', genbank.GenBankDna(id_range=(34549,34563),prefix='KM0')
def get_other_ebov_sources():
yield 'EBOV_1976', genbank.GenBankDna(ids=['AF272001','KC242801'])
yield 'EBOV_1995', genbank.GenBankDna(ids=['KC242796','KC242799'])
yield 'EBOV_2007', genbank.GenBankDna(id_range=(84,90),prefix='KC2427')
#原书中需要更新的代码
#这部分代码自己也不是太明白，反正目的是将序列的名字改成自己想要的格式
def gb_to_taxon(gb,taxon_namespace):
label = species + "_" + gb.accession
taxon = taxon_namespace.require_taxon(label=label)
return taxon
taxon_namespace = dendropy.TaxonNamespace()
sampled = open('sample.fasta','w')
for species, recs in get_other_ebolavirus_sources():
char_mat = recs.generate_char_matrix(taxon_namespace = taxon_namespace,gb_to_taxon_fn = gb_to_taxon)
char_mat.write_to_stream(sampled,'fasta')
def gb_to_taxon1(gb,taxon_namespace):
label = "EBOV_2014_" + gb.accession
taxon = taxon_namespace.require_taxon(label=label)
return taxon
for species, recs in get_ebov_2014_sources():
char_mat = recs.generate_char_matrix(taxon_namespace = taxon_namespace,gb_to_taxon_fn = gb_to_taxon1)
char_mat.write_to_stream(sampled,'fasta')
for species, rec in get_other_ebov_sources():
char_mat = recs.generate_char_matrix(taxon_namespace = taxon_namespace,gb_to_taxon_fn = gb_to_taxon1)
char_mat.write_to_stream(sampled,'fasta')
sampled.close()
接下来可以重复比对和建树了
