CN110211640B - 一种基于gpu并行计算的复杂疾病基因互作关联分析方法 - Google Patents

Abstract

本发明提供一种基于GPU并行计算的复杂疾病基因互作关联分析方法，包括如下步骤：步骤一、基于背景控制的主效+上位性检测算法研究，采用两阶段策略：第一阶段，在所有主效和上位性效应中，剔除无关变量，保留较少的变量；第二阶段，将通过第一阶段筛选的变量放入模型中进行最终确定变量选择方法；步骤二、建立基于GPU并行计算的加速算法，得到MRMLM‑GPU新方法；步骤三、MRMLM‑GPU新方法的性能分析。本发明是在MRMLM算法上进行背景控制的上位性拓展以及GPU加速改进，MRMLM算法上位性初筛阶段，SNP以及SNP对之间相互独立，PyCUDA提供了更为简单的编程方式实现GPU的并行。本发明为预测人类疾病的互作基因提供了可行的新方法。

Description

一种基于GPU并行计算的复杂疾病基因互作关联分析方法

技术领域

本发明涉及医学统计学技术领域，具体涉及一种基于GPU并行计算的复杂疾病基因互作关联分析方法。

背景技术

全基因组关联分析(genome-wide association study,GWAS)成功的揭示了数以万计的单核苷酸多态性(single-nucleotide polymorphisms,SNPs)与复杂疾病性状的关联。当前，高通量技术能够测定每个个体高达约500万个SNP，通过从参考数据中插补未测定的SNP可将其增加一倍以上。全基因组一次检测一个SNP的分析方法着重关注单个位点的边际效应，然而，对于大多数常见疾病，几乎所有的关联SNP都具有较小的效应，并且共同解释了表型可遗传变异很少的一部分。这种现象常常造成“缺失遗传率(missingheritability)”的难题。虽然很多因素如遗传变异、结构变异、表观遗传学、基因-环境相互作用可能导致缺失遗传率，但是由于生物***的复杂性，基因-基因相互作用(上位性)被认为是多因子疾病遗传学的重要组成部分。因此，复杂疾病基因互作的研究具有重要的意义。

在人类GWAS中，单核苷酸多态性(SNP)的数量大约为10⁶-10⁷数量级，因此，需要检测的SNP-SNP对的数量可以达到10¹²-10¹⁴的数量级，并且在一些大的研究队列中收集的个体有数十万人。用于数据分析的软件和计算策略需要重新设计以适应如此大量的信息。面对海量数据的挑战，一种可行的方式是通过并行算法设计提高运算速度。目前的图形处理单元(graphics processing units,GPUs)被设计为大规模并行处理器，其提供比中央处理器(central processing units,CPUs)更强的计算能力。GPU是一个强大的计算硬件而且价格可以接受，在笔记本电脑和台式电脑中使用具有计算能力的通用GPU得到广泛推广。因此，利用GPU强大的并行特性，进行基因-基因互作海量数据的分析，对提高运算效率具有实际意义。

由于需要进行大量的SNP对测试运算，基因-基因互作的计算非常耗时，给统计学算法和计算方法提出了新的挑战。例如，即使是由中等大小的数据集500000个SNPs组成，那么大约1250亿的两两相互作用要执行测试。对于中等大小的数据集，所有成对SNP组合的这种广泛分析可能需要几个小时或几天。大规模数据集甚至在标准计算平台上需要数周或数月时间。由于GWAS数据集的可用性和规模都在迅速增加，因此寻找更快的解决方案对于研究来说非常重要。

为了减轻计算负担，一种处理策略是使用两阶段方法进行SNP-SNP交互扫描。其中，首先筛选出具有显着主效应的SNP子集，然后在候选子集上分析SNP-SNP相互作用。然而，这样就忽视了具有不显著的主效而SNP-SNP交互显著的情况。另一种策略是利用CPU或者GPU多核特性及其集群，大大加快所有SNP-SNP互作搜索的速度。

随着GPU编程平台的发展，越来越多基于GPU并行的GWAS算法开发出。Hu等提出了处理风险上位性的SHEsisEpi算法，该算法能够控制由单位点边际效应引起的干扰，并且对来自WTCCC(Wellcome Trust Case Control Consortium)的50万SNPs躁郁症数据进行了分析。SHEsisEpi算法仅仅用了27h完成了所有SNPs互作的扫描，比基于CPU的分析快了300多倍。其研究结果表明神经通路相关基因ASTN2和SNAP29或ASTN2和PIK4CA或与这些对中的任何基因处于连锁不平衡中的其他基因之间的相互作用是高加索人和中国汉族人群中BPD的风险因子。

Wan等提出了检测基因-基因互作的快速方法BOOST算法，该方法利用快速初筛、对数线性模型检验的策略，穷举检测基因-基因交互作用。由于BOOST算法中在收集列联表时的海量内存操作，以及分析不同SNP对的独立性使其适合在GPU中实现。Yung等修改了BOOST算法输入数据结构并在筛选阶段对计算进行GPU并行化，提出了加速版本GBOOST，并对WTCCC的II型糖尿病数据进行了分析，较之前BOOST版本加速了40多倍。GBOOST没有考虑协变量的影响，这在某些情况下导致辛普森的悖论。Wang等在GBOOST算法的基础上进行了改进，提出了最多可以同时考虑5个协变量(如年纪、性别、吸烟史等)的新版本GBOOST 2.0，其具有更高的统计功效。

除了上述提到的方法，在处理病例-对照(case-control)表型方面还有很多GPU并行算法，如基于多因素降维方法(multifactor dimensionality reduction,MDR)并行化的MDR-GPU和GMDR-GPU算法，基于ROC(receiver operating characteristic)曲线无模型GWIS算法，以及基于皮尔逊相关系数差异和逻辑回归分析的两步策略EPIBLASTER算法等。

然而，在处理连续型性状(continuous trait)或者数量性状(quantitativetrait)方面，相关的GPU并行算法并不多。Hemani等提出了针对连续型性状的上位性检测快速算法EpiGPU算法，这种算法速度上有着令人瞩目的表现，但是其局限于基因型数值输入，不适用于任意实数的输入，如估算基因型数值。KamThong等提出了基于GPU的线性回归上位性检测GLIDE算法，该算法克服了EpiGPU算法数值输入的局限，不仅能够计算基因相互作用系数的统计得分，而且可以计算截距和单位点边际效应统计得分。Arkin等基于欧几里得空间转换和随机投影方法，提出了数量性状上位性检测EPIQ算法。

目前，已提出的大部分基于GPU并行运算的基因-基因互作检测方法是针对病例-对照表型的，在人类研究中像癌细胞增殖、工作记忆相关的大脑活动等连续型性状也很常见，然而，针对连续型表型的算法比较少。即使存在一些算法，如前面提到的EpiGPU算法、GLIDE算法和EPIQ算法，这些研究大都是基于单标记的算法，需要Bonferroni校正，处理微小效应的互作效果不佳。对于剔除单标记边际效应的影响；考虑性别、年龄、群体结构等协变量；效应值精确估计等的研究鲜有报道。

发明内容

本发明要解决的技术问题是提供一种基于GPU并行计算的复杂疾病基因互作关联分析方法，是一种基于多位点随机SNP效应混合线性模型MRMLM的改进算法，将模型拓展到主效+上位性模型，并对算法进行并行化，利用GPU加速，从而弥补和改进现有人类遗传中连续性表型上位性检测方法的不足，在保持高功效、低假阳率的同时，大大提高运算效率。利用新算法可发现更具功能的基因-基因互作致病途径，为清晰复杂疾病的致病机理，提供方法支持。

为解决上述技术问题，本发明的实施例提供一种基于GPU并行计算的复杂疾病基因互作关联分析方法，包括如下步骤：

步骤一、基于背景控制的主效+上位性检测算法研究，采用两阶段策略：

第一阶段，在所有主效和上位性效应中，剔除无关变量，保留少量的变量；

第二阶段，将通过第一阶段筛选的变量放入模型中进行最终确定变量选择方法；

步骤二、建立基于GPU并行计算的加速算法，得到MRMLM-GPU新方法；

步骤三、MRMLM-GPU新方法的性能分析。

其中，所述第一阶段的具体步骤包括：

1)第一阶段变量初筛：在单位点筛选的时候，引入背景控制的思想，将每个SNP标记当成随机效应，同时进行多基因背景控制，把遗传方差剖析为主效效应+主效多基因背景+误差方差三个方差分量，建立混合线性模型，然后在混合线性模型的基础上考虑上位性，则遗传方差剖析为主效效应+上位性效应+主效多基因背景+上位性多基因背景+误差方差五个方差分量，对MRMLM中的变量初筛进行背景控制，对包含上位性的五个方差分量进行快速求解；

2)第一阶段初筛阈值的选取：设定P＝0.01，并结合BIC准则确定。

其中，所述第二阶段的具体步骤为：将通过初筛的SNP标记及SNP-SNP对进行变量选择，用六种常用的变量选择方法进行对比，通过模拟数据，比较最终的统计功效、假阳率和效应值估计等指标，从而选择一种最优的变量选择方法。

其中，步骤一中，初筛阶段使用的模型为：

y＝1μ+Qv+G_mβ_m+G_lβ_l+Z_mu_m+Z_lu_l+ε (1)；

其中，y是n×1维数量性状表型向量，n表示样本个体数目；1表示n×1维的单位向量；μ表示表型均值；Q表示n×c维的固定效应矩阵，包括群体结构或者主成分组分、性别和年龄等；v表示c×1维固定效应向量，但不包括截距μ；G_m表示n×1维假定的数量性状QTN主效基因型，β_m是假定QTN的主效效应；G_l表示n×1维假定的数量性状两个位点QTN_i和QTN_j互作基因型，即QTN_i和QTN_j基因型的叉乘，β_l是假定互作基因型效应；Z_mu_m是主效多基因背景，Z_m＝(z_ij)_n×s是主效多基因背景效应u_m对应的设计矩阵，s是主效标记个数，

是主效多基因背景方差，K_m是主效亲缘系数矩阵；Z_lu_l是互作的多基因背景，Z_l＝(z_i#z_j)_n×q是互作多基因背景效应u_l对应的设计矩阵，“#”表示两个列向量的对应元素分别相乘，

是互作标记个数，

是互作多基因背景方差，K_l是互作亲缘系数矩阵；

表示n×1维剩余误差向量，

是剩余误差方差，I_n表示n×n维的单位阵；

若群体结构的影响存在，通过数量性状表型观察值与群体结构的回归分析，可剔除群体结构的效应影响，校正后的模型为：

基于模型(2)，y_-Q的方差可表示为：

通过Kang等提出的EMMA算法，可以获得λ_m和λ_l的值；B_m为半正定矩阵，可以对其进行特征分解，以简化运算；

其中，Q_m为正交矩阵，Λ_mr是具有正特征值的对角阵，mr＝Rank(B_m)，即为矩阵B_m的秩；Q_m1和Q_m2分别为Q_m的子块矩阵，其维度分别为n×mr和n×(n-mr)，0是对应的零块矩阵，同理，可得

其中，第二阶段中确定初筛阈值以及变量选择方法的步骤为：先将初筛阈值设定为0.01，即P<0.01的主效和上位性通过筛选，然后，通过BIC准则来确定初筛保留哪些变量；下面提到的筛选针对主效和上位性分别进行，运用的法则是相同的，如果P<0.01时的变量个数大于3000，初筛分别选择1000，2000和3000，变量进入第二阶段的变量选择，分别计算三种情况下的BIC值，记为BIC1、BIC2和BIC3，选取BIC＝min{BIC1,BIC2,BIC3}时为初筛变量的个数；如果P<0.01时的变量个数大于2000、小于3000，记为num(P<0.01)，则分别筛选1000，2000和num(P<0.01)变量进入第二阶段的变量筛选，与P<0.01变量个数大于3000时类似，分别计算三种情况下的BIC值，选取BIC值最小时的变量个数为初筛变量个数；如果P<0.01时的变量个数大于1000、小于2000，记为num(P<0.01)，则分别筛选1000和num(P<0.01)变量进入第二阶段的变量筛选，分别计算两种情况下的BIC值，记为BIC1和BIC2，选取BIC＝{BIC1,BIC2}时为初筛变量的个数；如果P<0.01时的变量个数小于1000，记为num(P<0.01)，则进入第二阶段筛选的变量个数为num(P<0.01)；

第二阶段变量选择方法从六种方法里面挑选，分别为LASSO，GroupLASSO，SCAD，SIS，LARS和BayesianLASSO，设置四组模拟，第一组只有上位性效应，第二组包含主效效应+上位性效应，第三组包含上位性效应+多基因背景效应，第四组包含主效效应+上位性效应+多基因背景效应；四组模拟中设置20个上位性位点，其中，遗传率为1％的有10个，遗传率为3％的有10个，第二组和第四组设置10个主效位点，遗传率为1％；第三组和第四组的多基因背景效应遗传率分别设置为5％和10％；通过对模拟数据统计功效和假阳率的综合比较，以及对不同情况的稳定性，来确定第二阶段变量选择方法。

其中，步骤二中，建立基于GPU并行计算的加速算法，得到MRMLM-GPU新方法的具体步骤为：

①准备一台具有CUDA运算能力的NVIDIA显卡，搭建Python运行环境，Windows***下可以使用Anaconda、Visual Studio金和Ecplise中的一种，安装pycuda包，导入GPU函数pycuda.autoinit、pycuda.compiler、pycuda.cumath、pycuda.curandom、pycuda.driver、pycuda.gpuarray和pycuda.tools；

②利用pycuda.driver.In,pycuda.driver.Out以及pycuda.driver.InOut参数进行简化内存和显存之间的数据拷贝；

③通过SourceModule()函数，在函数内用__global__标识符设置线程数量以及线程组织形式，以期将SNP标记批量的分配到线程中，调用syncthreads()函数，使同一块中的线程同步；

④在单标记初筛过程中，并行的难点是每个SNP标记都会执行似然函数的最优值的调用，直接在GPU上重新编写执行最优化函数optim，即用GPU函数调用GPU函数；遇到半正定矩阵求逆的问题，如(4)所示，进行近似运算，利用特征值分解，大大加快运算速度。单标记筛选是加速的核心过程，每个标记之间执行相同的操作过程，相对独立，适合进行GPU并行。

⑤初筛过后，进入多位点模型的标记个数小于等于3000，然后，多位点模型，标记之间不是相互独立的，进行局部GPU并行计算，涉及到矩阵运算的操作移植到GPU上运行；

⑥软件包的交互式界面通过tkinter包来实现，通过模块化设计，界面与实现方法分别完成，然后整合成一个软件包。

其中，步骤三中，MRMLM-GPU新方法的性能分析的步骤为：

首先，要保证一个算法的准确性，将基于背景控制并改进后能进行上位性检测的MRMLM串行算法用一套规模较小的数据运行一下结果；然后，将改写成GPU并行的新方法MRMLM-GPU进行对比，检验结果是否一致；其次，新方法与串行的MRMLM方法进行不同数量级SNP的互作时间比较，做出时间比图形；最后，将新方法与EpiGPU算法、GLIDE算法和EPIQ算法用四套模拟数据进行比较，主要比较统计功效、假阳率和计算时间三个方面；

新方法与其它三种方法运行WTCCC的真实数据，比较一下每种方法找到的SNP-SNP上位性位点的个数，确定新方法与其它三种方法找到的哪些互作位点是一致的。

本发明的上述技术方案的有益效果如下：本发明是在MRMLM算法上进行背景控制的上位性拓展以及GPU加速改进，MRMLM算法上位性初筛阶段，SNP以及SNP对之间相互独立，PyCUDA提供了更为简单的编程方式实现GPU的并行。由于计算能力和运算时间的限制，许多人类GWAS数据只是进行了单个位点的扫描，本发明为预测人类疾病的互作基因提供了可行的新方法。

附图说明

图1为本发明的流程框图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明提供一种基于GPU并行计算的复杂疾病基因互作关联分析方法，其技术路线框图如图1所示，包括如下步骤：

步骤一、基于背景控制的主效+上位性检测算法研究

上位性是造成遗传率缺失的重要因素，目前测序技术的发展，使得SNP标记的个数远远大于样本的个数，尤其是考虑上位性的情况，同时放入一个模型中，模型的维度会极度增加。这种超饱和模型，在存储空间和计算时间上都给变量选择造成极大困难。常常采用两阶段策略：

第一阶段，在所有主效和上位性效应中，剔除无关变量，保留少量的变量；

其中，所述第一阶段的具体步骤包括：

1)第一阶段变量初筛：现在变量初筛的方法大都是基于快速检测单个SNP对的方法。Wen等提出的FASTmrEMMA(2017，Briefings in Bioinformatics)算法在单位点筛选的时候，引入了背景控制的思想，将每个SNP标记当成随机效应，同时进行多基因背景控制，把遗传方差剖析为主效效应+主效多基因背景+误差方差三个方差分量，建立混合线性模型。如果在此基础上考虑上位性，那么遗传方差剖析为主效效应+上位性效应+主效多基因背景+上位性多基因背景+误差方差五个方差分量，对MRMLM(Wang等,2017)中的变量初筛进行背景控制并加以推广，研究如何对包含上位性的五个方差分量进行快速求解。

2)第一阶段初筛阈值的选取：对于单位点变量筛选方法，一般需要进行Bonferroni校正，即0.01/m或者0.05/m(m为标记个数)，这样的阈值过于严格，存在一些具有致病基因位点但计算出的P值略大于校正阈值，还存在一些小效应位点和稀有等位基因位点，这些位点往往不能通过Bonferroni校正的严格筛选，造成重要位点的遗漏。如果选取一个固定的P值作为阈值，如1e-8或者其它，则比较机械化，对于不同SNP数目不同样本大小的数据影响较大。如果进行Permutation选择合适的P值，这样会消耗大量时间，不利于算法加速。研究如何选取初筛阈值，既能保留重要位点，又能够保证较快的运算速度，对于算法最终结果影响重大。因此，初筛阈值选取时，设定P＝0.01，并结合BIC准则确定。

第二阶段，将通过第一阶段筛选的变量放入模型中进行最终确定变量选择方法，是将通过初筛的SNP标记及SNP-SNP对进行变量选择，用六种常用的变量选择方法进行对比，通过模拟数据，比较最终的统计功效、假阳率和效应值估计等指标，从而选择一种最优的变量选择方法。

具体步骤为：先将初筛阈值设定为0.01，即P<0.01的主效和上位性通过筛选，然后，通过BIC准则来确定初筛保留哪些变量；下面提到的筛选针对主效和上位性分别进行，运用的法则是相同的，如果P<0.01时的变量个数大于3000，初筛分别选择1000，2000和3000，变量进入第二阶段的变量选择，分别计算三种情况下的BIC值，记为BIC1、BIC2和BIC3，选取BIC＝min{BIC1,BIC2,BIC3}时为初筛变量的个数；如果P<0.01时的变量个数大于2000、小于3000，记为num(P<0.01)，则分别筛选1000，2000和num(P<0.01)变量进入第二阶段的变量筛选，与P<0.01变量个数大于3000时类似，分别计算三种情况下的BIC值，选取BIC值最小时的变量个数为初筛变量个数；如果P<0.01时的变量个数大于1000、小于2000，记为num(P<0.01)，则分别筛选1000和num(P<0.01)变量进入第二阶段的变量筛选，分别计算两种情况下的BIC值，记为BIC1和BIC2，选取BIC＝{BIC1,BIC2}时为初筛变量的个数；如果P<0.01时的变量个数小于1000，记为num(P<0.01)，则进入第二阶段筛选的变量个数为num(P<0.01)。

第二阶段变量选择方法从六种方法里面挑选，分别为LASSO，GroupLASSO，SCAD，SIS，LARS和BayesianLASSO，设置四组模拟，第一组只有上位性效应，第二组包含主效效应+上位性效应，第三组包含上位性效应+多基因背景效应，第四组包含主效效应+上位性效应+多基因背景效应；四组模拟中设置20个上位性位点，其中，遗传率为1％的有10个，遗传率为3％的有10个，第二组和第四组设置10个主效位点，遗传率为1％；第三组和第四组的多基因背景效应遗传率分别设置为5％和10％；通过对模拟数据统计功效和假阳率的综合比较，以及对不同情况的稳定性，来确定第二阶段变量选择方法。

步骤二、建立基于GPU并行计算的加速算法，得到MRMLM-GPU新方法，具体步骤为：

①准备一台具有CUDA运算能力的NVIDIA显卡，搭建Python运行环境，Windows***下可以使用Anaconda、Visual Studio金和Ecplise中的一种，安装pycuda包，导入GPU函数pycuda.autoinit、pycuda.compiler、pycuda.cumath、pycuda.curandom、pycuda.driver、pycuda.gpuarray和pycuda.tools；

②利用pycuda.driver.In,pycuda.driver.Out以及pycuda.driver.InOut参数进行简化内存和显存之间的数据拷贝；

③通过SourceModule()函数，在函数内用__global__标识符设置线程数量以及线程组织形式，以期将SNP标记批量的分配到线程中，调用syncthreads()函数，使同一块中的线程同步；

④在单标记初筛过程中，并行的难点是每个SNP标记都会执行似然函数的最优值的调用，直接在GPU上重新编写执行最优化函数optim，即用GPU函数调用GPU函数；遇到半正定矩阵求逆的问题，如(4)所示，进行近似运算，利用特征值分解，大大加快运算速度。单标记筛选是加速的核心过程，每个标记之间执行相同的操作过程，相对独立，适合进行GPU并行。

⑤初筛过后，进入多位点模型的标记个数小于等于3000，然后，多位点模型，标记之间不是相互独立的，进行局部GPU并行计算，涉及到矩阵运算的操作移植到GPU上运行；

⑥软件包的交互式界面通过tkinter包来实现，通过模块化设计，界面与实现方法分别完成，然后整合成一个软件包。

OpenCL是异构***并行编程的标准，它既可以在GPU上运行，也支持多核CPU的编程，但是编程难度很大。NVIDIA公司推出了支持GPU编程的CUDA平台，它可以看作C语言的拓展，编程难度适中。基于CUDA的思想，原有的算法要进行改写，适合GPU并行。如果单位点初筛，不进行背景控制，无论是SNP或者是SNP-SNP对都是独立的，借鉴前人的GPU程序，容易改写成GPU并行算法。第一阶段的初筛是最耗时的过程，解决好这一阶段的并行加速，对于整个算法的提速至关重要。

第二阶段的变量选择方法采用局部的GPU改写，能够通过初筛进入变量选择的SNP标记个数，应该不会很多。虽然GPU具有更多的内核，更强的并行计算能力但是，数据在内存与显存之间交换需要额外花费时间。因此，变量个数不多的时候，进行有选择的局部代码GPU改写。

在进行GPU改写的时候，由于需要算法设计要满足并行特性，有些地方可能会进行近似的处理。如果做了近似，需要对GPU加速版本跟原始的串行程序进行比较，至少保证检测到的主效和上位性的位点一致。

步骤三、MRMLM-GPU新方法的性能分析；

其中，MRMLM-GPU新方法的性能分析的步骤为：

首先，要保证一个算法的准确性，将基于背景控制并改进后能进行上位性检测的MRMLM串行算法用一套规模较小的数据运行一下结果；然后，将改写成GPU并行的新方法MRMLM-GPU进行对比，检验结果是否一致；其次，新方法与串行的MRMLM方法进行不同数量级SNP的互作时间比较，做出时间比图形；最后，将新方法与EpiGPU算法、GLIDE算法和EPIQ算法用四套模拟数据进行比较，主要比较统计功效、假阳率和计算时间三个方面；

新方法与其它三种方法运行WTCCC的真实数据，比较一下每种方法找到的SNP-SNP上位性位点的个数，确定新方法与其它三种方法找到的哪些互作位点是一致的。

本发明的步骤一中，

初筛阶段使用的模型为：

y＝1μ+Qv+G_mβ_m+G_lβ_l+Z_mu_m+Z_lu_l+ε (1)；

其中，y是n×1维数量性状表型向量，n表示样本个体数目；1表示n×1维的单位向量；μ表示表型均值；Q表示n×c维的固定效应矩阵，包括群体结构或者主成分组分、性别和年龄等；v表示c×1维固定效应向量，但不包括截距μ；G_m表示n×1维假定的数量性状QTN主效基因型，β_m是假定QTN的主效固定效应；G_l表示n×1维假定的数量性状两个位点QTN_i和QTN_j互作基因型，即QTN_i和QTN_j基因型的叉乘，β_l是假定互作基因型的固定效应；Z_mu_m是主效多基因背景，Z_m＝(z_ij)_n×s是主效多基因背景效应u_m对应的设计矩阵，s是主效标记个数，

是主效多基因背景方差，K_m是主效亲缘系数矩阵；Z_lu_l是互作的多基因背景，Z_l＝(z_i#z_j)_n×q是互作多基因背景效应u_l对应的设计矩阵，

是互作标记个数，

是互作多基因背景方差，K_l是互作亲缘系数矩阵；

表示n×1维剩余误差向量，

是剩余误差方差，I_n表示n×n维的单位阵。

若群体结构的影响存在，通过数量性状表型观察值与群体结构的回归分析，可剔除群体结构的效应影响，校正后的模型为：

基于模型(2)，y_-Q的方差可表示为：

通过Kang等提出的EMMA算法，可以获得λ_m和λ_l的值；B_m为半正定矩阵，可以对其进行特征分解，以简化运算；

其中，Q_m为正交矩阵，Λ_mr是具有正特征值的对角阵，mr＝Rank(B_m)，即为矩阵B_m的秩；Q_m1和Q_m2分别为Q_m的子块矩阵，其维度分别为n×mr和n×(n-mr)，0是对应的零块矩阵，同理，可得

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (5)

1.一种基于GPU并行计算的复杂疾病基因互作关联分析方法，其特征在于，包括如下步骤：

步骤一、基于背景控制的主效+上位性检测算法研究，采用两阶段策略：

第一阶段，在所有主效和上位性效应中，剔除无关变量，保留少量变量；

第二阶段，将通过第一阶段筛选的变量放入模型中进行最终确定变量选择方法；

步骤二、建立基于GPU并行计算的加速算法，得到MRMLM-GPU新方法；

步骤二中，建立基于GPU并行计算的加速算法，得到MRMLM-GPU新方法的具体步骤为：

①准备一台具有CUDA运算能力的NVIDIA显卡，搭建Python运行环境，Windows***下可以使用Anaconda、Visual Studio和Ecplise中的一种，安装pycuda包，导入GPU函数pycuda.autoinit、pycuda.compiler、pycuda.cumath、pycuda.curandom、pycuda.driver、pycuda.gpuarray和pycuda.tools；

②利用pycuda.driver.In,pycuda.driver.Out以及pycuda.driver.InOut参数进行简化内存和显存之间的数据拷贝；

③通过SourceModule()函数，在函数内用__global__标识符设置线程数量以及线程组织形式，以期将SNP标记批量的分配到线程中，调用syncthreads()函数，使同一块中的线程同步；

④直接在GPU上重新编写执行最优化函数optim，即用GPU函数调用GPU函数；

⑤初筛过后，进入多位点模型的标记个数小于等于3000，然后，多位点模型，标记之间不是相互独立的，进行局部GPU并行计算，涉及到矩阵运算的操作移植到GPU上运行。

2.根据权利要求1所述的基于GPU并行计算的复杂疾病基因互作关联分析方法，其特征在于，所述第一阶段的具体步骤包括：

1)第一阶段变量初筛：在单位点筛选的时候，引入背景控制的思想，将每个SNP标记当成随机效应，同时进行多基因背景控制，把遗传方差剖析为主效效应+主效多基因背景+误差方差三个方差分量，建立混合线性模型，然后在混合线性模型的基础上考虑上位性，则遗传方差剖析为主效效应+上位性效应+主效多基因背景+上位性多基因背景+误差方差五个方差分量，对MRMLM中的变量初筛进行背景控制，对包含上位性的五个方差分量进行快速求解；

2)第一阶段初筛阈值的选取：设定P＝0.01，并结合BIC准则确定。

3.根据权利要求1所述的基于GPU并行计算的复杂疾病基因互作关联分析方法，其特征在于，所述第二阶段的具体步骤为：将通过初筛的SNP标记及SNP-SNP对进行变量选择，用六种常用的变量选择方法进行对比，通过模拟数据，比较最终的统计功效、假阳率和效应值估计指标，从而选择一种最优的变量选择方法。

4.根据权利要求1所述的基于GPU并行计算的复杂疾病基因互作关联分析方法，其特征在于，

步骤一中，初筛阶段使用的模型为：

y＝1μ+Qv+G_mβ_m+G_lβ_l+Z_mu_m+Z_lu_l+ε (1)；

其中，y是n×1维数量性状表型向量，n表示样本个体数目；1表示n×1维的单位向量；μ表示表型均值；Q表示n×c维的固定效应矩阵，包括群体结构或者主成分组分、性别和年龄；v表示c×1维固定效应向量，但不包括截距μ；G_m表示n×1维假定的数量性状QTN主效基因型，β_m是假定QTN的主效效应；G_l表示n×1维假定的数量性状两个位点QTN_i和QTN_j互作基因型，即QTN_i和QTN_j基因型的叉乘，β_l是假定互作基因型效应；Z_mu_m是主效多基因背景，Z_m＝(z_ij)_n×s是主效多基因背景效应u_m对应的设计矩阵，s是主效标记个数，

是主效多基因背景方差，K_m是主效亲缘系数矩阵；Z_lu_l是互作的多基因背景，Z_l＝(z_i#z_j)_n×q是互作多基因背景效应u_l对应的设计矩阵，

是互作标记个数，

是互作多基因背景方差，K_l是互作亲缘系数矩阵；

表示n×1维剩余误差向量，

是剩余误差方差，I_n表示n×n维的单位阵；

若群体结构的影响存在，通过数量性状表型观察值与群体结构的回归分析，可剔除群体结构的效应影响，校正后的模型为：

基于模型(2)，y_-Q的方差可表示为：

通过EMMA算法，可以获得λ_m和λ_l的值；B_m为半正定矩阵；

其中，Q_m为正交矩阵，Λ_mr是具有正特征值的对角阵，mr＝Rank(B_m)，即为矩阵B_m的秩；Q_m1和Q_m2分别为Q_m的子块矩阵，其维度分别为n×mr和n×(n-mr)，0是对应的零块矩阵，同理，可得

5.根据权利要求1或3所述的基于GPU并行计算的复杂疾病基因互作关联分析方法，其特征在于，第二阶段中确定初筛阈值以及变量选择方法的步骤为：先将初筛阈值设定为0.01，即P<0.01的主效和上位性通过筛选，然后，通过BIC准则来确定初筛保留哪些变量；下面提到的筛选针对主效和上位性分别进行，运用的法则是相同的，如果P<0.01时的变量个数大于3000，初筛分别选择1000，2000和3000，变量进入第二阶段的变量选择，分别计算三种情况下的BIC值，记为BIC1、BIC2和BIC3，选取BIC＝min{BIC1,BIC2,BIC3}时为初筛变量的个数；如果P<0.01时的变量个数大于2000、小于3000，记为num(P<0.01)，则分别筛选1000，2000和num(P<0.01)变量进入第二阶段的变量筛选，与P<0.01变量个数大于3000时类似，分别计算三种情况下的BIC值，选取BIC值最小时的变量个数为初筛变量个数；如果P<0.01时的变量个数大于1000、小于2000，记为num(P<0.01)，则分别筛选1000和num(P<0.01)变量进入第二阶段的变量筛选，分别计算两种情况下的BIC值，记为BIC1和BIC2，选取BIC＝{BIC1,BIC2}时为初筛变量的个数；如果P<0.01时的变量个数小于1000，记为num(P<0.01)，则进入第二阶段筛选的变量个数为num(P<0.01)；

第二阶段变量选择方法从六种方法里面挑选，分别为LASSO，GroupLASSO，SCAD，SIS，LARS和BayesianLASSO，设置四组模拟，第一组只有上位性效应，第二组包含主效效应+上位性效应，第三组包含上位性效应+多基因背景效应，第四组包含主效效应+上位性效应+多基因背景效应；四组模拟中设置20个上位性位点，其中，遗传率为1％的有10个，遗传率为3％的有10个，第二组和第四组设置10个主效位点，遗传率为1％；第三组和第四组的多基因背景效应遗传率分别设置为5％和10％；通过对模拟数据统计功效和假阳率的综合比较，以及对不同情况的稳定性，来确定第二阶段变量选择方法。

Images