????此前我們已經(jīng)講過了如何使用R語言計算兩列數(shù)據(jù)相關(guān)性的分析方法坟冲,今天,我們來看一種檢驗兩個矩陣相關(guān)關(guān)系的方法——Mantel test方法曲横,這種分析方法可用在植物微生物群落與生態(tài)環(huán)境之間相關(guān)性的檢驗以及人體腸道微生物群落與人體疾病相關(guān)性檢驗等領(lǐng)域增炭。話不多說,直接邁入正文墓卦!
加載包
1)設(shè)置工作目錄
rm(list=ls())#clear Global Environment
setwd('D:\\桌面\\mantel test分析')#設(shè)置工作路徑
2)安裝、加載R包
#安裝包
# install.packages("ggplot2")
# install.packages("vegan")
# install.packages("dplyr")
# install.packages("devtools")
# devtools::install_github("houyunhuang/ggcor")
#加載包
library(vegan)
library(dplyr)
library(ggcor)
library(ggplot2)
加載數(shù)據(jù)
#OTU表格
df <- read.table("otu.txt",sep="\t",header = T,row.names = 1,check.names = F)
df <-data.frame(t(df))
#環(huán)境因子數(shù)據(jù)
env <- read.table("env.txt",sep="\t",header = T,row.names = 1,check.names = F)
head(df)
head(env)
環(huán)境因子的相關(guān)性分析及展示
quickcor(env, type = "lower",method = "spearman") +
geom_square()+
scale_fill_gradient2( high = 'orange', mid = 'white',low = 'navyblue') #顏色設(shè)置
Mantel test分析
1)計算OTU與環(huán)境因子之間的mantel test的r值和p值
df_mantel <- mantel_test(df, env, mantel.fun = 'mantel',
spec.dist.method = 'bray',
env.dist.method = 'euclidean',
spec.select = list(A = 1:4,
B = c(6,8,9,10),
C = c(5,7,11,12)))#將群落數(shù)據(jù)按組進行分開
2)定義標簽
df_mantel <- df_mantel %>%
mutate(df_r = cut(r, breaks = c(-Inf, 0.1, 0.2, 0.4, Inf),
labels = c("< 0.1", "0.1 - 0.2", "0.2 - 0.4", ">= 0.4")),#定義Mantel的R值范圍標簽
df_p = cut(p.value, breaks = c(-Inf, 0.01, 0.05, Inf),
labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))#定義Mantel的P值范圍標簽
繪圖
quickcor(env,method = "spearman", type = "upper", cor.test = T, cluster.type = "all") +#環(huán)境因子之間的相關(guān)性熱圖
geom_square() +#相關(guān)性顯示形式
geom_mark(r = NA,sig.thres = 0.05, size = 3.5, colour = "black")+#顯著性標簽
scale_fill_gradient2( high = 'red', mid = 'white',low = 'blue') + #顏色設(shè)置
anno_link(df_mantel, aes(color = df_p,
size = df_r))+
scale_size_manual(values = c(0.5, 1, 1.5, 2))+#連線粗細設(shè)置
scale_color_manual(values = c("green","blue","orange"))+#線條顏色設(shè)置
guides(fill = guide_colorbar(title = "correlation", order = 1),#圖例相關(guān)設(shè)置
size = guide_legend(title = "Mantel's r",order = 2),
color = guide_legend(title = "Mantel's p", order = 3),
linetype = "none")
參考:
1)https://zhuanlan.zhihu.com/p/507384776
2)https://zhuanlan.zhihu.com/p/110501058