微生物生态数据用ggtree进行外环注释方案一（进化树版本）

写在前面

首先感谢Y叔的ggtree。很有用且很顺手的R包。这部分的进化树加外环的方式是在Y叔实验室开发用用图层叠加方法绘制环形进化树; https://mp.weixin.qq.com/s/Il8yZqUoBVCvND7U7HKxxA 之前做的，这篇文章当时投到了SEL（中国微生物生态快报），但当时我完善了这部分的工作，目的是觉得这些操作可能在后续的工作中还会遇到。当然Y数开发的ggtreeExtra也是极好的，后面我会演示如何在微生物组数据中使用ggtreeExtra，尤其是在微生物生态过程中往往用的更多的是物种分类树，而不是物种进化树，（扩增子测序分辨率太低，进化层面上达不到，也没必要）的使用案例。

导入本次需要的R包

#-------------tax-----------Graphlan------Rraphlan----microbiomeTree------TreemicroTax----------------
library("tidyverse")
library(ggtree)
library(phyloseq)
#---导入提取otu和tax表格功能函数#----
vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}
vegan_tax <- function(physeq){
tax <- tax_table(physeq)
return(as(tax,"matrix"))
}

准备数据-从微生物组数据直接到树

基于phyloseq独享我们和容易过滤和剪枝，得到一个合适大小的进化树。

#导入数据#-----
ps = readRDS("../../data//ps_liu.rds")
# otu数量很多，所以选择一部分展示，一般树展示200个作用最为合适
ps = filter_taxa(ps, function(x) sum(x ) > 1200 , TRUE)
ps
ps = transform_sample_counts(ps, function(x) x / sum(x) )
#---提取进化树#---------
tree = phy_tree(ps)
ggtree(tree)

提取OTU表格和注释文件

#-提取otu表和tax表格#--------
otu = as.data.frame(t(vegan_otu(ps )))
tax = as.data.frame(vegan_tax(ps ))
#合并otu表格和tax表格#---------
otu_tax = merge(otu,tax,by = "row.names",all = F)
head(otu_tax)

这里使用比较老的方式将数据加入到树中

将微生物的分类门水平分类信息加入到树中，然后作为颜色设置的分组

#--枝条分组信息先不做了--------统一绘图tax树的仪一起做。#---构造分组信息文件#--------
colnames(otu_tax)[1] = c("id")
groupInfo <- split(row.names(otu_tax), otu_tax$Phylum) # OTU and phylum for group
tree <- groupOTU(tree, groupInfo)

p<- ggtree(tree,layout="fan", branch.length = "none", ladderize = FALSE,aes(color=Phylum))%<+% otu_tax +
# geom_point(aes(fill = Phylum),pch = 21,size = 3)
theme(legend.position = "right")

p = p + xlim(-10,NA)
p

进化树并不是一个暗箱，提取数据

直接从图形中提取数据这一点适合任何ggplot对象。这里我们过滤一下，只用叶节点的数据来进一步美化

#-提取可视化进化树中的数据#----
df <- p$data
# 去掉枝节点
df <- df[df$isTip, ]
head(df)
# df$label

合并OTU的平均丰度到数据中，使用ggplot图层注释到进化树上

#-----------为了更好的注释叶节点我需要构建phyloseq文件
mean = rowMeans(otu)
df = cbind(df,mean )
colnames(df)
#----叶节点控制
#----x坐标即为叶节点的坐标，我们通过调整x值达到控制整个图形的目的
p1 = p + geom_point(data = df, aes(x = x, y,size =mean,fill = Phylum ),
inherit.aes = FALSE,pch = 21,color = "black")
p1

如何拥有更多的点形状

ggsymbol是Y叔实验室的作品，这里可以用于填充更多的形状

library(ggplot2)
library(ggsymbol)
d <- data.frame(p=c(0:127),f=c(rep("g",26), rep("s", 7), rep("g", 95)))
d$f <- factor(d$f, levels=c("g", "s"))
p_sy <- ggplot() +
geom_symbol(data=d,
mapping=aes(x=p%%16, y=p%/%16, symbol=p, fill=f),
size=4, stroke=0.5, show.legend=FALSE) +
geom_text(data=d,
mapping=aes(x=p%%16, y=p%/%16+0.25, label=p),
size=3) +
scale_symbol_identity() +
scale_fill_manual(values=c("red", "blue")) +
xlab(NULL) +
ylab(NULL) +
theme(axis.ticks=element_blank(),
axis.text=element_blank())
p_sy

head(df)
df$y
df$mean
df$Phylum
p1 = p + geom_symbol(data = df, aes(x = x, y,size =mean,fill = Phylum ),
symbol = 27)
p1

用点的形状和大小双映射展示丰度信息

head(df)
# df$y
# df$mean
# df$Phylum
p1 = p + geom_symbol(data = df, aes(x = x, y,size =mean,fill = Phylum, symbol = Phylum )
)
p1

进化树叠加热图

这里我直接将OTU表格映射到了外环上，所以密密麻麻的

#-------下面的工作为优化注释#---------
#-------模仿graphlan使用形状填充节点
p$layers[[3]] = NULL
p1 = p + geom_point(data = df, aes(x = x, y,,size =mean,color = Phylum,shape = Phylum ),
inherit.aes = FALSE)
p1

#--------添加热图
mean = as.data.frame(mean)
otu
p2 = gheatmap(p1,otu, offset = 1, width=0.3, font.size=3, hjust=-.1,low = "white",high = "black")
p2

# #-按照分组添加热图-似乎热图只能使用同一套颜色银映射系统--所以后面谈添加热图只能使用geom_tile
# install.packages("ggnewscale")
library(ggnewscale)
p21 <- p2 + new_scale_fill()
p23 = gheatmap(p21,mean, offset = 3, width=0.3, font.size=3, hjust=-.1,low = "black",high = "white")
p23

控制x坐标 用于不同元素外环的信息添加

#---------外圈标记形状
p3 = p2 + geom_point(data = df, aes(x = x+10, y ),inherit.aes = FALSE,pch = 21,size = 2,fill = "red",color = "black")
p3

这里使用geom_tile填充 丰度映射演示width的用法

# library(ggcor) # + geom_star(data = df, aes(x = x+12, y ),inherit.aes = FALSE,fill = "blue",color = "black")
p4 = p3
p4
# p4 + geom_point(data = df, aes(x = x+12, y,shape = Phylum,color=Phylum ),inherit.aes = FALSE)
# 最外圈，添加矩形块，用于辨别丰度
df$width= runif(length(df$x), 0, 4)
p5 = p4 + geom_tile(data = df, aes(x = x+20, y ,width = width),color = "black",fill = "blue")
p5

# 最外圈，添加矩形块，用于辨别丰度
df$width= runif(length(df$x), 1, 5)
p5 = p4 + geom_tile(data = df, aes(x = x+15, y ,width = width),color = "black",fill = "blue")
p5

添加分割线

# 添加一条间隔线，可以设置宽度
p6 = p5 + geom_bar(data = df, aes(x+17, y =y,color = "black"),inherit.aes = FALSE,stat = "identity",width = 0.3)
p6

如何在外环添加柱状图

#geom_segment不能太小了，负责看上去怪异
p7 = p6 + geom_segment(data = df, aes(x+18, y =y,xend = x+18 +width, yend = y),color = "red",size = 5)
p7

添加标签

# 添加标签
p8 = p7 + geom_tiplab(data = df, aes(x +5, y =y,label= Phylum ,color = Phylum,angle = angle ),size=5, offset=20)
p8
#-

ggsave("./cs6.pdf",p8,width = 20,height = 18,limitsize = FALSE)

赞 (0)

相关推荐