MicrobiomeStatPlot | 冗余分析教程dbRDA analysis

# 基于CRAN安装R包，检测没有则安装p_list = c("vegan", "ggplot2", "ggpubr", "ggrepel", "rdacca.hp", "vegan", "psych", "reshape2")for(p in p_list){if (!requireNamespace(p)){install.packages(p)}    library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)}
# 加载R包 Load the packagesuppressWarnings(suppressMessages(library(vegan)))suppressWarnings(suppressMessages(library(ggplot2)))suppressWarnings(suppressMessages(library(ggpubr)))suppressWarnings(suppressMessages(library(ggrepel)))suppressWarnings(suppressMessages(library(rdacca.hp)))

# load data# 读入物种数据，以细菌 OTU 水平丰度表为例otu = read.csv('data/species_data.csv', head = T, row.names=1)otu <- data.frame(t(otu))# 分组数据# group datamatadata <- read.table(paste("data/group_data.txt",sep=""), header=T, row.names=1, sep="\t", comment.char="")otu = otu[rownames(otu) %in% rownames(matadata), ]# 读取环境数据# confounding factorsenv = read.csv('data/c_index_species_new4.csv', header = T, row.names=1)env = na.omit(env)env = env[rownames(env) %in% rownames(matadata), ]rownames = rownames(env)rownames = as.data.frame(rownames)otu$rownames = rownames(otu)otu = merge(otu, rownames, by = "rownames")rownames(otu) = otu$rownamesotu = otu[, -1]# 根据原理一步步计算 db-RDA# Calculate db-RDA step by step according to the principle# 计算样方距离，以 Bray-curtis 距离为例，详情 ?vegdist# Calculate the sample distance, taking Bray-curtis distance as an example, details ?vegdist# dis_bray <- vegdist(otu, method = 'bray')distance_mat = read.table(paste("data/Species_beta_diversity2.txt",sep=""), header=T, row.names=1, sep="\t", comment.char="")metadata2 = t(matadata)distance_mat2 = distance_mat[rownames(distance_mat) %in% rownames(matadata), ]distance_mat3 = distance_mat2[, colnames(distance_mat2) %in% colnames(metadata2)]distance_mat = distance_mat3# 或者直接使用现有的距离矩阵，这里同样为 Bray-curtis 距离# Or directly use the existing distance matrix, which is also the Bray-curtis distancedis_bray <- as.dist(distance_mat)# PCoA 排序，这里通过 add = TRUE校正负特征值，详情 ?cmdscale# PCoA sorting, here add = TRUE is used to correct negative eigenvalues, details ?cmdscalepcoa <- cmdscale(dis_bray, k = nrow(otu) - 1, eig = TRUE, add = TRUE)# 提取 PCoA 样方得分（坐标）# Extract PCoA sample scores (coordinates)pcoa_site <- pcoa$point# db-RDA，环境变量与 PCoA 轴的多元回归# db-RDA, multiple regression of environmental variables and PCoA axes# 通过 vegan 包的 RDA 函数 rda() 执行，详情 ?rda# Execute via the RDA function rda() of the vegan package, details ?rdadb_rda <- rda(pcoa_site, env, scale = FALSE)# summary(db_rda)# 被动拟合物种得分# Passive fitting species scorev.eig <- t(otu) %*% db_rda$CCA$u/sqrt(nrow(otu) - 1)db_rda$CCA$v <- decostand(v.eig, 'normalize', MARGIN = 2)v.eig <- t(otu) %*% db_rda$CA$u/sqrt(nrow(otu) - 1)db_rda$CA$v <- decostand(v.eig, 'normalize', MARGIN = 2)# 提取画图所需要的数据# Extract the data needed for drawingscore = scores(db_rda)# score$sitesdb_rda$CCA$biplot#>               RDA1        RDA2       RDA3#> status  0.76440409 -0.48475784  0.4250838#> sex     0.98352818  0.17911294 -0.0243079#> age    -0.01213581 -0.06874258 -0.9975606# score$speciesCAP1 = score$sites[,1]CAP2 = score$sites[,2]seg = as.data.frame(db_rda$CCA$biplot)CPA_data = as.data.frame(score$sites)CPA_data$group = rownames(CPA_data)CPA_data$group = as.character(CPA_data$group)CPA_data$group = gsub("[0-9]","", CPA_data$group)# 将上述提取的数据整合成画图所需的表格# Integrate the above extracted data into the table required for drawingplotdata = data.frame(rownames(score$sites), CAP1, CAP2, CPA_data$group)colnames(plotdata) = c('sample','CAP1','CAP2','Group') #为其加上列名(add column name)# write.csv(plotdata,'data/data5/dbRDA_npc_KO9.csv')# 计算第一、二主轴的解释度# Calculate the explanation of the first and second principal axesCAP1_exp = round(db_rda$CCA$eig[1]/sum(db_rda$CCA$eig)*100,2)CAP2_exp = round(db_rda$CCA$eig[2]/sum(db_rda$CCA$eig)*100,2)p1 = ggplot(plotdata, aes(CAP1, CAP2)) +  geom_point(aes(fill = Group, color = Group),size = 1.3) +   scale_fill_manual(values = c("#74add1","#a60026"))+  scale_color_manual(values = c("#74add1","#a60026"))+  #添加置信区间，两种选一个就行  #add confidence intervals  #stat_ellipse(linetype = 1,level = 0.95,aes(group = Treatment, color = Treatment)) +  #stat_chull(geom = 'polygon', aes(group = Treatment, color = Treatment, fill = Treatment), alpha = 0.1) +  xlab(paste('dbRDA1 ( ',CAP1_exp,'%',' )', sep = '')) +   ylab(paste('dbRDA2 ( ',CAP2_exp,'%',' )', sep = '')) +  geom_segment(data = seg, aes(x = 0, y = 0, xend = seg[,1], yend = seg[,2]),               colour = "red", size = 0.3,               arrow = arrow(angle = 30, length = unit(0.4, 'cm'))) +  geom_text_repel(data = seg, segment.colour = 'black',                  aes(x = seg[,1], y = seg[,2],                       label = rownames(seg)),size = 3) +  geom_vline(aes(xintercept = 0), linetype = 'dotted') +  geom_hline(aes(yintercept = 0), linetype = 'dotted') +  theme_bw()+  theme(text = element_text(family = 'sans', size = 12),        panel.grid.major = element_blank(),        panel.grid.minor = element_blank()        ,legend.position = 'none'        )+   theme(axis.text = element_text(colour = 'black',size = 12))+  stat_ellipse(#data=plotdata,               #geom = "polygon",               level=0.95,               linetype = 2,size=0.7,               aes(color=Group),               alpha=0.8)+  coord_equal(ratio=0.9)ggsave(paste("results/dbRDA_Species01.pdf",".pdf", sep=""), p1, width=149 * 1.5, height=80 * 1.5, unit='mm')#p1# 置换检验# Permutation testenv$status = as.numeric(env$status)env$sex = as.vector(env$sex)# env$region = as.vector(env$region)env$age = as.vector(env$age)envfit <- envfit(db_rda, env, permutations  = 999)r <- as.matrix(envfit$vectors$r)p <- as.matrix(envfit$vectors$pvals)env.p <- cbind(r,p)colnames(env.p) <- c("r2","p-value")KK <- as.data.frame(env.p)KK$p.adj = p.adjust(KK$`p-value`, method = 'BH')# KK# write.csv(as.data.frame(env.p),file="data/rdaenvfit_HP_new_species_01.csv")# 利用rdacca.hp计算每个环境因子的效应# Use rdacca.hp to calculate the effect of each environmental factordis_bray <- as.dist(distance_mat)bray = dis_braycap.hp = rdacca.hp(bray, env, method = 'dbRDA', type = 'R2', scale = FALSE)cap.hp$Total_explained_variation#> [1] 0.204cap.hp$Hier.part#>        Unique Average.share Individual I.perc(%)#> status 0.0760       -0.0009     0.0751     36.81#> sex    0.0788        0.0026     0.0814     39.90#> age    0.0539       -0.0063     0.0476     23.33# 导出每个环境因子的解释度# Export the explanation of each environmental factor# write.csv(cap.hp$Hier.part, 'data/env_effect_HP_new_Species_02.csv')# Unique:每个解释变量单独解释的总变化量的比例，即与其它解释变量不存在共存解释的部分# Unique: The proportion of the total variation explained by each explanatory variable alone, that is, the part that does not coexist with other explanatory variables# Average.share: 对每个解释变量与其它解释变量共同解释部分的分割# Average.share: Split each explanatory variable into its common explanation with other explanatory variables# Individual: 每个解释变量站总变化量的比例Individual等于Unique和Average.share的总和# Individual: The proportion of the total variation of each explanatory variable station is equal to the sum of Unique and Average.share# I.prec(%): 每个解释变量占总被解释变化量的比例，即各解释变量的Individual占Individual总和的百分比，Individual总和等同db-RDA的总解释方差占比# I.prec(%): The proportion of each explanatory variable to the total explained variance, that is, the percentage of each explanatory variable's individual to the individual sum, and the individual sum is equivalent to the total explained variance proportion of db-RDA

# load data# 读取数据# 读入物种数据，以细菌 OTU 水平丰度表为例otu = read.csv('data/pathway_unstrafied4.csv', head = T, row.names=1)           #读取物种相对丰度表otu <- data.frame(t(otu))# 分组数据# group datamatadata <- read.table(paste("data/group_data.txt",sep=""), header=T, row.names=1, sep="\t", comment.char="")otu = otu[rownames(otu) %in% rownames(matadata), ]# 读取环境数据# confounding factorsenv = read.csv('data/c_index_species_new4.csv', header = T, row.names=1)env = na.omit(env)env = env[rownames(env) %in% rownames(matadata), ]rownames = rownames(env)rownames = as.data.frame(rownames)otu$rownames = rownames(otu)otu = merge(otu, rownames, by = "rownames")rownames(otu) = otu$rownamesotu = otu[, -1]# Bray-curtis distance# Pathwaysdistance_mat = read.table(paste("data/bray_curtis_path3.txt",sep=""), header=T, row.names=1, sep="\t", comment.char="")metadata2 = t(matadata)distance_mat2 = distance_mat[rownames(distance_mat) %in% rownames(matadata), ]distance_mat3 = distance_mat2[, colnames(distance_mat2) %in% colnames(metadata2)]distance_mat = distance_mat3# 或者直接使用现有的距离矩阵，这里同样为 Bray-curtis 距离# Or directly use the existing distance matrix, which is also the Bray-curtis distancedis_bray <- as.dist(distance_mat)# PCoA 排序，这里通过 add = TRUE校正负特征值，详情 ?cmdscale# PCoA sorting, here add = TRUE is used to correct negative eigenvalues, details ?cmdscalepcoa <- cmdscale(dis_bray, k = nrow(otu) - 1, eig = TRUE, add = TRUE)# 提取 PCoA 样方得分（坐标）# Extract PCoA sample scores (coordinates)pcoa_site <- pcoa$point# db-RDA，环境变量与 PCoA 轴的多元回归# db-RDA, multiple regression of environmental variables and PCoA axes# 通过 vegan 包的 RDA 函数 rda() 执行，详情 ?rda# Execute via the RDA function rda() of the vegan package, details ?rdadb_rda <- rda(pcoa_site, env, scale = FALSE)# summary(db_rda)# 被动拟合物种得分# Passive fitting species scorev.eig <- t(otu) %*% db_rda$CCA$u/sqrt(nrow(otu) - 1)db_rda$CCA$v <- decostand(v.eig, 'normalize', MARGIN = 2)v.eig <- t(otu) %*% db_rda$CA$u/sqrt(nrow(otu) - 1)db_rda$CA$v <- decostand(v.eig, 'normalize', MARGIN = 2)# 提取画图所需要的数据# Extract the data needed for drawingscore = scores(db_rda)# score$sitesdb_rda$CCA$biplot#>             RDA1        RDA2        RDA3#> status 0.8611765  0.47808765 -0.17264755#> sex    0.6074361 -0.01168604 -0.79428262#> age    0.1040567 -0.99010545  0.09414559# score$speciesCAP1 = score$sites[,1]CAP2 = score$sites[,2]seg = as.data.frame(db_rda$CCA$biplot)CPA_data = as.data.frame(score$sites)CPA_data$group = rownames(CPA_data)CPA_data$group = as.character(CPA_data$group)CPA_data$group = gsub("[0-9]","", CPA_data$group)# 将上述提取的数据整合成画图所需的表格# Integrate the above extracted data into the table required for drawingplotdata = data.frame(rownames(score$sites), CAP1, CAP2, CPA_data$group)colnames(plotdata) = c('sample','CAP1','CAP2','Group') #为其加上列名(add column name)# 计算第一、二主轴的解释度# Calculate the explanation of the first and second principal axesCAP1_exp = round(db_rda$CCA$eig[1]/sum(db_rda$CCA$eig)*100,2)CAP2_exp = round(db_rda$CCA$eig[2]/sum(db_rda$CCA$eig)*100,2)p2 = ggplot(plotdata, aes(CAP1, CAP2)) +  geom_point(aes(fill = Group, color = Group),size = 1.3) +   scale_fill_manual(values = c("#74add1","#a60026"))+  scale_color_manual(values = c("#74add1","#a60026"))+  #添加置信区间，两种选一个就行  #add confidence intervals  #stat_ellipse(linetype = 1,level = 0.95,aes(group = Treatment, color = Treatment)) +  #stat_chull(geom = 'polygon', aes(group = Treatment, color = Treatment, fill = Treatment), alpha = 0.1) +  xlab(paste('dbRDA1 ( ',CAP1_exp,'%',' )', sep = '')) +   ylab(paste('dbRDA2 ( ',CAP2_exp,'%',' )', sep = '')) +  geom_segment(data = seg, aes(x = 0, y = 0, xend = seg[,1], yend = seg[,2]),               colour = "red", size = 0.3,               arrow = arrow(angle = 30, length = unit(0.4, 'cm'))) +  geom_text_repel(data = seg, segment.colour = 'black',                  aes(x = seg[,1], y = seg[,2],                       label = rownames(seg)),size = 3) +  geom_vline(aes(xintercept = 0), linetype = 'dotted') +  geom_hline(aes(yintercept = 0), linetype = 'dotted') +  theme_bw()+  theme(text = element_text(family = 'sans', size = 12),        panel.grid.major = element_blank(),        panel.grid.minor = element_blank()        ,legend.position = 'none'        )+   theme(axis.text = element_text(colour = 'black',size = 12))+  stat_ellipse(#data=plotdata,               #geom = "polygon",               level=0.95,               linetype = 2,size=0.7,               aes(color=Group),               alpha=0.8)+  coord_equal(ratio=0.4)ggsave(paste("results/dbRDA_PATHs01.pdf",".pdf", sep=""), p2, width=149 * 1.5, height=80 * 1.5, unit='mm')#p2# 置换检验# Permutation testenv$status = as.numeric(env$status)env$sex = as.vector(env$sex)env$age = as.vector(env$age)envfit <- envfit(db_rda, env, permutations  = 999)r <- as.matrix(envfit$vectors$r)p <- as.matrix(envfit$vectors$pvals)env.p <- cbind(r,p)colnames(env.p) <- c("r2","p-value")KK <- as.data.frame(env.p)KK$p.adj = p.adjust(KK$`p-value`, method = 'BH')# KK# write.csv(as.data.frame(env.p),file="data/rdaenvfit_HP_new_species_01.csv")# db-RDA，环境变量与 PCoA 轴的多元回归# db-RDA, multiple regression of environmental variables and PCoA axes# 通过 vegan 包的 RDA 函数 rda() 执行，详情 ?rda# Execute the RDA function rda() of the vegan package, see ?rda for details# db_rda <- rda(pcoa_site~., env, scale = FALSE)# summary(db_rda)# anova(db_rda, by = "term", permutations=999)# RsquareAdj(rda(pcoa_site ~ status, data=env))$adj.r.squared# RsquareAdj(rda(pcoa_site ~ status+sex, data=env))$adj.r.squared# RsquareAdj(rda(pcoa_site ~ status+sex+age, data=env))$adj.r.squared# 利用rdacca.hp计算每个环境因子的效应# Use rdacca.hp to calculate the effect of each environmental factor# bray = vegdist(otu, method = 'bray')# bray = dis_bray#dis_bray <- as.dist(distance_mat)bray = dis_braycap.hp = rdacca.hp(bray, env, method = 'dbRDA', type = 'R2', scale = FALSE)cap.hp$Total_explained_variation#> [1] 0.08cap.hp$Hier.part#>        Unique Average.share Individual I.perc(%)#> status 0.0381        0.0046     0.0427     53.38#> sex    0.0078        0.0082     0.0160     20.00#> age    0.0254       -0.0037     0.0217     27.12# 导出每个环境因子的解释度# Export the explanation of each environmental factorwrite.csv(cap.hp$Hier.part, 'results/env_effect_HP_new_paths_04.csv')

# Load dataallnew <- read.csv('data/data_heat_4.csv', header = T)
allnew[which(allnew$p<0.001),'sig'] <- '***'allnew[which(allnew$p<0.01 & allnew$p>=0.001),'sig'] <- '**'allnew[which(allnew$p<0.05 & allnew$p>=0.01),'sig'] <- '*'
allnew = allnew %>%  mutate(variable = ordered(variable,                         levels=c("Species", "Pathways")))
# Plotp3 <- ggplot(allnew, aes(variable,names, fill = deviance_explained)) +  geom_tile(aes(width = 1, height = 1), size = 12) +  geom_text(aes(label=deviance_explained), color="black", size=4) +   geom_text(aes(label=sig), color="black", size=8,vjust = 1.6) +   scale_fill_gradient2(low = "#afaa80", mid = "#e4cac6", high = "#5196d5",limit=c(0,1), breaks = seq(0, 0.25, 1)) +  theme(    panel.grid = element_blank(),    panel.background = element_blank(),    legend.title = element_text(angle = -90, hjust = 0.2),    axis.text.y = element_text(size=12, colour="black", family = "sans", angle = 0),    axis.text.x = element_text(size=12, colour="black", family = "sans", angle = 0, hjust = 0)  )+  theme(text = element_text(family = 'sans', size = 12),        panel.grid.major = element_blank(),        panel.grid.minor = element_blank()        #,legend.position = 'none'        )+  scale_color_discrete(breaks = c(0, 1.0, 0.25))+  theme(legend.direction='vertical')+  labs(x=NULL,y=NULL) +   guides(fill=guide_colorbar(barheight = 18,title = "Variance explained (%)", title.position = "right"))+  coord_equal(ratio = 1.7)ggsave(paste("results/dbRDA_HP_heat_n01_only_sex_age_new0318_3",".pdf", sep=""), p3, width=79 * 1.5, height=100 * 1.5, unit='mm')#p3

library(cowplot)p0 = plot_grid(p1, p2, p3, labels = c("A", "B", "C"), ncol = 3)ggsave("results/dbRDA_plot02.pdf", p0, width=179 * 1.5, height=100 * 1.5, units = "mm")