CN109493916A - 一种基于稀疏性因子分析的基因-基因交互作用识别方法 - Google Patents

Abstract

本发明公开一种基于稀疏性因子分析的基因‑基因交互作用识别方法(Sparse Factor Analysis for Epistasis,EPISFA)，其中包含以下步骤。1)输入基因型原始数据，并根据基因之间的连锁不平衡系数进行筛选；2)数据随机分块3)将数据按患病状态分为患病和非患病人群，分别计算两组的相关系数矩阵，利用Fisher变换扣除两人群的相关系数矩阵基因位点相关性；4)使用稀疏性因子分析方法学习模型权重；5)交叉验证，选择模型参数并识别相应的基因‑基因交互作用。试验表明此方法的统计功效和计算效率均较高，具有良好应用前景。

Description

一种基于稀疏性因子分析的基因-基因交互作用识别方法

1.技术领域

本发明涉及遗传学领域，具体涉及一种基于稀疏性因子分析的基因-基因交互作用识别方法。

2.背景技术

研究复杂疾病的遗传易感性一直是遗传学领域的重要问题。尽管近年来的全基因组关联研究发 现了大量与疾病存在关联的多态性位点，但仅纳入多态性位点的一维信息远未能解释复杂疾病在人 群中的遗传度。基因-基因交互作用是导致此类遗传度缺失的主要原因之一。

全基因组关联研究时代的遗传学研究往往一次检测大量的多态性位点，因而使用传统假设检验 方法难以克服“维度膨胀”带来的问题。为此，人们已经提出了很多基于机器学习的算法来识别基 因-基因交互作用。用于识别基因-基因交互作用的机器学习方法一般由降维算法和机器学习模型两 个部分组成。目前，常用降维算法有多因子降维法、组合优化法和稀疏学习方法，而常用的机器学 习模型则包括了Logic回归模型、随机森林模型、支持向量机、神经网络模型和贝叶斯网络模型等。

这些机器学习方法在一定程度上克服了传统方法在高维遗传学数据中的弱点，但在统计功效和 计算效率上依然存在不足，还难以满足在复杂疾病中进行基因-基因交互作用研究的需求。而近年 来，非监督式机器学习方法的发展逐渐受到重视，并在高维数据分析领域取得了一些成功，是基因 -基因交互作用方法学研究中一种较有前景的方向。

3.发明内容

本发明的目的是为了解决现有算法在处理高维情境中基因-基因交互作用识别能力不足的问 题。本发明基于稀疏性因子分析模型，提供了一种速度快，统计功效高，假阳性率低的方法。

为了达到上述目的，本发明采用的技术方案主要包括了以下步骤。

1)输入基因型原始数据，并根据基因之间的连锁不平衡系数进行筛选；

2)数据随机平均分为K块，其中K-1用于模型训练，剩余1块用于交叉验证；

3)将数据按患病状态分为患病和非患病人群，分别计算两组的相关系数矩阵，利用Fisher变换扣 除两组相关系数矩阵基因位点的相关性；

4)使用稀疏性因子分析方法学习模型权重，在训练集中先根据不同的模型参数拟合一组不同的稀疏 性因子模型，然后提取模型拟合权重运用于验证集并计算预测的相关系数矩阵。继之计算预测的 相关系数矩阵与数据中实际计算的相关系数矩阵之间的Kullback-Leibler距离作为预测值与观 察值之间的误差度量，并根据该距离选择相应的参数。。

5)重复步骤4)进行交叉验证，并识别数据存在中的基因-基因交互作用。

4.附图说明

附图1为本发明实施例的流程图

附图2为本发明模型结构示意图

5.具体实施方式

假设K＝{1,2,…,k}为一组SNP位点，编码x_k＝{-1,0,1},k∈K；y＝{0,1}为一个二分类质量性状， 定义M＝{1,2,…,m}，m<k为一组隐变量，编码为z_m,m∈M；n×k矩阵X是标准化后的基因型编码矩 阵，n×m矩阵Z是隐变量矩阵。定义一个维度为k×m的线性变换W，满足Z＝XW和X′＝ZW^T，并定义残 差矩阵为Ψ＝X-X′。

模型结构如图2所示，基因型编码矩阵X通过线性映射W，投影到隐变量矩阵Z上，再通过线性变换ZW^T还原到X′，并使得误差项Ψ最小。其中，根据稀疏性假设，Z的维度m＜＜k，并且W的绝大多数元素为0。

假设数据X和X′的误差函数为l，则模型可以表述为：

其中，ρ和γ均为调节参数，γ趋近于+∞时，模型趋近于LASSO，而γ趋近于1时，模型趋近于和MDR法类 似的阈值模型。通过尝试调节ρ和γ的不同取值，可以获得不同稀疏程度的W权重。其中，当ρ越大或γ越小时， 权重W则越稀疏。

公式(1)中的误差函数l(.)存在多种形式，为了能直接使用Hirose等的GEM方法，这里选择其在稀疏性因 子分析中使用的似然损失函数。即假设基因型编码，样本量为n，近似满足一个均值为0，方差为Σ的多元正态 分布，那么基于似然函数的损失函数l(.)可以表示为：

虽然SNP位点的基因型是一个离散型变量，并不满足多元正态分布，但公式(2)其实可以看作两个协方差矩 阵的Kullback-Leibler距离，因而也可以适用于不满足多元正态分布假设的情况。

设w_ij代表第i个基因型编码变量在第j个隐变量上的权重。如果一组之间存在两两相关， 那么存在一个隐变量z_j，使得权重w_.j同时不为0。如果权重向量w_.j不全为0，则非零权重的个数一定大于2。 这一性质确保了稀疏性因子分析方法总能发现2阶及以上的基因-基因交互作用。

当SNP位点在未患病个体中处于连锁平衡时，存在基因-基因交互作用的SNP位点在患病个体中就会存在相 关性。这一规律就是本方法的核心原理：当患病个体中，第i个隐变量上的权重同时不为0时， 则SNP位点{i₁,i₂,…,i_k}存在基因-基因交互作用。

假设一组家系数据中，患病个体的样本量为n^Aff，未患病个体的样本量为n^Un，患病个体中第i个和第j个位 点基因型的相关系数为未患病个体中第i个和第j个位点基因型的相关系数为利用Fisher’s Z变换 的原理，可以得到调整了连锁不平衡的相关系数和等效的样本量n^Adj。

其中，公式(3)中的tanh为双曲正切函数，将转换后的相关系数矩阵∑^Adj和等效样本量 n^Adj，重新带入EPISFA算法中，就可以得到扣除了位点之间连锁不平衡和人群分层效应的结果。 本研究中将这个补充算法命名为EPISFA-LD算法，具体步骤如下。

第1步，基因型编码和筛选

分别将在家系数据的患病和未患病个体，根据SNP位点的基因型进行编码，并根据基因型 编码的均数和标准差对编码矩阵分别进行标准化，得到患病和未患病个体的基因型编码矩阵X^Aff和X^Un，维度分别为n^Aff×p和n^Un×p，计算X^Un的相关系数矩阵∑^Un。设置隐变量的个数 m，1≤m＜p。

将SNP位点按照缺失率由高到低排序，计算未患病个体中各位点的Pearson相关系数，设 置相关系数阈值θ，按缺失率由高到低将相关系数超过θ的SNP进行修剪，直到SNP位点两两之间 相关性均小于θ。

第2步，随机分块。

将X^Aff按行平均分为k份，然后以为验证集，将剩余部分为训练集，分别计算验证集和训练集的相关 系数矩阵Σ^validate和Σ^Un，以及样本量n^validate＝n/k和n^train＝n×(k-1)/k。

第3步，相关性扣除。

将(Σ^train,n^train,Σ^Un,n^un)带入公式(3)和公式(4)，计算得到调整后的相关系数矩阵Σ^Adj和等效样本 量n^Adj。

第4步，稀疏性模型建立和权重估计

对调整后的相关系数矩阵Σ^Adj使用稀疏性因子分析模型，并令样本量为n^Adj。通过设置一系列不同的γ和ρ 参数组合，获得不同的权重矩阵W，然后在验证集中计算模型性能指标，选择其中性能最优的参数组合(γ^*,ρ^*)， 并获得在这个参数下的权重矩阵W。

第5步，交叉验证

提取不同列的权重向量w_.j，得到非零权重组合I₁₁,…,I_1m。重复第3-4步k次，得到k个不同的非零权重 I₁，并计算交叉验证一致性指标。交叉验证一致性指标为k次试验中，相同组合I_ij出现的次数。记录交叉验证 一次性指数最高的组合，并从{1,2,…,p}中将对应的SNP位点删除，在剩余位点中重新计算X^Aff，X^Un和Σ^Un。 重复第3-5步，直到算法中止条件满足，或权重矩阵W的所有元素均为0。

上述内容均采用了使用R语言进行实施。

本发明所采用的稀疏性因子分析技术，是一种将稀疏方法和非监督式机器学习方法进行结合的 方法。既能利用非监督式机器学习方法中统计功效高的优点，又能良好地控制假阳性率。此外，稀 疏方法还能显著提高算法的计算效率从而实现在较短时间内完成高维基因型数据的搜索。

Claims (4)

1.本发明根据基因分型所获得的基因型原始数据，计算和识别多个基因位点之间影响疾病风险的交互作用。相对于既有算法，可以在较短的计算时间内对全基因组关联研究的数据进行搜索，并在保持一类错误率不变的情况下，显著提高发现基因-基因交互作用的统计功效。本发明主要包括以下流程：

1)输入基因型原始数据，并根据基因位点之间的连锁不平衡系数进行筛选；

2)数据随机平均分为K块，其中K-1用于模型训练，剩余1块用于交叉验证；

3)将数据按患病状态分为患病和非患病人群，分别计算两组的相关系数矩阵，利用Fisher变换扣除两组相关系数矩阵基因位点的相关性；

4)使用稀疏性因子分析方法学习模型权重，在训练集中先根据不同的模型参数拟合一组不同的稀疏性因子模型，然后提取模型拟合权重运用于验证集并计算预测的相关系数矩阵。继之计算预测的相关系数矩阵与数据中实际计算的相关系数矩阵之间的Kullback-Leibler距离作为预测值与观察值之间的误差度量，并根据该距离选择相应的参数；

5)重复步骤4)进行交叉验证，并识别数据存在中的基因-基因交互作用。

2.根据权利1所述的基因-基因交互作用识别新算法，其特征在于：所述步骤1中，设置一个连锁不平衡系数阈值θ，并使用基因型原始数据计算位点之间的连锁不平衡系数矩阵，从第一个基因位点开始，依次排除超过连锁不平衡阈值的基因位点，直到连锁不平衡系数矩阵的所有非对角线元素均低于阈值。

3.根据权利1所述的基因-基因交互作用识别新算法，其特征在于：所述步骤4中，建立稀疏性因子分析模型，设置参数ρ和γ，并将非零因子载荷作为模型的权重。

4.根据权利1所述的基因-基因交互作用识别新算法，其特征在于：所述步骤5中，通过交叉验证的方式对参数ρ和γ进行选择，并在最优参数组合下训练算法，将此参数组合下所识别的非零因子载荷组合，作为存在交互作用的基因组合。