最近在學(xué)習(xí) Bioinformatics with python cookbook 這本書(shū)第六章 Phylogenetics 的內(nèi)容稍计，了解到python中與系統(tǒng)發(fā)育相關(guān)的兩個(gè)模塊 Dendropy和 ete3 (A Python framework for the analysis and visualization of trees)，瀏覽ete3的文檔的時(shí)候發(fā)現(xiàn)了很多非常漂亮的圖片吸耿，第一感覺(jué)是和R語(yǔ)言里的ggtree功能很相似卿堂，所以覺(jué)得還是有必要學(xué)習(xí)一下睡腿。以下內(nèi)容記錄自己重復(fù)ete3文檔中漂亮圖片的過(guò)程买羞。（題外話：個(gè)人感覺(jué)python繪圖系統(tǒng)的默認(rèn)配色比R語(yǔ)言的配色漂亮一點(diǎn)）

• 第一步 安裝
  自己 windows 的電腦按住了Anaconda3,直接在DOS命令行下使用easy_install即可安裝相應(yīng)的python模塊.(正常應(yīng)該使用pip install安裝也是可以的，但是自己嘗試的時(shí)候遇到了報(bào)錯(cuò)，沒(méi)有搞清楚是什么原因)

easy_install ete3

• 第一個(gè)簡(jiǎn)單的小例子
  讀入樹(shù)文件缓屠，查看奇昙，然后保存為pdf文件

from ete3 import Tree
t = Tree("../../Desktop/Malus.output.fasta.treefile")
t.show()

25.PNG

有點(diǎn)像figtree

未完待續(xù)......

更新

將讀入的樹(shù)文件寫(xiě)入到新文件中

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") #寫(xiě)入到文件中

文檔的內(nèi)容有些枯燥，還是先從重復(fù)美圖開(kāi)始吧
t.show()函數(shù)運(yùn)行后會(huì)跳出來(lái)ETE Tree Browser窗口敌完，將樹(shù)顯示到桌面上
t.render()函數(shù)可以將樹(shù)輸出到圖片里储耐，可以生成png,pdf,svg格式
一個(gè)簡(jiǎn)單的小例子

from ete3 import Tree, TreeStyle
t = Tree()
t.render("mytree.png",w=183,units="mm")

mytree.png

• 第二個(gè)簡(jiǎn)單的小例子

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、第三個(gè)簡(jiǎn)單的小例子

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個(gè)小例子

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

• 第五個(gè)小例子

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 模塊的內(nèi)容

比對(duì)格式之間的轉(zhuǎn)化滨溉，比較常用的比如從fasta格式轉(zhuǎn)換成newick格式什湘，或者newick轉(zhuǎn)換成nexus格式，自己之前遇到此類問(wèn)題一直使用的是在線工具 http://sing.ei.uvigo.es/ALTER/ 晦攒。今天瀏覽dendropy文檔時(shí)發(fā)現(xiàn)這個(gè)模塊也可以實(shí)現(xiàn)格式轉(zhuǎn)換闽撤，多了一種選擇，簡(jiǎn)單記錄。（具體都可以轉(zhuǎn)換那些格式自己還不是很清楚，自己目前知道的是fasta,newick,nexus,phylip）使用到的示例文件
https://pan.baidu.com/s/1chchsxMjP2fM-ghKaOaArQ

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利用上一步的比對(duì)文件建一棵樹(shù)哑蔫，導(dǎo)出為newick格式布疙，然后利用dendropy模塊轉(zhuǎn)化為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")

查看樹(shù)（viewing and displaying trees）
兩種方式

• print_plot()可以查看拓?fù)浣Y(jié)構(gòu)
• as_string()可以查看文本形式的樹(shù)

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數(shù)據(jù)庫(kù)下載fasta格式的數(shù)據(jù)（這部分是重復(fù)Bioinformatics with python cookbook 這本書(shū)第六章 Phylogenetics 的內(nèi)容第一步：下載誒博拉病毒的基因組數(shù)據(jù)慧脱，之前嘗試了好多次一直沒(méi)有看懂書(shū)中的代碼，嘗試原封不動(dòng)的重復(fù)一直遇到錯(cuò)誤，今天瀏覽dendropy的文檔的過(guò)程中找到了一直遇到報(bào)錯(cuò)的原因：dendropy的部分代碼已經(jīng)更新，書(shū)中提到的部分代碼已經(jīng)不再使用）
先重復(fù)文檔中的兩個(gè)小例子

import dendropy
from dendropy.interop import genbank
gb_dna = genbank.GenBankDna(ids = ['EU105474','EU105475'])
#如果序列號(hào)之間是連續(xù)的舍沙，還可以換一種寫(xiě)法
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"))
#寫(xiě)到文件里
fw = open("dendropy_practice_1.fasta","w")
char_mat.write_to_stream(fw,'fasta')
fw.close()

接下來(lái)重復(fù)書(shū)中下載序列用到的的部分代碼（書(shū)中的內(nèi)容還涉及到了 yield 函數(shù)，自己還沒(méi)有太搞懂這個(gè)函數(shù)的用法 剔宪，可以參考 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'])


#原書(shū)中需要更新的代碼
#這部分代碼自己也不是太明白拂铡，反正目的是將序列的名字改成自己想要的格式

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()

下載所有序列用到的完整代碼(小插曲：第一次試運(yùn)行遇到了報(bào)錯(cuò)，仔細(xì)檢查才發(fā)現(xiàn)把序列號(hào)中的數(shù)字0錯(cuò)看成了字母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')
    

    
    
    
#原書(shū)中需要更新的代碼
#這部分代碼自己也不是太明白葱绒，反正目的是將序列的名字改成自己想要的格式

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()

接下來(lái)可以重復(fù)比對(duì)和建樹(shù)了