跟着Molecular Ecology学数据分析：R语言lfmm包做环境数据和变异数据的关联分析

用户7010445

发布于 2023-11-27 17:05:59

3670

发布于 2023-11-27 17:05:59

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

论文

Genomic signals of local adaptation and hybridization in Asian white birch

https://onlinelibrary.wiley.com/doi/abs/10.1111/mec.16788

论文中提供的数据和代码的链接

https://github.com/GabrieleNocchi/betula-platyphylla-local-adaptation

我下载下来在电脑上存储的文件夹是20231024

一个小知识点，linux系统解压rar文件

mamba install unrar
unrar x birch.rar

之前有一篇推文介绍了LEA这个R包做环境数据和变异数据的关联分析，如果数据量比较大的话运行速度是非常慢的。之前推文的链接是跟着Nature Communications学数据分析：R语言LEA包做变异位点和环境数据的关联分析

https://bookdown.org/hhwagner1/LandGenCourse_book/WE_11.html

这个链接里提到了这个问题。因为LEA这个R包里用到的模型是MCMC

The previous version of LFMM (v1.5, implemented in the LEA package) uses an MCMC (Markov chain Monte Carlo) algorithm to identify GEAs while correcting for confounding. MCMC made it (very!) time-intensive for large data sets. LFMM v.2 computes LFMMs for GEA using a least-squares estimation method that is substantially faster than v1.5 (Caye et al., 2019).

https://popgen.nescent.org/2018-03-27_RDA_GEA.html

还有这个链接做参考，和上面的链接应该是同一个人写的

读取环境数据

X <- read.table("environmental_data.txt", h = T, stringsAsFactors = F)
X <- as.matrix(X[, c(4:14)])

读取基因型数据

Y <- data.table::fread(paste("birch", ".lfmm",
                                 sep = ""), header = F)

模型拟合

tictoc这个R包可以计算R语言的命令运行的时间，环境数据一次性全部放进去，不用一个一个算

library(tictoc)
tictoc()
K=3
mod.lfmm <- lfmm::lfmm_ridge(Y = Y,
                                 X = X,
                                 K = K)
pv <- lfmm::lfmm_test(Y = Y, X = X, lfmm = mod.lfmm,
                          calibrate = "gif")

运算时间就200秒

总共是1380000个位点，这也太快了点吧

pv$calibrated.pvalue

存放了11个环境变量的p值

计算qvalue

pvalues <- pv$calibrated.pvalue
my_qvalue <- function(x) {
      q <- qvalue::qvalue(x)
      q <- q$qvalues
      return(q)
    }

qvalues <- apply(pvalues, 2, my_qvalue)

这里用到的数据和代码都可以在https://github.com/GabrieleNocchi/betula-platyphylla-local-adaptation 链接下载

欢迎大家关注我的公众号

小明的数据分析笔记本