此前文章中小编已经给大家介绍了微生物群落组成多样性分析的三大主要板块——Alpha多样性分析、Beta多样性分析及群落结构组成分析。然而，在微生物高通量测序数据分析过程中，分析环境因子是否对微生物的群落结构造成一定影响也是其中很重要的一个板块。所以，小编今天要给大家讲的内容就是使用RDA分析方法分析环境因子对微生物群落的影响。
RDA，即冗余分析（Redundancy Analysis），是一种响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析，说的通俗点其实就是约束化的PCA分析。通过RDA分析可以实现将样本和与其对应的环境因素同时反映在同一个图中，以此揭示出环境因子是否对样本分布或则样本微生物群落结构等造成影响。
1、加载包
rm(list=ls())#clear Global Environment
setwd('D:\\桌面\\RDA')#设置工作路径
#加载包
library(vegan)
library(ggplot2)
2、加载数据
#OTU表格
df <- read.table("otu.txt",sep="\t",header = T,row.names = 1,check.names = F)
head(df)
#环境因子数据
env <- read.table("env.txt",sep="\t",header = T,row.names = 1,check.names = F)
head(env)
image.png
image.png
3、RDA分析
1）使用decorana函数检测我们的数据是否符合做RDA的要求
print(decorana(t(df)))
#根据DCA1的Axis Lengths值进行选择，如果>4.0选CCA；如果在3.0-4.0之间，选RDA和CCA都可以；如果<3.0, 选择RDA分析即可。
image.png
2）vegan包rda函数支持RDA分析
RDA <- rda(t(df),env,scale = T)
3）提取数据
#提取数据
df_rda <- data.frame(RDA$CCA$u[,1:2],rownames(env))
colnames(df_rda)=c("RDA1","RDA2","samples")
# 提取物种得分
df_rda_score <- data.frame(RDA$CCA$v[,1:2])
#计算轴标签数据（=轴特征值/sum(所有轴的特征值)）
RDA1 =round(RDA$CCA$eig[1]/sum(RDA$CCA$eig)*100,2)
RDA2 =round(RDA$CCA$eig[2]/sum(RDA$CCA$eig)*100,2)
4、绘制RDA图
#读入分组文件
group <- read.table("group.txt", sep='\t', header=T)
#修改列名
colnames(group) <- c("samples","group")
#将绘图数据和分组合并
df_rda <- merge(df_rda,group,by="samples")
color=c("#1597A5","#FFC24B","#FEB3AE")#颜色变量
p1<-ggplot(data=df_rda,aes(x=RDA1,y=RDA2,
color=group))+#指定数据、X轴、Y轴，颜色
theme_bw()+#主题设置
geom_point(size=3,shape=16)+#绘制点图并设定大小
theme(panel.grid = element_blank())+
geom_vline(xintercept = 0,lty="dashed")+
geom_hline(yintercept = 0,lty="dashed")+#图中虚线
geom_text(aes(label=samples, y=RDA2+0.03,x=RDA1+0.03, vjust=0),size=3)+#添加数据点的标签
# guides(color=guide_legend(title=NULL))+#去除图例标题
labs(x=paste0("RDA1 (",RDA1,"%)"),
y=paste0("RDA2 (",RDA2,"%)"))+#将x、y轴标题改为贡献度
stat_ellipse(data=df_rda,
level=0.95,
linetype = 2,size=0.8,
show.legend = T)+
scale_color_manual(values = color) +#点的颜色设置
scale_fill_manual(values = c("#1597A5","#FFC24B","#FEB3AE"))+
theme(axis.title.x=element_text(size=12),#修改X轴标题文本
axis.title.y=element_text(size=12,angle=90),#修改y轴标题文本
axis.text.y=element_text(size=10),#修改x轴刻度标签文本
axis.text.x=element_text(size=10),#修改y轴刻度标签文本
panel.grid=element_blank())#隐藏网格线
p1
image.png
5、添加环境因子数据
#提取环境因子得分
df_rda_env <- RDA$CCA$biplot[,1:2]
df_rda_env <- as.data.frame(df_rda_env)
head(df_rda_env)
# 添加环境因子数据
p1+geom_segment(data=df_rda_env,aes(x=0,y=0,xend=df_rda_env[,1],yend=df_rda_env[,2]),
color="black",size=0.8,
arrow=arrow(angle = 35,length=unit(0.3,"cm")))+
geom_text(data=df_rda_env,aes(x=df_rda_env[,1],y=df_rda_env[,2],
label=rownames(df_rda_env)),size=3.5,
color="blue",
hjust="inward",
vjust=0.5*(1-sign(df_rda_env[,2])))+
theme(legend.position = "top")
image.png
6、AI美化
image.png
源码及作图数据在公众号【科研后花园】后台回复即可获取！！！
