論文
Reference genome assemblies reveal the origin and evolution of allohexaploid oat
https://www.nature.com/articles/s41588-022-01127-7#Sec31
論文中公開的數(shù)據(jù)處理流程
論文里還公布了所有圖的原始數(shù)據(jù)赡突，我們可以試著用論文中的原始數(shù)據(jù)來模仿出論文中的圖
今天的推文我們來重復(fù)一下論文中的Figure3b 中的第一個樹狀圖
ggtree所有樹的布局
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
我現(xiàn)在用的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)
但是有些細(xì)節(jié)調(diào)整起來我還不清楚劫乱，最終出圖效果如下
目前能想到的辦法是
把作圖數(shù)據(jù)單獨(dú)提取出來，然后用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
這里添加色塊用到的函數(shù)是ggforce包中的geom_mark_hull()函數(shù)牡整，這里比較麻煩的是還需要自己手動計算色塊的邊界坐標(biāo)枫慷，算這些坐標(biāo)還挺費(fèi)時間的膛腐，還有一個問題是如何給色塊添加漸變色
拼圖
library(patchwork)
p1+p1+theme_void()
示例數(shù)據(jù)和代碼可以自己到論文中獲取遗增，或者給本篇推文點(diǎn)贊张咳，點(diǎn)擊在看帝洪，然后留言獲取
歡迎大家關(guān)注我的公眾號
小明的數(shù)據(jù)分析筆記本
小明的數(shù)據(jù)分析筆記本 公眾號 主要分享：1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子脚猾；2葱峡、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)龙助、群體遺傳學(xué)文獻(xiàn)閱讀筆記砰奕；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記提鸟！
處理論文中進(jìn)化樹文件遇到的報錯
論文中提供的數(shù)據(jù)是excel存儲的军援，首先把進(jìn)化樹的內(nèi)容復(fù)制到一個文本文件里
讀取樹文件
library(ggtree)
read.tree("data/20220725/tree01.nwk")
遇到報錯
Error in nchar(tree) : invalid multibyte string, element 1
查了一下報錯，有人提到可能是字符編碼問題
使用readLines（）函數(shù)讀取的話可以 看到有一個字符是亂碼
readLines("data/20220725/tree01.nwk")
對應(yīng)著到進(jìn)化樹中去看發(fā)現(xiàn)其是用減號鏈接的字符沽一，這樣可能不行盖溺？
對應(yīng)著把這個減號改成下劃線就好了,如果確實(shí)想在圖上體現(xiàn)這個減號，可以出圖后再編輯
然后讀取
read.tree("data/20220725/tree01.nwk")
又遇到報錯
Error in FUN(X[[i]], ...) :
numbers of left and right parentheses in Newick string not equal
報錯的意思就是進(jìn)化樹里的半括號數(shù)不匹配
搜索找到 參考
https://github.com/YuLab-SMU/ggtree/issues/432
有人說可能是進(jìn)化樹的文本標(biāo)簽 里有分號铣缠，我搜了一下我的沒有
統(tǒng)計一下半括號的數(shù)量
readLines("data/20220725/Figure3b_1.txt") %>%
str_count("\\(")
readLines("data/20220725/Figure3b_1.txt") %>%
str_count("\\)")
一個13 一個14 確實(shí)不匹配
暫時想不到如何用代碼去找是哪個括號多了烘嘱，只能一個一個看了昆禽，還好進(jìn)化樹文件不多
把枝長信息去掉
readLines("data/20220725/Figure3b_1.txt") %>%
str_replace_all(":0.[0-9]+","")
把結(jié)果
((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寫代碼的地方，遇到逗號就換行蝇庭，就能夠找到多的那個右括號
但實(shí)際應(yīng)該是少了一個左括號醉鳖，在文件的最左邊添加上就可以了
可能是在將樹文件復(fù)制到excel的時候少選了一個左邊的括號？
