冗余分析（redundancy analysis, RDA）自己之前也听过，好像是生态学研究中用的比较多，主要是用来探索环境和一些样本指标之间的关系。最近自己在看一些群体遗传相关的内容，发现RDA也可以用在群体遗传方面 ，比如这个参考链接 https://popgen.nescent.org/2018-03-27_RDA_GEA.html 就介绍了这个分析，主要研究内容自己还没有看明白：大体好像是利用芯片技术测了一些狼的基因型，同时采集了狼生活地点的环境数据，利用RDA同时分析基因型数据和环境数据。这个看的还有些模棱两可，还需要仔细看看。这个链接对应的两篇论文
- https://onlinelibrary.wiley.com/doi/abs/10.1111/mec.13364
- https://onlinelibrary.wiley.com/doi/full/10.1111/mec.14584
找资料的时候还找到了另外一篇论文
- https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.12906
image.png
论文对应的数据和代码
https://datadryad.org/stash/landing/show?id=doi:10.5061/dryad.1s7v5
https://github.com/Capblancq/RDA-genome-scan
今天的推文重复一下这个论文里的冗余分析的代码
首先是读入数据
sim1.csv这个数据集1:14列是环境数据，后面都是基因型数据
geno<-read.csv("sim1.csv")[,-c(1:14)]
env<-read.csv("sim1.csv")[,c(1:14)]
geno[1:6,1:6]
head(env)对基因型数据进行过滤
这里又涉及到了最小等位基因频率这个概念
代码语言：javascript复制MAF <- 0.05
frequencies <- colSums(geno)/(2*nrow(geno))
maf <- which(frequencies > MAF & frequencies < (1-MAF))
geno <- geno[,maf]接下来就是RDA分析了
代码语言：javascript复制library(vegan)
RDA <- rda(geno ~ env$envir1 env$envir2 env$envir3 env$envir4 env$envir5 env$envir6 env$envir7 env$envir8 env$envir9 env$envir10, env)
library(ggplot2)
p1<-ggplot()
geom_line(aes(x=c(1:length(RDA$CCA$eig)), y=as.vector(RDA$CCA$eig)), linetype="dotted", size = 1.5, color="darkgrey")
geom_point(aes(x=c(1:length(RDA$CCA$eig)), y=as.vector(RDA$CCA$eig)), size = 3, color="darkgrey")
scale_x_discrete(name = "Ordination axes", limits=c(1:9))
ylab("Inertia")
theme_bw()
#library(robustbase)
#install.packages("robust")
# 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))
}
res_rdadapt<-rdadapt(RDA, 5)
p2<-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()
which(res_rdadapt[,2] < 0.1)
p3<-ggplot()
geom_point(aes(x=RDA$CCA$v[,1], y=RDA$CCA$v[,2]), col = "gray86")
geom_point(aes(x=RDA$CCA$v[which(res_rdadapt[,2] < 0.1),1], y=RDA$CCA$v[which(res_rdadapt[,2] < 0.1),2]), col = "orange")
geom_segment(aes(xend=RDA$CCA$biplot[,1]/10, yend=RDA$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$CCA$biplot[,1]/10, y=1.2*RDA$CCA$biplot[,2]/10, label = colnames(env[,2:11])))
xlab("RDA 1") ylab("RDA 2")
theme_bw()
theme(legend.position="none")
library(patchwork)
p1/(p2 p3)
image.png
