R包vegan的群落去趋势对应分析（DCA）

在对应分析（Correspondence analysis，CA）中，有时可能会产生“弓形效应（arch effect）”，对排序的精度产生影响，此时可通过去趋势对应分析（Detrended Correspondence Analysis，DCA）来解决。

弓形特征和DCA的基本原理描述，可参考前文。

本文结合微生物群落研究中的16S扩增子测序数据，简介R包vegan中的DCA实现方法。

示例文件、R脚本的百度盘链接：

https://pan.baidu.com/s/1qMAvJmmpv-vGug9EPtBjdA

CA中产生弓形效应

网盘附件“otu_table.txt”为OTU水平的物种丰度表。

在一开始，期望通过CA排序，比较数据中各样方的群落组成差异。然而，发现排序图中弓形效应严重。

library(vegan)

#读取 OTU 丰度表
otu <- read.delim('otu_table.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))

#CA 排序
ca_otu <- cca(otu)

ordiplot(ca_otu, display = 'site', type = 'n', choices = c(1, 2), scaling = 1)
points(ca_otu, display = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 10), rep('orange', 10), rep('green3', 10)), col = NA)
text(ca_otu, display = 'site', choices = c(1, 2), scaling = 1, cex = 0.8)

DCA排序

然后期望通过DCA，消除CA中的弓形效应。

DCA执行

vegan中，可通过decorana()函数执行DCA。

decorana()中，ira = 0，即执行去趋势对应分析（DCA）；若ira = 1，则执行常规的对应分析（CA）。不过不推荐通过decorana()作CA，因为无法选择标尺。

#继续使用上述的 OTU 丰度表
#DCA 排序，详情 ?decorana
dca_otu <- decorana(otu, ira = 0) #注，ira = 0 执行 DCA，ira = 1 则执行 CA

dca_otu

summary(dca_otu)
#注，summary() 中默认展示的样方和物种得分，均为标准化后的值

直接输入排序结果“dca_otu”回车后，屏幕输出如下。

Eigenvalues为主要DCA轴的特征值，默认展示前四轴。

Axis lengths反映了各DCA轴的梯度长度（gradient length）信息。所谓梯度长度，前篇提到，在DCA中，坐标轴缩放为物种更替标准差单位（standard deviation units of species turnover，简称SD），可用于表征环境异质性或物种β多样性特征。通常，SD≤3可表明物种沿排序轴更有可能是线性分布；SD≥4表明物种沿排序轴更有可能是单峰分布。

我们还可借助DCA对环境梯度的评估结果，确定选择PCA还是CA。DCA排序默认产生4个轴，我们取其中的最大值。一般经验，对应上文描述，该值小于3推荐选PCA（代表线性模型），大于4推荐选CA（代表单峰模型），介于3-4之间选PCA或CA均可。

提取主要内容

提取DCA轴特征值，计算各轴解释率，以及获取样方排序坐标等。

#各 DCA 轴特征值
dca_otu$evals
#各 DCA 轴解释量
dca_otu$evals/sum(dca_otu$evals)

#样方得分
dca_site <- dca_otu$rproj #原始坐标，不同于 summary() 显示的坐标
head(dca_site)

dca_site_scale <- scale(dca_otu$rproj, scale = FALSE) #标准化后的坐标，等同于 summary() 显示的坐标
head(dca_site_scale)
#write.table(dca_site_scale, 'dca_site_scale.txt', col.names = NA, sep = '\t', quote = FALSE)

#物种得分
dca_sp <- dca_otu$cproj #原始坐标，不同于 summary() 显示的坐标
head(dca_sp)

dca_sp_scale <- scale(dca_otu$cproj, scale = FALSE) #标准化后的坐标，等同于 summary() 显示的坐标
head(dca_sp_scale)
#write.table(dca_sp_scale, 'dca_sp_scale.txt', col.names = NA, sep = '\t', quote = FALSE)

#提取各轴“梯度长度”
#各 DCA 轴的梯度长度，即通过该轴中，所有样方坐标的最大值减去最小值获得
apply(dca_otu$rproj, 2, max) - apply(dca_otu$rproj, 2, min)

注意，“dca_otu$rproj”和“dca_otu$cproj”中存储的坐标数据，与上述通过summary(dca_otu)查看的坐标数值是不一致的。原因在于，“dca_otu$rproj”和“dca_otu$cproj”中的是未经过标准化之前的数值，而summary(dca_otu)默认展示的是标准化后的数值。若想让二者保持一致，可在提取“dca_otu$rproj”或“dca_otu$cproj”中数值时通过scale()执行标准化处理。

对比先前的PCA或CA，可以发现此处DCA中没有I型或II型标尺选择一说。猜测这个标准化可能就类似标尺缩放的过程吧。

排序图

根据计算好的排序坐标，作图就可以了。

#简洁版排序图
par(mfrow = c(1, 2))
plot(dca_otu, display = 'sites', main = '仅样方')
plot(dca_otu, main = '双序图')    #物种太多时，还是算了
#ggplot2 作图
#提取样方坐标，前两轴
dca_site_scale <- data.frame(scale(dca_otu$rproj, scale = FALSE))[1:2]

#添加分组信息
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)

dca_site_scale$sample <- rownames(dca_site_scale)
dca_site_scale <- merge(dca_site_scale, group, by = 'sample')

#ggplot2 作图
library(ggplot2)
library(ggrepel) #用于 geom_text_repel() 添加标签

p <- ggplot(dca_site_scale, aes(DCA1, DCA2)) +
geom_point(aes(color = group1)) +
scale_color_manual(values = c('green3', 'red', 'orange')) +
theme(panel.grid.major = element_line(color = 'gray', size = 0.2), panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent'), legend.title = element_blank()) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_text_repel(aes(label = sample, color = group1), size = 3, box.padding = unit(0.3, 'lines'), show.legend = FALSE)

p

#ggsave('DCA.pdf', p, width = 6, height = 5)
ggsave('DCA.png', p, width = 6, height = 5)

比较CA与DCA的差异

排序图出来后即可得知，相对于一开始的CA结果，这里DCA明显消除了弓形效应。我们将两种排序结果放在一起比较下。

#上述已经计算好了 CA 和 DCA，直接做图
par(mfrow = c(1, 2))

#CA
ordiplot(ca_otu, display = 'site', choices = c(1, 2), scaling = 1, type = 'n', main = 'CA')
points(ca_otu, display = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 10), rep('orange', 10), rep('green3', 10)), col = NA)
text(ca_otu, display = 'site', choices = c(1, 2), scaling = 1, cex = 0.8)

#DCA
ordiplot(dca_otu, display = 'sites', choices = c(1, 2), type = 'n', main = 'DCA')
points(dca_otu, display = 'site', choices = c(1, 2), pch = 21, bg = c(rep('red', 10), rep('orange', 10), rep('green3', 10)), col = NA)
text(dca_otu, display = 'site', choices = c(1, 2), cex = 0.8)

友情链接