CN109448794B - 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法 - Google Patents

Abstract

本发明公开了一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，包括：1、将基因型数据转换为二进制表示的布尔型数据；2、利用逻辑与操作快速地计算任意SNP位点对与表型间条件互信息，取出top‑N节点对，构建包含SNP位点的初始网络图；3、基于初始网络个体，通过随机增加边、删除边、逆转边生成新的个体，直到网络个体数量达到种群大小规模；4、通过遗传算法的三种操作与贝叶斯网络的打分机制，对贝叶斯网络结构进行演化，找到网络结构的最优解，快速准确的获取到影响表型性状的上位性基因位点。本发明可以帮助生物学研究者得到影响特定表型性状的上位性基因位点，进而辅助基因功能挖掘，以及为不同物种的复杂数量性状的遗传基础解析提供借鉴。

Description

一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法

技术领域

本发明涉及生物信息技术领域，尤其涉及一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法。

背景技术

随着人们生活水平和医疗环境的不断提高和改善，那些仅仅由环境因素决定的疾病(比如传染病、营养不良等)基本得到了控制，而复杂疾病和孟德尔遗传病成为了目前影响人类健康的主要疾病。孟德尔遗传病是一种单基因疾病，其遗传过程遵循孟德尔遗传定律，目前研究者利用定位克隆的方法确定了相关遗传基因，基本阐明了其遗传方式。复杂疾病占人类疾病的大约80％以上，对人类健康造成了极大的伤害。哮喘、癌症、糖尿病、高血压、老年痴呆症、类风湿性关节炎、精神分裂症、心脏病、心血管疾病、肥胖、肿瘤等常见慢性疾病，统称为复杂疾病。复杂疾病的病因非常复杂，涉及到环境、基因以及它们之间的相互作用等多种因素。因此，急需阐明复杂疾病的致病原因及遗传机制，给复杂疾病的诊断和治疗提供科学依据，为人类健康提供保障，也具有重要的研究意义。

从生物遗传学的角度看，决定生物复杂性状的遗传因素主要包括三个方面：基因主效应、基因与基因之间的相互作用和基因与环境之间的相互作用。通过生物学大量实验研究发现，控制生物复杂性状的主要原因是基因与基因之间的相互作用。基因与基因之间的相互作用，又称为上位性(Epistasis)，它主要表现为SNP之间的相互作用。同时，随着高通量技术的迅速发展，目前产生了海量的生物数据。利用全基因组关联研究(Genome-wideAssociation Study,GWAS)方法从基因组范围内的数据中筛选出和疾病显著关联的SNPs，从而阐释复杂疾病的遗传机制是当前生物信息学研究的一个热点问题。GWAS方法主要侧重于主效基因的检测，在前期研究尽管利用该方法找到了很多与表型相关的位点，但也只能解释极少数的遗传变异。其中一个最重要的原因就是这些研究忽略了基因与基因之间的相互作用，即上位性。可见，进行上位性位点挖掘是目前解释复杂疾病遗传机制的主要手段。然而，目前上位性检测方法仍然存在计算困难、算法复杂度高、效率低下以及假阳性率高等问题，导致不能准确高效地检测出与疾病相关联的SNP位点及其组合。因此，在全基因组范围内提出更有效、更准确的上位性检测算法具有十分重要的研究意义，也对复杂疾病致病机理的发现、诊断、治疗和预防有着非常重要的作用。

发明内容

本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法。

本发明解决其技术问题所采用的技术方案是：

本发明提供一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，包括以下步骤：

步骤1、对SNP基因型数据，将基因型数据表示为0、1、2形式的数据，0表示纯合子常见基因型，1表示杂合子，2表示纯合子少见基因型；获取待挖掘的基因样本，以样本数为单位分为0、1、2三组，将基因型数据转换为二进制形式表示的布尔型数据0、1；

步骤2、基于信息熵理论，计算任意SNP位点对与表型性状间条件互信息，根据计算的互信息大小对节点对进行排序，取出top-N节点对，构建包含SNP位点对的初始网络图；

步骤3、在不生成环的前提下，对初始网络个体通过随机增加边、删除边、逆转边操作生成下一个网络个体，然后在下一个网络个体的基础上再生成新的网络个体；重复以上生成新网络个体的操作，直到网络个体数量达到初始种群规模大小；

步骤4、通过禁忌搜索优化的遗传算法的三种操作，包括选择、交叉和变异，以及贝叶斯网络的打分机制，对步骤3得到的初始网络种群进行演化，初始网络种群为包括SNP位点的贝叶斯网络，找到网络结构的最优解，从而获取到影响表型性状的上位性基因位点。

进一步地，本发明的该方法还包括对构建的网络进行判断的方法：

步骤5、采用适应度函数作为评判网络个体优劣的标准，采用BIC打分的方法对网络的优劣进行判断。

进一步地，本发明的步骤2和步骤5中，将基因型数据转换为二进制形式表示的布尔型数据，直接利用逻辑与运算对二进制数据进行操作，进而快速的进行节点间条件互信息和贝叶斯网络的BIC打分计算。

进一步地，本发明的步骤2中构建包含SNP位点对的初始网络图的具体方法为：

步骤2.1、设待挖掘的上位性基因位点个数nlocus，对所有位点中nlocus个位点进行排列组合，基于信息熵理论，利用逻辑与操作快速地计算不同组合的nlocus个位点与表型性状间条件互信息；

步骤2.2、根据计算的条件互信息大小对不同的节点对进行排序，取出top-N节点对，其中N的大小根据实验结果进行确定；对于未包含在top-N节点对中SNP位点，选择其第一次出现的节点对，将其插入到top-N节点对中；

步骤2.3、将所有的基因SNP位点看作网络中节点，根据步骤2.2得到的top-N节点对，将不同节点对相应的边插入到网络图中，构建初始网络图。

进一步地，本发明的步骤4中进行演化的具体方法为：

步骤4.1、选择操作；利用贝叶斯网络的评分方法对网络进行打分，将打分最高的最优贝叶斯网络个体放在种群的初始位置，采用轮盘赌选择方法选择网络个体进入下一代；

步骤4.2、禁忌交叉操作；采用多列交叉方法对两个网络进行演化，并进行生成环判断；为了避免普通交叉操作产生早熟现象，利用禁忌搜索具有的记忆功能，进行交叉操作后把产生的子代网络与禁忌表中的个体进行比较；如果不属于禁忌列表，将这个子代网络个体进入到下一代，并将其存储到禁忌表中；如果该个体已经属于禁忌表，则抛弃这个子代个体，重新进行禁忌交叉操作，直到产生的子代不属于禁忌表为止；

步骤4.3、禁忌变异操作；对网络个体以一定的变异概率进行增加边、删除边、逆转边操作，选择使网络评分增加最多的变异，从而得到优化网络结构；利用禁忌搜索具有的记忆功能，将变异产生的可改进当前适应值的劣解存入禁忌表中。

进一步地，本发明的步骤2.1中的具体计算方法为：

当挖掘影响表型Class的k个上位性SNP位点时，I(Class|SNP₁,...SNP_k)表示k个上位性SNP位点与表型Class间的条件互信息，其计算的公式为：

I(Class|SNP₁,...SNP_k)＝H(Class)+H(SNP₁,...SNP_k)-H(Class,SNP₁,...SNP_k)

计算Class的信息熵H(Class)的公式为：

计算k个SNP位点的信息熵H(SNP₁,…,SNP_k)的公式为：

进一步地，本发明的计算BIC评分的具体方法为：

用D表示样本数据，G表示贝叶斯网络结构，根据贝叶斯公式可得：

P(G|D)＝P(D|G)P(G)/P(D)

其中，P(G)表示网络结构的先验知识；

用θ_G表示网络结构的参数，通过边缘积分对上式展开可得：

P(D|G)＝∫P(D|G,θ_G)P(θ_G|G)dθ_G

进而得到贝叶斯网络的BIC评分方法：

其中m表示样本的总数量，n表示变量的个数，r_i表示第i个变量的取值个数，q_i表示第i个变量的父变量的组合个数，m_ijk表示第i个变量取第k个值，其父变量取第j个组合的样本个数。

本发明产生的有益效果是：本发明的基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，首先将基因型数据转换为二进制表示的布尔型数据，利用逻辑与操作快速地计算任意SNP位点对与表型间条件互信息。根据计算的互信息大小，在对SNP位点对进行排序的基础上，取出top-N节点对，构建包含SNP位点的初始网络图。然后，在不生成环的前提下，基于初始网络个体通过随机增加边、删除边、逆转边三个操作生成新的个体，直到达到种群大小规模。将禁忌搜索算法的记忆思想用于遗传算法的进化操作中，利用遗传算法的三种操作(选择、交叉和变异)与贝叶斯网络的BIC打分方法，对贝叶斯网络结构进行演化，找到最优的网络结构，快速准确的获取到影响表型性状的上位性基因位点，辅助基因功能挖掘。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明具体实施的流程示意图；

图2基因型数据表示图；

图3二进制布尔型数据表示图；

图4贝叶斯网络结构矩阵编码表示图；

图5避免环结构产生的处理过程；

图6交叉操作没有产生新的个体示意图；

图7贝叶斯网络变异操作；

图8不同方法2位点上位性检测准确率比较；

图9不同方法2位点上位性检测效率比较；

图10不同方法3位点上位性检测准确率比较。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

1、用0，1，2形式的数据表示基因型数据，如对SNP基因型为AT的数据表示如下：AA用0表示，TT用2表示，AT/TA用1表示。图2所示为11个SNP(SNP_A～SNP_K)对应的4个样本的基因型数据，最后一列Class表示表型性状，其中Class＝1表示case(患病)，Class＝0表示control(对照)。为了提高后续节点间条件互信息以及贝叶斯网络打分的计算效率，将基因型数据转换为二进制形式表示的布尔型数据0/1。将图2中SNP_A～SNP_D的基因型数据转变为二进制形式的布尔型数据格式如图3所示。

图3中，前4列为基因型为0的二进制表示，中间4列为基因型为1的二进制表示，最后4列为基因型为2的二进制表示。当特定SNP在某样本的基因型为0时，则将其前4列对应位置的二进制用1表示，如图2和图3方框中对应数据所示。

2、基于二进制布尔型形式表示的基因型数据，利用信息熵理论计算SNP位点对与表型性状间条件互信息，对节点对进行排序并取出top-N节点对，构建包含SNP位点的初始网络图。

(1)当挖掘影响表型Class的k个上位性SNP位点时，采用式(1)计算k个SNP位点与Class间的条件互信息。采用式(2)计算式(1)中Class的信息熵H(Class)，采用式(3)计算式(1)中k个SNP位点的信息熵H(SNP₁,…,SNP_k)。

I(Class|SNP₁,...SNP_k)＝H(Class)+H(SNP₁,...SNP_k)-H(Class,SNP₁,...SNP_k) (1)

基于二进制形式表示的基因型数据，通过逻辑与操作可以快速的计算节点间条件互信息。例如，采用式(4)计算图3的M_bit中条件互信息I(Class|SNP_B,SNP_C)。

I(Class|SNP_B,SNP_C)＝H(Class)+H(SNP_B,SNP_C)-H(Class,SNP_B,SNP_C) (4)

如采用式(5)计算式(4)中H(SNP_B,SNP_C)。

当对式(5)进行计算时，可以直接利用二进制逻辑与操作，通过对二进制字符串中字符1进行计数求解。如图3的M_bit中下划线部分的二进制值所示，采用式(6)计算式(5)中p(1,1)，

(2)根据上述计算的k个上位性SNP位点与Class间的条件互信息大小，由高到低对不同的节点对进行排序，取出top-N节点对，其中N的大小可以根据实验结果进行确定。对于未包含在top-N节点对中SNP节点，在排好序的节点对中选择其第一次出现的节点对，将其插入到top-N节点对中。

(3)将所有的基因SNP位点看作网络中节点，对上述计算得到的top-N节点对进行循环，将节点对相应的边插入到网络图中，以此类推，构建初始网络图。

3、在贝叶斯网络结构的空间中进行搜索时，本发明采用的遗传禁忌算法中的每个个体对应一个贝叶斯网络结构。用邻接矩阵对贝叶斯网络个体进行表示，设网络中SNP位点个数为n，每个个体表示为n×n的邻接矩阵C。采用0/1的编码方案，如果节点i是节点j的父节点，C_ij＝1，否则，C_ij＝0。如图4所示。

遗传算法初始种群中的贝叶斯网络个体由SNP节点和节点间的边组成。在不生成环的前提下，对初始网络个体通过随机增加边、删除边、逆转边操作生成下一个网络个体。然后在下一个网络个体的基础上再生成新的网络个体，依此类推直到网络个体数量达到初始种群规模大小。遗传算法以初始种群作为起始点开始迭代。

4、在初始网络种群的基础上，通过遗传操作(选择、交叉和变异)和贝叶斯网络的打分机制，对包括SNP位点的贝叶斯网络进行演化，找到网络结构的最优解，进而快速准确的获取到影响表型性状的上位性基因位点。同时，为了增强种群的多样性和获得全局最优解，将禁忌搜索策略应用于遗传算法的交叉和变异进化操作中，也加速了算法的收敛。

(1)选择操作

利用选择操作从当前群体中选出优良的个体，使优良个体作为父代为下一代繁殖后代，选择操作的原则是适应度越强的个体被选中的概率越大。利用贝叶斯网络的评分方法对网络进行打分，对打分较高的网络进行选择，将打分最高的最优贝叶斯网络个体放在种群的初始位置，采用轮盘赌选择方法选择较优的网络个体进入下一代。设第i个网络的适应度值为f_i，则i被选中的概率P_i如Eq.(7)所示，其中N表示种群的大小。

(2)禁忌交叉操作

通过交叉操作可以得到新一代中较好的个体，新个体继承了其父辈个体的特性。为了加快收敛速度，采用多列交叉方法。当交叉使得网络发生变化时，进行生成环判断。假设种群中两个网络个体Individual₁与Individual₂，随机选出Individual₁两列f₁,f₂和Individual₂两列s₁,s₂。将Individual₁的f₁列与Individual₂的s₁列进行交换，将Individual₁的f₂列与Individual₂的s₂列进行交换。即：

进而得到：Individual₁[...s₁...s₂...]，Individual₂[...f₁...f₂...]。进行交叉操作的时候要判断网络中是否有环生成，当Individual₁和Individual₂都没有环结构时，则认为Individual₁和Individual₂为新的子代个体。若交换操作产生环结构，则跳过此位点的交换操作，继续判断下一个，直到两列全部交换完毕。如图5所示。图5中执行过程如下：

<1>Individual₁的第二列(I1.f2)与Individual₂的第一列(I2.s1)逐行进行交换，如图中“交换第一列”部分所示，交叉成功，从而得到交叉完第一列的两个个体。

<2>Individual₁的第三列(I1.f3)与Individual₂的第三列(I2.s3)逐行交换，如图中“交换第二列”部分所示。第一行中，若Individual₁第三列的第一行对应值0与Individual₂第三列的第一行对应值1交换，将导致Individual₁产生环结构，如图中红圈部分所示。算法则跳过第一行，从第二行开始进行交换。最后，交换完毕得到最终两个个体。

在不同的种群子代中，普通的交叉操作可能会产生相同的后代，导致群体中的染色体具有局部相似性，从而使搜索停滞不前，容易产生“早熟”现象。如图6中，在第一次迭代中，Individual₁的I1.f1和Individual₂的I2.s2以及Individual₁的I1.f2和Individual₂的I2.s3进行交叉操作，得到第二次迭代中所示的三个新的子代网络个体。在第二次迭代中，Individual₁的I1.f2和Individual₃的I3.p2以及Individual₁的I1.f3和Individual₃的I3.p3进行交叉操作。得到第三次迭代中的两个子代个体I1、I3与父代重复，这次交叉操作没有产生新的子代个体。

普通交叉操作容易产生早熟现象，为了解决此问题，利用禁忌搜索具有的记忆功能，进行交叉操作后把产生的子代网络与禁忌表中的个体一一进行比较。如果不属于禁忌列表，将这个子代网络个体进入到下一代，并将其存储到禁忌表中。如果该个体已经属于禁忌列表，则抛弃这个子代个体，重新进行禁忌交叉操作，直到产生的子代不属于禁忌表为止。

(3)禁忌变异操作

变异操作首先在群体中随机选择一个个体，对于选中的网络个体以一定的变异概率P_m进行增加边、删除边、逆转边操作，从而增加种群的多样性。如图7中，矩阵中变异操作对应网络中删除节点A与节点C之间的边。在保证不生成环的情况下，选择使网络评分增加最多的变异，从而得到较优的网络结构。普通的变异操作具有较强的随机性且容易破坏适应度较优的网络个体，存在爬山能力差、容易陷入局部最优等问题。本发明利用禁忌搜索具有的记忆功能，在变异产生劣解时，将该解存入禁忌表中，然后在当前基础上进行搜索。该方法可以避免迂回搜索，跳出局部最优，改进变异操作的爬山能力，进而有助于快速查找到较优的网络个体。

5、适应度函数评价。

适应度函数是用于评判网络个体好坏的标准，它决定哪些优秀的个体被保留，哪些较差的个体被淘汰。本发明中遗传禁忌算法是对贝叶斯网络进行打分来计算适应度，进而作为进化搜索的依据，采用常用的BIC打分方法对网络的优劣进行判断。

贝叶斯评分主要是在给定的先验知识和样本数据的情况下，选择后验概率最大的贝叶斯网络结构。用D表示样本数据，G表示贝叶斯网络结构，根据贝叶斯公式可得式(7)。其中，P(G)表示网络结构的先验知识。

P(G|D)＝P(D|G)P(G)/P(D) (7)

用θ_G表示网络结构的参数，通过边缘积分对式(7)展开可得式(8).

P(D|G)＝∫P(D|G,θ_G)P(θ_G|G)dθ_G (8)

进而得到贝叶斯网络的BIC评分方法，如式(9)所示。

其中m表示样本的总数量，n表示变量的个数，r_i表示第i个变量的取值个数，q_i表示第i个变量的父变量的组合个数，m_ijk表示第i个变量取第k个值，其父变量取第j个组合的样本个数。同理，基于二进制形式表示的布尔型数据，利用逻辑与操作进行BIC打分计算，可以实现快速计算。

6、算法终止

当达到最大迭代次数或者最优个体经过k代评分值都保持不变时，本发明中算法结束。否则，用经过选择、交叉、变异后的新一代个体取代上一代个体，并返回继续迭代执行。

7、通过实验来说明基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法的高效性，分别比较二节点和三节点上位性检测的准确率和效率。

下面是应用本发明的方法在GAMETES软件生成数据集上进行上位性位点挖掘的实施例，通过相关的实验来详细说明本发明方法挖掘上位性位点的高效性。GAMETES软件是一款业界常用的用于生成Epistasis模拟数据的软件^[4]，该软件可以快速准确地生成Epistasis模拟数据，通过改变不同的参数生成特定的两位点甚至多位点Epistasis模型。可以设置的参数包括：SNP位点的个数、遗传率(heritability,h²)、最小等位基因频率(MAF)以及患病率(prevalence)等。生成模拟数据的文件中第1行为位点名称，最后1列为Class标签，1表示患病，0表示对照。基因型数据用0，1，2表示，0表示纯合子常见基因型，1表示杂合子，2表示纯合子少见基因型。

将本发明中遗传禁忌优化的贝叶斯网络上位性位点挖掘方法记为Epi-GTBN，实验比较的上位性检测方法包括以下几种：BEAM,AntEpiSeeker,SNPRuler,MDR,BOOST和贝叶斯网络学习方法hill-climbing。通过设置不同的遗传率h²(0.025,0.05,0.1,0.2,0.3,0.4)和最小等位基因频率MAF(0.1,0.2,0.3,0.4)，采用GAMETES软件生成了不同的数据集，每个参数设置下的数据集包括100个文件。采用式(10)计算上位性位点挖掘的准确率。其中Num_edge表示能检测到目标上位性位点的数据集个数。

实验1.2位点上位性检测准确率和效率对比

在本实验中，比较了设置不同遗传率h²和最小等位基因频率MAF情况下上位性位点挖掘的准确率。图8和图9给出了不同方法2位点上位性挖掘的准确率和效率比较。

通过图8可见，在不同的h²和MAF取值情况下，BEAM和hill-climbing方法的2位点上位性检测准确率要远低于其他5种方法。Epi-GTBN，MDR，BOOST和AntEpiSeeker方法的检测准确率最大，基本上是100％，SNPRuler方法的检测准确性要稍微低于其他4种方法。

通过图9可见，AntEpiSeeker方法的2位点上位性检测时间最多，要远多于其他6种方法。BEAM，BOOST和SNPRuler三种方法的检测时间最少，MDR，hill-climbing和Epi-GTBN三种方法所用时间居中。其中，Epi-GTBN所用时间要少于MDR和hill-climbing。在Epi-GTBN方法中，将基因型数据转换为二进制形式表示的布尔型数据，可以利用逻辑与操作快速的计算节点间条件互信息和对网络进行BIC打分。在构建初始网络和对网络进行打分计算时，这可以节省大量的计算时间。

总之，MDR，BOOST和AntEpiSeeker方法的二位点上位性检测准确率与Epi-GTBN方法基本一致，基本上是100％。但MDR和AntEpiSeeker两种方法所用的检测时间明显多于Epi-GTBN方法。另外，AntEpiSeeker方法的参数设置比较复杂，其结果与参数设置紧密相关。BOOST方法仅能用于2位点上位性检测，不能用于多个位点的上位性检测。可见，本发明中Epi-GTBN方法在不影响上位性检测效率的前提下，具有较好的检测准确率。

实验2.3位点上位性检测准确率对比

在本实验中，比较了不同方法在设置不同遗传率h²和最小等位基因频率MAF情况下的3位点上位性挖掘的准确率，如图10所示。

由实验结果可见，图10中3位点上位性检测准确性与图8中2位点上位性检测结果基本一致。BEAM和hill-climbing方法的上位性检测准确率最低。MDR与Epi-GTBN的检测准确率最高，基本是100％，高于SNPRuler方法。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (7)

1.一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，其特征在于，包括以下步骤：

步骤1、对SNP基因型数据，将基因型数据表示为0、1、2形式的数据，0表示纯合子常见基因型，1表示杂合子，2表示纯合子少见基因型；获取待挖掘的基因样本，以样本数为单位分为0、1、2三组，将基因型数据转换为二进制形式表示的布尔型数据0、1；

步骤2、基于信息熵理论，计算任意SNP位点对与表型性状间条件互信息，根据计算的互信息大小对节点对进行排序，取出top-N节点对，构建包含SNP位点对的初始网络图；

步骤3、在不生成环的前提下，对初始网络个体通过随机增加边、删除边、逆转边操作生成下一个网络个体，然后在下一个网络个体的基础上再生成新的网络个体；重复以上生成新网络个体的操作，直到网络个体数量达到初始种群规模大小；

步骤4、通过禁忌搜索优化的遗传算法的三种操作，包括选择、交叉和变异，以及贝叶斯网络的打分机制，对步骤3得到的初始网络种群进行演化，初始网络种群为包括SNP位点的贝叶斯网络，找到网络结构的最优解，从而获取到影响表型性状的上位性基因位点。

2.根据权利要求1所述的基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，其特征在于，该方法还包括对构建的网络进行判断的方法：

步骤5、采用适应度函数作为评判网络个体优劣的标准，采用BIC打分的方法对网络的优劣进行判断。

3.根据权利要求2所述的基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，其特征在于，步骤2和步骤5中，将基因型数据转换为二进制形式表示的布尔型数据，直接利用逻辑与运算对二进制数据进行操作，进而快速的进行节点间条件互信息和贝叶斯网络的BIC打分计算。

4.根据权利要求1所述的基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，其特征在于，步骤2中构建包含SNP位点对的初始网络图的具体方法为：

步骤2.1、设待挖掘的上位性基因位点个数nlocus，对所有位点中nlocus个位点进行排列组合，基于信息熵理论，利用逻辑与操作快速地计算不同组合的nlocus个位点与表型性状间条件互信息；

步骤2.2、根据计算的条件互信息大小对不同的节点对进行排序，取出top-N节点对，其中N的大小根据实验结果进行确定；对于未包含在top-N节点对中SNP位点，选择其第一次出现的节点对，将其插入到top-N节点对中；

步骤2.3、将所有的基因SNP位点看作网络中节点，根据步骤2.2得到的top-N节点对，将不同节点对相应的边插入到网络图中，构建初始网络图。

5.根据权利要求1所述的基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，其特征在于，步骤4中进行演化的具体方法为：

步骤4.1、选择操作；利用贝叶斯网络的评分方法对网络进行打分，将打分最高的最优贝叶斯网络个体放在种群的初始位置，采用轮盘赌选择方法选择网络个体进入下一代；

步骤4.2、禁忌交叉操作；采用多列交叉方法对两个网络进行演化，并进行生成环判断；为了避免普通交叉操作产生早熟现象，利用禁忌搜索具有的记忆功能，进行交叉操作后把产生的子代网络与禁忌表中的个体进行比较；如果不属于禁忌列表，将这个子代网络个体进入到下一代，并将其存储到禁忌表中；如果该个体已经属于禁忌表，则抛弃这个子代个体，重新进行禁忌交叉操作，直到产生的子代不属于禁忌表为止；

步骤4.3、禁忌变异操作；对网络个体以一定的变异概率进行增加边、删除边、逆转边操作，选择使网络评分增加最多的变异，从而得到优化网络结构；利用禁忌搜索具有的记忆功能，将变异产生的可改进当前适应值的劣解存入禁忌表中。

6.根据权利要求4所述的基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，其特征在于，步骤2.1中的具体计算方法为：

当挖掘影响表型Class的k个上位性SNP位点时，I(Class|SNP₁,...SNP_k)表示k个上位性SNP位点与表型Class间的条件互信息，其计算的公式为：

I(Class|SNP₁,...SNP_k)＝H(Class)+H(SNP₁,...SNP_k)-H(Class,SNP₁,...SNP_k)

计算Class的信息熵H(Class)的公式为：

计算k个SNP位点的信息熵H(SNP₁,…,SNP_k)的公式为：

7.根据权利要求2所述的基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法，其特征在于，计算BIC评分的具体方法为：

用D表示样本数据，G表示贝叶斯网络结构，根据贝叶斯公式可得：

P(G|D)＝P(D|G)P(G)/P(D)

其中，P(G)表示网络结构的先验知识；

用θ_G表示网络结构的参数，通过边缘积分对上式展开可得：

P(D|G)＝∫P(D|G,θ_G)P(θ_G|G)dθ_G

进而得到贝叶斯网络的BIC评分方法：

其中m表示样本的总数量，n表示变量的个数，r_i表示第i个变量的取值个数，q_i表示第i个变量的父变量的组合个数，m_ijk表示第i个变量取第k个值，其父变量取第j个组合的样本个数。

Images