论文
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文件
代码语言：javascript复制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
还有这个链接做参考，和上面的链接应该是同一个人写的
读取环境数据
代码语言：javascript复制X <- read.table("environmental_data.txt", h = T, stringsAsFactors = F)
X <- as.matrix(X[, c(4:14)])
读取基因型数据
代码语言：javascript复制Y <- data.table::fread(paste("birch", ".lfmm",
sep = ""), header = F)
模型拟合
tictoc这个R包可以计算R语言的命令运行的时间，环境数据一次性全部放进去，不用一个一个算
代码语言：javascript复制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个位点，这也太快了点吧
代码语言：javascript复制pv$calibrated.pvalue
存放了11个环境变量的p值
计算qvalue
代码语言：javascript复制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 链接下载
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
