一种基于大样本的基因-环境互作关联分析方法

摘要

本发明公开了一种基于大样本的基因‑环境互作关联分析方法，通过采集目标个体的基因序列、环境信息、表型信息，构建基因‑环境检测线性混合模型，并基于预处理共轭梯度法、矩估计方法对基因‑环境检测线性混合模型进行求解，构建服从于卡方分布的基因‑环境互作得分测试统计量，然后应用基因‑环境互作得分测试统计量，完成对基因序列与环境信息是否存在互作效应的判断，该方法可以快速识别多个环境信息对基因‑环境互作位点的基因的影响，适用于大样本的基因‑环境互作关联快速检测。

说明书

技术领域

本发明涉及遗传算法领域，具体涉及一种基于大样本的基因-环境互作关联分析方法。

背景技术

人类很多基因受到环境因子的调控，环境因子包括饮食、身体活动和其他生活方式协变量，其中基因-环境相互作用(GEI)值得研究。目前，拥有数百万参与者的基因组数据集正在以前所未有的规模被收集，这些数据集包括All of Us Research Program(约100万参与者)、UK BioBank(约50万参与者)和Million Veteran Program(多于100万参与者)。虽然如此大的基因组数据为基因-环境互作(GEI)研究提供了极好的机会，但GEI分析工具的研发的滞后一直保持不变，完全无法适应基因-环境相互作用调查的需要。现有的多数方法没有针对数十万或数百万个人的分析进行优化，计算资源要求过高，并且缺乏多线程能力来减少运行时间。QUICKTEST方法、PLINK2方法、CGEN方法和GxEScan方法虽然可以在大样本中进行有效分析，但不能进行稳健的推理和多重交互项计算；SPAGE方法目前仅支持二元性状的分析，且缺少多线程计算能力。StructLMM方法是一种计算高效的方法，可以识别基因与多个环境相互作用，并且能够应用于数十万人的大样本分析，它不能处理相关个体，忽略相关个体会造成统计功效的损失。因此，有必要开发一种适用于相关大样本的基因-环境互作关联分析方法。

发明内容

本发明目的：在于提供一种基于大样本的基因-环境互作关联分析方法，实现快速识别多个环境信息对基因-环境互作位点的基因的影响，适用于大样本的基因-环境互作关联快速检测。

为实现以上功能，本发明设计一种基于大样本的基因-环境互作关联分析方法，执行如下步骤S1-步骤S5，获得基因-环境互作得分测试统计量，然后应用基因-环境互作得分测试统计量，完成对基因序列与环境信息是否存在互作效应的判断；

S1:采集目标个体的基因序列、环境信息、表型信息，并基于目标个体的基因序列、环境信息、表型信息构建目标个体的基因-环境样本；

S2:基于目标个体的基因-环境样本中的基因序列、环境信息，以基因序列中预设基因-环境互作位点的基因为焦点变体，针对该焦点变体，构建焦点变体所对应的基因型向量，基于目标个体的环境信息，构建环境信息所对应的环境矩阵；

S3:基于步骤S2所获得的基因型向量、环境矩阵，通过矩阵乘法，构建基因-环境交互设计矩阵，并根据对效应的假定，构建基因型向量所对应的焦点变体固定效应向量、环境矩阵所对应的环境随机效应向量、基因-环境交互设计矩阵所对应的基因-环境交互随机效应向量；

基于目标个体的表型信息，构建数量性状表型向量；

基于基因型向量、焦点变体固定效应向量、环境矩阵、环境随机效应向量、基因-环境交互设计矩阵、基因-环境交互随机效应向量，以及数量性状表型向量，构建基因-环境检测线性混合模型；

S4:基于步骤S3所构建的基因-环境检测线性混合模型，采用预处理共轭梯度法、矩估计方法对基因-环境检测线性混合模型进行求解，构建服从于卡方分布的基因-环境互作得分测试统计量，所述基因-环境互作得分测试统计量用于判断所述环境信息是否对目标个体的预设所有基因-环境互作位点的基因产生影响；

S5:基于步骤S4所构建的基因-环境互作得分测试统计量，通过对基因序列中各基因-环境互作位点的随机抽样，计算基因-环境互作得分测试统计量所对应的卡方分布的系数a值，进而获得所有基因-环境互作位点的卡方分布的P值，并预设基因-环境互作阈值，当卡方分布P值小于基因-环境互作阈值，则判定环境信息对该基因-环境互作位点的基因产生影响，即存在基因-环境互作效应，否则判定环境信息对该基因-环境互作位点的基因不产生影响，即不存在基因-环境互作效应。

作为本发明的一种优选技术方案：步骤S3中构建基因-环境交互设计矩阵如下式：

S＝G⊙E＝diag(G)E

式中，S为基因-环境交互设计矩阵，其形式为N×Q维矩阵，N为样本量大小，Q为环境信息个数，⊙表示哈达玛积矩阵乘法计算，G为焦点变体所对应的基因型向量，E为环境信息所对应的环境矩阵，其形式为N×Q维的矩阵；

步骤S3中所构建的基因-环境检测线性混合模型如下式：

Y＝Xβ

式中，Y为数量性状表型向量，其形式为N×1维向量，N为样本量大小，X为包括截距的固定效应设计矩阵，其形式为N×P维的协变量矩阵，β

作为本发明的一种优选技术方案：步骤S4中采用预处理共轭梯度法、矩估计方法对基因-环境检测线性混合模型进行求解，构建服从于卡方分布的基因-环境互作得分测试统计量的具体步骤如下：

S41:分别针对环境矩阵所对应的环境随机效应向量β

式中，diag(G)E＝G⊙E＝S，∑

S42:基于步骤S41所获得的多元正态分布数量性状表型向量Y，构建服从于卡方分布的基因-环境互作得分测试统计量T如下式：

式中：

K＝diag(G)∑

＝[diag(G)E][diag(G)E]

基因-环境互作得分测试统计量T表示为如下形式：

式中：

S＝G⊙E＝diag(G)E

式中：

其中H

式中，式中

S43:针对各方差组分

式中，YY

S44:采用矩估计方法对方差组分进行求解如下式：

式中，V表示为如下形式：

V＝I

式中，D为一个2×2维矩阵，具体表示为如下形式：

其中，Tr表示矩阵的迹；

式中，b为一个2×1维列向量，具体表示为如下形式：

式中，c为一个2×1维列向量，具体表示为如下形式：

S45:针对步骤S44所获得的矩阵D，通过近似迭代法，对其进行预设次数的迭代，并采用Hutchinson估计器获得矩阵D中各元素的迹。

作为本发明的一种优选技术方案：步骤S5中通过对基因序列中各基因-环境互作位点的随机抽样，计算基因-环境互作得分测试统计量T所对应的卡方分布的系数a值，进而获得所有基因-环境互作位点的卡方分布的P值,并基于预设基因-环境互作阈值，对环境信息是否对该基因-环境互作位点的基因产生影响进行判定的具体步骤如下：

S51：基因-环境互作得分测试统计量T服从卡方分布，具体表示为如下形式：

式中，a为卡方分布的系数，

S52：对目标个体的基因序列中各基因-环境互作位点进行n次随机抽样，分别针对每次随机抽取的基因-环境互作位点的基因，根据S

S53：将步骤S52获得卡方分布的系数a的估计值

基于卡方分布表，获得基因-环境互作得分测试统计量T所对应的卡方分布P值。

本发明还设计一种基于大样本的基因-环境互作关联分析系统，包括：

一个或多个处理器；

存储器，存储可被操作的指令，所述指令在通过所述一个或多个处理器执行时使得所述一个或多个处理器执行操作，通过以下步骤获得基因-环境互作得分测试统计量，然后应用基因-环境互作得分测试统计量，完成对基因序列与环境信息是否存在互作效应的判断：

S1:采集目标个体的基因序列、环境信息、表型信息，并基于目标个体的基因序列、环境信息、表型信息构建目标个体的基因-环境样本；

S2:基于目标个体的基因-环境样本中的基因序列、环境信息，以基因序列中预设基因-环境互作位点的基因为焦点变体，针对该焦点变体，构建焦点变体所对应的基因型向量，基于目标个体的环境信息，构建环境信息所对应的环境矩阵；

S3:基于步骤S2所获得的基因型向量、环境矩阵，通过矩阵乘法，构建基因-环境交互设计矩阵，并根据对效应的假定，构建基因型向量所对应的焦点变体固定效应向量、环境矩阵所对应的环境随机效应向量、基因-环境交互设计矩阵所对应的基因-环境交互随机效应向量；

基于目标个体的表型信息，构建数量性状表型向量；

基于基因型向量、焦点变体固定效应向量、环境矩阵、环境随机效应向量、基因-环境交互设计矩阵、基因-环境交互随机效应向量，以及数量性状表型向量，构建基因-环境检测线性混合模型；

S4:基于步骤S3所构建的基因-环境检测线性混合模型，采用预处理共轭梯度法、矩估计方法对基因-环境检测线性混合模型进行求解，构建服从于卡方分布的基因-环境互作得分测试统计量，所述基因-环境互作得分测试统计量用于判断所述环境信息是否对目标个体的预设所有基因-环境互作位点的基因产生影响；

S5:基于步骤S4所构建的基因-环境互作得分测试统计量，通过对基因序列中各基因-环境互作位点的随机抽样，计算基因-环境互作得分测试统计量所对应的卡方分布的系数a值，进而获得所有基因-环境互作位点的卡方分布的P值，并预设基因-环境互作阈值，当卡方分布P值小于基因-环境互作阈值，则判定环境信息对该基因-环境互作位点的基因产生影响，即存在基因-环境互作效应，否则判定环境信息对该基因-环境互作位点的基因不产生影响，即不存在基因-环境互作效应。

本发明还设计一种存储软件的计算机可读取介质，所述可读取介质包括能通过一个或多个计算机执行的指令，所述指令在被所述一个或多个计算机执行时，执行如所述权利要求1-4中任意一项所述一种基于大样本的基因-环境互作关联分析方法的操作。

有益效果：相对于现有技术，本发明的优点包括：

1.本发明设计了一种基于大样本的基因-环境互作关联分析方法，该方法可同时处理一个基因与几个甚至上百个环境的互作效应，并考虑个体之间的相关性，弥补了StructLMM方法仅能分析无关个体的不足，能够充分利用所有样本的信息，可以高效的分析几十万的大样本遗传数据集；同时，本发明的方法可以进行稳健的推理计算，支持多线程并行运算，既可以分析数量性状，又可以分析二元性状，弥补了QUICKTEST方法和SPAGE等方法存在的不足。利用该方法可发现更具生物学意义的基因-环境互作位点，为清晰环境信息对基因的影响的机理提供了方法支持。

2.通过Monte Carlo模拟实验分析，本发明设计的方法由于考虑了相关个体，充分利用所有样本信息，其基因-环境互作的检测功效更高；由于本方法采用了多线程并行计算技术，其运算速度大幅提高，CPU的核数越多加速越明显，与传统方法比较，加速比可达2～50倍；由于对巨型矩阵进行了自适应分块，本方法可在一般内存的电脑(32GB及以上)上运行，很大程度提高了本方法的适用范围；由于本方法利用稳健的标准误差方法来处理违反线性模型假设的情况，针对不同数据，相比现有技术，本方法更稳健。

附图说明

图1是根据本发明实施例提供的基于大样本的基因-环境互作关联分析方法的流程图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

参照图1，本发明实施例提供的一种基于大样本的基因-环境互作关联分析方法，执行如下步骤S1-步骤S5，获得基因-环境互作得分测试统计量，然后应用基因-环境互作得分测试统计量，完成对基因序列与环境信息是否存在互作效应的判断；

S1:采集目标个体的基因序列、环境信息、表型信息，并基于目标个体的基因序列、环境信息、表型信息构建目标个体的基因-环境样本；

S2:基于目标个体的基因-环境样本中的基因序列、环境信息，以基因序列中预设基因-环境互作位点的基因为焦点变体，针对该焦点变体，构建焦点变体所对应的基因型向量，基于目标个体的环境信息，构建环境信息所对应的环境矩阵；

S3:基于步骤S2所获得的基因型向量、环境矩阵，通过矩阵乘法，构建基因-环境交互设计矩阵，并根据对效应的假定，构建基因型向量所对应的焦点变体固定效应向量、环境矩阵所对应的环境随机效应向量、基因-环境交互设计矩阵所对应的基因-环境交互随机效应向量；

基于目标个体的表型信息，构建数量性状表型向量；

基于基因型向量、焦点变体固定效应向量、环境矩阵、环境随机效应向量、基因-环境交互设计矩阵、基因-环境交互随机效应向量，以及数量性状表型向量，构建基因-环境检测线性混合模型；

步骤S3中构建基因-环境交互设计矩阵如下式：

S＝G⊙E＝diag(G)E

式中，S为基因-环境交互设计矩阵，其形式为N×Q维矩阵，N为样本量大小，Q为环境信息个数，其中N>0，Q>0，⊙表示哈达玛积矩阵乘法计算，G为焦点变体所对应的基因型向量，E为环境信息所对应的环境矩阵，其形式为N×Q维的矩阵；

步骤S3中所构建的基因-环境检测线性混合模型如下式：

Y＝Xβ

式中，Y为数量性状表型向量，例如：目标个体身高、BMI(身体质量指数)和一秒呼气量等表型，其形式为N×1维向量，N为样本量大小，X为包括截距的固定效应设计矩阵，其形式为N×P维的协变量矩阵，其中P>0，β

S4:基于步骤S3所构建的基因-环境检测线性混合模型，采用预处理共轭梯度法、矩估计方法对基因-环境检测线性混合模型进行求解，构建服从于卡方分布的基因-环境互作得分测试统计量，所述基因-环境互作得分测试统计量用于判断所述环境信息是否对目标个体的预设所有基因-环境互作位点的基因产生影响；

步骤S4中采用预处理共轭梯度法、矩估计方法对基因-环境检测线性混合模型进行求解，构建服从于卡方分布的基因-环境互作得分测试统计量的具体步骤如下：

S41:分别针对环境矩阵所对应的环境随机效应向量β

式中，diag(G)E＝G⊙E＝S，∑

S42:基于步骤S41所获得的多元正态分布数量性状表型向量Y，构建服从于卡方分布的基因-环境互作得分测试统计量T如下式：

式中：

K＝diag(G)∑

＝[diag(G)E][diag(G)E]

基因-环境互作得分测试统计量T表示为如下形式：

式中：

S＝G⊙E＝diag(G)E

式中：

其中H

式中，式中

对于大样本基因-互作位点的矩阵而言，基因-环境互作得分测试统计量无法直接构建，其实现过程需要将巨型矩阵进行分块，具体方法如下：

S43:针对各方差组分

式中，YY

S44:采用矩估计方法对方差组分进行求解如下式：

式中，V表示为如下形式：

V＝I

式中，D为一个2×2维矩阵，具体表示为如下形式：

其中，Tr表示矩阵的迹；

式中，b为一个2×1维列向量，具体表示为如下形式：

式中，c为一个2×1维列向量，具体表示为如下形式：

S45:针对步骤S44所获得的矩阵D，由于D中含有求解高维矩阵相乘的迹，直接对其求解时间复杂度非常高，因此通过近似迭代法，对其进行预设次数的迭代，并采用Hutchinson估计器获得矩阵D中各元素的迹。

S5:基于步骤S4所构建的基因-环境互作得分测试统计量，通过对基因序列中各基因-环境互作位点的随机抽样，计算基因-环境互作得分测试统计量所对应的卡方分布的系数a值，进而获得所有基因-环境互作位点的卡方分布的P值，并预设基因-环境互作阈值，当卡方分布P值小于基因-环境互作阈值，则判定环境信息对该基因-环境互作位点的基因产生影响，即存在基因-环境互作效应，否则判定环境信息对该基因-环境互作位点的基因不产生影响，即不存在基因-环境互作效应。

步骤S5中通过对基因序列中各基因-环境互作位点的随机抽样，计算基因-环境互作得分测试统计量T所对应的卡方分布的系数a值，进而获得所有基因-环境互作位点的卡方分布的P值,并基于预设基因-环境互作阈值，对环境信息是否对该基因-环境互作位点的基因产生影响进行判定的具体步骤如下：

S51：基因-环境互作得分测试统计量T服从卡方分布，具体表示为如下形式：

式中，a为卡方分布的系数，

若对于大样本的基因-环境互作位点求得卡方分布的系数a的精确值，其计算负担无疑是巨大的，因此在一个实施例中，通过对基因序列中各基因-环境互作位点的随机抽样，计算基因-环境互作得分测试统计量所对应的卡方分布的系数a值的估计值；

S52：对目标个体的基因序列中各基因-环境互作位点进行n次随机抽样，分别针对每次随机抽取的基因-环境互作位点的基因，根据S

在一个实施例中，此处随机抽取的n个基因-环境互作位点的基因，n＝30；

在一个实施例中，此处随机抽取的n个基因-环境互作位点的基因，n＝50；

S53：将步骤S52获得卡方分布的系数a的估计值

基于卡方分布表，获得基因-环境互作得分测试统计量T所对应的卡方分布P值。

本发明实施例提供的一种基于大样本的基因-环境互作关联分析系统，包括：

一个或多个处理器；

存储器，存储可被操作的指令，所述指令在通过所述一个或多个处理器执行时使得所述一个或多个处理器执行操作，通过以下步骤获得基因-环境互作得分测试统计量，然后应用基因-环境互作得分测试统计量，完成对基因序列与环境信息是否存在互作效应的判断：

S1:采集目标个体的基因序列、环境信息、表型信息，并基于目标个体的基因序列、环境信息、表型信息构建目标个体的基因-环境样本；

S2:基于目标个体的基因-环境样本中的基因序列、环境信息，以基因序列中预设基因-环境互作位点的基因为焦点变体，针对该焦点变体，构建焦点变体所对应的基因型向量，基于目标个体的环境信息，构建环境信息所对应的环境矩阵；

S3:基于步骤S2所获得的基因型向量、环境矩阵，通过矩阵乘法，构建基因-环境交互设计矩阵，并根据对效应的假定，构建基因型向量所对应的焦点变体固定效应向量、环境矩阵所对应的环境随机效应向量、基因-环境交互设计矩阵所对应的基因-环境交互随机效应向量；

基于目标个体的表型信息，构建数量性状表型向量；

基于基因型向量、焦点变体固定效应向量、环境矩阵、环境随机效应向量、基因-环境交互设计矩阵、基因-环境交互随机效应向量，以及数量性状表型向量，构建基因-环境检测线性混合模型；

S4:基于步骤S3所构建的基因-环境检测线性混合模型，采用预处理共轭梯度法、矩估计方法对基因-环境检测线性混合模型进行求解，构建服从于卡方分布的基因-环境互作得分测试统计量，所述基因-环境互作得分测试统计量用于判断所述环境信息是否对目标个体的预设所有基因-环境互作位点的基因产生影响；

S5:基于步骤S4所构建的基因-环境互作得分测试统计量，通过对基因序列中各基因-环境互作位点的随机抽样，计算基因-环境互作得分测试统计量所对应的卡方分布的系数a值，进而获得所有基因-环境互作位点的卡方分布的P值，并预设基因-环境互作阈值，当卡方分布P值小于基因-环境互作阈值，则判定环境信息对该基因-环境互作位点的基因产生影响，即存在基因-环境互作效应，否则判定环境信息对该基因-环境互作位点的基因不产生影响，即不存在基因-环境互作效应。

本发明实施例提供的一种存储软件的计算机可读取介质，所述可读取介质包括能通过一个或多个计算机执行的指令，所述指令在被所述一个或多个计算机执行时，执行如所述权利要求1-4中任意一项所述一种基于大样本的基因-环境互作关联分析方法的操作。

基于本发明实施例提供的基于大样本的基因-环境互作关联分析方法，利用现代CPU多核多线程并行计算，是常用的技术，例如Intel MKL、OpenMP和MPI等；

本发明采用Intel MKL技术实现，在一个实施例中是通过C++语言，链接到IntelMKL优化后的LAPACK和BLAS库等，实现运算速度的提高；

在另一个实施例中是通过Microsoft R Open，直接利用微软版本的R软件，实现并行加速，基于本发明实施例提供的基于大样本的基因-环境互作关联分析方法所对应的程序可同时实现数据并行和任务并行。

针对本发明实施例提供的基于大样本的基因-环境互作关联分析方法，通过模拟实验验证本发明所设计的方法的有效性，具体包括统计功效、均方误差、运算时间等参数；在一个实施例中，模拟实验验证包括Monte Carlo模拟实验分析，以及UK BioBank真实数据分析。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。