溫馨提示×

溫馨提示×

您好，登錄后才能下訂單哦！

密碼登錄×
登錄注冊×
其他方式登錄
點(diǎn)擊 登錄注冊 即表示同意《億速云用戶服務(wù)條款》

如何用R語言做冗余分析

發(fā)布時間：2021-11-22 15:44:31 來源：億速云 閱讀：651 作者：柒染 欄目：大數(shù)據(jù)

今天就跟大家聊聊有關(guān)如何用R語言做冗余分析，可能很多人都不太了解，為了讓大家更加了解，小編給大家總結(jié)了以下內(nèi)容，希望大家根據(jù)這篇文章可以有所收獲。

冗余分析（redundancy analysis, RDA）自己之前也聽過，好像是生態(tài)學(xué)研究中用的比較多，主要是用來探索環(huán)境和一些樣本指標(biāo)之間的關(guān)系。最近自己在看一些群體遺傳相關(guān)的內(nèi)容，發(fā)現(xiàn)RDA也可以用在群體遺傳方面 ，比如這個參考鏈接 https://popgen.nescent.org/2018-03-27_RDA_GEA.html 就介紹了這個分析，主要研究內(nèi)容自己還沒有看明白：大體好像是利用芯片技術(shù)測了一些狼的基因型，同時采集了狼生活地點(diǎn)的環(huán)境數(shù)據(jù)，利用RDA同時分析基因型數(shù)據(jù)和環(huán)境數(shù)據(jù)。這個看的還有些模棱兩可，還需要仔細(xì)看看。這個鏈接對應(yīng)的兩篇論文

  • https://onlinelibrary.wiley.com/doi/abs/10.1111/mec.13364

如何用R語言做冗余分析

  • https://onlinelibrary.wiley.com/doi/full/10.1111/mec.14584

如何用R語言做冗余分析

找資料的時候還找到了另外一篇論文

  • https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.12906

如何用R語言做冗余分析  

論文對應(yīng)的數(shù)據(jù)和代碼 


如何用R語言做冗余分析


https://github.com/Capblancq/RDA-genome-scan


如何用R語言做冗余分析

今天的推文重復(fù)一下這個論文里的冗余分析的代碼

 首先是讀入數(shù)據(jù)

sim1.csv這個數(shù)據(jù)集1:14列是環(huán)境數(shù)據(jù)，后面都是基因型數(shù)據(jù)

geno<-read.csv("sim1.csv")[,-c(1:14)]
env<-read.csv("sim1.csv")[,c(1:14)]
geno[1:6,1:6]
head(env)
   對基因型數(shù)據(jù)進(jìn)行過濾

這里又涉及到了最小等位基因頻率這個概念

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)
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語言做冗余分析  

看完上述內(nèi)容，你們對如何用R語言做冗余分析有進(jìn)一步的了解嗎？如果還想了解更多知識或者相關(guān)內(nèi)容，請關(guān)注億速云行業(yè)資訊頻道，感謝大家的支持。

向AI問一下細(xì)節(jié)

免責(zé)聲明：本站發(fā)布的內(nèi)容（圖片、視頻和文字）以原創(chuàng)、轉(zhuǎn)載和分享為主，文章觀點(diǎn)不代表本網(wǎng)站立場，如果涉及侵權(quán)請聯(lián)系站長郵箱：is@yisu.com進(jìn)行舉報(bào)，并提供相關(guān)證據(jù)，一經(jīng)查實(shí)，將立刻刪除涉嫌侵權(quán)內(nèi)容。

AI
助
手