Python用Lasso改进线性混合模型Linear Mixed Model分析拟南芥和小鼠复杂性状遗传机制多标记表型预测可视化

全文链接：https://tecdat.cn/?p=38800

原文出处：拓端数据部落公众号

在生物医学领域，探究可遗传性状的遗传基础是关键挑战之一。对于受多基因位点多因素控制的性状，准确检测其关联存在诸多困难，且易受群体结构等混杂因素影响产生假阳性结果。本文帮助客户建立Lasso线性混合模型，它能实现多位点定位并校正混杂效应，具有无需调参、可有效控制群体结构、适用于全基因组数据集等优点，能同时发现可能的因果变体并进行基于基因型的多标记表型预测。通过在拟南芥和小鼠相关研究中的应用，验证了其有效性，为解析复杂性状的遗传机制提供了有力工具。

引言

人类、动植物中诸多数量性状虽具遗传性，但人们对其潜在遗传结构的全面认识仍不足。像全基因组关联研究和连锁图谱分析虽已揭示出部分控制性状变异的因果变体，可单基因变体往往只能解释表型变异的小部分，个体效应量小。样本间的关联以及群体结构会导致假关联模式出现，影响对复杂性状真实遗传结构的理解。在分析遗传数据时，传统单独评估单个位点显著性的方法存在局限，而多元回归虽能一定程度解决问题，但相关方法也存在不足。为此，人们一直在探索更好的模型与方法来应对这些挑战，以准确捕捉基因位点对表型的影响。

多元线性混合模型

（一）模型构建基础

我们的方法基于多元统计，通过个体遗传效应和随机混杂变量的总和来解释表型变异。简单来说，对于 m 个样本的表型 y（包含 y1 到 ym 等元素），可以表示为 n 个单核苷酸多态性（SNPs）的总和、混杂因素以及噪声的叠加，如下式：

1.   # y 表示样本表型，S 表示基因型数据（包含多个 SNPs 信息），beta 表示遗传效应系数，u 表示混杂影响，t 表示噪声
2.   y = sum(beta_j * s_j for j in range(n)) + u + t

这里的 t 代表观测噪声，u 是混杂影响，在遗传图谱分析中，u 通常不能直接观测，但很多情况下其高斯协方差 K 可从观测数据中估计出来，比如利用能体现样本间整体遗传相似性的关系矩阵等方式来估计，以此来考虑群体结构等造成的混杂情况。以往多是针对单个候选 SNPs 考虑这种混合模型，不过这样在复杂遗传结构下可能存在问题，比如一些遗传因素相互掩盖等情况。所以我们提出一种更有效的联合推断方法，也就是 Lasso线性混合模型 方法，它能同时考虑所有 SNPs 及其相互依赖关系，还能结合经典 Lasso 的优点，实现稀疏回归，将解释的表型方差分解为个体 SNP 效应和群体结构效应两部分。

线性混合模型 Lasso 细节

设 S 为 m 个个体的 n 个 SNPs 构成的 m×n 矩阵，sj 就是代表 SNP j 的 m×1 向量。我们把 m 个个体的表型 y 建模为 SNPs 的遗传效应和混杂影响 u 的总和（见前面提到的公式），遗传效应当作固定效应，混杂影响当作随机效应，并且通过拉普拉斯收缩先验让大部分 SNPs 的效应大小为 0。随机变量 u 假设服从高斯分布，其协方差为 K，观测噪声 t 也假设服从高斯分布，通过对随机变量 u 边缘化处理，可以得到权重向量 b 的条件后验分布，这里涉及到一些参数如拉普拉斯先验的稀疏性超参数、残差噪声方差以及随机效应分量的方差等。

参数推断

联合学习超参数和权重是个较难的非凸优化问题。我们采用先在零模型下拟合部分参数（排除个体 SNP 效应），再转化为标准 Lasso 回归问题的方法。

• 零模型拟合：为得到实用且可扩展的算法，先在零模型下通过最大似然法优化随机效应方差和噪声方差相关的参数，通过一些计算技巧来高效优化它们之间的比率，比如先对协方差 K 进行特征分解来转换数据，使正态分布的协方差矩阵各向同性，然后对边际似然进行一维数值优化等。

1.   """
2.   训练随机效应模型
3.   :param y: 表型数据
4.   :param K: 亲缘关系矩阵
5.   :param numintervals: 用于 delta 线搜索的区间数量
6.   :param ldeltamin: 最小 delta 值（对数空间）
7.   :param ldeltamax: 最大 delta 值（对数空间）
8.   :param debug: 是否开启调试模式
9.   :return: 相关参数及监控信息
10.   """
11.    
12.    
13.   n_s = y.shape[0]
14.    
15.   # 旋转数据
16.   S,U = LA.eigh(K)
17.   Uy = SP.dot(U.T,y)
18.    
19.   # 网格搜索
20.   nllgrid = SP.ones(numintervals + 1) * SP.inf
21.   ldeltagrid = SP.arange(numintervals + 1) / (numintervals * 1.0) * (ldeltamax - ldeltamin) + ldeltamin
22.   nllmin = SP.inf
23.   for i in SP.arange(numintervals + 1):
24.   nllgrid[i] = nLLeval(ldeltagrid[i], Uy, S)

• 转化为标准 Lasso 问题：固定相关比率后，再次利用 K 的特征分解转换数据，使得确定权重的任务等价于 Lasso 回归模型，通过最小化相应的目标函数来求解，合适的稀疏性超参数可以通过交叉验证等方法来确定，不过在多元情况下，完全的交叉验证优化成本较高，所以实验中多在零模型上优化相关参数。

表型预测

给定基于一组基因型和表型训练好的 Lasso线性混合模型，就能预测测试个体未观测到的表型。预测分布可通过对所有个体的联合分布基于训练个体进行条件化推导得出，其均值预测是 Lasso 部分和随机效应部分贡献的总和，类似最佳线性无偏预测，这里涉及到不同个体间协方差矩阵等的相关计算。

1.   预测表型
2.    
3.   输入:
4.   y_t: 训练表型数据，形状为 n_train x 1
5.   X_t: 训练 SNP 矩阵，形状为 n_train x n_f
6.   X_v: 验证 SNP 矩阵，形状为 n_val x n_f
7.   K_tt: 训练集亲缘关系矩阵，形状为 n_train x n_train
8.   K_vt: 验证集与训练集的亲缘关系矩阵，形状为 n_val x n_train
9.   ldelta: 核参数
10.   w: Lasso 权重，形状为 n_f x 1
11.    
12.   输出:
13.   y_v: 预测的表型数据，形状为 n_val x 1
14.   """
15.    
16.    
17.    
18.   [n_train,n_f] = X_t.shape
19.   n_test = X_v.shape[0]
20.    
21.   if y_t.ndim==1:
22.   y_t = SP.reshape(y_t,(n_train,1))
23.   if w.ndim==1:
24.   w = SP.reshape(w,(n_f,1))
25.    
26.   delta = SP.exp(ldelta)
27.   idx = w.nonzero()[0]
28.    
29.   if idx.shape[0]==0:
30.   return SP.dot(K_vt,LA.solve(K_tt + delta*SP.eye(n_train),y_t))
31.    

考虑群体结构的随机效应协方差选择

随机效应协方差 K 的选择有多种方式，比如基于血缘关系矩阵、状态相同矩阵、实现关系矩阵（RRM）等，后续实验中我们采用的是实现关系矩阵。从贝叶斯角度看，使用 RRM 作为协方差矩阵相当于在特定模型下对权重设置独立高斯先验并进行积分，这样的选择可以看作是对因群体结构混杂或微小累加效应的遗传效应进行建模，而效应足够大的单个 SNPs 则直接纳入 Lasso 模型。

方法与材料

（一）拟南芥研究

从相关文献获取了多达 199 份拟南芥的基因型和表型数据，每份基因型包含大量 SNPs，我们着重研究与植物开花时间相关的表型，排除了部分样本量小的表型数据，最终考虑 20 种开花表型，拟南芥个体间相关性范围广，群体结构复杂。

（二）小鼠近交群体研究

获取了 1940 只来自多亲本近交群体小鼠的基因型和表型数据，每只小鼠基因型包含一定数量的 SNPs，这些小鼠来自八个近交品系杂交产生的异质群体，其表型涵盖从生化到行为等多种不同测量指标，我们聚焦于有数值或二元值的 273 种表型。

（三）半经验数据

为评估变量选择替代方法的准确性，基于扩展的拟南芥数据集构建了半经验示例，包含真实表型数据以体现受群体结构影响的背景信号，还添加了不同效应大小和遗传模型复杂度的模拟关联，具体模拟过程和不同方法恢复关联的评估详见补充材料。

（四）数据预处理

对 SNP 数据进行标准化，使得效应大小的先验依赖于次要等位基因频率；对表型数据进行 Box–Cox 变换后再标准化。

（五）模型选择

对于 Lasso 方法模型复杂度的变化，可通过选择活跃 SNPs 数量或者改变相关超参数来实现，我们选择改变活跃 SNPs 数量。在进行表型预测时，采用 10 折交叉验证，将数据随机分成 10 份，轮流选一份作测试集，其余作训练集，选能使测试集解释方差最大的模型；在进行变量选择评估个体特征显著性时，采用稳定性选择，固定活跃 SNPs 数量，多次随机抽取部分数据，选在多数重复中出现的 SNPs，并根据其选择频率推导显著性估计等。

稳定性选择函数 ：

1.    
2.   运行稳定性选择
3.    
4.   输入:
5.   X: SNP 矩阵，形状为 n_s x n_f
6.   y: 表型数据，形状为 n_s x 1
7.   K: 亲缘关系矩阵，形状为 n_s x n_s
8.   mu: L1 惩罚项
9.    
10.   n_reps: 重复次数
11.   f_subset: 用于创建一个自助抽样数据集的数据集占比
12.    
13.   输出:
14.   所有 SNPs 的选择频率，形状为 n_f x 1
15.   """
16.   time_start = time.time()
17.   [n_s,n_f] = X.shape
18.   n_subsample = SP.ceil(f_subset * n_s)
19.   freq = SP.zeros(n_f)
20.    
21.   for i in range(n_reps):
22.   print('Iteration %d' % i)
23.   idx = SP.random.permutation(n_s)[:n_subsample]
24.   res = train(X[idx],K[idx][:,idx],y[idx],mu,**kwargs)
25.   snp_idx = (res['weights']!=0).flatten()
26.   freq[snp_idx] += 1.

模型训练

1.   # 训练零模型
2.   S,U,ldelta0,monitor_nm = tnullmodel(y,K,numintervals,ldeltamin,ldeltamax,debug=debug)
3.    
4.   # 在残差上训练 Lasso
5.   delta0 = SP.exp(ldelta0)
6.   Sdi = 1./(S+delta0)
7.   Sdi_sqrt = SP.sqrt(Sdi)
8.   SUX = SP.dot(U.T,X)

结果

（一）半经验设置下的验证

在模拟的半经验数据集中评估 Lasso线性混合模型 恢复真实基因型 - 表型关联的能力，与标准 Lasso、单变量线性混合模型（LMM）以及标准单变量线性模型（LM）等现有技术对比。Lasso线性混合模型在恢复有真实模拟关联的 SNPs 准确性上表现更好，尤其在考虑群体结构校正方面优势明显，在不同模拟设置下，如改变群体结构强度、遗传模型复杂度、信号与噪声比等情况下，Lasso线性混合模型 大多都优于其他方法，说明在存在较多弱关联等情况下，考虑能体现弱效应的遗传协方差是有优势的。

图 1 说明：展示了在模拟半经验全基因组关联研究（GWAS）数据集上，使用不同方法恢复模拟因果 SNPs 的精确率 - 召回率曲线（a）以及接收者操作特征曲线（ROC）（b），可以看出不同方法的效果对比。

图 2 说明：呈现了在不同模拟设置下，不同方法在半经验 GWAS 数据集上寻找真实模拟关联的精确率 - 召回率曲线下面积，体现各方法在不同情况的性能表现。

（二）Lasso线性混合模型在模型系统中对复杂性状遗传结构的解析

在拟南芥和小鼠中验证Lasso线性混合模型对基因型 - 表型映射的建模能力，聚焦拟南芥的 20 种开花时间表型和小鼠与人类健康相关的 273 种表型。

• 更准确的表型预测与更稀疏的遗传模型：对比不同方法解释基因型对表型联合效应的能力，通过评估可解释的表型变异分数来衡量预测能力，发现Lasso线性混合模型在多数拟南芥和小鼠表型中比 Lasso 更准确地预测表型，能解释更多由遗传因素导致的表型变异，且其选择的 SNPs 数量更少，意味着构建的遗传模型更稀疏，不过仍体现出多基因位点的复杂遗传控制。
  
  图 3 说明：展示了 Lasso 和 Lasso线性混合模型应用于模型系统中数量性状时，拟合遗传模型的预测能力和稀疏性，包括拟南芥开花表型和小鼠生化、生理表型相关情况。
• 区分个体遗传效应和群体结构效应：以拟南芥开花时间表型为例，研究 Lasso线性混合模型区分个体遗传效应和群体结构效应的能力，随着模型中包含的 SNPs 数量增加，Lasso线性混合模型 解释方差从主要受群体结构影响逐渐向个体 SNPs 效应转移，且对比 Lasso 有更好的解释方差表现。
  
  图 4 说明：呈现了将方差分解为个体 SNP 效应和由群体结构驱动的全局遗传背景的情况，展示了拟南芥开花表型在独立测试集上解释方差随活跃 SNPs 数量的变化。
• 与已知候选基因的关联富集：考虑不同方法检索到的关联在拟南芥中与已知开花相关候选基因的接近程度，发现 Lasso线性混合模型 找到与候选基因相连的 SNPs 数量在多数表型中多于 Lasso，且在不同选择阈值下大多表现更优，此外，Lasso线性混合模型 选择的与候选基因相连的 SNPs 中相邻成对等情况的比例也有一定特点，一定程度体现其检测遗传异质性的能力。

讨论

本文建立的Lasso线性混合模型结合了混合模型校正混杂效应和多标记模型考虑基因标记联合效应的优点，能更好地恢复真实遗传效应，在复杂遗传结构、个体标记弱效应或强混杂效应等具有挑战性的情况下也有较好表现。它不仅能提升从基因型预测表型的能力，还能区分不同来源的表型变异，检索到的遗传关联在已知候选基因方面有富集情况。虽然考虑群体结构和对遗传数据进行多元建模的概念本身不是全新的，但现有相关方法存在各种局限，比如有的不能考虑随机效应控制混杂，有的缺乏可扩展算法等，而 Lasso线性混合模型能明确将混杂解释为随机效应，有助于解决个体遗传效应和因群体结构导致的表型变异之间的模糊性问题，是解析基因型 - 表型关系的有用工具，不过对于 Lasso 方法在评估统计显著性方面仍存在挑战，有待未来进一步研究。

参考文献：

• Atwell, S., et al. (2010). Genome-wide association study of 107 phenotypes in Arabidopsis thaliana inbred lines. Nature, 465, 627-631.
• Bradley, J. K., et al. (2011). Parallel coordinate descent for l1-regularized loss minimization. ICML, 321-328.
• Bühlmann, P. (2012). Statistical significance in high-dimensional linear models. arXiv:1202.1377v2.
• Craddock, N., et al. (2010). Genome-wide association study of cnvs in 16,000 cases of eight common diseases and 3,000 shared controls. Nature, 464, 713-720.
• Flint, J., & Mackay, T. F. (2009). Genetic architecture of quantitative traits in mice, flies, and humans. Genome Res., 19, 723-733.
• Foster, S., et al. (2007). Incorporating lasso effects into a mixed model for quantitative trait loci detection. J. Agric. Biol. Environ. Stat., 12, 300-314.
• Fusi, N., et al. (2012). Joint modelling of confounding factors and prominent genetic regulators provides increased accuracy in genetical genomics studies. PLoS Comput. Biol., 8, e1002330.
• Goddard, M. E., et al. (2009). Estimating effects and making predictions from genome-wide marker data. Stat. Sci., 24, 517-529.
• Hastie, T., et al. (2003). The Elements of Statistical Learning. Corrected edition. Springer New York Inc., New York, NY, USA.
• Hayes, B. J., et al. (2009). Increased accuracy of artificial selection by using the realized relationship matrix. Genet. Res. (Camb.), 91, 47-60.
• Hoggart, C. J., et al. (2008). Simultaneous analysis of all SNPs in genome-wide and re-sequencing association studies. PLoS Genet., 4, e1000130.
• Horton, M. W., et al. (2012). Genome-wide patterns of genetic variation in worldwide Arabidopsis thaliana accessions from the RegMap panel. Nat. Genet., 44, 212-216.
• Kang, H. M., et al. (2010). Variance component model to account for sample structure in genome-wide association studies. Nat. Genet., 42, 348-354.
• Kang, H., et al. (2008). Efficient control of population structure in model organism association mapping. Genetics, 178, 1709-1723.
• Kim, S., & Xing, E. P. (2009). Statistical estimation of correlated genome associations to a quantitative trait network. PLoS Genet., 5, e1000587.
• Lee, S., & Xing, E. P. (2012). Leveraging input and output structures for joint mapping of epistatic and marginal eQTLs. Bioinformatics, 28, i137-i146.
• Lippert, C., et al. (2011). FaST linear mixed models for genome-wide association studies. Nat. Methods, 8, 833-835.
• Listgarten, J., et al. (2010). Correction for hidden confounders in the genetic analysis of gene expression. Proc. Natl Acad. Sci. USA, 107, 16465-16470.
• Li, J., et al. (2011). The bayesian lasso for genome-wide association studies. Bioinformatics, 27, 516-523.
• Mackay, T. F. C., et al. (2009). The genetics of quantitative traits: challenges and prospects. Nat. Rev. Genet., 10, 565-577.
• McCarthy, M., et al. (2008). Genome-wide association studies for complex traits: consensus, uncertainty and challenges. Nat. Rev. Genet., 9, 356-369.
• Meinshausen, N., & Bühlmann, P. (2010). Stability selection. J. R. Stat. Soc. Series B Stat. Methodol., 72, 417-473.
• Meinshausen, N., et al. (2009). P-values for high-dimensional regression. J. Am. Stat. Assoc., 104, 1671-1681.
• Newman, D., et al. (2001). The importance of genealogy in determining genetic associations with complex traits. Am. J. Hum. Genet., 69, 1146-1148.
• Ober, U., et al. (2012). Using whole-genome sequence data to predict quantitative trait phenotypes in Drosophila melanogaster. PLoS Genet., 8, e1002685.
• Platt, A., et al. (2010a). Conditions under which genome-wide association studies will be positively misleading. Genetics, 186, 1054-1052.
• Platt, A., et al. (2010b). The scale of population structure in Arabidopsis thaliana. PLoS Genet., 6, e1000843.