冗余分析（RDA）——R包vegan

从概念上讲，冗余分析（redundancy analysis, RDA）是响应变量矩阵与解释变量矩阵之间多元多重线性回归的拟合值矩阵的PCA分析，也是多响应变量（multi-response）回归分析的拓展。简单一点来说，RDA是通过线性回归分析结合主成分分析的排序方法，目的是寻找能最大程度解释响应变量矩阵变差的一些列的解释变量的线性组合，也就是环境对于样本的影响，因此RDA是被解释变量约束的排序。对于算法的详细介绍，就不在这里赘述了，可以参考以下链接：群落分析的冗余分析（RDA）概述 (qq.com)

下面介绍可以进行RDA的一个R包，生态学统计分析常用的vegan，官网：https://github.com/vegandevs/vegan

在基因组环境关联分析（GEA）中，常用到的是LFMM和RDA两种方法，但是这里要讲到的进化中的RDA和群落分析中的有一些区别，是RDA基因组扫描 RDA-based genome scan，论文：Evaluation of redundancy analysis to identify signatures of local adaptation。文章的脚本有上传在GitHub，链接：https://github.com/Capblancq/RDA-genome-scan/blob/master/Script_RDA.R

0.安装

> install.packages("remotes")
> remotes::install_github("vegandevs/vegan")
> library(vegan)

1.数据文件

RDA用到两个文件，响应变量（基因型数据、物种数据、细菌门水平丰度表等）和解释变量（环境数据：温度、湿度、海拔、有机质含量等）。

示例数据：Dryad Data -- Evaluation of redundancy analysis to identify signatures of local adaptation (datadryad.org)

sim1.csv文件中包含环境数据和基因型数据，直接读取：

> data<-read.csv("sim1.csv",header = T)

sim1.csv

第1列为分群信息，第2列到第11列为环境数据，第12列到14列是表型数据，14列之后是基因型数据。

1.1响应变量

在本贴中用到的响应数据是基因型数据，即sim1中14列之后的数据。

#提取响应数据
> geno<-data[,-c(1:14)]

geno

1.2解释变量

sim1.csv中的第2到11列为环境数据

> env<-data[,2:11]

env

2. RDA基因组扫描

> library("robustbase")
> library("robust")
> library("qvalue")
> rdadapt<-function(rda,K)
 {
   loadings<-rda$CCA$v[,1:as.numeric(K)]
   resscale <- apply(loadings, 2, scale)
   resmaha <- covRob(resscale, distance = TRUE, na.action= na.omit, estim="pairwiseGK")$dist
   lambda <- median(resmaha)/qchisq(0.5,df=K)
   reschi2test <- pchisq(resmaha/lambda,K,lower.tail=FALSE)
   qval <- qvalue(reschi2test)
   q.values_rdadapt<-qval$qvalues
   return(data.frame(p.values=reschi2test, q.values=q.values_rdadapt))
}

3.冗余分析

3.1使用全部环境变量

> rda_all <- rda(geno~., env, scale = FALSE)
# 保留RDA中四个约束轴
> res_rdadapt<-rdadapt(rda_all, 4)
> row.names(res_rdadapt)<-colnames(geno)
> head(res_rdadapt)
        p.values   q.values
X   0.0002819648 0.01566471
X.1 0.3924069078 0.95071092
X.2 0.1651672760 0.92863749
X.3 0.1539308576 0.92863749
X.4 0.3538912022 0.95071092
X.5 0.7169689587 0.98712367
# 输出SNP的p值和q值
> write.csv(res_rdadapt,"rda_pvalue.csv")

• geno：响应变量
• ~.：省略写法，表示将解释变量中所有的环境因子都作为解释变量，带入RDA
• scale：响应变量是否执行标准化。当有多个响应变量，且量纲不一致时，推荐使用标准化scale = TRUE。

结果解读：

> rda_all.scaling1 <- summary(rda_all, scaling = 1)
> rda_all.scaling1

如果对排序样方之间的距离更感兴趣，或者大多数解释变量为因子类型，则考虑I型标尺 scaling = 1；如果对变量之间的相关关系更感兴趣，则考虑II型标尺scaling = 2。

3.2使用部分环境变量

> rda_part <- rda(geno~envir1+envir3+envir5, env, scale = FALSE)

当只考虑部分环境因子时，可在~后指定环境因子，不同的因子间用+连接。

3.3考虑环境变量间的互作

> rda_inter <- rda(geno~envir1+envir3+envir5+envir3*envir5, env, scale = FALSE)

因子间的互作用*表示，envir3 * envir5即envir3和envir5两个因子间的互作。

4.结果绘制

> library("ggplot2")
> p1<- ggplot() +
      geom_line(aes(x=c(1:length(rda_all$CCA$eig)), y=as.vector(rda_all$CCA$eig)), linetype="dotted", size = 1.5, color="darkgrey") +
      geom_point(aes(x=c(1:length(rda_all$CCA$eig)), y=as.vector(rda_all$CCA$eig)), size = 3, color="darkgrey") +
      scale_x_discrete(name = "Ordination axes", limits=c(1:9)) +
      ylab("Inertia") +
      theme_bw()

4.1RDA排序图

根据是否展示物种向量，排序图可分为双序图（仅展示样方和环境变量二者关系）和三序图（展示样方、物种及环境变量三者关系）。

> plot(rda_all, scaling = 1, main = 'I 型标尺', display = c('wa', 'sp', 'cn'))

RDA三序图中将排序对象（样方），展示为“标签点”；响应变量（基因型）和解释变量（环境变量）均以向量表示。

display参数用于定义在排序图中展示哪些信息：

• wa：使用基因位点（物种）加权和计算的样方坐标，
• sp：基因位点（物种）；
• cn：约束成分（即解释变量）

RDA三序图
图形解释：群落分析的冗余分析（RDA）概述 (qq.com)

这个图该有的信息都有了，就是不够好看，美化一下。

> p2<- ggplot() +
   geom_point(aes(x=rda_all$CCA$v[,1], y=rda_all$CCA$v[,2]), col = "gray86") +
   geom_point(aes(x=rda_all$CCA$v[which(res_rdadapt[,2] < 0.1),1], y=rda_all$CCA$v[which(res_rdadapt[,2] < 0.1),2]), col = "orange") +
   geom_segment(aes(xend=rda_all$CCA$biplot[,1]/10, yend=rda_all$CCA$biplot[,2]/10, x=0, y=0), colour="black", size=0.5, linetype=1, arrow=arrow(length = unit(0.02, "npc"))) +
   geom_text(aes(x=1.2*rda_all$CCA$biplot[,1]/10, y=1.2*rda_all$CCA$biplot[,2]/10, label = colnames(env))) +
   xlab("RDA 1") + ylab("RDA 2") +
   theme_bw() +
   theme(legend.position="none")

RDA三序图

4.2Manhattan图

> p3<- ggplot() +
  geom_point(aes(x=c(1:length(res_rdadapt[,1])), y=-log10(res_rdadapt[,1])), col = "gray83") +
  geom_point(aes(x=c(1:length(res_rdadapt[,1]))[which(res_rdadapt[,2] < 0.1)], y=-log10(res_rdadapt[which(res_rdadapt[,2] < 0.1),1])), col = "orange") +
  xlab("SNPs") + ylab("-log10(p.values)") +
  theme_bw()
# 此处筛选标准为q<0.1为显著位点，可根据自己的数据进行调整

Manhattan plot

#筛选显著位点
> rda_sig_snp<-subset(res_rdadapt,res_rdadapt$q.values < 0.1)
> head(rda_sig_snp)
         p.values     q.values
X    2.819648e-04 1.566471e-02
X.10 1.497575e-06 1.069696e-04
X.20 2.065122e-08 2.065122e-06
X.30 4.986654e-04 2.624555e-02
X.40 3.810703e-05 2.381690e-03
X.50 5.882242e-05 3.460142e-03
> write.csv(rda_sig_snp,"rda_sig_snp.csv",row.names = T)

4.3图形排列

> library(ggpubr)
> ggarrange(ggarrange(p1, p2, ncol = 2, labels = c("A", "B")),
  ggarrange(p3,ncol = 1, labels = c("C")),
  nrow = 2)

引用转载请注明出处，如有错误敬请指出。