一步一步教你使用ggtree

SYSU星空

发布于 2022-05-05 12:58:23

8K0

发布于 2022-05-05 12:58:23

ggtree是R语言中一个强大的系统发育树可视化及注释软件包，在Bioconductor中发布，同时兼有ggplot2的优点。ggtree可以读取多种格式（包括newick，nexus，NHX，jplace和phylip）的系统发育树，并结合不同类型的相关数据进行注释分析。在R中ggtree的安装方法如下：

source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")

ggtree需要依赖Bioconductor中的treeio，以及ggplot2、ggstance、ape等软件包，如果安装失败，可能是没有预先安装依赖包。作为ggplot2的拓展包，ggtree可以充分利用ggplot2来进行系统发育树的注释和美化，做出更加丰富多彩的图形。

⑴系统发育树及其注释的可视化

常用的系统发育树为newick格式，在这里我们以FastTree创建的系统发育树为例。首先我们可以简单的将系统发育树可视化：

library(ggplot2)
library(ggtree) #加载需要的软件包
tree=read.tree("top50_OTUs_with_genus.tre")
ggtree(tree, layout="rectangular", size=0.8, col="deepskyblue3") #调整展示形状（矩形）设置树枝大小以及颜色

其中layout为发育树的展示形状，有矩形"rectangular"、斜八字形"slanted"、风扇形"fan"、环形"circular"、放射形"radial"、无根树"unrooted"；open.angle为当形状为风扇形的时候调整角度；branch.length="none"则分枝末端齐平；size可以调整树枝的宽度，col可以调整树枝的颜色。在深入分析之前，我们可以将tree转换为数据框列表来查看其内容，以方便后面脚本的理解：

data=as.data.frame(tree) #或者下面命令
data=fortify(tree)

上面所做的系统发育树仍十分简略，所含信息很少，接下来我们使用ggplot2函数对其图层进行扩展注释：

tree=read.tree("top50_OTUs_with_genus.tre")
data=fortify(tree)
tregraph=ggtree(tree, layout="rectangular", size=0.8, col="deepskyblue3") +
geom_tiplab(size=3, color="purple4") + #显示物种信息，并设置颜色大小
geom_tippoint(size=2, color="deepskyblue3") + #显示物种标识，并设置颜色大小
geom_text2(aes(subset=!isTip, label=node), hjust=-0.3, size=2, color="deepskyblue4") +#显示节点支持率，并设置其位置、大小以及颜色
geom_nodepoint(color="orange", alpha=1/4, size=4) + #显示节点标识及其颜色大小，alpha值为透明度
theme_tree2() + #显示坐标轴（绝对遗传距离）
xlim(NA, max(data$x)*1.2) #调节x轴范围，使得物种信息不超出边界
tregraph #查看图形

上面脚本中geom_tiplab和geom_tippoint控制显示物种及其标记，geom_nodepoint和geom_text2控制显示节点及其节点支持率，theme_tree2控制显示x轴，xlim则调节x轴的范围，通过脚本可以看出ggplot2的语法特征，图片元素通过图层叠加的方法来进行调整。

⑵系统发育树与其他数据整合展示

除了系统发育树内置数据的注释，ggtree还可以整合其他数据进行可视化注释，接下来我们使用facet_plot函数在发育树后面绘制每个物种的序列分布柱状图，完整脚本如下：

library(ggplot2)
library(ggtree)
library(ggstance)
tree=read.tree("top50_OTUs_with_genus.tre") 
data=fortify(tree)
tregraph=ggtree(tree, layout="rectangular", size=0.8, col="deepskyblue4") + 
  geom_tiplab(size=2.5, color="gray15", align=T, linetype=3, linesize=0.5, hjust=-0.02) + #align为添加水平比对线
  geom_tippoint(size=1.5, color="deepskyblue4") + 
  geom_text2(aes(subset=!isTip, label=node), hjust=-0.3, size=2, color="black") + 
  geom_nodepoint(color="orange", alpha=1/2, size=3) + 
  theme_tree2()
table=read.table(file="top50_otu_table.txt", header=T) #读取OTU序列数据，在这里top50_otu_table.txt即为与代表序列对应的count table文件
rownames(table)=table[,1]
ndata=data[1:50,]
rownames(ndata)=ndata[,"label"] 
count=table[rownames(ndata),] #根据tree的顺序重排OTU序列数据
count1=as.vector(t(as.matrix(count[,-1]))) #将矩阵转换为列向量
otu_count=data.frame(id=rep(count$OTU_ID, each=length(colnames(table))-1), value=count1, Group=rep(colnames(table)[-1])) #将数据矩阵转换为ggplot2作图所需要的格式
otu_count$Group=factor(otu_count$Group, levels=c("Ctrl", "DO1", "DO2", "DO3")) #强制样品排列顺序
graph=facet_plot(tregraph + xlim_tree(max(data$x)*1.4), panel="Count", data=otu_count, geom=geom_barh, mapping=aes(x=value, fill=Group), stat='identity') + 
  scale_fill_manual(values = c("orangered","cyan4","green3","blue","brown")) + #设置填充颜色
  theme(legend.position="right") #显示图例并调整其位置
graph #查看图形

除了柱状图，还可以选择其他的展示方式例如箱型图（geom_boxploth()）、线图（geom_linerangeh()）等。接下来我们还可以使用gheatmap在发育树后面绘制每个物种的序列分布热图，gheatmap支持矩阵作为输入数据，完整脚本如下：

library(ggplot2)
library(ggtree)
library(ggstance)
tree=read.tree("top50_OTUs_with_genus.tre")
data=fortify(tree)
tregraph=ggtree(tree, layout="rectangular", size=0.8, col="deepskyblue4") + 
  geom_tiplab(size=2.5, color="gray15", align=T, linetype=3, linesize=0.5, hjust=-0.02) + 
  geom_tippoint(size=1.5, color="deepskyblue4") + 
  geom_text2(aes(subset=!isTip, label=node), hjust=-0.3, size=2, color="deepskyblue4") + 
  geom_nodepoint(color="orange", alpha=1/2, size=3) + 
  theme_tree2() + 
  xlim_tree(max(data$x)*1.4)
table=read.table(file="top50_otu_table.txt", header=T) 
rownames(table)=table[,1]
ndata=data[1:50,]
rownames(ndata)=ndata[,"label"] 
count=table[rownames(ndata),] 
count=count[,-1]
graph=gheatmap(tregraph, count, offset=0.25, colnames=F) %>% scale_x_ggtree #创建热图并融合两边坐标轴
graph #查看图形

上面图形仍十分粗操，接下来对图形进行调整美化，调节展示方式、颜色范围、图例位置等，完整脚本如下：

library(ggplot2)
library(ggtree)
library(ggstance)
tree=read.tree("top50_OTUs_with_genus.tre")
data=fortify(tree)
tregraph=ggtree(tree, layout="circular", size=0.8, col="deepskyblue4", branch.length="none") + #发育树环形展示，树枝末端齐平
  geom_tiplab(size=3, color="black", hjust=-0.02, offset=5.5, aes(angle=angle+300)) + #设置大的offset值使物种信息展示在热图外围，并使字体原本角度+300度旋转
  geom_tippoint(size=1.5, color="deepskyblue4") + 
  geom_text2(aes(subset=!isTip, label=node), hjust=-0.3, size=2.5, color="black") + 
  geom_nodepoint(color="orange", alpha=1/2, size=4) + 
  xlim_tree(max(data$x)*1.35)
table=read.table(file="top50_otu_table.txt", header=T) 
rownames(table)=table[,1]
ndata=data[1:50,]
rownames(ndata)=ndata[,"label"] 
count=table[rownames(ndata),] 
count=count[,-1]
graph=gheatmap(tregraph, count, offset=0.1, colnames=T, width=0.3, colnames_offset_y=0.1, font.size=2.5, low="palegreen3", high="darkorange3", colnames_angle=-45) + 
  theme(legend.position=c(0.8,0.3)) #设置最低点和最高点颜色，并调整热图的宽度，字体大小，调整图例位置正好在环状开口处
open_tree(graph, 80) %>% rotate_tree(0) #使环状图开口80度以避免热图过于稀疏，并旋转0度

⑶系统发育树内插注释图形

ggtree软件包的inset函数可以实现系统发育树节点或末端内插注释图形，从而极大丰富系统发育树的展示内容，下面我们在系统发育树tip处添加序列分布饼图，完整脚本如下所示：

library(ggplot2)
library(ggtree)
library(ggstance)
tree=read.tree("top50_OTUs_with_genus.tre")
data=fortify(tree)
tregraph=ggtree(tree, layout="rectangular", size=0.8, col="deepskyblue4", branch.length="none", ladderize=F) + 
  geom_tiplab(size=3, color="gray15", offset=5, align=T, linetype=3, linesize=0.3) + 
  geom_tippoint(size=1, color="deepskyblue4") + 
  geom_text2(aes(subset=!isTip, label=node), hjust=-0.3, size=2, color="gray15") + 
  geom_nodepoint(color="orange", alpha=1/3, size=4) + 
  theme_tree2() + 
xlim(NA, max(data$x)*35) 
table=read.table(file="top50_otu_table.txt", header=T)
rownames(table)=table[,1]
ndata=data[1:50,]
rownames(ndata)=ndata[,"label"]
count=table[rownames(ndata),]
count1=as.vector(t(as.matrix(count[,-1])))
count2=count[,-1]
size=numeric(50)
for (i in 1:50) {
size[i]=sum(count2[i,])
} #限定x轴范围
n=length(colnames(count))-1
otu=lapply(1:50, function(x) {y=count1[(n*(x-1)+1):(n*x)]}) #将矩阵转换为list系列
bar=lapply(otu, function(y) {
  rcount=data.frame(Group=colnames(count2),y=y) #将列表转换为数据框
ggplot(data=rcount, mapping=aes(y=1, x=y, fill =Group)) + 
    geom_barh(stat='identity',  width=1) + 
xlim(NA, max(size)) +
    scale_fill_manual(values = c("orangered","cyan4","green3","blue","brown")) +
    theme_inset()
}) #对每一个list作图
names(bar)=1:50 #将一系列图命名
bar1=lapply(otu[1], function(y) {
rcount=data.frame(Group=colnames(count2),y=y)
ggplot(data=rcount, mapping=aes(y=1, x=y, fill =Group)) + 
    geom_barh(stat='identity',  width=1) + 
xlim(NA, max(size)) +
    scale_fill_manual(values = c("orangered","cyan4","green3","blue","brown")) +
    theme_inset(legend.position=c(2.2,-20))
})
names(bar1)=1 #bar1主要作用是添加图例
tr=inset(tregraph, bar, width=0.4, height=0.04, hjust=-2.5, vjust=0.2) #插入bar
inset(tr, bar1, width=0.4, height=0.04, hjust=-2.5, vjust=0.2)插入legend