TVP

# R语言做冗余分析（RDA）的一个简单小例子

3.8K0

• 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://github.com/Capblancq/RDA-genome-scan

###### 首先是读入数据

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

``````geno<-read.csv("sim1.csv")[,-c(1:14)]
geno[1:6,1:6]
###### 对基因型数据进行过滤

``````MAF <- 0.05
frequencies <- colSums(geno)/(2*nrow(geno))
maf <- which(frequencies > MAF & frequencies < (1-MAF))
geno <- geno[,maf]``````
###### 接下来就是RDA分析了
``````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)
{
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)
}

p2<-ggplot() +
xlab("SNPs") + ylab("-log10(p.values)") +
theme_bw()
p3<-ggplot() +
geom_point(aes(x=RDA\$CCA\$v[,1], y=RDA\$CCA\$v[,2]), col = "gray86") +
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

0 条评论

• 首先是读入数据
• 对基因型数据进行过滤
• 接下来就是RDA分析了