R包vegan的基于距离的冗余分析(db-RDA)
常规的冗余分析( RDA ),样方或物种的降维排序实质上反映了欧氏距离。有时我们可能期望关注非欧式距离样方或物种关系的 RDA ,就可使用基于距离的冗余分析( Distance-based redundancy analysis , db-RDA )来解决。 db-RDA 将主坐标分析( PCoA )计算的样方得分矩阵应用在 RDA 中,其好处是可以基于任意一种距离测度(例如 Bray-curtis 距离等,而不再仅仅局限于欧氏距离)进行 RDA 排序。 因此从概念上来看, db-RDA 是一种回归分析结合主坐标(而非主成分)分析的排序方法,但仍然反映了样方与环境的线性响应关系。 db-RDA 在群落分析中应用广泛。 有关RDA、db-RDA等的基本概念描述,可参考前文 。 前篇已经简介了R包vegan的 常规 RDA 方法。鉴于db-RDA在群落分析中非常重要,本篇单独阐述它。本篇以某16S扩增子测序所得细菌群落数据,简介R包vegan的db-RDA过程。 示例文件、 R 脚本等的百度盘链接: https://pan.baidu.com/s/1itLEWnV_sLpKohHleMwggA
某研究对三种环境中的土壤进行了采样(环境A、B、C,各环境下采集9个土壤样本),结合二代测序观测了土壤细菌群落物种组成,并对多种土壤理化指标进行了测量。获得两个数据集。
(1)文件“otu_table.txt”为通过16S测序所得的细菌OTU丰度表格,下文统称为“物种”,尽管OTU并不完全等于物种。
(2)文件“env_table.txt”记录了多种土壤理化指标信息,即环境数据。
接下来期望通过db-RDA分析这些数据,用以分析群落物种组成与环境梯度的关系。
vegan包中db-RDA的初步计算
首先展示db-RDA的初步计算过程,对于内容详解见下文。
手动分步计算
按照上述db-RDA原理,不妨先从头走一遍,以便加深理解。
先计算PCoA,再将结果作为RDA的输入,最后被动添加物种得分。方法大致如下。
library(vegan)
#读入物种数据,以细菌 OTU 水平丰度表为例
otu otu
#读取环境数据
env
##根据原理一步步计算 db-RDA
#计算样方距离,以 Bray-curtis 距离为例,详情 ?vegdist
dis_bray
#或者直接使用现有的距离矩阵,这里同样为 Bray-curtis 距离
dis_bray
#PCoA 排序,这里通过 add = TRUE校正负特征值,详情 ?cmdscale
pcoa
#提取 PCoA 样方得分(坐标)
pcoa_site
#db-RDA,环境变量与 PCoA 轴的多元回归
#通过 vegan 包的 RDA 函数 rda() 执行,详情 ?rda
db_rda
#被动拟合物种得分
v.eig db_rda$CCA$v v.eig db_rda$CA$v
#先作图展示下,详情 ?plot.cca
#样方展示为点,物种暂且展示为“+”,环境变量为向量
par(mfrow = c(1, 2))
plot(db_rda, type = 'n', display = c('wa', 'cn'), choices = 1:2, scaling = 1, main = 'I型标尺,双序图')
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
plot(db_rda, type = 'n', display = c('wa', 'sp', 'cn'), choices = 1:2, scaling = 1, main = 'I型标尺,三序图')
points(db_rda, choices = 1:2, scaling = 1, display = 'sp', pch = 3, col = 'gray', cex = 1)
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
capscale()函数
即上述分步过程相比,capscale()是vegan中打包好的db-RDA执行函数,可以直接一步出结果。平常说的“CAP(Canonical analysis of principal coordinates)”就是指它,就是db-RDA。
#或者,capscale() 提供了直接运行的方法,详情 ?capscale
#输入原始数据,指定距离类型,例如样方距离以 Bray-curtis 距离为例
#计算结果中包含物种得分
db_rda
#或者直接输入已经计算/或读入的距离测度
db_rda #但是这种情况下,物种得分丢失,需要被动添加
v.eig db_rda$CCA$v v.eig db_rda$CA$v
#先作图展示下
#样方展示为点,物种暂且展示为“+”,环境变量为向量
par(mfrow = c(1, 2))
plot(db_rda, type = 'n', display = c('wa', 'cn'), choices = 1:2, scalin