溫馨提示×
您好，登錄后才能下訂單哦！
點擊 登錄注冊 即表示同意《億速云用戶服務條款》
您好，登錄后才能下訂單哦！
今天就跟大家聊聊有關如何用R語言做冗余分析，可能很多人都不太了解，為了讓大家更加了解，小編給大家總結了以下內容，希望大家根據這篇文章可以有所收獲。
冗余分析（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://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)
這里又涉及到了最小等位基因頻率這個概念
MAF <- 0.05
frequencies <- colSums(geno)/(2*nrow(geno))
maf <- which(frequencies > MAF & frequencies < (1-MAF))
geno <- geno[,maf]
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)
看完上述內容，你們對如何用R語言做冗余分析有進一步的了解嗎？如果還想了解更多知識或者相關內容，請關注億速云行業資訊頻道，感謝大家的支持。
免責聲明：本站發布的內容（圖片、視頻和文字）以原創、轉載和分享為主，文章觀點不代表本網站立場，如果涉及侵權請聯系站長郵箱：is@yisu.com進行舉報，并提供相關證據，一經查實，將立刻刪除涉嫌侵權內容。
