之前有两篇推文介绍过这个内容
使用R语言的phytools包在进化树上标记自己测序取样的样本
如何使用R语言ggtree包在进化树上标记自己取样测序的样本
但是之前的两篇推文都没有达到自己想要的效果，这几天有点想明白了，记录一下自己的想法
首先使用ggtree读取进化树，然后画个进化树的图
library(ggtree)
read.tree("C:/Users/lenovo/Downloads/fig2.nwk") -> tree
ggtree(tree)+
geom_tiplab()
假如我想标记B H J 三个样本
首先是获取 每个tip （样本名就是tip）对应的node number，这个node number对应的和tree$tip.label样本出现的顺序，比如A 第一个出现 A对应的就是1 B是2 这样以此下来
ggtree(tree)+
geom_tiplab()+
geom_tiplab(aes(label=node),offset = 2)
phangorn::Ancestors(x=tree,node = 1,type = "all")
这个1节点对应的父节点是 20 19 18 17 16 15 14
14是根节点的node,可以用函数
ape::Ntip(tree)+1
获取这个节点的node
这个函数可以获取指定节点的所有父节点
ggtree(tree)+
geom_tiplab()+
geom_tiplab(aes(label=node),offset = 2)+
geom_nodelab(aes(label=node))
phytools::getDescendants(tree,16)
这个函数可以获取指定节点的所有子节点
如果用ggtree(tree) 画个最基本的进化树
ggtree(tree) -> tree.plot
ggplot_build(tree.plot)$data
这个data里有两个数据，第一个数据是用来画进化树横线的数据，第二个数据是用来画进化树竖线的数据
ggplot_build(tree.plot)$data[[1]] %>%
ggplot()+
geom_segment(aes(x=x,xend=xend,y=y,yend=yend)) -> p1
ggplot_build(tree.plot)$data[[2]] %>%
ggplot()+
geom_segment(aes(x=x,xend=xend,y=y,yend=yend)) -> p2
library(patchwork)
p1+p2
标记A样本的代码
ggtree(tree)+
geom_tiplab()+
geom_segment(data=phangorn::Ancestors(x=tree,node = 1,type = "all") %>%
as.data.frame() %>%
magrittr::set_colnames("node") %>%
bind_rows(data.frame(node=1)) %>%
left_join(ggplot_build(tree.plot)$data[[1]]),
aes(x=x,xend=xend,y=y,yend=yend),
color="red")+
geom_segment(data=phangorn::Ancestors(x=tree,node = 1,type = "all") %>%
as.data.frame() %>%
magrittr::set_colnames("node") %>%
bind_rows(data.frame(node=1)) %>%
left_join(ggplot_build(tree.plot)$data[[2]]),
aes(x=x,xend=xend,y=y,yend=yend),
color="red")
一次性标记多个样本
ggtree(tree)+
geom_tiplab()+
geom_segment(data = lapply(fortify(tree) %>%
filter(label%in%tips.new) %>%
pull(node), phangorn::Ancestors,x=tree,type="all") %>%
unlist() %>%
unique() %>%
as.data.frame() %>%
magrittr::set_colnames("id") %>%
bind_rows(data.frame(id=fortify(tree) %>%
filter(label%in%tips.new) %>%
pull(node))) %>%
left_join(ggplot_build(tree.plot)$data[[1]],
by=c("id"="node")),
aes(x=x,xend=xend,y=y,yend=yend),color="red")+
geom_segment(data = lapply(fortify(tree) %>%
filter(label%in%tips.new) %>%
pull(node), phangorn::Ancestors,x=tree,type="all") %>%
unlist() %>%
unique() %>%
as.data.frame() %>%
magrittr::set_colnames("id") %>%
bind_rows(data.frame(id=fortify(tree) %>%
filter(label%in%tips.new) %>%
pull(node))) %>%
left_join(ggplot_build(tree.plot)$data[[2]],
by=c("id"="node")),
aes(x=x,xend=xend,y=y,yend=yend),color="red")
fortify(tree)
可以把树转换成一个数据框，数据框里有每个tip对应的node number
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！