也做了挺多次RDA分析，自己现在小结一下RDA分析流程：
image.png
就我个人而言，虚线前面都是不太经历的步骤，我一般不会主动删去样品的环境信息，因为我接触的菌群这块本来就没有什么多余的环境信息-_-||，所以我的重点放在怎么去除多余OTU或菌群上面。
一般而言，我首先会做一次差异分析，挑选有差异的OTU或菌群进行展示（phyloseq推荐使用DESeq2和edgeR，详见Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible），这里不是重点不在赘述。
但是差异OTU或菌群还有可能太多，RDA呈现出来密密麻麻的，调也得调好久，最后还是好不美观。
偶然间，发现envfit不仅可以评估环境因子的显著性，也可以评估物种的相关性和显著性，这为我们进一步去取冗余物种提供了条件，值得记录下来学习。
示例：
library(BiodiversityR) # also loads vegan
library(ggplot2)
library(ggrepel)
data(dune)
data(dune.env)
attach(dune.env)
# script generated by the BiodiversityR GUI from the constrained ordination menu
dune.Hellinger <- disttransform(dune, method='hellinger')
#rda
model <- rda(dune.Hellinger ~ Management, data=dune.env, scaling="species")
model
# vif.cca(model)
# summary(model)
plot2 <- ordiplot(model, choices=c(1,2),type = "text")
默认图，中间点过于密集
#提取坐标，ggplot重绘
sites <- sites.long(plot2)
head(sites) #样本坐标
species <- species.long(plot2)
head(species) #物种坐标
centroids <- centroids.long(sites,Management,centroids.only = T)
centroids # 环境因子坐标
axis <- axis.long(model)
axis #坐标轴分数
ggplot() +
geom_vline(xintercept = c(0), color = "grey", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey", linetype = 2) +
xlab(axis$label[1]) +
ylab(axis$label[2]) +
geom_point(data=sites, aes(x=axis1, y=axis2, colour=Management, shape=Management), size=5) +
geom_point(data=species, aes(x=axis1, y=axis2))+
theme_bw()
ggplot2重绘，中间物种点还是密集
#去除冗余物种信息
(spec.envfit <- envfit(plot2, env=dune.Hellinger))#评估物种相关性和显著性
spec.data.envfit <- data.frame(r=spec.envfit$vectors$r, p=spec.envfit$vectors$pvals)#提取物种相关信息
species <- species.long(plot2, spec.data=spec.data.envfit)#提取坐标
species <- species[species$r >= 0.6 & species$p<0.05, ]#筛选
species
ggplot() +
geom_vline(xintercept = c(0), color = "grey", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey", linetype = 2) +
xlab(axis$label[1]) +
ylab(axis$label[2]) +
geom_point(data=sites, aes(x=axis1, y=axis2, colour=Management, shape=Management), size=5) +
geom_point(data=species, aes(x=axis1, y=axis2))+
theme_bw()
中间点明显变少
#加上物种点的信息，再优化一下图层
#同时添加一些可调节的细节设置
a=1
b=1
ellipse.type="norm"
ellipse.level=0.95
ggplot()+
geom_vline(xintercept = 0, color = 'grey', size = 0.6,linetype=2) +
geom_hline(yintercept = 0, color = 'grey', size = 0.6,linetype=2) +
#add site scores
geom_point(data=sites,aes(x=axis1,y=axis2,color=Management),size=3)+#Management是环境因子，实际上可以改为样本原先分组信息
# ggrepel::geom_text_repel(data=sites,size=3,aes(x = axis1, y = axis2, label = labels))+
# stat_ellipse(data=sites,type= ellipse.type,level= ellipse.level,aes(x=axis1,y=axis2,color=Management))+
scale_color_manual(values=ggpubr::get_palette("lancet",4))+
# add species scores
geom_point(data=species,aes(x=axis1*a,y=axis2*a),color=ggpubr::get_palette("npg",5)[5])+
# geom_segment(data = species,color=get_palette("npg",5)[5],size=0.8,aes(x = 0,y = 0, xend = axis1*a, yend =axis2*a),arrow = grid::arrow(length = grid::unit(0.2, "cm")))+ # add variables
ggrepel::geom_text_repel(data=species,size=3,aes(x = axis1*a, y = axis2*a, label = labels))+
# add environmental arrows
geom_segment(data = centroids,color=ggpubr::get_palette("jco",5)[1],size=0.8,aes(x = 0,y = 0, xend = axis1c*b, yend =axis2c*b),arrow = grid::arrow(length = grid::unit(0.2, "cm")))+ # add variables
ggrepel::geom_text_repel(data=centroids,size=4,color=ggpubr::get_palette("jco",5)[1],segment.alpha =0,aes(x = axis1c*b, y = axis2c*b, label = Centroid))+
xlab(axis$label[1]) +
ylab(axis$label[2]) +
ggtitle("RDA plot")+
theme_bw(base_size = 14)+
theme(aspect.ratio=0.8,
plot.title = element_text(size = 12,hjust = 0.5),
legend.title = element_text(size = 12,colour = "black"),
legend.text = element_text(size = 10,colour = "black"),
axis.title =element_text(size=12,colour = "black"),
axis.text = element_text(size=10,colour = "black"),
axis.line = element_line(color = "black"),
axis.ticks.length= unit(-0.25, 'cm'), #设置刻度向内
)
最终图像
