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

從概念上講，冗余分析（redundancy analysis, RDA）是響應(yīng)變量矩陣與解釋變量矩陣之間多元多重線性回歸的擬合值矩陣的PCA分析，也是多響應(yīng)變量（multi-response）回歸分析的拓展。簡單一點來說，RDA是通過線性回歸分析結(jié)合主成分分析的排序方法，目的是尋找能最大程度解釋響應(yīng)變量矩陣變差的一些列的解釋變量的線性組合，也就是環(huán)境對于樣本的影響，因此RDA是被解釋變量約束的排序。對于算法的詳細(xì)介紹，就不在這里贅述了，可以參考以下鏈接：群落分析的冗余分析（RDA）概述 (qq.com)

下面介紹可以進行RDA的一個R包，生態(tài)學(xué)統(tǒng)計分析常用的vegan，官網(wǎng)：https://github.com/vegandevs/vegan

在基因組環(huán)境關(guān)聯(lián)分析（GEA）中，常用到的是LFMM和RDA兩種方法，但是這里要講到的進化中的RDA和群落分析中的有一些區(qū)別，是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.數(shù)據(jù)文件

RDA用到兩個文件，響應(yīng)變量（基因型數(shù)據(jù)、物種數(shù)據(jù)、細(xì)菌門水平豐度表等）和解釋變量（環(huán)境數(shù)據(jù)：溫度、濕度、海拔、有機質(zhì)含量等）。

示例數(shù)據(jù)：Dryad Data -- Evaluation of redundancy analysis to identify signatures of local adaptation (datadryad.org)

sim1.csv文件中包含環(huán)境數(shù)據(jù)和基因型數(shù)據(jù)，直接讀?。?/p> 

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

第1列為分群信息，第2列到第11列為環(huán)境數(shù)據(jù)，第12列到14列是表型數(shù)據(jù)，14列之后是基因型數(shù)據(jù)。

1.1響應(yīng)變量

在本貼中用到的響應(yīng)數(shù)據(jù)是基因型數(shù)據(jù)，即sim1中14列之后的數(shù)據(jù)。

#提取響應(yīng)數(shù)據(jù)
> geno<-data[,-c(1:14)]
geno

1.2解釋變量

sim1.csv中的第2到11列為環(huán)境數(shù)據(jù)

> 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使用全部環(huán)境變量

> 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：響應(yīng)變量
  • ~.：省略寫法，表示將解釋變量中所有的環(huán)境因子都作為解釋變量，帶入RDA
  • scale：響應(yīng)變量是否執(zhí)行標(biāo)準(zhǔn)化。當(dāng)有多個響應(yīng)變量，且量綱不一致時，推薦使用標(biāo)準(zhǔn)化scale = TRUE。

結(jié)果解讀：

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

如果對排序樣方之間的距離更感興趣，或者大多數(shù)解釋變量為因子類型，則考慮I型標(biāo)尺 scaling = 1；如果對變量之間的相關(guān)關(guān)系更感興趣，則考慮II型標(biāo)尺scaling = 2。


3.2使用部分環(huán)境變量

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

當(dāng)只考慮部分環(huán)境因子時，可在~后指定環(huán)境因子，不同的因子間用+連接。

3.3考慮環(huán)境變量間的互作

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

因子間的互作用*表示，envir3 * envir5即envir3和envir5兩個因子間的互作。

4.結(jié)果繪制

> 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排序圖

根據(jù)是否展示物種向量，排序圖可分為雙序圖（僅展示樣方和環(huán)境變量二者關(guān)系）和三序圖（展示樣方、物種及環(huán)境變量三者關(guān)系）。

> plot(rda_all, scaling = 1, main = 'I 型標(biāo)尺', display = c('wa', 'sp', 'cn'))

RDA三序圖中將排序?qū)ο螅臃剑?，展示為“?biāo)簽點”；響應(yīng)變量（基因型）和解釋變量（環(huán)境變量）均以向量表示。

display參數(shù)用于定義在排序圖中展示哪些信息：

  • wa：使用基因位點（物種）加權(quán)和計算的樣方坐標(biāo)，
  • 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()
# 此處篩選標(biāo)準(zhǔn)為q<0.1為顯著位點，可根據(jù)自己的數(shù)據(jù)進行調(diào)整
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)

引用轉(zhuǎn)載請注明出處，如有錯誤敬請指出。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容