论文
Reference genome assemblies reveal the origin and evolution of allohexaploid oat
https://www.nature.com/articles/s41588-022-01127-7#Sec31
论文中公开的数据处理流程
论文里还公布了所有图的原始数据，我们可以试着用论文中的原始数据来模仿出论文中的图
今天的推文我们来重复一下论文中的Figure3b 中的第一个树状图
image.png
ggtree所有树的布局
image.png
https://yulab-smu.top/treedata-book/chapter4.html
和论文中比较像的布局是 dayight这个布局
使用ggtree作图的时候
ggtree(tree01,layout = "daylight")+
geom_tiplab()
使用daylight这个布局一直报错
Error: C stack usage 15924720 is too close to the limit
我现在用的R是4.0.3
换成4.1版本的R就没有这个问题
读取树文件
library(ggtree)
library(ggplot2)
library(ggforce)
vert.tree<-read.tree("data/20220725/tree01.nwk")
作图的时候最方便就是直接使用ggtree
ggtree(vert.tree)
ggtree(vert.tree,layout = "daylight")+
geom_tiplab() -> p
ggtree(vert.tree,layout = "daylight")+
geom_nodelab(aes(label=node))
p+geom_nodelab(aes(label=node))
p+geom_highlight(node=15,expand=0.01)
p+geom_highlight(node=15)+
xlim(-0.15,0.1)+ylim(-0.1,0.2)
但是有些细节调整起来我还不清楚，最终出图效果如下
image.png
目前能想到的办法是
把作图数据单独提取出来，然后用ggplot2操作
ggplot_build(p)$data[[1]] -> df1
ggplot_build(p)$data[[2]] -> df2
ggplot_build(p)$data[[3]] -> df3
作图代码
ggplot()+
geom_segment(data=df1,
aes(x=x,y=y,xend=xend,yend=yend))+
geom_text(data=df2,aes(x=x,y=y,
label=label,
angle=angle,
hjust=hjust,
vjust=vjust))+
geom_text(data=df3,aes(x=x,y=y,
label=label,
angle=angle,
hjust=hjust,
vjust=vjust))+
xlim(-0.15,0.1)+ylim(-0.1,0.2)+
geom_mark_hull(data=data.frame(x=c(0.025,0.1,0.1,0.07),
y=c(0.08,0.15,0.02,0.02)),
aes(x=x,y=y),
fill="red",
color="transparent",
#con.colour="white",
expand = unit(5,'mm'),
radius = unit(15,'mm'))+
geom_mark_hull(data=data.frame(x=c(-0.05,-0.12,-0.15,-0.08,-0.13),
y=c(0.08,0.15,0.08,-0.02,0.02)),
aes(x=x,y=y),
fill="blue",
color="transparent")+
geom_mark_hull(data=data.frame(x=c(-0.025,0.02,0.075),
y=c(-0.08,0.05,0)),
aes(x=x,y=y),
fill="orange",
color="transparent",
expand = unit(0,'mm'))+
geom_mark_hull(data=data.frame(x=c(0,0,0.06),
y=c(0.125,0.2,0.14)),
aes(x=x,y=y),
fill="darkgreen",
color="transparent",
expand = unit(3,'mm')) -> p1
p1
image.png
这里添加色块用到的函数是ggforce包中的geom_mark_hull()函数，这里比较麻烦的是还需要自己手动计算色块的边界坐标，算这些坐标还挺费时间的，还有一个问题是如何给色块添加渐变色
拼图
library(patchwork)
p1+p1+theme_void()
image.png
示例数据和代码可以自己到论文中获取，或者给本篇推文点赞，点击在看，然后留言获取
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
处理论文中进化树文件遇到的报错
论文中提供的数据是excel存储的，首先把进化树的内容复制到一个文本文件里
读取树文件
library(ggtree)
read.tree("data/20220725/tree01.nwk")
遇到报错
Error in nchar(tree) : invalid multibyte string, element 1
查了一下报错，有人提到可能是字符编码问题
使用readLines（）函数读取的话可以 看到有一个字符是乱码
readLines("data/20220725/tree01.nwk")
image.png
对应着到进化树中去看发现其是用减号链接的字符，这样可能不行？
image.png
对应着把这个减号改成下划线就好了,如果确实想在图上体现这个减号，可以出图后再编辑
然后读取
read.tree("data/20220725/tree01.nwk")
又遇到报错
Error in FUN(X[[i]], ...) :
numbers of left and right parentheses in Newick string not equal
报错的意思就是进化树里的半括号数不匹配
搜索找到 参考
https://github.com/YuLab-SMU/ggtree/issues/432
有人说可能是进化树的文本标签 里有分号，我搜了一下我的没有
统计一下半括号的数量
readLines("data/20220725/Figure3b_1.txt") %>%
str_count("\\(")
readLines("data/20220725/Figure3b_1.txt") %>%
str_count("\\)")
一个13 一个14 确实不匹配
暂时想不到如何用代码去找是哪个括号多了，只能一个一个看了，还好进化树文件不多
把枝长信息去掉
readLines("data/20220725/Figure3b_1.txt") %>%
str_replace_all(":0.[0-9]+","")
把结果
((AbarCN65538_AB,(((Asat_Ogle_ACD_A,Sanfensan_ACD_A),(AlonCN58139_Al,AlonCN58138_Al)),((((AinsINS_4_CD_D,AinsCN108634_CD_D),AmarCN57945_CD_D),AmurCN21989_CD_D),(AagaCN25869_AB,(AcanCN23017_Ac,AdamCN19457_Ad)))))),(AnudCN58062_As,AstrCN88610_As),AwieCN90217_As);
放到rstudio写代码的地方，遇到逗号就换行，就能够找到多的那个右括号
但实际应该是少了一个左括号，在文件的最左边添加上就可以了
可能是在将树文件复制到excel的时候少选了一个左边的括号？
网友评论