R语言在树状图的末端标注物种值

R语言数据分析指南

发布于 2024-03-02 08:28:00

1050

发布于 2024-03-02 08:28:00

文章被收录于专栏：R语言数据分析指南R语言数据分析指南

欢迎关注R语言数据分析指南

❝本节来分享一个进化树与棒棒糖图结合的案例来进行系统发育可视化展示，案例主要使用phytools包+基础绘图语法来进行展示，当然也可以使用ggplot语法来实现相同的功能。更多详细内容可参考官方文档 ❞

官方链接

❝http://blog.phytools.org/2024/02/plotting-species-values-at-tips-of-tree.html ❞

library(phytools) # 加载phytools包

other attached packages:
[1] phytools_2.1-1 maps_3.4.2     ape_5.7-1

data(eel.tree) # 加载鳗鱼树数据
eel.tree<-ladderize(eel.tree,right=FALSE) # 对鳗鱼树进行阶梯化处理，使树更加紧凑
data(eel.data) # 加载鳗鱼数据
eel_size<-setNames(eel.data$Max_TL_cm,rownames(eel.data)) # 将鳗鱼的最大全长数据设置为名称
h<-max(nodeHeights(eel.tree)) # 获取树的最大节点高度
plotTree(eel.tree,ftype="off",lwd=1,direction="upwards",ylim=c(0,2*h), # 绘制鳗鱼树
  mar=c(0.1,3.1,0.1,0.1))

pp <-get("last_plot.phylo",envir=.PlotPhyloEnv) # 获取最后一次绘制的树的信息

polygon(c(0,max(pp$xx)+1,max(pp$xx)+1,0), # 绘制一个多边形作为背景
  c(1.1*h,1.1*h,1.9*h,1.9*h),border=FALSE,
  col="#F2F2F2")

lines(range(pp$xx),rep(1.1*h,2),col="black",lty="dotted") # 绘制一条虚线
rel.eel_size<-eel_size/max(eel_size)*0.75*h # 计算相对大小
segments(x0=pp$xx[1:Ntip(eel.tree)],y0=rep(1.1*h,Ntip(eel.tree)), # 绘制线段表示鳗鱼大小
  x1=pp$xx[1:Ntip(eel.tree)],y1=rel.eel_size[eel.tree$tip.label]+1.1*h)
axis(2,at=1.1*h+max(rel.eel_size)/max(eel_size)*pretty(eel_size), # 添加y轴
  labels=pretty(eel_size),las=1,cex.axis=0.6)
cols<-setNames( # 设置颜色
  viridisLite::viridis(n=100)[ceiling(100*eel_size/max(eel_size))],
  names(eel_size))
points(pp$xx[1:Ntip(eel.tree)],rel.eel_size[eel.tree$tip.label]+1.1*h, # 绘制点表示鳗鱼
  pch=21,cex=1.5,bg=cols[eel.tree$tip])

data(anoletree) # 加载变色龙树数据
anole.tree<-as.phylo(anoletree) # 将变色龙数据转换为phylo对象
data(anole.data) # 加载变色龙数据
anole_resid<-phyl.resid(anole.tree,x=as.matrix(anole.data[,"SVL",drop=FALSE]), # 计算残差
                        Y=as.matrix(anole.data[,c(6,3)]))
anole_data<-cbind(anole_resid$resid,exp(anole.data[,"SVL",drop=FALSE])) # 组合数据
h<-max(nodeHeights(anole.tree)) # 获取树的最大节点高度
plotTree(anole.tree,ftype="off",direction="upwards",ylim=c(0,4*h), # 绘制变色龙树
         mar=c(0.1,5.1,0.1,0.1),lwd=1)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv) # 获取最后一次绘制的树的信息
titles<-c("relative log(TL)","relative log(HLL)","SVL") # 设置标题
for(i in ncol(anole_data):1){ # 遍历数据列
  d<-max(c(diff(range(anole_data[,i])),max(anole_data[,i]))) # 计算距离
  x<-setNames(anole_data[,i]/d*0.8*h, # 计算相对位置
              rownames(anole.data))
  polygon(c(0,max(pp$xx)+1,max(pp$xx)+1,0), # 绘制背景多边形
          (i+c(0.05,0.05,0.95,0.95))*h,border=FALSE,
          col="#F2F2F2")
  hh<-(i+0.1)*h-min(c(0,min(x))) # 计算高度
  lines(range(pp$xx),rep(hh,2),col="black",lty="dotted") # 绘制虚线
  segments(x0=pp$xx[1:Ntip(anole.tree)],y0=rep(hh,Ntip(anole.tree)), # 绘制线段表示数据
           x1=pp$xx[1:Ntip(anole.tree)],y1=x[anole.tree$tip.label]+hh)
  labs<-pretty(anole_data[,i]) # 设置标签
  labs[!(labs>max(anole_data[,i]))]->labs
  labs[!(labs<min(anole_data[,i]))]->labs
  axis(2,at=hh+max(x)/max(anole_data[,i])*labs, # 添加y轴
       labels=labs,las=1,cex.axis=0.6)
  cols<-setNames( # 设置颜色
    viridisLite::viridis(n=100)[ceiling(99*((x-min(x))/diff(range(x))))+1],
    names(x))
  points(pp$xx[1:Ntip(anole.tree)],x[anole.tree$tip.label]+hh, # 绘制点表示数据
         pch=21,cex=1.2,bg=cols[anole.tree$tip.label])
  mtext(titles[i],2,line=3,at=mean(hh+max(x)/max(anole_data[,i])*labs), # 添加标题
        cex=0.8)
}