CN109390032B - 一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的snp组合的方法 - Google Patents

Abstract

本发明公开了一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，包括如下步骤：步骤一、对群体和个体记录表进行初始化，并且对群体中的个体进行评价指标的计算；其中，所述评价指标包括：ce、gini、k2、g、cec、ginic、k2c、gc；步骤二、对所述评价指标进行排序融合；步骤三、判断群体的进化是否达成终止条件，如果达到了终止条件，则输出进化结果；步骤四、产生一个0～1之间的随机数，判断所述随机数是否大于所述探索概率，根据判断的结果决定用探索或利用的方式产生新个体；步骤五、调整所述新个体，计算调整后的新个体的评价指标，将其追加到个体记录中，判断新个体的八个评价指标是否都大于当前群体中维护的对应评价指标的最大值。

Description

一种基于进化算法在全基因组关联分析的数据中探索与疾病 相关的SNP组合的方法

技术领域

本发明涉及进化算法技术领域，具体涉及一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法。

背景技术

随着高通量基因分型技术的快速发展，越来越多的基于全基因组的单核苷酸多态性(single-nucleotide polymorphism SNP)的患病组-对照组(case-control)数据涌现出来，其规模往往包含数千个样本与数十万的SNP，研究者们期望利用统计、计算机与生物学的各种方法来分析这些数据，找出与疾病关联的SNP，从而进一步探索疾病潜在的患病模型，对疾病的解释、预防、治疗带来更多的指导，这个研究方向称之为全基因组关联分析(Genome-wide association study GWAS)。由于上位性(epistasis)的存在，某些SNP只有当与其他一些SNP组合在一起分析时才会表现出与疾病的相关性，为了对GWAS数据进行全面的分析，不得不考虑SNP组合与疾病的相关性，若只考虑两阶的上位性，即考虑每两个SNP与疾病的关系，其搜索空间也高达数百亿，因此设计出可以快速、准确地在GWAS中找出与疾病相关的SNP组合的方法有着重要的意义。目前，在GWAS数据上探索相关的SNP或SNP组合的方法主要分为五大类：穷举搜索，随机搜索，过滤搜索，基于模型搜索，基于进化算法搜索。穷举搜索，例如MDR算法，它需要计算每一个SNP组合与疾病的相关性，随着考虑的上位性的阶数的提高，其计算复杂性是难以承受的，这种算法往往无法在有效的时间内完成对GWAS数据的分析。随机搜索，例如BEAM，这类算法通过一个随机抽样的过程来衡量每个SNP与疾病的相关性，虽然有相对较高的效率，但往往精度很低。过滤搜索，例如BOOST、FashtChi，这类算法往往分为两个阶段，第一阶段，它们会利用一个简单、快速的评分来衡量每个SNP或每个SNP组合，只有通过检查的SNP或SNP组合可以进入第二阶段，这类算法旨在利用一些快速的指标来减小搜索空间，在第二阶段，只需要在很小的一个搜索空间上对每个SNP组合进行评估，但由于第一阶段所使用的评分往往无法达到令人满意的精度，导致很多重要的SNP无法进入第二阶段进行精密的分析，因此，这类算法的精度仍然无法令人满意。基于模型的搜索，例如ts-RF、AdaBoost，这类算法利用主流的分类器在GWAS数据上构建一个分类精确的模型，而后利用模型在构建完成时对变量的评分来衡量每个SNP与疾病的相关性，但由于模型在构建过程中倾向于选择拥有较高边际效应的变量，因此，这类算法在探索与疾病相关，但无边际效应的SNP时，表现较差。基于进化算法搜索，例如FHSA-SED、MACOED、CSE，这类算法以一种评价SNP组合和疾病之间的关联性的指标作为目标函数，利用遗传算法或蚁群算法这样的进化算法来搜索使目标函数最优的SNP组合，但是进化算法的设计和目标函数的选择成为研究的难点。综上所述，尽管当前存在很多在GWAS数据上探索与疾病相关的SNP或SNP组合的算法，但是归结起来它们的缺陷如下：1、精度上，对SNP组合和疾病的关联性的度量往往只选择一种度量，这种度量无法很好地找到不符合其对致病模型假设的SNP组合。2、内存上，很多算法在分析GWAS数据时对内存的要求很高，大多数计算平台无法运行。3、速度上，多数算法的计算复杂性偏高，甚至很多算法无法在合理的时间内完成对GWAS数据的分析。4、算法的运行结果无法在生物学上得到很好地解释，导致后继的研究者无法很好地利用其结果。

发明内容

本发明设计开发了一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，通过本发明设计的SEE(Sort,Exploration and Exploitation)算法，将融合了八种评价SNP组合的指标作为该进化算法的目标函数，克服了现有方法在精度、内存、速度、可解释性上的缺陷，不仅大大降低了计算所需要的时间和空间，也提高了结果的精度和可解释性。

本发明提供的技术方案为：

一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，包括如下步骤：

步骤一、对群体和个体记录表进行初始化，并且对群体中的个体进行评价指标的计算；其中，所述评价指标包括：ce、gini、k2、g、cec、ginic、k2c、gc；

步骤二、对所述评价指标进行排序融合，通过如下公式计算探索概率er：

式中，numSpecies是当前群体中SNP的去重数目，numPop是群体中个体的数目，l是每个个体的长度，pe为SEE算法的参数；

步骤三、如果群体的进化达到终止条件，产生结果；否则，产生一个0～1之间的随机数；如果所述随机数大于所述探索概率，以利用的方法产生新个体，如果所述随机数小于所述探索概率，以探索的方法产生新个体；

步骤四、计算所述新个体的评价指标，如果新个体的八个评价指标都大于当前群体中维护的对应评价指标的最大值，则认为新个体对于群体是无用的，则再次进行所述步骤三，如果新个体的八个评价指标中至少有一个小于当前群体中维护的对应评价指标的最大值，则认为新个体对于群体为有用的，用该新个体替换当前群体中最差的个体，则再次进行步骤二。

优选的是，在所述步骤一中，

cec的计算方法为：cec(Y,I)＝ce(Y|I)-ce(Y|E)；式中，ce(Y|I)为所述SNP组合的ce值，E为所述SNP组合中最小ce值的SNP，ce(Y|E)为E的ce值；

ginic的计算方法为：ginic(Y,I)＝gini(Y|I)-gini(Y|E)；式中，gini(Y|I)为所述SNP组合的gini值，E为所述SNP组合中最小gini值的SNP，gini(Y|E)为E的gini值；

k2c的计算方法为：

式中，k2(Y,I)为所述SNP组合的k2值，E为所述SNP组合中最小k2值的SNP，k2(Y,E)为E的k2值；

gc的计算方法为：

式中，g(Y,I)为所述SNP组合的g值，E为所述SNP组合中最小g值的SNP，g(Y|E)为E的g值。

优选的是，在所述步骤一中，

ce的计算方法为：

式中，I为所述SNP组合的一个个体；Y为表型、样本状态；mi为互信息；H为信息熵；C为所述SNP组合的可能取值集合；S为表型可取值的集合；p(c,s)为样本中SNP组合取值为c并且表型取值为s的样本比例；p(c)为样本中SNP组合取值为c的样本比例；

gini的计算方法为：

式中，p(s)为样本中表型取值为s的样本比例，C为SNP组合的可能取值集合，p(c)为样本中SNP组合取值为c的样本比例；S是表型可取值的集合；P(s|c)为所述SNP组合取值为c的样本中，表型取值为s的样本的比例；

k2的计算方法为：

式中，C为SNP组合的可能取值集合，SL为SNP组合的长度；C为SNP组合的可能取值集合，S为表型可取值的集合，m_c为样本中SNP组合取值为c的样本数，m_c,s为样本中SNP组合取值为c且表型取值为s的样本数目；

g的计算方法为：

式中，E_c,s为当所述SNP组合与表型独立时，样本中所述SNP组合取值为c且表型取值为s的样本数的期望值，C为所述SNP组合的可能取值集合，S为表型可取值的集合，m_c为样本中所述SNP组合取值为c的样本数，m_s为样本中表型取值为s的样本数，m_c,s为样本中SNP组合取值为c且表型取值为s的样本数目，pvaue_of是在卡方分布下计算变量的pvalue值，自由度为C的长度减1。

优选的是，在所述步骤一中，随机产生多个个体组成初始群体，并且令G＝0。

优选的是，在所述步骤二中，对所述评价指标进行排序融合包括：对所有个体的每个指标进行排序，每个个体会得到8个序号，分别代表在对应的指标下，比它小的个体的数目，对这些序号进行加权累加得到权值，对所述权值进行融合排序。

优选的是，所述权值初始设定默认值为1。

优选的是，在所述步骤三中，当G超过设定的上限时，达到所述终止条件。

优选的是，在所述步骤四中，以探索的方法产生新个体过程包括：

在区间[0,numPop)上随机产生两个整数，将其中较大的整数赋值给第一变量A；在区间[0,1)上随机产生一个实数，赋值给第二变量B；

如果第二变量小于第一变量和numPop的商，则随机产生一个个体，e即是SEE以探索的方法产生的新个体；否则，则将群体上第A个个体赋值给e，在区间[0,o)上随机产生一个整数第三变量C，在区间[0,n)上随机产生一个整数，赋值给e的第C个元素，e即是SEE以探索的方法产生的新个体；

其中，numPop为群体中个体的个数；o为每个个体包含的元素的个数；n为数据中SNP的数目。

优选的是，在所述步骤四中，以利用的方法产生新个体过程包括：

在区间[0,numPop)上随机产生两个随机整数，将较小的那个赋值给第四变量D，在区间[0,numPop)上随机产生两个随机整数，将较小的那个赋值给第五变量E，将群体pop的第D个元素赋值给e，在区间[0,o)上随机产生两个实数F和G，将e的第F个元素替换为群体pop中第E个个体的第G个元素，e即是SEE以利用的方法产生的新个体；

其中，numPop为群体中个体的个数；o为每个个体包含的元素的个数；n为数据中SNP的数目。

本发明与现有技术相比较所具有的有益效果：本发明针对当前算法的精度不足的问题，提出通过排序的方法融合8个评价指标，从而提高精度。针对当前算法的空间复杂性高的问题，使用类似BOOST的存储结构的同时，调整个体时所使用的算法的空间复杂性很小，降低了内存的需求，这使得SEE软件可以在大多数计算平台上稳定地运行。针对当前算法速度上的不足，需要设计新的进化算法，结合调整新个体的算法，加快了分析GWAS数据的速度。针对结果可解释性不足的问题，本发明将SNP映射到基因上，利用cytoscape软件绘制网络，增加了结果的可解释性。

附图说明

图1为SEE算法的主体流程。

图2为通过举例说明如何通过排序来融合8个评价指标。

图3(a)为通过举例说明如何对新个体进行调整。

图3(b)为通过举例说明如何对新个体进行调整。

图3(c)为通过举例说明如何对新个体进行调整。

图3(d)为通过举例说明如何对新个体进行调整。

图4为利用SEE算法在CD数据集上进行分析所得到的部分结果。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

如图1所示，本发明提供了一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，具体步骤如下：

步骤S110、初始化群体和个体记录表；

随机生成若干(numPop，用户指定的参数)个个体(SNP组合)组成初始的群体，令G＝0(在应用进化算法时，通常记录群体的进化次数G，初始群体时，G＝0，进化算法每产生一个新的个体，G加1，SEE算法运行时需指定参数maxIter，当G等于maxIter时，SEE运行终止，更大的maxIter会让SEE的结果更好，但也需要更多的运行时间)，初始化个体(SNP组合)记录表tracing为空，而后将群体中现有的所有个体追加到tracing中。

步骤S120、对群体中每个个体计算其评价指标；

为了评价一个SNP组合(个体)与疾病的关系，本发明定义了8个评价指标，其中前四个指标为该领域常用的用于评价一个SNP组合和疾病关系的度量，后四种为本发明中的设计，旨在衡量一个SNP组合和疾病的关系是否是由于其中某个SNP拥有超强的边际效应而导致的，这八个指标越小，所代表的SNP组合和疾病的关系越显著、明显，对于每个个体，这8个指标的值都将被逐一计算。

在本发明中，这8个指标分别为：

(1)、ce(Conditional Entropy)是来自于对mi(互信息)的推导，多年来，很多研究者使用mi作为衡量SNP组合和疾病关系的一种指标，由于在对一个GWAS数据进行分析的过程中，H(Y)始终保持不变，所以ce可以成为mi的等价替换，ce的计算方法如下：

mi(I,Y)＝H(Y)-ce(Y|I)

式中，I为一个个体(SNP组合)；Y为表型、样本状态；mi为互信息；H为信息熵；C为SNP组合(个体I)的可能取值集合；S为表型可取值的集合；p(c,s)样本中SNP组合取值为c并且表型取值为s的样本比例；p(c)为样本中SNP组合取值为c的样本比例；

(2)、gini是CART(分类回归树)所使用的一种衡量变量之间关系的方法，在该领域，它也被常用于衡量SNP组合和疾病之间的关系，gini的计算方法如下：

gini(Y|I)＝Gini(Y,I)；

式中，p(s)为样本中表型取值为s的样本比例；C为SNP组合的可能取值集合；p(c)为样本中SNP组合取值为c的样本比例；S为表型可取值的集合；P(s|c)为在所有SNP组合取值为c的样本中，表型取值为s的样本的比例；

(3)、k2是一种基于贝叶斯网络的评分准则，也经常被用做衡量SNP组合和疾病之间关系的指标，k2的计算方法如下：

式中，C为SNP组合的可能取值集合；SL为SNP组合的长度(包含多少个SNP)；C为SNP组合的可能取值集合；S为表型可取值的集合；m_c为样本中SNP组合取值为c的样本数；m_c,s为样本中SNP组合取值为c且表型取值为s的样本数目；

(4)、g是指G-test的pvalue值，G-test作为一种检验两个变量是否互相独立的统计方法广泛地应用在判断SNP组合与疾病之间的关系，g的计算方法如下：

g(Y,I)＝pvalue_of(G-statistic(Y,I))；

式中，E_c,s为当SNP组合与表型独立时，样本中SNP组合取值为c且表型取值为s的样本数的期望值；C为SNP组合的可能取值集合；S为表型可取值的集合；m_c为样本中SNP组合取值为c的样本数；m_s为样本中表型取值为s的样本数；m_c,s为样本中SNP组合取值为c且表型取值为s的样本数目；pvaue_of为在卡方分布下计算变量的pvalue值，自由度为C的长度减1；

(5)、cec是用来衡量在ce这种指标下SNP组合的分数比它拥有最小ce的SNP还要小多少，cec的计算方法如下：

cec(Y,I)＝ce(Y|I)-ce(Y|E)；

式中，ce(Y|I)为该SNP组合的ce值；E为该SNP组合中拥有最小ce值的SNP；ce(Y|E)为E的ce值；

(6)、ginic用来衡量在gini这种指标下SNP组合的分数比它拥有最小gini的SNP还要小多少，ginic的计算方法如下：

ginic(Y,I)＝gini(Y|I)-gini(Y|E)；

式中，gini(Y|I)为该SNP组合的gini值；E为该SNP组合中拥有最小gini值的SNP；gini(Y|E)为E的gini值；

(7)、k2c是用来衡量在k2这种指标下SNP组合的分数比它拥有最小k2的SNP还要小多少，k2c的计算方法如下：

式中，k2(Y,I)为该SNP组合的k2值；E为该SNP组合中拥有最小k2值的SNP；k2(Y,E)为E的k2值；

(8)、gc是用来衡量在g这种指标下SNP组合的分数比它拥有最小g的SNP还要小多少，g的计算方法如下：

式中，g(Y,I)为该SNP组合的g值；E为该SNP组合中拥有最小g值的SNP；g(Y|E)是E的g值。

步骤S130、排序群体，计算探索概率；

为了同时使用这8个指标，本发明中需要融合这8个指标为一个新的指标，然而这8个指标来自不同的领域，取值范围、分布非常不同，简单地相加或相乘无法有效地利用这8个指标，本发明设计了一种基于排序的方法来融合这8个评价指标，图2举例阐述了本发明所设计的方法，为了让SEE算法可以找到最优的解，进化的过程中必须同时考虑“S160、探索”和“S170、利用”这两个过程，它们分别代表丰富群体中SNP的数目以防止群体陷在局部最优的算法和充分利用群体中的SNP组合产生更好的解的方法，它们是互补、对立的两种方法，为了找到最优解，缺一不可。

当群体的SNP的数目很多时，SEE应该更偏向调用“S160、探索”，反之，应该更多地调用“S170、利用”；针对这种思想，本发明设计了探索概率(er)这个变量，其计算方法如下：

er＝1-coverage^pe

式中，numSpecies为当前群体中SNP的数目(去重)，numPop是群体中个体的数目，l是每个个体的长度，pe为SEE算法的其中一个参数，需要用户指定，pe越大，在整个进化的过程中，更偏向于用“S160、探索”来产生新个体，pe越小，则更偏向于“S170、利用”。

步骤S140、是否达成终止条件；

判断是否终止群体的进化，一般SEE会持续维持G变量，算法刚开始时，G＝0，随着SEE的运行，G会越来越大，当它达到本发明预先设定的上限maxIter(由用户指定的参数，maxIter越大，SEE运行时间越长，结果越好)时，SEE算法终止群体的进化。

步骤S150、产生一个0到1之间的随机数；

产生一个0到1之间的随机数，若其小于探索概率，执行“S160、探索”，否则，执行“S170、利用”。

步骤S160、以探索的方法产生新个体；

当群体中SNP的数目较少时，SEE偏向于调用该方法来产生新的个体，当由该方法产生的个体被群体收录时，一般会增加群体的SNP数目(去重)，该方法的详细描述如下：

算法1以探索(Exploration)的方法产生新个体：

定义：numPop是群体中个体的个数，由用户运行SEE时指定；o是每个个体包含的元素的个数，由用户运行SEE时指定；n是数据中SNP的数目，由SEE分析的数据的内容决定。

1.在区间[0,numPop)上随机产生两个整数，将其中较大的整数赋值给i1；

2.在区间[0,1)上随机产生一个实数，赋值给ran；

3.如果ran小于i1/numPop，执行4，否则，执行5；

4.随机产生一个个体e，即随机产生e的每一个元素，每个元素都在区间[0,n)上，执行8；

5.将群体pop上第i1个个体pop[i1]赋值给e；

6.在区间[0,o)上随机产生一个整数i2；

7.在区间[0,n)上随机产生一个整数，赋值给e的第i2个元素，即e[i2]；

8.e即是SEE以探索的方法产生的新个体。

该方法的伪代码如下：

步骤S170、以利用的方法产生新个体；

当群体中SNP的数目较多时，SEE偏向于调用该方法来产生新的个体，其思想主要参考于遗传算法中的交叉操作，该方法的详细描述如下：

算法2以利用(Exploitation)的方法产生新个体：

定义：pop，numPop，o的定义同算法1。

1.在区间[0,numPop)上随机产生两个随机整数，将较小的那个赋值给i1；

2.在区间[0,numPop)上随机产生两个随机整数，将较小的那个赋值给i2；

3.将群体pop的第i1个元素赋值给e，即e＝pop[i1]；

4.在区间[0,o)上随机产生两个实数i3和i4；

5.将e的第i3个元素替换为群体pop中第i2个个体的第i4个元素，即e[i3]＝pop[i2][i4]；

6.e即是SEE以利用的方法产生的新个体；

该方法的伪代码如下：

步骤S180、调整新个体；

通过探索或者利用产生的新个体未必可以直接作为一个可能的解添加到群体中，本发明需要对其进行调整，使它成为一个可能的解，并根据个体记录表的内容，避免重复计算相同的解的同时，充分利用重复的解，即当进化算法重复地访问同一个解时，与这个重复的解相邻的解会被优先选择成为调整后的个体，该思想有效地加快了进化算法对全局最优的搜索。图3举例说明调整新个体是如何进行的，调整的详细描述如下：

算法3调整个体：

定义：x是输入的、待调整的个体；stepInTable为用户指定的参数，用于决定调整一个个体时允许的最大距离；tracing是个体记录表，由SEE运行时维护，存储了所有曾经被计算过的个体；n，o的定义同算法1。

1.初始化两个空的队列，waiting0和waiting1；

2.如果x包含两个一样的元素，调整失败，返回，否则，执行3；

3.如果x不在tracing中，调整成功，返回x作为调整后的个体，否则，执行4；

4.将x添加到waiting0和waiting1，初始化i＝0；

5.如果waiting1不是空，从waiting1中取出一个个体赋值给x，执行6，否则执行10；

6.如果x不在tracing中，调整成功，x就是调整后的个体，清空waiting1和waiting0，返回，否则，执行7；

7.初始化j＝0；

8.将x赋值给xx，xx[j]＝xx[j]+1；

9.如果xx[j]在[0,n)之上，并且xx之上的所有元素满足由小到大的顺序，将xx添加到waiting1，j＝j+1，如果j<o，执行8，否则执行10；

10.如果waiting0不是空，从waiting0中取出一个个体赋值给x，执行11，否则执行15；

11.如果x不在tracing中，调整成功，x就是调整后的个体，清空waiting1和waiting0，返回，否则，执行12；

12.初始化j＝0；

13.将x赋值给xx，xx[j]＝xx[j]-1；

14.如果xx[j]在[0,n)之上，并且xx之上的所有元素满足由小到大的顺序，将xx添加到waiting0，j＝j+1，如果j<o，执行13，否则执行15；

15.i＝i+1，如果i<stepInTable，执行5，否则，调整失败，清空waiting1和waiting0，返回；

调整的伪代码如下：

S190、计算调整后的新个体，将其追加到个体记录表中；

计算调整后的新个体的8个评价指标，将这个个体追加到个体记录表中，记下这个个体已经被计算过了。

S200、新个体对于群体是否有用；

SEE算法在群体的进化过程中会持续维护群体中8个评价指标的最大值(SEE算法运行时，每个个体都有8个值，对应于8个指标，对于每个指标，取所有个体中最大的值，这8个最大值会随着SEE算法的运行、群体的进化而不断更新)，如果新个体的8个评价指标中的任何一个小于当前群体中对应指标的最大值，本发明认为这个新个体对于群体是有用的，否则，这个新个体是无用的。

S210、以新个体替换群体中最差的个体；

以新个体替换群体中最差的个体，由于群体是根据rankSum排序的，由小到大，最差的个体便是群体中最后一个个体，即拥有最大的rankSum的个体，以新个体替换群体中最后一个个体即可。

S220、产生结果。

SEE算法会为8个评价指标中的每个提供阈值，用户可以通过运行程序时指定参数来设置这些阈值，当算法运行结束时，只有8个评价指标的值都不大于其对应阈值的个体才会作为结果返回；此外，鉴于仅以“与疾病相关的SNP组合”作为结果无法在生物学上得到很好地解释，导致后继的研究者无法很好地利用其结果，因此，本发明利用NCBI数据库将SNP转化为基因，将SNP组合转化为基因组合，其结果有更好的解释性，方便后继研究者的利用，图4展示了本发明利用SEE算法在CD数据集上进行分析所得到的部分结果。

实施例

如图4所示，通过本发明中的SEE算法分析了WTCCC1七种疾病数据集中的CD(Crohn's Disease)数据集，探索与CD相关的2阶的SNP组合，最后将结果中的SNP根据NCBI转换为基因标识，得到一系列基因对，此图为根据这些基因对绘制的网络，其中每一条边代表在相连的这两个基因上，存在至少4个SNP组合和CD是相关的，由图4可以清晰地看出SEE算法判定对于CD这种疾病，LDB2、LOC107986262、RRP15、SMG1P5等基因可能起着一些关键的作用，值得进一步研究。

本发明的技术方案针对当前算法的精度不足的问题，提出利用排序的方法融合8个不同的评价SNP组合与疾病关系的指标，其中前4个指标ce、gini、k2、g来源于已经在该领域广泛应用的指标，cec、ginic、k2c、gc是本发明设计用于评价SNP组合上位性的新指标，旨在衡量一个SNP组合和疾病的关系是否是由于其中某个SNP拥有超强的边际效应而导致的，这8个指标越小，SNP组合和疾病的关系越强烈，为了融合这8个指标，本发明在群体上分别排序这8个指标，再针对每个个体取所有指标序号的加权和，得到最终的融合后的指标rankSum，以此作为评价SNP组合的标准，融合多个指标不仅提高了结果的精度，也使进化算法在运行的过程中群体的信息量更大，对加快算法的运行也起到了一定的作用。

本发明的技术方案针对当前算法往往需要超大内存的问题，使用了类似BOOST的存储结构，利用“位”来存储GWAS数据，大大降低了所需的内存空间，对于规模约为5000个样本、500000个SNP的GWAS数据，SEE软件在运行时所需内存约为1GB，这使得大多数计算平台都可以利用SEE软件进行GWAS数据的分析。

本发明的技术方案针对当前算法的速度不足的问题，提出了一种新的进化算法，通过“探索”和“利用”这两个过程不断进化群体，而通过群体中SNP的数目来决定调用这两个过程产生新个体的概率，SEE算法利用个体记录表对新个体进行“调整”，不仅不会重复计算同一个SNP组合，也充分利用了群体的信息，即当进化算法反复访问同一个个体时，与这个个体相邻的个体也应该被优先访问。此外，由于本发明用“位”来存储数据，计算8个指标的过程也充分利用了“位操作”，这也使SEE软件的速度得到了巨大的提升。

本发明的技术方案针对当前算法的可解释性不足的问题，提出了将SNP映射到基因上来展示结果，在找到和疾病相关的SNP组合之后，本发明利用NCBI数据库将SNP映射到基因上，将结果转化为基因的组合，接着，本发明利用cytoscape软件绘制了基因与基因的网络，从中可以看到哪些基因在这个疾病中有较大的影响，提高了算法的可解释性。

如图2所示，在另一种实施例中，在步骤S130中，通过排序来融合8个评价指标。首先本发明对所有个体的每个指标进行排序，每个个体会得到8个序号，分别代表在对应的指标下，比它小的个体的数目，对这些序号进行加权累加得到的值，便是本发明融合了8个评价指标所产生的指标，这些权值控制着SEE算法更重视哪些指标或舍弃哪些指标，本发明推荐所有的权值都设为默认值1。

如图3所示，在另一种实施例中，在步骤S180中，调整新个体过程具体实施过程包括：假设个体的长度为2，并且由利用或探索产生的新个体是[3,4]。在图3(b)，3(c)，3(d)中，F代表这个解是不可行的，因为本发明要求个体中的SNP是增序且不相同，x代表这个个体已经被计算过了，存储在个体记录表中，o代表即将会被调整算法返回的个体，调整结束后将返回标记为o的个体中的一个；图3(a)展示了每一个可行解与[3,4]之间的曼哈顿距离，在图3(b)中，当个体[3,4]未被计算过，返回[3,4]作为调整后的个体，在图3(c)中，当个体[3,4]已经被计算过了，返回距离[3,4]最近的且没有被计算过的个体[3,5]作为调整后的个体，在图3(d)中，当[3,4]以及与它曼哈顿距离为1的个体都被计算过，返回与它距离为2并且未被计算过的个体作为调整后的个体，在SEE算法中，stepInTable是用来控制调整深度的参数，当不存在曼哈顿距离小于stepInTable且未被计算过的个体，返回调整失败。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (4)

1.一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，其特征在于，包括如下步骤：

步骤一、对群体和个体记录表进行初始化，并且对群体中的个体进行评价指标的计算；其中，所述评价指标包括：ce、gini、k2、g、cec、ginic、k2c、gc；

步骤二、对所述评价指标进行排序融合，通过如下公式计算探索概率er：

式中，numSpecies是当前群体中SNP的去重数目，numPop是群体中个体的数目，l是每个个体的长度，pe为参数；

步骤三、如果群体的进化达到终止条件，产生结果；否则，产生一个0～1之间的随机数；如果所述随机数大于所述探索概率，以利用的方法产生新个体，如果所述随机数小于所述探索概率，以探索的方法产生新个体；

步骤四、计算所述新个体的评价指标，如果新个体的八个评价指标都大于当前群体中维护的对应评价指标的最大值，则认为新个体对于群体是无用的，则再次进行所述步骤三，如果新个体的八个评价指标中至少有一个小于当前群体中维护的对应评价指标的最大值，则认为新个体对于群体为有用的，用该新个体替换当前群体中最差的个体，则再次进行步骤二；

在所述步骤一中，

ce的计算方法为：

式中，I为所述SNP组合的一个个体；Y为表型、样本状态；mi为互信息；H为信息熵；C为所述SNP组合的可能取值集合；S为表型可取值的集合；p(c,s)为样本中SNP组合取值为c并且表型取值为s的样本比例；p(c)为样本中SNP组合取值为c的样本比例；

gini的计算方法为：

式中，p(s)为样本中表型取值为s的样本比例，C为SNP组合的可能取值集合，p(c)为样本中SNP组合取值为c的样本比例；S是表型可取值的集合；P(s|c)为所述SNP组合取值为c的样本中，表型取值为s的样本的比例；

k2的计算方法为：

式中，C为SNP组合的可能取值集合，SL为SNP组合的长度；C为SNP组合的可能取值集合，S为表型可取值的集合，m_c为样本中SNP组合取值为c的样本数，m_c,s为样本中SNP组合取值为c且表型取值为s的样本数目；

g的计算方法为：

式中，E_c,s为当所述SNP组合与表型独立时，样本中所述SNP组合取值为c且表型取值为s的样本数的期望值，C为所述SNP组合的可能取值集合，S为表型可取值的集合，m_c为样本中所述SNP组合取值为c的样本数，m_s为样本中表型取值为s的样本数，m_c,s为样本中SNP组合取值为c且表型取值为s的样本数目，pvalue_of是在卡方分布下计算变量的pvalue值，自由度为C的长度减1，G为群体的进化次数；

cec的计算方法为：cec(Y,I)＝ce(Y|I)-ce(Y|E)；式中，ce(Y|I)为所述SNP组合的ce值，E为所述SNP组合中最小ce值的SNP，ce(Y|E)为E的ce值；

ginic的计算方法为：ginic(Y,I)＝gini(Y|I)-gini(Y|E)；式中，gini(Y|I)为所述SNP组合的gini值，E为所述SNP组合中最小gini值的SNP，gini(Y|E)为E的gini值；

k2c的计算方法为：

式中，k2(Y,I)为所述SNP组合的k2值，E为所述SNP组合中最小k2值的SNP，k2(Y,E)为E的k2值；

gc的计算方法为：

式中，g(Y,I)为所述SNP组合的g值，E为所述SNP组合中最小g值的SNP，g(Y|E)为E的g值；

在所述步骤二中，对所述评价指标进行排序融合包括：对所有个体的每个指标进行排序，每个个体会得到8个序号，分别代表在对应的指标下，比它小的个体的数目，对这些序号进行加权累加得到权值，对所述权值进行融合排序；

在所述步骤三中，以探索的方法产生新个体过程包括：

在区间[0,numPop)上随机产生两个整数，将其中较大的整数赋值给第一变量A；在区间[0,1)上随机产生一个实数，赋值给第二变量B；

如果第二变量小于第一变量和numPop的商，则随机产生一个个体e，e即是以探索的方法产生的新个体；否则，则将群体上第A个个体赋值给e，在区间[0,o)上随机产生一个整数第三变量C，在区间[0,n)上随机产生一个整数，赋值给e的第C个元素，e即是以探索的方法产生的新个体；

其中，numPop为群体中个体的个数；o为每个个体包含的元素的个数；n为数据中SNP的数目；

以利用的方法产生新个体过程包括：

在区间[0,numPop)上随机产生两个随机整数，将较小的那个赋值给第四变量D，在区间[0,numPop)上随机产生两个随机整数，将较小的那个赋值给第五变量E，将群体pop的第D个元素赋值给e，在区间[0,o)上随机产生两个实数F和G，将e的第F个元素替换为群体pop中第E个个体的第G个元素，e即是以利用的方法产生的新个体；

其中，numPop为群体中个体的个数；o为每个个体包含的元素的个数；n为数据中SNP的数目。

2.如权利要求1所述的基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，其特征在于，在所述步骤一中，随机产生多个个体组成初始群体，并且令G＝0。

3.如权利要求2所述的基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，其特征在于，所述权值初始设定默认值为1。

4.如权利要求1所述的基于进化算法在全基因组关联分析的数据中探索与疾病相关的SNP组合的方法，其特征在于，在所述步骤三中，当G超过设定的上限时，达到所述终止条件。

Images