一种在生态学领域玩转门特尔检验的技法
写在前面
对于生态学来讲，各种指标只有多，没有少，我们在这众多的指标中寻找对我们目标变量影响最大的指标。例如：
环境指标对微生物群落的影响？
生态学指标对生物固氮的影响？
环境指标对温室气体排放的影响？
人体血液等指标对肠道微生物稳态的影响？
植物指标对植物抵御病害能力的影响？
这些指标之间往往有相互作用，排除其他指标，寻找目标指标对目标变量的影响。
例子
下图是一个例子，最外圈是各种环境因子，四个圈代表对这些环境因子分为四类，去除每组分类指标后分别对其他指标对目标指标做相关分析。我们今天也教给大家如何绘制这样一个图形。
实战
本次数据在此：（链接失效后台留言）
链接：https://pan.baidu.com/s/1HmkeO6cqvnU2yEkgx1oXlw
提取码：qon8
library(vegan)
library(tidyverse)
library(phyloseq)
library(readxl)
data <- read_xlsx("./env.xlsx",sheet = 1) %>% as.data.frame()
head(data)
dim(data)
data$ID = NULL
env.std <- data[-30]
#-分组#--------
group = read_xlsx("./env.xlsx",sheet = 2) %>% as.data.frame()
head(group)
colnames(group) = c("Group","ID")
num = unique(group$Group) %>% length()
BC.beta = vegan::vegdist((data["env30"]), method="bray")
不去除指标计算每个环境因子目标矩阵的相关
#--不做去除
report =c()
for(i in 1:ncol(env.std)){ ##
envdis = dist(env.std)
# pmantel.BC = mantel.partial(BC.beta, envdis, envdis2, na.rm=T)
pmantel.BC = vegan::mantel(envdis, BC.beta, na.rm=T)
# pmantel.JC = mantel.partial(JC.beta, envdis, envdis2, na.rm=T)
report = rbind(report,c(colnames(env.std)[i], pmantel.BC$statistic, pmantel.BC$signif))
}
colnames(report)<- c("Envs",paste(rep(c("r","p")),rep(c("BC","JC")),sep = "."))
report = as.matrix(report)
report = as.data.frame(report)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.character)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.numeric)
report1 <- report
去除土壤组变量的矩阵计算环境变量和目标变量相关
#---仅仅去除土壤
j = 4
rem <- group$ID[group$Group %in% unique(group$Group)[2:j]]
# JC.beta = vegan::vegdist(data["CUE"], method="jaccard",binary=T)
report =c()
for(i in 1:ncol(env.std)){
##
envdis = dist(env.std)
envdis2 = dist(env.std[,rem])
pmantel.BC = mantel.partial(BC.beta, envdis, envdis2, na.rm=T)
# pmantel.JC = mantel.partial(JC.beta, envdis, envdis2, na.rm=T)
report = rbind(report,c(colnames(env.std)[i], pmantel.BC$statistic, pmantel.BC$signif))
}
colnames(report)<- c("Envs",paste(rep(c("r","p")),rep(c("BC","JC")),sep = "."))
report = as.matrix(report)
report = as.data.frame(report)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.character)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.numeric)
report2 <- report
report2
去除土壤 + 植物组变量的矩阵计算环境变量和目标变量相关
#---去除土壤+ "plant"
rem <- group$ID[group$Group %in% unique(group$Group)[3:j]]
# JC.beta = vegan::vegdist(data["CUE"], method="jaccard",binary=T)
report =c()
for(i in 1:ncol(env.std)){ ##
envdis = dist(env.std)
envdis2 = dist(env.std[,rem])
pmantel.BC = mantel.partial(BC.beta, envdis, envdis2, na.rm=T)
# pmantel.JC = mantel.partial(JC.beta, envdis, envdis2, na.rm=T)
report = rbind(report,c(colnames(env.std)[i], pmantel.BC$statistic, pmantel.BC$signif))
}
colnames(report)<- c("Envs",paste(rep(c("r","p")),rep(c("BC","JC")),sep = "."))
report = as.matrix(report)
report = as.data.frame(report)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.character)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.numeric)
report3 <- report
report3
去除土壤 + 植物 + 酶活组变量的矩阵计算环境变量和目标变量相关
#---去除土壤+ "plant" + "enzyme
rem <- group$ID[group$Group %in% unique(group$Group)[4:j]]
# JC.beta = vegan::vegdist(data["CUE"], method="jaccard",binary=T)
report =c()
for(i in 1:ncol(env.std)){ ##
envdis = dist(env.std)
envdis2 = dist(env.std[,rem])
pmantel.BC = mantel.partial(BC.beta, envdis, envdis2, na.rm=T)
# pmantel.JC = mantel.partial(JC.beta, envdis, envdis2, na.rm=T)
report = rbind(report,c(colnames(env.std)[i], pmantel.BC$statistic, pmantel.BC$signif))
}
colnames(report)<- c("Envs",paste(rep(c("r","p")),rep(c("BC","JC")),sep = "."))
report = as.matrix(report)
report = as.data.frame(report)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.character)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.numeric)
report4 <- report
report4
去除该指标外全部指标后计算环境变量和目标变量相关
#---其他全部
report =c()
for(i in 1:ncol(env.std)){ ##
envdis = dist(env.std)
envdis2 = dist(env.std[,-i])
pmantel.BC = mantel.partial(BC.beta, envdis, envdis2, na.rm=T)
# pmantel.JC = mantel.partial(JC.beta, envdis, envdis2, na.rm=T)
report = rbind(report,c(colnames(env.std)[i], pmantel.BC$statistic, pmantel.BC$signif))
}
colnames(report)<- c("Envs",paste(rep(c("r","p")),rep(c("BC","JC")),sep = "."))
report = as.matrix(report)
report = as.data.frame(report)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.character)
report[,2:dim(report)[2]]<-lapply(report[,2:dim(report)[2]],as.numeric)
report5 <- report
report5
整理数据准备出图
cor = rbind(report1$r.BC,report2$r.BC,report3$r.BC,report4$r.BC,report5$r.BC)
colnames(cor) = report1$Envs
row.names(cor) = c("zone","soil","plant","enzyme","all")
corp = rbind(report1$p.JC,report2$p.JC,report3$p.JC,report4$p.JC,report5$p.JC)
colnames(corp) = report1$Envs
row.names(corp) = c("zone","soil","plant","enzyme","all")
#-----------开始绘图#--------------
occor.r = cor
occor.p = corp
# occor.r[occor.p > p.threshold|abs(occor.r) < r.threshold] = 0
head(occor.r)
occor.r = as.data.frame(occor.r)
data <- cor %>% as.data.frame()
data$id = row.names(data)
pcm = reshape2::melt(data, id = c("id"))
pcm1 = pcm
head(group)
group2 = group %>% filter(Group != "mubiao")
head(group2)
colname = colnames(data)
rowname = row.names(data)
df0 = NULL
for (i in 1:4) {
b = group2$ID[group2$Group == unique(group2$Group)[i]]
pcmsub <- pcm %>% filter(variable %in% b)
pcmsub1 = pcmsub
if (i == 1) {
pcmsub$id = factor(pcmsub$id,levels = unique(pcmsub$id)) %>% as.numeric()
pcmsub$variable = factor(as.character(pcmsub$variable),levels = unique(pcmsub$variable)) %>% as.numeric()
# df = pcmsub
colnames(pcmsub) = c("x","y","r")
df0 = pcmsub
#---y轴标签#----
laby.data0 = data.frame(x = length(rowname),y = 1:max(df0$y),
laby = unique(pcmsub1$variable))
} else {
pcmsub$id = factor(pcmsub$id,levels = unique(pcmsub$id)) %>% as.numeric()
pcmsub$variable = factor((pcmsub$variable),levels = unique(pcmsub$variable)) %>% as.numeric()
pcmsub$variable = pcmsub$variable + 1 + max(df0$y)
colnames(pcmsub) = c("x","y","r")
df0 = rbind(df0,pcmsub)
laby.data1 = data.frame(x = length(rowname),y = min(pcmsub$y):max(pcmsub$y),laby = unique(pcmsub1$variable))
laby.data0 = rbind(laby.data0,laby.data1)
}
}
df = df0
laby.data = laby.data0
head(df)
p = ggplot() +
geom_tile(aes(x = x, y = y,fill=r),data = df,color='gray',size=0.5)+
# geom_point(aes(x = x, y = y,fill=r),df, shape=22,color='white')+
scale_size(range = c(1, 8))+
scale_fill_distiller(palette="RdYlBu") +
geom_text(aes(x = x,y = y,label = round(r,2)),data = df)
pp = ggplot() +
geom_tile(aes(x = x, y = y,fill=r),data = df,color='gray',size=0.5)+
# geom_point(aes(x = x, y = y,fill=r),df, shape=22,color='white')+
scale_size(range = c(1, 8))+
scale_fill_distiller(palette="RdYlBu") +
geom_text(aes(x = x,y = y,label = round(r,2)),data = df)
p
ggsave("cs.pdf",p,width = 8,height = 12)
p + geom_tile(aes(x = max(df$x) + 1, y = df$y ),color='black',fill = 'white',size=0.5)
p +coord_fixed(ratio = 1)
#------设置标签
data$id = NULL
n =dim(data)
# colname = colnames(data)
rowname = row.names(data)
#--x轴标签#----
labx.data = data.frame(x = 1:n[1],y = 0,labx = rowname)
p + geom_text(aes(x = x + 1 ,y = y,label = laby),data = laby.data,hjust = 0)
p = p + geom_text(aes(x = x +1,y = y,label = laby),data = laby.data,hjust = 0) +
geom_text(aes(x = x,y = y + 11,label = labx),data = labx.data) +
geom_text(aes(x = x ,y = y + 5,label = labx),data = labx.data) +
theme_void()
p
p + coord_polar(theta = "y",start = pi / 3) +
xlim(c(-15,25)) + theme_void()
### 设定角度
head(df)
ID = paste(1:max(laby.data$y))
angle1 = 90 - 360 * ( as.numeric(ID) - 0.5) /length(ID)
angle1 = angle1[laby.data$y]
library(ggtree)
df$angle = sort(rep(angle1 -90,5), decreasing =T )
p <- ggplot() +
geom_tile(aes(x = x, y = y,fill=r),data = df,color='white',size=1.5)+
# geom_point(aes(x = x, y = y,fill=r),df, shape=22,color='white')+
scale_size(range = c(1, 8))+
scale_fill_distiller(palette="RdYlBu") +
# geom_text(aes(x = x,y = y,label = round(r,2)),data = df , angle = df$angle) +
coord_polar(theta = "y",start = 0) +
xlim(c(-10,8)) + theme_void() +
geom_text(aes(x = x,y = y,label = labx),data = labx.data,hjust = 0) +
# geom_text(aes(x = x,y = y + 11,label = labx),data = labx.data,angle = -10) +
# geom_text(aes(x = x ,y = y + 5,label = labx),data = labx.data,angle = 90,hjust = 0) +
geom_text(aes(x = x + 3,y = y,label = laby),data = laby.data,hjust = 0.5, angle = angle1+ 90)
p
occor.r[occor.p > p.threshold] = 0
head(occor.r)
occor.r = as.data.frame(occor.r)
data <- occor.r
data$id = row.names(data)
pcm = reshape2::melt(data, id = c("id"))
pcm1 = pcm
colname = colnames(data)
rowname = row.names(data)
df0 = NULL
for (i in 1:4) {
b = group2$ID[group2$Group == unique(group2$Group)[i]]
pcmsub <- pcm %>% filter(variable %in% b)
pcmsub1 = pcmsub
if (i == 1) {
pcmsub$id = factor(pcmsub$id,levels = unique(pcmsub$id)) %>% as.numeric()
pcmsub$variable = factor(as.character(pcmsub$variable),levels = unique(pcmsub$variable)) %>% as.numeric()
# df = pcmsub
colnames(pcmsub) = c("x","y","r")
df0 = pcmsub
#---y轴标签#----
laby.data0 = data.frame(x = length(rowname),y = 1:max(df0$y),
laby = unique(pcmsub1$variable))
} else {
pcmsub$id = factor(pcmsub$id,levels = unique(pcmsub$id)) %>% as.numeric()
pcmsub$variable = factor((pcmsub$variable),levels = unique(pcmsub$variable)) %>% as.numeric()
pcmsub$variable = pcmsub$variable + 1 + max(df0$y)
colnames(pcmsub) = c("x","y","r")
df0 = rbind(df0,pcmsub)
laby.data1 = data.frame(x = length(rowname),y = min(pcmsub$y):max(pcmsub$y),laby = unique(pcmsub1$variable))
laby.data0 = rbind(laby.data0,laby.data1)
}
}
df0
df1 = df0
df2 <- df1 %>% filter(r != 0)
p = p + geom_tile(aes(x = x, y = y,fill=r),data = df2,color='black',size=0.5)
p <- p +
geom_text(aes(x = x,y = y,label = round(r,2)),data = df , angle = df$angle) +
geom_tile(aes(x = max(df$x)+ 2, y = df$y ),color='black',fill = 'white',size=0.5)
p
ggsave("cs5.pdf",p,width = 12,height = 12)
根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。