CN117649876A - 基于gwo算法在gwas数据上检测与复杂疾病相关snp组合的方法 - Google Patents
基于gwo算法在gwas数据上检测与复杂疾病相关snp组合的方法 Download PDFInfo
- Publication number
- CN117649876A CN117649876A CN202410118311.1A CN202410118311A CN117649876A CN 117649876 A CN117649876 A CN 117649876A CN 202410118311 A CN202410118311 A CN 202410118311A CN 117649876 A CN117649876 A CN 117649876A
- Authority
- CN
- China
- Prior art keywords
- wolf
- snp
- wolves
- algorithm
- samples
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 151
- 201000010099 disease Diseases 0.000 title claims abstract description 111
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 title claims abstract description 111
- 238000000034 method Methods 0.000 title claims abstract description 59
- 241000282461 Canis lupus Species 0.000 claims abstract description 120
- 238000001514 detection method Methods 0.000 claims abstract description 22
- 238000012098 association analyses Methods 0.000 claims abstract description 6
- 241000282421 Canidae Species 0.000 claims description 90
- 230000006870 function Effects 0.000 claims description 59
- 238000004364 calculation method Methods 0.000 claims description 26
- 239000013598 vector Substances 0.000 claims description 15
- 230000035772 mutation Effects 0.000 claims description 12
- 238000012360 testing method Methods 0.000 claims description 9
- 239000002773 nucleotide Substances 0.000 claims description 6
- 238000001162 G-test Methods 0.000 claims description 5
- 238000003491 array Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 abstract description 21
- 238000011160 research Methods 0.000 abstract description 6
- 230000008506 pathogenesis Effects 0.000 abstract description 4
- 238000005516 engineering process Methods 0.000 abstract description 2
- 230000002265 prevention Effects 0.000 abstract description 2
- 238000005065 mining Methods 0.000 abstract 1
- 230000001737 promoting effect Effects 0.000 abstract 1
- 230000001717 pathogenic effect Effects 0.000 description 20
- 230000008569 process Effects 0.000 description 14
- 230000000875 corresponding effect Effects 0.000 description 13
- 238000011156 evaluation Methods 0.000 description 11
- 230000008901 benefit Effects 0.000 description 8
- 230000000694 effects Effects 0.000 description 8
- 108090000623 proteins and genes Proteins 0.000 description 7
- 238000002474 experimental method Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 5
- 238000009825 accumulation Methods 0.000 description 3
- 244000144992 flock Species 0.000 description 3
- 230000002596 correlated effect Effects 0.000 description 2
- 238000010845 search algorithm Methods 0.000 description 2
- 239000000243 solution Substances 0.000 description 2
- 206010020772 Hypertension Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 206010012601 diabetes mellitus Diseases 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000008571 general function Effects 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 208000019622 heart disease Diseases 0.000 description 1
- 238000002347 injection Methods 0.000 description 1
- 239000007924 injection Substances 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000009191 jumping Effects 0.000 description 1
- 238000007477 logistic regression Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 125000003729 nucleotide group Chemical group 0.000 description 1
- 230000003234 polygenic effect Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000013179 statistical model Methods 0.000 description 1
Landscapes
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，属于计算机科学技术与生物学的交叉领域；包括一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法。本发明提出了一种新的方法，用于挖掘数据中蕴含的疾病信息。全基因组关联分析是一种常用的探索复杂疾病的方法，提出高效且准确的检测方法一直是该领域的研究热点之一。本发明提出的方法基于灰狼优化算法实现，能够有效地在该数据上检测与复杂疾病相关的SNP组合。相比现有的其它算法，该方法具有更高的检测能力。本发明的研究成果将有助于进一步探索复杂疾病的致病机理，并有助于推进复杂疾病的解释、预防和治疗。
Description
技术领域
本发明属于计算机科学技术与生物学的交叉领域，尤其涉及一种基于灰狼优化(Grey Wolf Optimizer ，GWO)算法在全基因组关联分析数据上检测与复杂疾病相关SNP组合的方法。
背景技术
复杂疾病包含心脏病、高血压、糖尿病等常见疾病，它们是由多基因、多因素共同引起的疾病，它们很难被预防，一旦发现，也很难治愈，很多甚至无法治愈，只能依靠药物来维持，影响着患者的生活。全基因组关联分析（Genome-wide association study ，GWAS）是一种遗传学研究方法，用于揭示基因与复杂疾病之间的关系，GWAS一般针对某一种复杂疾病，收集数千的正常人与患者，通过检测他们DNA上几十万的SNP（single-nucleotidepolymorphism）的基因型，探索SNP与某种复杂疾病之间的相关性，其结果对于复杂疾病致病机理的探索，复杂疾病的预防与治疗都有着重要的意义。然而，GWAS数据上普遍存在的上位性现象为我们检测与复杂疾病相关的SNP带来了巨大的困难。上位性的原意指的是一种基因相互作用的类型，其中一个基因能够掩盖或支配另一个基因的表现。当在GWAS数据上使用统计的方法检测与疾病相关的SNP时，上位性表现为，很多SNP在单独分析时，其与疾病没有显著的相关性，但当它们与其他的一些SNP组合在一起进行分析时，它们与复杂疾病之间产生了非常显著的相关性，这种现象为检测与疾病相关的SNP带来了巨大的问题。其原因是，复杂疾病是由多基因、多因素共同引起的，往往SNP以组合的方式共同与疾病的发生有关系，上位性普遍存在于GWAS数据中，为了充分地分析复杂疾病的致病机理，就需要去分析数据中每两个、每三个（甚至更多）SNP与疾病之间的相关性，其计算量随着上位性阶数的升高成指数级增长。因此，提出高效、精确的在GWAS数据上检测与复杂疾病相关SNP组合的方法是当前的研究热点之一。
近年来，针对这一问题，提出了非常多的算法，AntEpiSeeker（Ant代表蚂蚁，因为它是一种蚁群算法，Epi代表上位性，因为它是在数据中找上位性SNP组合，Seeker是表示它是一种搜索算法）是一种基于蚁群优化算法检测SNP的方法，通过蚁群算法探索与疾病显著相关的SNP组合，再对蚁群算法的结果做进一步的穷举搜索，获取结果，它提出了用蚁群优化算法检测SNP的方法，证明了基于群智能优化算法检测与复杂疾病相关的SNP的可行性。MACOED（a multi-objective ant colony optimization algorithm for snp epistasisdetection）是一种多目标启发式优化方法，它将逻辑回归和贝叶斯网络方法结合作为蚁群优化算法的目标函数，其实验结果证明了在群智能优化算法中引入多个目标函数可以达到互补的效果，从而提升群智能算法检测致病SNP组合的效率。BOOST（Boolean Operation-based Screening and Testing）是一种快速检测二阶SNP组合与疾病关系的方法，它提出了一种基于统计模型检测显著与疾病相关的SNP组合的下界函数，这个下界函数可以快速地用穷举的方法检测每一对SNP是否有可能与疾病具有显著的相关性，只需要对下界函数的值高于指定阈值的SNP组合进行分析就可以在考虑二阶上位性的前提下检测到与疾病相关的SNP，BOOST还提出了一种在计算机中GWAS数据的二进制存储与计算的方法，可以提高算法在GWAS数据上的计算速度。
现有的在GWAS数据上检测与复杂疾病相关的SNP组合的方法主要有以下缺点：
1.多数基于群智能优化算法的检测方法都是基于一个目标函数设计的，而复杂疾病的潜在致病模型往往非常复杂，而一个目标函数无法充分、合理地度量SNP组合与疾病之间的相关性；
2.绝大多数算法都需要用户指定考虑的上位性的阶数，而复杂疾病的致病机理非常复杂，用户无法预判一个真实的复杂疾病潜在的上位性阶数，为用户的使用带来了不便；
3.绝大多数算法在设计的时候着重考虑了二阶上位性的情况，而忽略了高阶（>=3）上位性SNP组合，针对算法的实验也只在蕴含二阶上位性SNP致病关系的模拟数据上进行，因此绝大多数算法在检测高阶上位性SNP组合的能力上表现欠佳，从而导致无法检测出致病的高阶SNP组合；
4.该领域普遍使用的K2、CE、Gini等用于评测SNP组合与疾病关系的目标函数，普遍存在“样本分的份数越多，计算结果越相关”的问题，导致评测函数上的不公，阻碍了算法对相关SNP组合的检测。
发明内容
本发明目的在于提供一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，以解决在GWAS数据上检测与疾病相关的SNP需要非常大的计算量的技术问题，本发明的宗旨是在同等计算量下，追求检测的精度。
为实现上述目的，本发明的一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法的具体技术方案如下：
复杂疾病是由多基因、多因素共同引起的常见疾病，全基因组关联分析（GWAS，Genome-wide association study）数据蕴含着单核苷酸多态性（SNP， single-nucleotidepolymorphism）与复杂疾病之间的相关性，但由于上位性现象的存在，在GWAS数据上检测与疾病相关的SNP需要非常大的计算量，提出高效、精确的检测算法是当前的研究热点之一。
本发明提出了一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，针对现有方法存在的缺陷，主要的解决方案如下：
1.本发明提出的方法是一种基于多个目标函数实现的灰狼优化算法，同时使用K2、CE、Gini三个函数作为算法优化的目标，并最后使用G-test检测SNP组合与疾病是否显著相关；
2.本发明提出了一种自动根据数据中的样本数目推测计算过程中最大上位性阶数的方法，而不需要用户指定上位性阶数，算法考虑了小于或等于最大阶数的所有的上位性；
3.本发明提出了一种基于K2和G-test在高阶SNP组合上检测显著相关的SNP组合的方法，该方法使算法可以检测到与GWAS数据中的疾病显著相关的所有小于或等于最大上位性阶数的SNP组合；
4.针对K2、CE、Gini等函数评测不公的问题，本发明提出了一种对计算过程中的列联表进行合并的方案，可以有效地消除函数计算中的不公现象。
一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，包括以下步骤，且以下步骤顺次进行：
步骤S1：载入全基因组关联分析Genome-wide association study ，GWAS数据到内存
使用与Boolean Operation-based Screening and Testing，BOOST算法类似的基于二进制的存储方式；
步骤S2：根据数据中的样本数计算最大上位性阶数和列联表最大长度
步骤S3：初始化狼群
初始化狼群中的狼，狼的数目由用户通过参数指定，每一只狼是一个长度为mo的整型向量，向量上的每一个数字代表GWAS数据中某一个single-nucleotide polymorphism，SNP的索引下标，而每一只狼对应的就是一个SNP组合；
步骤S4：更新变异率
变异率的计算方式如式（3）所示，它是0到最大变异率mr之间的一个值；
（3）
其中，vr是变异率，ns是狼群中所有SNP去重之后的数目，代表了狼群的复杂性，nw是狼群中狼的数目，mr是算法可以接受的最大的变异率，mr变量的值由用户通过参数指定；
步骤S5：计算狼群中所有狼的K2、CE、Gini目标函数的值
其中，X代表一个SNP组合，Y代表疾病状态，k2(X,Y)是SNP组合X与疾病状态Y之间的K2值，XG是SNP组合的组合基因型集合，YG代表疾病状态的集合，对于GWAS研究，通常样本的状态只有患病与正常，因此YG通常为{0,1}，一般用0代表正常样本，1代表患病样本，mx是样本中组合基因型为x的样本的数目，mx,y是样本中组合基因型为x并且样本状态为y的样本数目；ce(X,Y)是X与Y之间的CE值，p(x,y)是SNP组合基因型为x并且样本状态为y的样本的数目与样本总数的比例，p(x)是SNP组合基因型为x的数目与样本总数的比例；gini(X,Y)是X与Y之间的Gini值，p(y|x)是在所有SNP组合基因型为x的样本中，样本状态为y的样本的占比；
步骤S6：根据K2、CE、Gini目标函数的值以及每匹狼的头狼次数选取三匹头狼，并更新头狼次数；
步骤S7：在三匹头狼上检测与疾病相关的SNP组合
在每一次狼群中选出三匹头狼之后，针对每一匹头狼，基于K2值反复移除头狼中的噪声SNP，直到不存在噪声，如果最后剩余的SNP的数目大于1，则算法找到了一个SNP组合，并且这个组合中的SNP彼此联系且与疾病状态相关；
步骤S8：狼群向三匹头狼移动
狼是一个长度为mo的整型向量，每一个整数代表了GWAS数据中一个SNP的下标，检测完三匹头狼之后，对于狼群中其它的狼，向三匹头狼移动，以达到狼群寻优的目的；
步骤S9：移动狼群中的每一匹头狼之外的狼，从而完成了一次狼群的移动，而后判断灰狼算法的迭代次数，如果迭代次数到达最大迭代次数，则算法执行完成，将算法记录的结果输出到结果文件，如果未达到最大迭代次数，更新变异率，进入下一次循环。
进一步，所述步骤S1中，在内存中使用两个三维数组分别存储患病样本与正常样本的基因型数据，对于每一个数组，第一维度代表数据中不同的SNP，第二维度代表对应的SNP不同的基因型，第三个维度代表对应SNP与其对应的基因型下有哪些样本。
进一步，所述步骤S8具体包括以下步骤：
对于一匹狼a和一匹头狼b，a向b移动后的狼为c，它们都是长度为mo的整型向量，在生成c中的每一个元素的时候，依变异率随机从GWAS数据的所有SNP中随机选择一个填入c中。
进一步，所述步骤S2中，最大上位性阶数的计算方法如式（1）所示：
(1)
其中，mo是最大上位性阶数，log是以自然常数为底的对数函数，m0是GWAS数据中的正常样本数，m1是GWAS中的患病样本数，min是取两个数最小值的函数；
列联表最大长度的计算如式（2）所示：
(2)
其中，ml是列联表最大长度，min、m0、m1的定义同公式1。进一步，所述步骤S6中，根据式（7）的计算结果进行头狼的选择：
（7）
其中，X是一个SNP组合，即狼群中的一匹狼，ls(X)是选择头狼时依据的值，在狼群中计算每一匹狼的ls值，选择值最小的狼作为头狼，lw(X)是X已经做过头狼的次数，在算法开始运行时，将所有SNP组合的头狼次数置0，当任意一匹狼X被选择为头狼时，对应的lw(X)的计数加一，f(X)是X的K2、CE、Gini值之一，在狼群初始化之后，或每一次向三匹头狼进行移动之后，针对每一个目标函数，都根据公式7在狼群中选择一匹头狼，选出三匹头狼之后，狼群参照三匹头狼的位置进行移动。
进一步，所述步骤S7中，根据公式8计算该SNP组合与疾病状态之间关系的显著性，如果其显著性超过用户指定的阈值，保存这个组合作为算法的输出结果之一；
其中，g(X,Y)是SNP组合X与样本状态Y进行G-test独立性检验的p-value值，F(X,Y)是独立性检验的自由度，mx是SNP组合基因型为x的样本的数目，my是样本状态为y的样本数目，m为样本总数，E(x,y)是当假设X与Y之间独立时，SNP组合基因型为x且样本状态为y的期望样本数，S(X,Y)是一个统计量，用于综合度量不同组合基因型下真实样本数与期望样本数之间的差距，根据独立向检验的知识，当X与Y独立时，S(X,Y)服从自由度为F(X,Y)的卡方分布，pvalueOfG(X,Y)是根据统计量S(X,Y)计算自由度为F(X,Y)的卡方分布下的pvalue值的函数，g(X,Y)的值越低，代表X与Y越不可能是互相独立的，其他符号的定义与前文相同。
本发明的一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法具有以下优点：
1.本发明提出的方法是一种基于多个目标函数实现的灰狼优化算法，以不同的角度度量与疾病状态相关的SNP组合，针对GWAS数据具有更好的检测能力；
2.本发明提出的方法不需要用户在运行算法时指定上位性阶数，可以检测到与疾病状态显著相关的任意的小于或等于最大上位性阶数的SNP组合，更加方便用户使用的同时，也更加地符合GWAS研究的实际情况；
3.本发明提出了一种在度量SNP组合与疾病状态相关性时，合并列联表的方法，解决了K2、CE、Gini等函数对列联表长度的偏性问题，从而使算法中的目标函数的计算更加的合理，提升了算法检测的能力。
附图说明
图1为本发明一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法的算法流程图；
图2为本发明实施例中的不同算法在DME 100数据集上的结果对比图；
图3为本发明实施例中的不同算法在DNME 100数据集上的结果对比图；
图4为本发明实施例中的不同算法在DNME3 100数据集上的结果对比图。
具体实施方式
为了更好地了解本发明的目的、结构及功能，下面结合附图，对本发明一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法做进一步详细的描述。
如图1所示，本发明提出了一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，命名为DE-GWO算法，其主要流程图如图 1所示，接下来将对流程图中的每个部分做详细地介绍。
一、载入GWAS数据到内存
DE-GWO算法基于Java编程语言实现，将GWAS数据载入内存时，使用与BOOST算法类似的基于二进制的存储方式，在内存中使用两个三维数组分别存储患病样本与正常样本的基因型数据，对于每一个数组，其第一维度代表数据中不同的SNP，第二维度代表对应的SNP不同的基因型，第三个维度代表对应SNP与其对应的基因型下有哪些样本。DE-GWO通过这种二进制的方式存储GWAS数据并提供了对应的二进制计算方式，大幅加快了算法中计算K2、CE、Gini等函数的速度。
二、根据数据中的样本数计算最大上位性阶数和列联表最大长度
最大上位性阶数（mo）用于指定DE-GWO算法可以检测到的与疾病相关的SNP组合的最大长度，它受到样本数的限制，因为在计算SNP组合与疾病相关性的时候，对于任何一种评价函数，如果SNP组合的长度过长，每一种SNP组合的基因型下只有非常少的样本，则会导致评价函数的失效。本发明使用的最大上位性阶数的计算方法如公式1所示，其目的是限制DE-GWO算法执行过程中基于各种评价函数计算的SNP组合的长度，保证在计算过程中每种组合基因型下的平均样本数为自然常数（e），避免因参与计算的SNP组合长度过长而导致的评价函数失效。
(1)
其中，mo是最大上位性阶数，log是以自然常数为底的对数函数，m0是GWAS数据中的正常样本数，m1是GWAS中的患病样本数，min是取两个数最小值的函数。
列联表最大长度（ml）是用来限制评价函数计算的过程中产生的列联表的最大长度的，在基于评价函数计算SNP组合与疾病的相关性的时候，即使对于同样长度的SNP组合，在生成列联表时，由于某些组合基因型下可能没有样本等现象，非零的列联表的长度并不统一，而几乎所有的评价函数在计算上对于不同长度的列联表是有偏性的，越长的非零列联表可能导致SNP组合与疾病越大的相关性，因此，在DE-GWO中计算评价函数时，引入了对列联表长度的限制方法，以公平地评价SNP组合与疾病的相关性。列联表最大长度的计算如公式2，其基本思想是使列联表中的每一个单元的平均样本数为10。
（2）
其中，ml是列联表最大长度，min、m0、m1的定义同公式1；注：“” 是向下取整，公式的完整意思是，取m0和m1的最小值，然后除以l0，得到的数，向下取整，得到ml必须是一个整数。
三、初始化狼群
初始化狼群中的狼，狼的数目（nw）由用户通过参数指定，默认值是40，每一只狼是一个长度为mo的整型向量，向量上的每一个数字代表GWAS数据中某一个SNP的索引下标，而每一只狼对应的就是一个SNP组合，灰狼优化算法通过使群体不断地进化，优化狼群，检测与疾病相关的SNP组合。
四、更新变异率
变异率（vr）用于在DE-GWO算法的执行过程中度量狼群中个体的复杂性。随着灰狼优化算法求优过程的进行，狼群会陷入到局部的最优解上，导致狼群中的狼彼此之间的相似度极高，为了解决这一问题，DE-GWO中引入了变异率，用以基于狼群个体之间的相似性控制狼群的进化方式，从而达到使狼群跳出局部最优解。变异率的计算方式如公式3所示，它是0到mr（最大变异率）之间的一个值。
（3）
其中，vr是变异率，ns是狼群中所有SNP去重之后的数目，代表了狼群的复杂性，nw是狼群中狼的数目，mr是算法可以接受的最大的变异率，该变量的值由用户通过参数指定。
五、计算狼群中所有狼三个目标函数的值
本发明提出的DE-GWO算法使用三个度量SNP组合和疾病之间相关性的函数作为群体进化的目标，这三个函数分别是K2、CE、Gini，它们被广泛地应用于GWAS数据的分析之中，具体的计算方法如公式4-6所示，这三个函数从不同的角度度量一个SNP组合与疾病之间的相关性，其值越低，对应SNP组合与疾病之间的关系越强。
在上面的公式中，X代表一个SNP组合，Y代表疾病状态，k2(X,Y)是SNP组合X与疾病状态Y之间的K2值，XG是SNP组合的组合基因型集合，如对于一个二阶SNP组合，通常情况下，每个SNP的基因型可以是0、1或2，此时XG是{(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2)}，YG代表疾病状态的集合，对于GWAS研究，通常样本的状态只有患病与正常，因此YG通常为{0,1}，一般用0代表正常样本，1代表患病样本，mx是样本中组合基因型为x的样本的数目，mx,y是样本中组合基因型为x并且样本状态为y的样本数目。ce(X,Y)是X与Y之间的CE值，p(x,y)是SNP组合基因型为x并且样本状态为y的样本的数目与样本总数的比例，p(x)是SNP组合基因型为x的数目与样本总数的比例。gini(X,Y)是X与Y之间的Gini值，p(y|x)是在所有SNP组合基因型为x的样本中，样本状态为y的样本的占比。
虽然以上的三个函数被广泛地应用于度量SNP组合与疾病之间的关系，然而它们都存在一个普遍的问题，其值对于SNP组合基因型数目有很大的偏性，即在计算的过程中，SNP组合基因型将样本分的份数越多，上面的三个函数的值越小，三个函数计算的结果表明SNP组合与疾病的关系越大，这种现象为灰狼算法的优化过程带来了巨大的阻碍。本发明提出了一种基于列联表最大长度（ml）的通用的函数优化方案，分析公式4-6，可以看出每一个函数都是通过基于XG的累加获得最后的函数值的（其中，K2的函数可以通过log的形式将累乘转为累加），而三个函数都是值越小越相关，因此在函数的计算过程中，DE-GWO列联表上对SNP组合基因型进行排序，只保留ml-1个值最小的单元，其余的单元在列联表上进行样本的合并，从而保证列联表的最大长度始终小于或等于ml，从而解决了三个函数在计算上对于SNP组合基因型数目上的偏性。
六、根据三个目标函数的值以及每匹狼的头狼次数选取三匹头狼，并更新头狼次数。
在DE-GWO算法的运行过程中，维护了三匹头狼用于指引算法优化的方向，三匹头狼分别是狼群中三个目标函数值最小的狼，并保证三匹头狼是当前狼群中不同的狼，但与传统的灰狼算法不同，DE-GWO是一个检测SNP组合的方法，在GWAS数据中，理论上同时存在多个与疾病相关的SNP组合，因此，DE-GWO必须具有较强地跳出局部极小点的能力，不能被某一匹很强的头狼一直引导，因此需要在选择头狼时综合考虑目标函数的值与头狼被选择的次数。综合考虑，在选择头狼时，DE-GWO根据公式7的计算结果进行选择。
（7）
其中，X是一个SNP组合，即狼群中的一匹狼，ls(X)是选择头狼时依据的值，在狼群中计算每一匹狼的ls值，选择值最小的狼作为头狼，lw(X)是X已经做过头狼的次数，在算法开始运行时，将所有SNP组合的头狼次数置0，当任意一匹狼X被选择为头狼时，对应的lw(X)的计数加一，f(X)是X的K2、CE、Gini值之一，在狼群初始化之后，或每一次向三匹头狼进行移动之后，针对每一个目标函数，都会根据公式7在狼群中选择一匹头狼，选出三匹头狼之后，狼群参照三匹头狼的位置进行移动，以求在灰狼算法的求优过程中从不同的角度充分考虑狼群中蕴含的信息。
七、在三匹头狼上检测与疾病相关的SNP组合
DE-GWO算法使用了一种基于K2函数检测与疾病相关的SNP组合的方法，其基本思想是针对一个SNP组合，已经根据公式4计算出了它与疾病状态之间的K2值k，当将一个SNPx从这个组合中移除时，根据公式4计算剩余的SNP与疾病状态之间的K2值kx，如果x是一个与疾病状态相关的SNP，或者x是一个与组合内其他SNP共同影响疾病状态的SNP（上位性），则kx>k，反之，如果x是一个噪声SNP，则kx<k。
基于上述思想，在DE-GWO算法求优的过程中，在每一次狼群中选出三匹头狼之后，针对每一匹头狼，基于K2值反复移除头狼中的噪声SNP，直到不存在噪声，而如果最后剩余的SNP的数目大于1，则算法找到了一个SNP组合，并且这个组合中的SNP彼此联系且与疾病状态相关，最后根据公式8计算该SNP组合与疾病状态之间关系的显著性，如果其显著性超过用户指定的阈值（默认值0.05），DE-GWO保存这个组合作为算法的输出结果之一。
公式8中，g(X,Y)是SNP组合X与样本状态Y进行G-test独立性检验的p-value值，F(X,Y)是独立性检验的自由度，mx是SNP组合基因型为x的样本的数目，my是样本状态为y的样本数目，m为样本总数，E(x,y)是SNP组合基因型为x且样本状态为y的期望样本数，pvalueOfG是根据统计量计算其在卡方分布下的pvalue值的函数，其他符号的定义与前文相同。
八、狼群向三匹头狼移动
在DE-GWO算法检测完三匹头狼之后，对于狼群中其它的狼，需要向三匹头狼移动，以达到狼群寻优的目的，与传统的灰狼算法不同，在DE-GWO中的狼是一个长度为mo的整型向量，每一个整数代表了GWAS数据中一个SNP的下标，基于这些向量的基础数学运算无法达到狼向头狼移动的目的，为此，DE-GWO算法提出了一种新的狼移动的算法，在考虑GWAS数据特殊性的前提下，完成一匹狼向三匹头狼移动的过程。
具体的说，对于一匹狼a和一匹头狼b，a向b移动后的狼为c，它们都是长度为mo的整型向量，c向量中一半的数字来自于a，另一半的数字来自于b，通过这样的方式完成a向b的移动，同时，考虑到狼群中的个体的复杂性，避免狼群陷在局部极小点，在生成c中的每一个元素的时候，依变异率（vr）随机从GWAS数据的所有SNP中随机选择一个填入c中，通过这样的方式DE-GWO算法完成了一匹狼向一匹头狼的移动。在DE-GWO算法中存在三匹头狼，因此在一匹狼的移动过程中，会同时存在三个向量c1、c2、c3，它们分别代表着一匹狼向三匹头狼移动的结果，为了充分在狼群的移动过程中考虑三匹头狼的信息，DE-GWO算法在c1、c2、c3中分别抽取三分之一的数字，组成了新的移动后的狼，以上就是狼群中一匹狼的移动过程。
依照上述方法移动狼群中的每一匹头狼之外的狼，从而完成了一次狼群的移动，而后判断灰狼算法的迭代次数，如果迭代次数到达最大迭代次数（用户指定的参数），则算法执行完成，将算法记录的结果输出到结果文件，如果未达到最大迭代次数，更新变异率，进入下一次循环。
具体实施例
1. 模拟数据集
为了验证本发明提出的DE-GWO算法的检测能力，本发明在共22个模拟数据集上进行了与其他检测算法的对比实验，这些模拟数据集都是使用GAMETES_2.1软件基于不同的致病模型生成的，其简介如表 1所示。
表 1 模拟数据集介绍
DME 100数据集：包含8个DME数据集，这8个数据集由8个文件夹构成，每个文件夹对应一个不同的DME（disease loci with marginal effects)模型，其下有100个文件，这些文件是由该文件夹对应的一个DME模型生成的，每个文件的内容是一个模拟的GWAS数据，其中包含1600个样本（800个患病，800个正常），100个SNP（2个致病，98个无关）。在该模拟数据集上的实验可以验证算法检测具有边际效应的与疾病相关的SNP组合的能力。
DNME 100数据集：包含8个DNME数据集，这8个数据集由8个文件夹构成，每个文件夹对应一个不同的DNME（disease loci without marginal effects)模型，其他内容与DME100数据集类似。DNME数据集中模拟的致病模型相对于DME数据集缺少了边际效应，因此在检测致病SNP组合时，更加的困难。在该模拟数据集上的实验可以验证算法检测无边际效应的与疾病相关的SNP组合的能力。
DNME3 100数据集：包含6个DNME数据集，这6个数据集由6个文件夹构成，每个文件夹对应一个不同的DNME模型（3阶），其内容与DNME 100数据集类似，不同之处在于DNME3100模型模拟的是3阶致病SNP组合模型，因此，数据集中的每个GWAS模拟数据文件包含100个SNP（3个致病，97个无关）。在该模拟数据集上的实验可以进一步验证算法在检测高阶无边际效应的与疾病相关的SNP组合上的能力。
2. 评价指标
为了比较不同算法在模拟数据集上检测与疾病相关的SNP组合的能力，本发明使用F-measure和power作为评价指标用于衡量算法的检测能力，它们都广泛应用于该领域中，其计算方法如公式9所示。
其中power和F-measure是用于评价算法检测能力的函数，针对数据集中的一个模型，有100个模拟数据文件，算法在模型上的power是使用该算法分析这个模型对应的100个文件之后，正确检测出文件中的致病模型的比例。而针对每一个模拟文件，当一个算法执行完毕之后，会检测到若干个SNP组合，算法认为它们是与疾病显著相关的SNP组合，此时，TP是算法检测到的致病SNP组合的数目，FN是算法未能检测出的致病SNP组合的数目，FP是算法检测出的SNP组合中非致病的组合数目，recall是召回率，precision是精确率，而针对数据集中的每一个模型，算法执行完毕后会获得100个F-measure，对这些值取平均，作为算法在模型上的F-measure。
3. 实验结果
为了验证本发明提出的DE-GWO算法的检测能力，本发明分别使用AntEpiSeeker、DECMDR、SNPHarvester、HSMMGKG、SEE、DE-GWO算法基于表 2中的参数分析了大量的模拟数据集，其中AntEpiSeeker和SNPHarvester无法用于检测3阶SNP组合。
如表 2所示，为了公平地比较不同算法在模拟数据集上检测SNP组合的能力，实验中的参数设置遵循各算法文献中的最优建议，并且，由于这些算法多数是基于群智能算法检测致病SNP组合的，因此不同的算法在相同的模拟数据集上进行分析时，群智能算法中的群体规模和迭代次数都是一样的。
表 2 算法在模拟数据集上的执行参数
在DME 100上的实验结果如表 3所示，其中p和f分别代表power和F-measure，其计算方法如公式9所示，本发明使用不同的算法分析了该数据集上的8个模拟数据集，这8个数据集由8个不同的模拟的致病模型构建，在大多数模型上，DE-GWO算法的检测能力都有明显的优势。平均上说，DE-GWO算法在DME 100数据集上的power和F-measure分别是0.93与0.80，而在其他算法中，最高的平均power和F-measure分别是0.75与0.65。
实验结果表明，DE-GWO算法在检测具有边际效应的致病SNP组合时，相对于当前的其他算法有明显的优势。图 2展示了不同算法在DME 100数据上的结果对比图。
表 3 不同算法在DME 100数据集上的对比结果
在DNME 100数据集上的执行结果如表 4所示，在该数据集上的8个不同的模拟致病模型上，DE-GWO算法仍然有很大的优势。平均上说，DE-GWO算法在该数据集上的power和F-measure分别是0.95与0.84，而在其他算法中，最高的平均power和F-measure分别是0.90与0.75，虽然AntEpiSeeker的power是0.90，只是稍低于DE-GWO算法，但其F-measure仅为0.11，与DE-GWO算法有很大的差距，可见AntEpiSeeker算法虽然可以有效地检测出致病的SNP组合，但相对于DE-GWO算法，其假阳性过高，检测结果中存在大量的噪声。
实验结果表明，DE-GWO算法在检测无边际效应的致病SNP组合时，相对于当前的其他算法仍然有明显的优势。图 3展示了不同算法在DNME 100数据上的结果对比图。
表 4 不同算法在DNME 100数据集上的对比结果
在DNME3 100数据集上的执行结果如表 5所示，在该数据集上的6个不同的模拟致病模型上，DE-GWO算法具有更大的优势。平均上说，DE-GWO算法在该数据集上的power和F-measure分别是0.96与0.76，而在其他算法中，最高的平均power和F-measure分别是0.02与0.02，而这种明显的优势主要是由于当前的其他算法在设计时，并没有充分考虑对于高阶致病SNP组合的检测，这些算法缺少这方面的能力，而复杂疾病的往往是由多个基因共同作用引起的常见疾病，在GWAS数据上蕴含的致病模型也可能非常的复杂，DE-GWO算法更加地适合应用于高阶致病SNP组合的检测任务。
实验结果表明，DE-GWO算法在检测高阶致病SNP组合时，相对于当前的其他算法有明显的优势。图 4展示了不同算法在DNME3 100数据上的结果对比图。
表 5 不同算法在DNME3 100数据集上的对比结果
可以理解，本发明是通过一些实施例进行描述的，本领域技术人员知悉的，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行各种改变或等效替换。另外，在本发明的教导下，可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本申请的权利要求范围内的实施例都属于本发明所保护的范围内。
Claims (6)
1.一种基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，其特征在于，包括以下步骤，且以下步骤顺次进行：
步骤S1：载入全基因组关联分析Genome-wide association study ，GWAS数据到内存
使用与Boolean Operation-based Screening and Testing，BOOST算法类似的基于二进制的存储方式；
步骤S2：根据数据中的样本数计算最大上位性阶数和列联表最大长度
步骤S3：初始化狼群
初始化狼群中的狼，狼的数目由用户通过参数指定，每一只狼是一个长度为mo的整型向量，向量上的每一个数字代表GWAS数据中某一个single-nucleotide polymorphism ，SNP的索引下标，而每一只狼对应的就是一个SNP组合；
步骤S4：更新变异率
变异率的计算方式如式（3）所示，它是0到最大变异率mr之间的一个值；
（3）
其中，vr是变异率，ns是狼群中所有SNP去重之后的数目，代表了狼群的复杂性，nw是狼群中狼的数目，mr是算法可以接受的最大的变异率，mr变量的值由用户通过参数指定；
步骤S5：计算狼群中所有狼的K2、CE、Gini目标函数的值
（4）
（5）
（6）
其中，X代表一个SNP组合，Y代表疾病状态，k2(X,Y)是SNP组合X与疾病状态Y之间的K2值，XG是SNP组合的组合基因型集合，YG代表疾病状态的集合，对于GWAS研究，通常样本的状态只有患病与正常，因此YG通常为{0,1}，一般用0代表正常样本，1代表患病样本，mx是样本中组合基因型为x的样本的数目，mx,y是样本中组合基因型为x并且样本状态为y的样本数目；ce(X,Y)是X与Y之间的CE值，p(x,y)是SNP组合基因型为x并且样本状态为y的样本的数目与样本总数的比例，p(x)是SNP组合基因型为x的数目与样本总数的比例；gini(X,Y)是X与Y之间的Gini值，p(y|x)是在所有SNP组合基因型为x的样本中，样本状态为y的样本的占比；
步骤S6：根据K2、CE、Gini目标函数的值以及每匹狼的头狼次数选取三匹头狼，并更新头狼次数；
步骤S7：在三匹头狼上检测与疾病相关的SNP组合
在每一次狼群中选出三匹头狼之后，针对每一匹头狼，基于K2值反复移除头狼中的噪声SNP，直到不存在噪声，如果最后剩余的SNP的数目大于1，则算法找到了一个SNP组合，并且这个组合中的SNP彼此联系且与疾病状态相关；
步骤S8：狼群向三匹头狼移动
狼是一个长度为mo的整型向量，每一个整数代表了GWAS数据中一个SNP的下标，检测完三匹头狼之后，对于狼群中其它的狼，向三匹头狼移动，以达到狼群寻优的目的；
步骤S9：移动狼群中的每一匹头狼之外的狼，从而完成了一次狼群的移动，而后判断灰狼算法的迭代次数，如果迭代次数到达最大迭代次数，则算法执行完成，将算法记录的结果输出到结果文件，如果未达到最大迭代次数，更新变异率，进入下一次循环。
2.根据权利要求1所述的基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，其特征在于，所述步骤S1中，在内存中使用两个三维数组分别存储患病样本与正常样本的基因型数据，对于每一个数组，第一维度代表数据中不同的SNP，第二维度代表对应的SNP不同的基因型，第三个维度代表对应SNP与其对应的基因型下有哪些样本。
3.根据权利要求1所述的基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，其特征在于，所述步骤S8具体包括以下步骤：
对于一匹狼a和一匹头狼b，a向b移动后的狼为c，它们都是长度为mo的整型向量，在生成c中的每一个元素的时候，依变异率随机从GWAS数据的所有SNP中随机选择一个填入c中。
4.根据权利要求1所述的基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，其特征在于，所述步骤S2中，最大上位性阶数的计算方法如式（1）所示：
(1)
其中，mo是最大上位性阶数，log是以自然常数为底的对数函数，m0是GWAS数据中的正常样本数，m1是GWAS中的患病样本数，min是取两个数最小值的函数；
列联表最大长度的计算如式（2）所示：
(2)
其中，ml是列联表最大长度。
5.根据权利要求1所述的基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，其特征在于，所述步骤S6中，根据式（7）的计算结果进行头狼的选择：
(7)
其中，X是一个SNP组合，即狼群中的一匹狼，ls(X)是选择头狼时依据的值，在狼群中计算每一匹狼的ls值，选择值最小的狼作为头狼，lw(X)是X已经做过头狼的次数，在算法开始运行时，将所有SNP组合的头狼次数置0，当任意一匹狼X被选择为头狼时，对应的lw(X)的计数加一，f(X)是X的K2、CE、Gini值之一，在狼群初始化之后，或每一次向三匹头狼进行移动之后，针对每一个目标函数，都根据公式7在狼群中选择一匹头狼，选出三匹头狼之后，狼群参照三匹头狼的位置进行移动。
6.根据权利要求1所述的基于GWO算法在GWAS数据上检测与复杂疾病相关SNP组合的方法，其特征在于，所述步骤S7中，根据公式8计算该SNP组合与疾病状态之间关系的显著性，如果其显著性超过用户指定的阈值，保存这个组合作为算法的输出结果之一；
(8)
其中，g(X,Y)是SNP组合X与样本状态Y进行G-test独立性检验的p-value值，F(X,Y)是独立性检验的自由度，mx是SNP组合基因型为x的样本的数目，my是样本状态为y的样本数目，m为样本总数，E(x,y)是当假设X与Y之间独立时，SNP组合基因型为x且样本状态为y的期望样本数，S(X,Y)是一个统计量，用于综合度量不同组合基因型下真实样本数与期望样本数之间的差距，根据独立向检验的知识，当X与Y独立时，S(X,Y)服从自由度为F(X,Y)的卡方分布，pvalueOfG(X,Y)是根据统计量S(X,Y)计算自由度为F(X,Y)的卡方分布下的pvalue值的函数，g(X,Y)的值越低，代表X与Y越不可能是互相独立的。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410118311.1A CN117649876B (zh) | 2024-01-29 | 2024-01-29 | 基于gwo算法在gwas数据上检测与复杂疾病相关snp组合的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410118311.1A CN117649876B (zh) | 2024-01-29 | 2024-01-29 | 基于gwo算法在gwas数据上检测与复杂疾病相关snp组合的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117649876A true CN117649876A (zh) | 2024-03-05 |
| CN117649876B CN117649876B (zh) | 2024-04-12 |
Family
ID=90045425
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410118311.1A Active CN117649876B (zh) | 2024-01-29 | 2024-01-29 | 基于gwo算法在gwas数据上检测与复杂疾病相关snp组合的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117649876B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107832830A (zh) * | 2017-11-17 | 2018-03-23 | 湖北工业大学 | 基于改进型灰狼优化算法的入侵检测系统特征选择方法 |
| CN109390032A (zh) * | 2018-11-02 | 2019-02-26 | 吉林大学 | 一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的snp组合的方法 |
| CN113241122A (zh) * | 2021-06-11 | 2021-08-10 | 长春工业大学 | 自适应弹性网与深度神经网络融合的基因数据变量选择及分类方法 |
| WO2022141213A1 (zh) * | 2020-12-30 | 2022-07-07 | 中南大学 | 一种智慧城市智轨车辆故障基因预测方法及系统 |
-
2024
- 2024-01-29 CN CN202410118311.1A patent/CN117649876B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107832830A (zh) * | 2017-11-17 | 2018-03-23 | 湖北工业大学 | 基于改进型灰狼优化算法的入侵检测系统特征选择方法 |
| CN109390032A (zh) * | 2018-11-02 | 2019-02-26 | 吉林大学 | 一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的snp组合的方法 |
| WO2022141213A1 (zh) * | 2020-12-30 | 2022-07-07 | 中南大学 | 一种智慧城市智轨车辆故障基因预测方法及系统 |
| CN113241122A (zh) * | 2021-06-11 | 2021-08-10 | 长春工业大学 | 自适应弹性网与深度神经网络融合的基因数据变量选择及分类方法 |
Non-Patent Citations (2)
| Title |
|---|
| ZHILIN DONG ETC.: "Time-Shift Multi-ScaleWeighted Permutation Entropy and GWO-SVM Based Fault Diagnosis Approach for Rolling Bearing", ENTROPY, 25 June 2019 (2019-06-25) * |
| 刘二辉;姚锡凡;刘敏;金鸿;: "基于改进灰狼优化算法的自动导引小车路径规划及其实现原型平台", 计算机集成制造系统, no. 11, 15 November 2018 (2018-11-15) * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117649876B (zh) | 2024-04-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wang et al. | Deep learning for plant genomics and crop improvement | |
| Chen et al. | Random forests for genomic data analysis | |
| Schliep et al. | Using hidden Markov models to analyze gene expression time course data | |
| Sinha | On counting position weight matrix matches in a sequence, with application to discriminative motif finding | |
| CN112466404B (zh) | 一种宏基因组重叠群无监督聚类方法及系统 | |
| CN107679367B (zh) | 一种基于网络节点关联度的共调控网络功能模块识别方法及系统 | |
| Linder et al. | Predicting RNA-seq coverage from DNA sequence as a unifying model of gene regulation | |
| CN109390032B (zh) | 一种基于进化算法在全基因组关联分析的数据中探索与疾病相关的snp组合的方法 | |
| Whitehouse et al. | Timesweeper: accurately identifying selective sweeps using population genomic time series | |
| Wang et al. | An efficient gene bigdata analysis using machine learning algorithms | |
| CN117649876B (zh) | 基于gwo算法在gwas数据上检测与复杂疾病相关snp组合的方法 | |
| Orzechowski et al. | Propagation-based biclustering algorithm for extracting inclusion-maximal motifs | |
| Kwarciak et al. | Tabu search algorithm for DNA sequencing by hybridization with multiplicity information available | |
| Robinson et al. | Improving computational predictions of cis-regulatory binding sites | |
| CN108897990B (zh) | 面向大规模高维序列数据的交互特征并行选择方法 | |
| CN115769300A (zh) | 变体致病性评分和分类及其用途 | |
| Vukusic et al. | Applying genetic programming to the prediction of alternative mRNA splice variants | |
| Sitarčík et al. | epiBAT: Multi-objective bat algorithm for detection of epistatic interactions | |
| KR101701168B1 (ko) | 유전자 패스웨이 활성지수의 세부적 정량화를 위한 유전자 프로파일 방법 | |
| Bonomo et al. | Prediction of Disease–lncRNA Associations via Machine Learning and Big Data Approaches | |
| Tan et al. | DITGOssi: a two-stage invasive tumor growth optimization algorithm for the detection of SNPSNP interactions | |
| CN113571134B (zh) | 基于骨干粒子群算法的基因数据特征选择方法及装置 | |
| Mezlini | Machine Learning Methodologies for Uncovering the Mechanisms of Complex Diseases | |
| Acar et al. | Adaptive updates for MAP configurations with applications to bioinformatics | |
| Anderson et al. | A path integral approach for allele frequency dynamics under polygenic selection |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |