示例数据
提取码：szjy
1. 数据准备
rm(list = ls()) #清空当前工作环境
library('vegan')
grp = read.table('grp.txt', row.names = 1, header = T, sep = '\t')
env = grp[,-c(1:2)]
otu=read.table("OTU.txt",header=T,row.names = 1, sep = '\t') #读取因变量矩阵
library('dplyr')
otu = otu[,rownames(grp)] %>% .[rowSums(.) !=0,] %>% t() #去除均为0的行，对列排序，转置
otu.helli=decostand(otu,method = "hellinger") #标准化
2. 模型选择
decorana(otu.helli) #执行DCA去趋势分析
根据看分析结果中Axis Lengths的第一轴的大小
- 如果大于4.0,就应选CCA（基于单峰模型，典范对应分析）
- 如果在3.0-4.0之间，选RDA和CCA均可
- 如果小于3.0, RDA的结果会更合理（基于线性模型，冗余分析）
3. 构建模型
fit.0 <- rda (otu.helli ~ 1, data = env) #因变量矩阵和指示变量逐个进行RDA分析
fit.all <- rda (otu.helli ~ ., data = env) #因变量矩阵与所有环境变量进行RDA分析
fit <- ordiR2step (fit.0, scope = formula (fit.all)) #执行前向选择
若是cca()，则使用ordistep()
4. 模型检验
vif.cca(fit) # 除去值大于10的环境因子，因其与其他变量存在高度线性相关
(total = summary(fit)) # 查看RDA分析结果
summary(eigenvals(fit)) # 查看每个排序轴的特征值和能够解释的变异
#查看模型统计特征
anova.cca(fit) # 查看RDA模型显著性
anova.cca(fit, by="term", parallel=2) # 查看解释变量的显著性
anova.cca(fit, by="axis") # 查看排序轴的显著性
5. 求各个环境因子对整体变异的解释率（贡献率）
基本模型是 rda(x,y,z)或者rda(x~y+condition(z))
-
x为响应矩阵，通常为群落矩阵，必须有 -
y为限制性矩阵，通常为环境变量的矩阵，可无，则rda(x)相当于prcomp(x)执行主成分分析 -
z为条件矩阵，通常也是环境变量，意味着排除（partialed out）该条件矩阵的影响，可无，则为常规rda()或者cca()分析
# 挑选出保留下来的环境因子
subenv = env[,c("AN", "AP", "TP", "AK", "TN")]
# 新建空矩阵，用来保存结果
res = matrix(0,nrow = 1, ncol = ncol(subenv)+4,
dimnames = list(c('Inertia'),c('total','Unconstrained','constrained',names(subenv),'joint')))
res[,1] = total$tot.chi #总体方差
res[,2] = total$unconst.chi # 未解释方差
res[,3] = total$constr.chi # 当前变量所能解释的方差
for (f in names(subenv)){
of = setdiff(names(subenv),f)
temp <- rda(otu.helli, subenv[f], subenv[,of])
temp$CCA$eig -> res[,f] #各个环境因子单独解释的方差
}
res[,ncol(res)] = total$constr.chi - sum(res[,-c(1:3)]) # 多个环境因子联合解释的方差
res = res/res[,1]
# 饼图展示效果
# jpeg('pie.jpg',width = 2000, height = 2000, res = 300)
pie(res[,-1],init.angle = 60,radius = 0.5,
labels = paste(colnames(res)[-1]," ", round((as.numeric(res[,-1])),3)*100,'%'),
main = "Variance explained by each factor",
col=rainbow(8))
# dev.off()
严格来说，其中存在包含关系constrained = AN + AP + AK + TN + TP + Joint.
这里使用饼图只是为了方便展示
5. 基本绘图
plot(fit, type="n") #绘制RDA图
text(fit, dis="cn") #绘制RDA图
points(fit, pch=21, col="blue", bg="blue", cex=2) #绘制RDA图
text(fit, "sites", col='red', cex=0.5, adj = 0.5, pos=2) #绘制RDA图
可以在图窗口点击export将图复制为meta file到Power Point或者Visio中修改
6. 使用ggplot2绘图
#获取样点坐标
sites = total$sites[,1:2] %>% data.frame() %>% merge(grp[,1:2],by = 'row.names')
#获取指示变量坐标
biplot =total$biplot[,1:2] %>% data.frame()
biplot$env = rownames(biplot)
#定义样点组合在图例中的出场顺序
sites$treatment = factor(sites$treatment, levels = c('S','SH','H'))
#重命名，须注意与上面的对应关系
levels(sites$treatment) = c('水田','水旱轮作','旱田')
library('ggplot2')
library('ggrepel')
p <- ggplot(sites, aes(x = RDA1, y =RDA2, color = treatment)) +
theme_bw() +
labs( x= "RDA1 (23.42%)", y = "RDA2 (11.84%)", color = '') +
guides(color = guide_legend(override.aes = list(size=5)))+ #增加图例中点的大小
geom_hline(yintercept=0, linetype=2,color='grey') + #添加经过坐标原点的横线
geom_vline(xintercept=0, linetype=2,color='grey') + #添加经过坐标原点的纵线
geom_point(size = 3) +
stat_ellipse(show.legend = F) + #添加置信区间
geom_segment(data = biplot,
aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(length = unit(1/2, 'picas')), lwd = 1,
colour = "blue") + #添加箭头
geom_text_repel(data = biplot,
aes(x=RDA1,y=RDA2,label=env),
size= 5, fontface='bold',color='black')+ #添加指示变量文本
theme(legend.position = c(0.9,0.2),
legend.background = element_blank(),
legend.text = element_text(face = 'bold',color='black',size=12),
axis.title = element_text(face = 'bold',color='black',size=14),
axis.text = element_text(face = 'bold',color='black',size=12),
panel.grid = element_blank())
