冗余分析(redundancy analysis, RDA)自己之前也听过，好像是生态学研究中用的比较多，主要是用来探索环境和一些样本指标之间的关系。最近自己在看一些群体遗传相关的内容，发现RDA也可以用在群体遗传方面

，比如这个参考链接 https://popgen.nescent.org/2018-03-27_RDA_GEA.html 就介绍了这个分析，主要研究内容自己还没有看明白：大体好像是利用芯片技术测了一些狼的基因型，同时采集了狼生活地点的环境数据，利用RDA同时分析基因型数据和环境数据。这个看的还有些模棱两可，还需要仔细看看。这个链接对应的两篇论文

找资料的时候还找到了另外一篇论文

image.png

image.png

image.png

今天的推文重复一下这个论文里的冗余分析的代码

首先是读入数据

sim1.csv这个数据集1:14列是环境数据，后面都是基因型数据

geno

env

geno[1:6,1:6]

head(env)

对基因型数据进行过滤

这里又涉及到了最小等位基因频率这个概念

MAF

frequencies

maf MAF & frequencies < (1-MAF))

geno

接下来就是RDA分析了

library(vegan)

RDA

library(ggplot2)

p1

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

{

loadings

resscale

resmaha

lambda

reschi2test

qval

q.values_rdadapt

return(data.frame(p.values=reschi2test, q.values=q.values_rdadapt))

}

res_rdadapt

p2

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

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

代码能够跑完，但是具体是什么意思还不太明白，仔细看看这些论文

欢迎大家关注我的公众号

小明的数据分析笔记本