昨晚看到了Y叔发布的ggtreeExtra-今天我和phylsoeq无缝对接-绘制微生物组的大圈图
写在前面
之前我绘制圈图已经有很长时间的历史了，用过ggtree，也用过graphlan，虽然graphlan非常强大，但是确实很难搞。Y叔的这个包还是很强大的。
实战
R包导入
这里尽量安装最新的R包，需要注意的地方我都进行了注释，大家参考一下还是非常容易安装成功的。
library(ggplot2)
library(dplyr)
library(magrittr)
library(RColorBrewer)
# devtools::install_github("YuLab-SMU/treeio") # 安装1.7以上版本的才能支持MicrobiotaProcess
# devtools::install_github("YuLab-SMU/MicrobiotaProcess")
# BiocManager::install("treeio")
library(MicrobiotaProcess)
library(tibble)
library(ggsignif)
library(ggtree)
library(ggtreeExtra)
library(ggplot2)
library(ggstar)
library(MicrobiotaProcess)
library(ggnewscale)
library(grid)
这里我直接使用ggClusterNet中的phyloseq对象做实例。自己可以通过readRDS导入phyloseq对象。
library(ggClusterNet)
library(phyloseq)
data(ps)
基本出图-圈图
这里注意ps对象直接选择丰度最高的前100个OTU进行展示。函数converttotreedata将注释表格转化为树文件。注意转化后名称发生的变化，增加了前缀：“st”。所以在后续的数据添加过程中，我也增加了这个标签。
alltax = ps %>%
filter_OTU_ps(100) %>%
vegan_tax() %>%
as.data.frame()
alltax$OTU = row.names(alltax)
head(alltax)
trda <- convert_to_treedata(alltax)
p <- ggtree(trda, layout="inward_circular", size=0.2, xlim=c(30,NA))
p
添加数据准备映射
这里我们模仿Y树的数据，定义的OTU标签颜色的表格为文件tippoint。使用门水平作为颜色填充。%<+% 符号定义了对ggplot数据集中添加数数据，要注意公共列相同。
这里我仿照Y叔的代码，颜色映射使用了自定义，所以如果你的数据数量不对的话，要修改，或者干脆点，直接去掉：scale_color_manual。
tax = ps %>%
filter_OTU_ps(100) %>%
vegan_tax() %>%
as.data.frame()
tippoint = data.frame(OTU = row.names(tax),Taxa = tax$Phylum,Level = "Phylum")
tippoint$OTU = paste("st__",tippoint$OTU,sep = "")
tippoint$names <- gsub("st__ASV_","",tippoint$OTU)
head(tippoint)
p <- p %<+% tippoint
p$data
tippoint$Taxa %>% unique()
p1 <- p +
geom_tippoint(
mapping=aes(
color=Taxa,
shape=Level
),
size=1,
alpha=0.8
) +
scale_color_manual(values=c("#EF3B2C", "#1D91C0", "#FEB24C", "grey60",
"#7FBC41", "#4D9221", "#276419"),
guide=guide_legend(
keywidth=0.5,
keyheight=0.5,
order=2,
override.aes=list(shape=c("Actinobacteria"=20,
"Proteobacteria" =20,
"Bacteroidetes" =20,
"Chloroflexi" =20,
"Firmicutes" =18,
"Unassigned" =18
),
size=2
),
na.translate=TRUE
)
) +
scale_shape_manual(values=c("Phylum"=20, "Class"=18), guide="none" )
p1
ggsave("cs1.pdf",p1,width = 4,height = 4)
添加环内OTU相关
这部分内环相关结合进化树的图形有一些，但是不多，所以学一学，装13还是没问题的。corMicro函数用于计算相关，甚至可以计算sparcc相关，这个函数在ggLusterNet中。这里我使用了tidyfst包的函数将矩阵转化为列表。
#------对OTU之间的关系进行连线
otu = ps %>%
filter_OTU_ps(100) %>%
vegan_otu() %>%
t() %>%
as.data.frame()
head(otu)
result = corMicro (ps = pssub ,N = 0,r.threshold=0.8,p.threshold=0.05,method = "pearson")
#--提取相关矩阵
cor = result[[1]]
diag(cor) = 0
# BiocManager::install("tidyfst")
# install.packages('Rcpp')# error: function 'Rcpp_precious_remove' not provided by package 'Rcpp'
library(tidyfst)
linktab = mat_df(cor) %>% filter(value != 0) %>%
mutate(direct = ifelse(value> 0, "positive", "nagetive"))
colnames(linktab) = c("Inhibitor","Sensitive","Interaction","direct")
head(linktab)
linktab$Inhibitor = paste("st__",linktab$Inhibitor,sep = "")
linktab$Sensitive = paste("st__",linktab$Sensitive,sep = "")
head(linktab)
p2 <- p1 +
new_scale_color() +
geom_taxalink(
data=linktab,
mapping=aes(
taxa1=Inhibitor,
taxa2=Sensitive,
color=direct
),
alpha=0.6,
offset=0.1,
size=0.15,
ncp=10,
hratio=1,
arrow=grid::arrow(length = unit(0.005, "npc"))
) +
scale_colour_manual(values=c("chocolate2", "#3690C0", "#009E73"),
guide=guide_legend(
keywidth=0.8, keyheight=0.5,
order=1, override.aes=list(alpha=1, size=0.5)
)
)
p2
ggsave("cs2.pdf",p2,width = 5,height = 5)
添加外环单个样本丰度值
提取OTU表格，然后数据宽变长，然后映射到图片上即可。
这里选择了log转化
otu$id = row.names(otu)
ringdat = otu %>% longer_dt(id)
ringdat$value = log2(ringdat$value+1)
ringdat$id = paste("st__",ringdat$id,sep = "")
head(ringdat)
library(RColorBrewer)
p3 <- p2 +
geom_fruit(
data=ringdat,
geom=geom_star,
mapping=aes(
y=id,
x=name,
size=value,
fill= name
),
starshape = 13,
starstroke = 0,
offset=-0.9,
pwidth=0.8,
grid.params=list(linetype=3)
) +
scale_size_continuous(range=c(0, 2),
limits=c(sort(ringdat$value)[2], max(ringdat$value)),
breaks=c(1, 2, 3),
name=bquote(paste(Log[2],"(",.("Count+1"), ")")),
guide=guide_legend(keywidth = 0.4, keyheight = 0.4, order=4,
override.aes = list(starstroke=0.3))
) +
scale_fill_manual(
values=c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3",brewer.pal(11,"Spectral")),
guide=guide_legend(keywidth = 0.4, keyheight = 0.4, order=3)
)
p3
ggsave("cs3.pdf",p3,width = 8,height = 8)
外环标签添加
这里是OTU的标签，我去除了ASV_这个全部都一样的字符，只留下了数字。
这部分数据就在p$data中，直接添加label即可。
p4 <- p3 +
geom_tiplab(
mapping=aes(
label=names
),
align=TRUE,
size=2,
linetype=NA,
offset=16
)
p4
ggsave("cs4.pdf",p4,width = 6,height = 6)
最外环堆叠柱状图展示分组丰度
这里我将计算的每个分组的的平均丰度进行添加。
head(otu)
dat <- pssub %>%
vegan_otu() %>% as.data.frame()
bartab = cbind(dat,sample_data(pssub)[,1:3]) %>%
group_by(Group) %>%
summarise_if(is.numeric,mean) %>%
as.data.frame()%>% longer_dt(Group)
head(bartab)
bartab$name = paste("st__",bartab$name,sep = "")
head(weighttab)
p5 <- p4 +
new_scale_fill() +
geom_fruit(
data=bartab,
geom=geom_bar,
mapping=aes(
x=value,
y=name,
fill= Group
),
stat="identity",
orientation="y",
offset=0.48,
pwidth=1.5,
axis.params=list(
axis = "x",
text.angle = -45,
hjust = 0,
vjust = 0.5,
nbreak = 4
)
) +
scale_fill_manual(
name = "Number of interactions",
values=c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3"),
guide=guide_legend(keywidth=0.5,keyheight=0.5,order=5)
) +
theme(
legend.background=element_rect(fill=NA),
legend.title=element_text(size=6.5),
legend.text=element_text(size=5),
legend.spacing.y = unit(0.02, "cm"),
legend.margin=margin(0.1, 0.9, 0.1,-0.9, unit="cm"),
legend.box.margin=margin(0.1, 0.9, 0.1, -0.9, unit="cm"),
plot.margin = unit(c(-1.2, -1.2, -1.2, 0.1),"cm")
)
p5
ggsave("cs5.pdf",p5,width = 7,height = 7)
根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。