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 <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))
 
#读取环境数据
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
 
##根据原理一步步计算 db-RDA
#计算样方距离，以 Bray-curtis 距离为例，详情 ?vegdist
dis_bray <- vegdist(otu, method = 'bray')
 
#或者直接使用现有的距离矩阵，这里同样为 Bray-curtis 距离
dis_bray <- as.dist(read.delim('bray_distance.txt', row.names = 1, sep = '\t', check.names = FALSE))
 
#PCoA 排序，这里通过 add = TRUE校正负特征值，详情 ?cmdscale
pcoa <- cmdscale(dis_bray, k = nrow(otu) - 1, eig = TRUE, add = TRUE)
 
#提取 PCoA 样方得分（坐标）
pcoa_site <- pcoa$point
 
#db-RDA，环境变量与 PCoA 轴的多元回归
#通过 vegan 包的 RDA 函数 rda() 执行，详情 ?rda
db_rda <- rda(pcoa_site, env, scale = FALSE)
 
#被动拟合物种得分
v.eig <- t(otu) %*% db_rda$CCA$u/sqrt(nrow(otu) - 1)
db_rda$CCA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
v.eig <- t(otu) %*% db_rda$CA$u/sqrt(nrow(otu) - 1)
db_rda$CA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
 
#先作图展示下，详情 ?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 <- capscale(otu~., env, distance = 'bray', add = TRUE)
 
#或者直接输入已经计算/或读入的距离测度
db_rda <- capscale(dis_bray~., env, add = TRUE)
#但是这种情况下，物种得分丢失，需要被动添加
v.eig <- t(otu) %*% db_rda$CCA$u/sqrt(nrow(otu) - 1)
db_rda$CCA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
v.eig <- t(otu) %*% db_rda$CA$u/sqrt(nrow(otu) - 1)
db_rda$CA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
 
#先作图展示下
#样方展示为点，物种暂且展示为“+”，环境变量为向量
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)

dbrda()函数

dbrda()是vegan中打包好的db-RDA另一执行函数，直接一步出结果。dbrda()和capscale()在计算方法上存在小差异，但返回值基本是一致的。

#以及 dbrda()，详情 ?dbrda
 
#输入原始数据，指定距离类型，例如样方距离以 Bray-curtis 距离为例
db_rda <- dbrda(otu~., env, distance = 'bray', add = TRUE)
 
#或者直接输入已经计算/或读入的距离测度
db_rda <- dbrda(dis_bray~., env, add = TRUE)
 
#两种情况下，dbrda() 均默认不计算物种得分，可通过 sppscores() 添加，详情 ?sppscores
sppscores(db_rda) <- otu
 
#先作图展示下
#样方展示为点，物种暂且展示为“+”，环境变量为向量
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)

关于在db-RDA中使用欧氏距离或卡方距离

上述均以群落分析中常用的Bray-curtis距离作演示。

基于欧式距离的db-RDA与常规RDA的比较

由概念还可知，若在db-RDA中使用欧式距离，将得到和常规RDA一样的结果。

#使用原始物种多度数据执行 RDA，推荐先作个 Hellinger 转化
hel <- decostand(otu, method = 'hellinger')
 
#分别执行基于 Hellinger 转化后物种数据的 tb-RDA 以及与使用欧氏距离的 db-RDA
tb_rda_test <- rda(hel~., env, scale = FALSE)
db_rda_test <- capscale(hel~., env, distance = 'euclidean')
 
par(mfrow = c(1, 2))
plot(tb_rda_test, scaling = 1, main = 'tb-RDA三序图')
plot(db_rda_test, scaling = 1, main = 'db-RDA三序图')

可看到二者结果是一致的。

注：不要在意两张图“方向相反”，排序图的方向是没有任何意义的，以排序空间中样方、物种和环境变量的相对位置为准。

基于卡方距离的db-RDA与CCA的比较

顺便再比较一下使用卡方（χ²）距离的db-RDA和CCA的区别。

#但是基于卡方距离的 db-RDA，和 CCA 并不完全相同
cca_test <- cca(otu~., env, scale = FALSE)
 
chi <- decostand(otu, method = 'chi.square')  #先做卡方标准化，再计算欧氏距离，即可得卡方距离
db_rda_test <- capscale(chi~., env, distance = 'euclidean')
 
par(mfrow = c(2, 2))
plot(cca_test, scaling = 1, main = 'db-RDA三序图')
plot(db_rda_test, scaling = 1, main = 'CCA三序图')
plot(cca_test, display = c('wa', 'cn'), scaling = 1, main = 'db-RDA双序图')
plot(db_rda_test, display = c('wa', 'cn'), scaling = 1, main = 'CCA双序图')

二者的样方得分和环境变量得分是一致的，区别在于物种得分。

从计算过程来看，CCA是CA与多元回归的结合，因此它和基于χ²距离的db-RDA（基于χ²距离的PCoA与多元回归的结合）结果的相似性是可以预料到的。对于物种得分的不同，CCA在计算过程中用加权平均法对物种排序，而db-RDA中物种排序是在排序结束后被动加入的，方法的不同带来了差异。

db-RDA结果的初步解读

上述计算结果的返回值，存储结构同rda()的结构，解读方法也同常规RDA。对于RDA的基本概念，如变差（方差）解释程度、排序图的详解、I型或II型标尺的概念等，可参考前文。

其中可能capscale()（如上所述，俗称CAP）最为常用，下文就以它的结果为例作简介。

db-RDA结果总概

summary()的解读方式和常规RDA一致。

#db-RDA 结果解读，以 capscale() 函数结果为例，简介
db_rda <- capscale(otu~., env, distance = 'bray', add = TRUE)
 
#查看统计结果信息，以 I 型标尺为例
db_rda.scaling1 <- summary(db_rda, scaling = 1)
db_rda.scaling1

db-RDA排序图

一些简略的排序图展示，便于观测数据，排序图解读方式可参考前文。

关于更美观的可视化方案，可参考上文示例，或者文末的ggplot2方法。

#作图查看排序结果，详情 ?plot.cca
#三序图，包含 I 型标尺和 II 型标尺，样方坐标展示为使用物种加权计算的样方得分
par(mfrow = c(1, 2))
 
plot(db_rda, scaling = 1, main = 'I 型标尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling1 <- scores(db_rda, choices = 1:2, scaling = 1, display = 'sp')
#arrows(0, 0, rda_sp.scaling1[ ,1], rda_sp.scaling1[ ,2], length =  0, lty = 1, col = 'red')
 
plot(db_rda, scaling = 2, main = 'II 型标尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling2 <- scores(db_rda, choices = 1:2, scaling = 2, display = 'sp')
#arrows(0, 0, rda_sp.scaling2[ ,1], rda_sp.scaling2[ ,2], length =  0, lty = 1, col = 'red')
#隐藏物种，以 I 型标尺为例展示双序图
#比较分别使用物种加权计算的样方坐标以及拟合的样方坐标的差异
par(mfrow = c(1, 2))
plot(db_rda, scaling = 1, main = 'I 型标尺，加权', display = c('wa', 'cn'))
plot(db_rda, scaling = 1, main = 'I 型标尺，拟合', display = c('lc', 'cn'))

主要结果内容提取

若有需要，考虑在结果中将需要的信息提取出来。

#scores() 提取排序得分（坐标），以 I 型标尺为例，前四轴为例
#使用物种加权和计算的样方得分
db_rda_site.scaling1 <- scores(db_rda, choices = 1:4, scaling = 1, display = 'wa')        
#物种变量（响应变量）得分
db_rda_sp.scaling1 <- scores(db_rda, choices = 1:4, scaling = 1, display = 'sp')
#环境变量（解释变量）得分
db_rda_env.scaling1 <- scores(db_rda, choices = 1:4, scaling = 1, display = 'bp')
 
#或者在 summary() 后提取，以 I 型标尺为例，前四轴为例
db_rda.scaling1 <- summary(db_rda, scaling = 1)
#使用物种加权和计算的样方得分
db_rda_site.scaling1 <- db_rda.scaling1$site[ ,1:4]
#物种
db_rda_sp.scaling1 <- db_rda.scaling1$species[ ,1:4]
#环境
db_rda_env.scaling1 <- db_rda.scaling1$biplot[ ,1:4]
 
#若需要输出在本地
#样方
write.table(data.frame(db_rda_site.scaling1), 'db_rda_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物种
write.table(data.frame(db_rda_sp.scaling1), 'db_rda_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#环境
write.table(data.frame(db_rda_env.scaling1), 'db_rda_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
 
#不建议直接在原始数据集中提取，因为这样提取的坐标数值未经标尺缩放处理，不利于反映生物学问题
#db_rda$CCA$u[ ,1:4]
#db_rda$CCA$v[ ,1:4]
#db_rda$CCA$biplot[ ,1:4]
#RDA 基于多元回归，若对典范特征系数（即每个解释变量与每个约束轴之间的回归系数）感兴趣
#coef() 提取典范系数
rda_coef <- coef(db_rda)

db-RDA的R²校正和约束轴的置换检验

关于RDA的总变差（即响应变量矩阵的总变差能够被解释变量解释的部分，如果用比例表示，即等同于多元回归的R²，也就是约束轴的总解释量）以及各约束轴解释量等，不建议直接在原始结果中获取。初始所得的RDA模型需要经过校正和检验后才可使用，否则可能会带来较大的偏差，原因可参考前文。

校正R²及约束轴特征值

方法同常规RDA，对于校正前后R²的提取，可使用RsquareAdj()完成。

R²经过校正后通常会降低，这是必然的。

#RsquareAdj() 提取 R2，详情 ?RsquareAdj() 
r2 <- RsquareAdj(db_rda)
db_rda_noadj <- r2$r.squared #原始 R2
db_rda_adj <- r2$adj.r.squared       #校正后的 R2

对于示例数据的db-RDA排序结果，通过上文，我们已经得知了响应变量矩阵的总变差为4.111，未校正前的约束轴解释的变差总和为2.476（占比0.6024，即未校正前的R²）；那么根据校正后的R²（0.39186），可还原校正R²后的约束轴解释的变差总和约为“4.111×0.39186=1.6109”。

对于各约束轴校正后的解释量。例如，已知RDA1的解释变差占约束轴总解释变差的0.4427，由于校正后的R²为0.39186，那么可推断校正R²后的RDA1解释率为“0.4427×0.39186=0.1735”；同理，校正R²后RDA2解释率“0.2532×0.39186=0.0992”；其余约束轴的解释量以此类推。

#概念如上所述，如下为代码实现
#关于约束轴承载的特征值或解释率，应当在 R2 校正后重新计算
db_rda_exp_adj <- db_rda_adj * db_rda$CCA$eig/sum(db_rda$CCA$eig)
db_rda_eig_adj <- db_rda_exp_adj * db_rda$tot.chi

约束轴的置换检验及p值校正

我们看到本示例的约束轴解释率还算可以（以前两轴为例，解释量0.1735+0.0992=0.2727，不算很低），那么，这些约束轴是否有效呢？

对于置换检验，首先检验全模型，通过后再检验各约束轴。

#置换检验
#所有约束轴的置换检验，即全局检验，基于 999 次置换，详情 ?anova.cca
db_rda_test <- anova.cca(db_rda, permutations = 999)
 
#各约束轴逐一检验，基于 999 次置换
db_rda_test_axis <- anova.cca(db_rda, by = 'axis', permutations = 999)
 
#p 值校正（Bonferroni 为例）
db_rda_test_axis$`Pr(>F)` <- p.adjust(db_rda_test_axis$`Pr(>F)`, method = 'bonferroni')

本示例中仅前两轴显著，也就是说，第三及之后的约束轴理论上不能再用作分析了。

注：如果检验未通过，在考虑更换其它方法之前，不妨先试试剔除一些环境变量，因为也有可能是某些环境变量与群落特征之间缺乏较好的线性关系所致。

非约束残差轴的重要性确定

此外，当db-RDA约束轴承载的变差比例较低时，可能需要审视非约束残差轴。

#断棍模型和 Kaiser-Guttman 准则帮助确定重要的残差轴
pcoa_eig <- db_rda$CA$eig
n <- length(pcoa_eig)
bsm <- data.frame(j=seq(1:n), p = 0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
 
barplot(t(cbind(100 * pcoa_eig/sum(pcoa_eig), bsm$p[n:1])), beside = TRUE, main = '% 变差', col = c('orange', 'bisque'), las = 2)
abline(h = mean(100 * pcoa_eig/sum(pcoa_eig)), col = 'red')
legend('topright', c('% 特征根', '断棍模型', '平均特征根'), pch = c(15, 15, NA), col = c('orange', 'bisque', 'red'), lwd = c(NA, NA, 1), bty = 'n')

结果指示可能还有1-2个残差非约束轴是比较重要的。若有必要，可绘制前1-2个非约束轴的排序图，查看为何还有相当一部分变差未能解释。

db-RDA的变量选择（前向选择为例）

RDA中，有时解释变量太多，需要想办法减少解释变量的数量。减少解释变量的数量可能有两个并不冲突的原因：第一是寻求简约的模型，利于我们对模型的解读；第二是当解释变量过多时会导致模型混乱，例如有些解释变量之间可能存在较强的线性相关，即共线性问题，可能会造成回归系数不稳定。

下文以前向选择为例展示db-RDA的解释变量选择过程，前向选择概念可参考前文。

方差膨胀因子（Variance Inflation Factor，VIF）显示环境变量间共线性很严重。有几个变量的VIF值特别大（超过10甚至20），所以有必要剔除一些变量。

#计算方差膨胀因子，详情 ?vif.cca
vif.cca(db_rda)

db-RDA的前向选择在vegan中的具体实现过程，和常规RDA的方法一致。

#前向选择，以 ordiR2step() 的方法为例，基于 999 次置换检验，详情 ?ordiR2step
db_rda_forward_pr <- ordiR2step(capscale(otu~1, env, distance = 'bray', add = TRUE), scope = formula(db_rda), R2scope = TRUE, direction = 'forward', permutations = 999)
 
#以 db_rda 和 db_rda_forward_pr 为例，简要绘制双序图比较变量选择前后结果
par(mfrow = c(1, 2))
plot(db_rda, scaling = 1, main = '原始模型，I 型标尺', display = c('wa', 'cn'))
plot(db_rda_forward_pr, scaling = 1, main = '前向选择后，I 型标尺', display = c('wa', 'cn'))
 
#细节部分查看
summary(db_rda_forward_pr, scaling = 1)
#比较选择前后校正后 R2 的差异，详情 ?RsquareAdj
#可以看到变量选择后，尽管去除了很多环境变量，但总 R2 并未损失很多
RsquareAdj(db_rda)$adj.r.squared
RsquareAdj(db_rda_forward_pr)$adj.r.squared
 
#所有约束轴的全局检验，999 次置换，详情 ?anova.cca
db_rda_forward_pr_test <- anova.cca(db_rda_forward_pr, permutations = 999)
 
#各约束轴逐一检验，999 次置换
db_rda_forward_pr_test_axis <- anova.cca(db_rda_forward_pr, by = 'axis', permutations = 999)
 
#p 值校正（Bonferroni 为例）
db_rda_forward_pr_test_axis$`Pr(>F)` <- p.adjust(db_rda_forward_pr_test_axis$`Pr(>F)`, method = 'bonferroni')

结果显示前两轴显著，我们考虑将简约db-RDA模型的前两轴坐标提取出，以便用于后续分析或可视化。

#提取或输出变量选择后的排序坐标，以 I 型标尺为例，前两轴为例
#提取方式可参考上文，这里通过 scores() 提取
 
#使用物种加权和计算的样方得分
db_rda_forward_pr_site.scaling1 <- scores(db_rda_forward_pr, choices = 1:2, scaling = 1, display = 'wa')
write.table(data.frame(db_rda_forward_pr_site.scaling1), 'db_rda_forward_pr_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物种变量（响应变量）得分
db_rda_forward_pr_sp.scaling1 <- scores(db_rda_forward_pr, choices = 1:2, scaling = 1, display = 'sp')
write.table(data.frame(db_rda_forward_pr_sp.scaling1), 'db_rda_forward_pr_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#环境变量（解释变量）得分
db_rda_forward_pr_env.scaling1 <- scores(db_rda_forward_pr, choices = 1:2, scaling = 1, display = 'bp')
write.table(data.frame(db_rda_forward_pr_env.scaling1), 'db_rda_forward_pr_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)

此外，如前所述，尽管变量选择策略具有明显的优势，但一定切记不应盲目依赖自动选择程序在回归模型中选择相关的环境变量，即仅通过统计学手段得到的最优变量集合可能并非具有很好的生物学意义。

附带一个ggplot2可视化示例

上文已经展示了通过内置的plot()函数，可绘制排序图。

同样地，也可在导出坐标后，通过其它作图包（如ggplot2）绘制。

#以前向选择后的简约模型 db_rda_forward_pr 为例作图展示前两轴
#plot(db_rda_forward_pr) 方法可参考上文
#下面是 ggplot2
 
#提取样方和环境因子排序坐标，前两轴，I 型标尺
db_rda_forward_pr.scaling1 <- summary(db_rda_forward_pr, scaling = 1)
db_rda_forward_pr.site <- data.frame(db_rda_forward_pr.scaling1$sites)[1:2]
db_rda_forward_pr.env <- data.frame(db_rda_forward_pr.scaling1$biplot)[1:2]
 
#手动添加分组
db_rda_forward_pr.env$name <- rownames(db_rda_forward_pr.env)
db_rda_forward_pr.site$name <- rownames(db_rda_forward_pr.site)
db_rda_forward_pr.site$group <- c(rep('A', 9), rep('B', 9), rep('C', 9))
 
#计算校正 R2 后的约束轴解释率
exp_adj <- RsquareAdj(db_rda_forward_pr)$adj.r.squared * db_rda_forward_pr$CCA$eig/sum(db_rda_forward_pr$CCA$eig)
rda1_exp <- paste('RDA1:', round(exp_adj[1]*100, 2), '%')
rda2_exp <- paste('RDA1:', round(exp_adj[2]*100, 2), '%')
 
#ggplot2 作图
library(ggplot2)
 
p <- ggplot(db_rda_forward_pr.site, aes(CAP1, CAP2)) +
geom_point(aes(color = group)) +
scale_color_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.title = element_blank(), legend.key = element_rect(fill = 'transparent')) + 
labs(x = rda1_exp, y = rda2_exp) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) + 
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = db_rda_forward_pr.env, aes(x = 0, y = 0, xend = CAP1, yend = CAP2), arrow = arrow(length = unit(0.2, 'cm')), size = 0.3, color = 'blue') +
geom_text(data = db_rda_forward_pr.env, aes(CAP1 * 1.1, CAP2 * 1.1, label = name), color = 'blue', size = 3)
 
p

友情链接