社区首页 >专栏 >跟着Nature ecology&evolution学python：toytree模块画进化树

跟着Nature ecology&evolution学python：toytree模块画进化树

用户7010445

发布于 2023-01-06 11:13:18

63100

代码可运行

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

运行总次数：0

代码可运行

论文

Replicated radiation of a plant clade along a cloud forest archipelago

https://www.nature.com/articles/s41559-022-01823-x#data-availability

数据代码链接

https://zenodo.org/record/5504439#.YtpqTXZBzic

https://github.com/eaton-lab/Oreinotinus-phylogeny

论文里提供了几乎所有数据和代码，要好好学习一下这篇论文

今天推文的主要内容是重复其中 Extended Data Fig. 1 ，不对绘图代码进行仔细研究了，知道python里有一个 toytree的模块可以展示进化树，看了代码，应该是不如ggtree好用帮助文档链接 https://toytree.readthedocs.io/en/latest/3-installation.html

github主页链接https://github.com/eaton-lab/toytree

还有一个模块是ipyrad，用来做群体遗传分析，帮助文档链接 https://ipyrad.readthedocs.io/en/master/

完整的作图代码

import toytree
import pandas as pd
fulldata = pd.read_csv("../Oreinotinus-phylogeny-main/Raw_data/oreinotinus_samples_database.csv")
colors = pd.read_csv("../Oreinotinus-phylogeny-main/Raw_data/oreinotinus_color_codes.csv")
sdata = fulldata[["NameInAssembly","Lastest_SP_name", "Num_for_Publication", "UltimateName"]]

namedict = {}
for i in range(sdata.shape[0]):
    if sdata.iloc[i, 2]:
        number = " (" + str(sdata.iloc[i, 2]) + ")"
    else:
        number = ""
    namedict[sdata.iloc[i, 0]] = f"V. {sdata.iloc[i, 1]}{number}"
        

# load color data and put in a dictionary
colordata = colors[["Species","Color"]]
colordict = {colordata.iloc[i, 0]: str(colordata.iloc[i, 1]) for i in range(colordata.shape[0])}

treeFile = "../Oreinotinus-phylogeny-main/Analyses/Individual-level-tree-reconstruction/analysis-raxml/RAxML_bipartitions.fulldataset_withAyava_10scaff_mcov025_rmcov01_mar2021"
tre = toytree.tree(treeFile)

rtre = tre.root(wildcard="dentatum")

# do some rotations to fit with geo
for i in [309,262,251,252,239,233]:
    rtre.idx_dict[i].children.reverse()
    rtre._coords.update()

# set new names
labels_updated = [namedict[i] for i in rtre.get_tip_labels()]
color_labels = []

# set color base on leaf form
for i in labels_updated:
    result = "Black"
    for key, item in colordict.items():
        if i.find(key) > -1:
            result = item
    color_labels.append(result)

# collapse weak supported nodes
# rtre = rtre.collapse_nodes(min_support=75)

# define threshold
support_value_threshold = 84

# plot the tree
canvas, axes, marks = rtre.draw(
    height=1800, width=600, 
    use_edge_lengths=True,
    tip_labels_align=True,
    tip_labels_style={"font-size": "12px"},
    tip_labels=labels_updated,
#     tip_labels_colors=color_labels,
    node_sizes=[5 if i else 0 for i in rtre.get_node_values()],
    node_colors=['black' if (i and int(i) > support_value_threshold) else 'white' for i in rtre.get_node_values('support', 1, 1)],
#     node_colors=colors,
    node_style={"stroke": "black", "stroke-width": 1},
#     node_labels="support"
    node_labels=['' if (i and int(i) > support_value_threshold) else i for i in rtre.get_node_values('support', 1, 0)],
    node_labels_style= {
        "-toyplot-anchor-shift": "10px",
        "baseline-shift": "0px",
        "text-shadow": "0.5px 0.5px #fff, -0.5px 0.5px #fff, 0.5px -0.5px #fff, -0.5px -0.5px #fff",
        "fill": "#000",
        "font-size": 8,
    },
#     node_labels="idx",
);
import toyplot.pdf
toyplot.svg.render(canvas, "./RAxML_bipartitions.fulldataset_withAyava_10scaff_mcov025_rmcov01_mar2021.svg")