在自然群体（区别于强人工选择）中，如果我们感兴趣的数量性状表现出与特定的地理环境变量有高度的关联性，随着环境变量的改变而变化，则这些环境变量往往反映了环境作用于个体表型的选择性压力，并最终反映在群体水平的遗传统计量，进而可以探究连续环境变量梯度下不同的选择压力以及群体的适应模式。
单变量潜在因素线性混合模型（LFMM）可用于基因组中环境适应特征的筛选，可以理解为另一种类型的全基因组关联分析（GWAS），通过鉴定与环境变量高度相关的SNP位点来研究环境适应的模式。LFMM可以检测环境变量（或表型值）与等位基因频率之间的显著性关系，从而确定可能受气候因子（或表型）影响的SNP位点，专注于整体环境的评估，而不是温度、湿度、纬度等单一的变量。
基于C语言的LFMM1.5和R包LEA都是基于贝叶斯自举方法实现潜在因子混合模型，最新的LFMM则是通过最小二乘估计来实现的，即R包lfmm。“lfmm”R包对于较大的数据量计算速度更快，有时比贝叶斯版本更准确，可运用于全基因组（genome-wide，GW），基因组-环境（genome-environment，GE），全表观基因组（epigenome-wide，EW）的关联分析中。
LFMM官网：http://membres-timc.imag.fr/Olivier.Francois/lfmm/index.htm
R包lfmm：https://bcm-uga.github.io/lfmm/index.html
说明书：http://membres-timc.imag.fr/Olivier.Francois/lfmm/files/manual.pdf
R包LEA：https://bioconductor.org/packages/release/bioc/html/LEA.html
0.安装
lfmm的安装需要devtools，没有devtools的需要先安装，已安装的可忽略第一步。
> install.packages("devtools")
> devtools::install_github("bcm-uga/lfmm")
> library(lfmm)
1.数据文件
lfmm用到两个数据文件，分别为分子数据（可以是SNP基因型、基因表达水平、标准甲基化数据等）和环境数据（表型值等），要注意两个数据文件都只支持数值。
1.1分子数据
在说明书中提供了几种基因型数据可以选择，这里只介绍最方便的ped格式，plink直接将vcf转出ped，无需手动修改数据，可以避免人工修改数据的错误，格式转换详细参考以下链接：
plink安装及基础用法 - 简书 (jianshu.com)
genotype data
随后删除ped文件的前六列，只保留基因型数据：
> Y<-read.table("branch.hic.id.ped")
> Y<-Y[,-c(1:6)]
1.2环境数据
环境数据中包含n行和D列，每一行为一个个体，每一列为一个环境，[2,3]即表示第二个个体在第三个环境下的表型值。这里要注意环境数据与基因型数据中的个体顺序必须一致。
environmental data
> X<-read.csv("phe.csv",header = F,stringsAsFactors = F)
> head(X)
V1 V2
1 68.33 78.83
2 55.33 66.40
3 71.20 83.17
4 55.83 64.33
5 76.50 85.17
2.参数估计
lfmm主要估计两个参数，λ和K。
2.1 λ
λ是基因组膨胀因子（genomic inflation factor），用来量化总体膨胀程度和检测结果过高的假阳性率，相当于一个平衡参数，在实际运行过程中，可以对拟合的LFMM可以进行λ的校正，获得标准化P值。在说明书中，默认λ=1e-05。
lfmm_ridge()
2.2 K
K是潜在因素的数目，有三种方法可以估计：
-
a.可以通过祖先成分推断的结果来确定，即使用STRUCTURE和Admixture等软件估计出的祖先成分；
5. GWAS：群体结构——Admixture - 简书 (jianshu.com) -
b.基于PCA方法选择合理K值；
6. GWAS：主成分分析——GCTA - 简书 (jianshu.com) -
c. 通过交叉验证的方法，设定备选的K值列表和lambda值列表，通过lfmm_ridge_CV函数返回的prediction选择合适的K值和lambda值，该方法需要较长的运算时间和较多的计算资源。
> errs <- lfmm_ridge_CV(Y = Y,
X = X,
n.fold.row = 5,
n.fold.col = 5,
lambdas = c(1e-10, 1,1e+10),
Ks = c(1,2,3,4,5,6))
=== params
lambda K
1 1e-10 1
=== params
lambda K
2 1 1
=== params
lambda K
3 1e+10 1
……
Y：分子数据（基因型）
X：环境数据（表型）
n.fold.row：行交叉验证折叠数
n.fold.col：列交叉验证折叠数
lambdas：备选λ的范围值
Ks：回归模型中潜在因素数目
这里要注意，准确的λ值并不是我们要关注的重点，要点是控制好潜在的混合效应，得到正确的P值，因此不用过于纠结。
3.LFMM泛岭估计
3.1 lfmm_ridge()
lfmm_ridge()计算估计B、U、V，分别为SNP位点效应的大小、潜在变量得分矩阵和潜在变量加载矩阵：
> mod.lfmm <- lfmm_ridge(Y = Y,
+ X = X,
+ K = 3)
Y：分子数据（基因型）
X：环境数据（表型）
K：admixture/structure/PCA得到的潜在因素数目，示例数据中K=3
3.2 lfmm_test()
随后使用拟合模型进行关联测试
> pv <- lfmm_test(Y = Y,
X = X,
lfmm = mod.lfmm,
calibrate = "gif")
3.3 QQ图与曼哈顿图展示显著性
> pvalues <- pv$calibrated.pvalue
> qqplot(rexp(length(pvalues), rate = log(10)),
-log10(pvalues), xlab = "Expected quantile",
pch = 19, cex = .4)
> abline(0,1)
QQ plot
引用转载请注明出处，如有错误敬请指出。
网友评论