CN110570909A - 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法 - Google Patents

一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法 Download PDF

Info

Publication number
CN110570909A
CN110570909A CN201910859461.7A CN201910859461A CN110570909A CN 110570909 A CN110570909 A CN 110570909A CN 201910859461 A CN201910859461 A CN 201910859461A CN 110570909 A CN110570909 A CN 110570909A
Authority
CN
China
Prior art keywords
network structure
snp
bee
node
network
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
Application number
CN201910859461.7A
Other languages
English (en)
Other versions
CN110570909B (zh
Inventor
刘建晓
杨晨
高辉
杨轩
黄穗宇
阚裕隆
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong Agricultural University
Original Assignee
Huazhong Agricultural University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Huazhong Agricultural University filed Critical Huazhong Agricultural University
Priority to CN201910859461.7A priority Critical patent/CN110570909B/zh
Publication of CN110570909A publication Critical patent/CN110570909A/zh
Application granted granted Critical
Publication of CN110570909B publication Critical patent/CN110570909B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/29—Graphical models, e.g. Bayesian networks
    • G06F18/295—Markov models or related models, e.g. semi-Markov models; Markov random fields; Networks embedding Markov models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Biophysics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Medical Informatics (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Biotechnology (AREA)
  • Molecular Biology (AREA)
  • Mathematical Physics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Computing Systems (AREA)
  • Analytical Chemistry (AREA)
  • Computational Linguistics (AREA)
  • Biomedical Technology (AREA)
  • Genetics & Genomics (AREA)
  • Bioethics (AREA)
  • Chemical & Material Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Epidemiology (AREA)
  • Public Health (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明涉及生物信息技术领域，提供了一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，包括S1‑S4四个步骤。本发明的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，首先利用扩张、收缩、对称性检测三个阶段，通过条件互信息计算得到节点的马尔可夫毯，构建初始蜜源网络结构，然后，基于初始蜜源通过随机加边、减边、逆转边3种操作生成新的蜜源，直到达到最大初始蜜源数目，利用人工蜂群算法的三种操作(采蜜蜂、观察蜂和侦察蜂)与贝叶斯网络的BIC与MIT打分方法，对贝叶斯网络结构进行演化，找到最优的网络结构，快速准确的获取到影响表型性状的上位性基因位点，辅助基因功能挖掘。

Description

一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法
技术领域
本发明涉及生物信息技术领域，具体为一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法。
背景技术
随着人们生活质量的提高、医疗卫生环境的改善和相关生物技术的迅速发展，人类高发病率的疾病类型发生了巨大的改变。营养不良、传染病等主要受环境影响的疾病得到了有效地控制，而复杂疾病和孟德尔遗传病成为了目前困扰人类的主要疾病，孟德尔遗传病是一种单基因疾病，其遗传过程遵循孟德尔遗传定律，目前研究者利用定位克隆的方法确定了相关遗传基因，基本阐明了其遗传方式。复杂疾病是由多个基因，环境因素共同作用的结果，其病因非常复杂。其中复杂疾病占人类疾病的80％以上，主要包括癌症、2型糖尿病、高血压、心脏病等，因此，探索复杂疾病其中的遗传机制具有重要的研究意义。
从生物遗传学的角度看，决定生物复杂性状的遗传因素主要包括三个方面：基因主效应、基因与基因之间的相互作用和基因与环境之间的相互作用。全基因组关联分析(Genome-wide Association Studies,GWAS)是用人类全基因组范围内数以百万计的单核苷酸多态性(Single Nucleotide Polymorphism,SNP)作为遗传标记，通过病例/对照分析，来进一步探索出复杂疾病的相关遗传机制。但是早期开展的全基因组关联研究，主要侧重于主效基因的检测上，忽略了基因与基因之间的相互作用，即上位性(Epistasis)。而对于复杂疾病来说，SNP之间的上位性作用的影响往往更为重要，因此在全基因组范围内进行上位性检测对于探寻复杂疾病的致病病因有着重要的意义。
近年来，研究人员提出了多种上位性检测方法，但是这些方法仍然存在计算困难、效率低下、算法复杂度高以及假阳性率高等问题，导致不能准确高效地检测出与疾病相关联的SNP位点及其组合。
因此，在全基因组范围内提出更有效、更准确的上位性检测算法具有十分重要的研究意义，也对复杂疾病致病机理的发现、诊断、治疗和预防有着非常重要的作用。
发明内容
本发明的目的在于提供一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，通过人工蜂群算法优化贝叶斯网络构建包括SNP位点和表型性状的网络结构，进而挖掘出上位性位点。
为实现上述目的，本发明实施例提供如下技术方案：一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，包括如下步骤：
S1，利用扩张、收缩、对称性检测三个阶段，将目标节点分为表型Class节点和非Class节点两种情况，通过条件互信息计算获取不同节点的马尔可夫毯；
S2，根据得到不同节点的马尔可夫毯，构建包括SNP位点和表型Class的网络，得到初始蜜源，即初始网络结构；
S3，在不生成环的前提下，基于初始蜜源网络个体，通过随机增加边、删除边、逆转边的操作生成新的网络个体，直到网络个体数量达到最大初始蜜源数目；
S4，在步骤3所得初始蜜源的基础上，通过人工蜂群优化算法的三种蜜蜂行为与贝叶斯网络的BIC和MIT两种打分函数，对包括SNP位点和表型Class的贝叶斯网络进行演化，得到较优的贝叶斯网络结构，从而准确快速地获取到影响表型性状的上位性基因位点，其中三种蜜蜂分别为采蜜蜂、观察蜂和侦察蜂。
进一步，在所述S1步骤前，将基因型数据表示为0，1，2形式的数据，对于特定SNP基因型数据，例如AT，进行如下表示：AA用0表示，TT用2表示，AT/TA用1表示，其中0表示纯合子常见基因型，1表示杂合子，2表示纯合子少见基因型，表型Class取值为0或1，0表示为未患病，1表示为患病。
进一步，在所述S1步骤中，扩张阶段具体为：在构建非Class节点的马尔可夫毯时，假设表型Class节点已经加入到目标SNP节点的马尔可夫毯中，计算任意一个SNP节点与Class在目标SNP节点条件下的条件互信息，并将所述条件互信息小于阈值的SNP位点加入目标节点的马尔可夫毯中；在构建Class节点的马尔可夫毯时，假设任意一个SNP节点已经加入表型Class节点的马尔可夫毯中，计算另外任意一个SNP节点与已加入马尔可夫毯的SNP节点在Class条件下的条件互信息，并将小于阈值的SNP位点加入Class节点的马尔可夫毯中。
进一步，在所述S1步骤中，收缩阶段具体为：去除多余节点；对于任意目标节点SNP1或者Class节点，计算在去除SNP2的马尔可夫毯的条件下目标节点SNP1与SNP2的条件互信息，若所述条件互信息大于阈值，则删除SNP2并更新节点SNP1的马尔可夫毯。
进一步，在所述S1步骤中，对称性检测具体为：由于本发明中的方法是基于条件独立性检测，因此任意目标节点马尔可夫毯中的节点之间应该相互依赖，假设SNP1存在于SNP2的马尔可夫毯中，SNP1与SNP2应该相互依赖，与之对应，SNP2应该存在于SNP1的马尔可夫毯中，假如两个节点的马尔可夫毯并不对称，那么分别在两个节点的马尔可夫毯中删除另外一节点。
进一步，在所述S4步骤中，采蜜蜂阶段具体为：
S40，首先采蜜蜂对初始蜜源的邻域进行搜索，即通过加边、减边、逆转边三种操作对网络结构进行更新；
S41，然后利用BIC和MIT两种评分函数同时对新的网络结构进行打分，并与初始网络结构的分数进行比较，若新的网络结构优于初始网络结构，则将该网络结构进行保存，并记录该网络结构提高的分数；
S42，最后，根据提高的分数使用轮盘赌的方式随机选择新的网络结构，此时采蜜蜂对其进行采蜜；
以上的网络结构即蜜源。
进一步，在所述采蜜蜂阶段的方法中，使网络分数提高越大的操作被选择的概率越大，这不仅可以保证网络演变的正确性，又不至于陷入局部最优，假如轮盘赌方法选择放弃本次加边、减边或者逆转边操作，则说明该蜜源邻域已无更优网络结构，表示该网络结构已是最优的网络结构，此时采蜜蜂对其进行采蜜。
进一步，在所述S4步骤中，观察蜂阶段具体为：观察蜂阶段和采蜜蜂阶段交替进行，每当采蜜蜂搜索一次邻域得到较优的网络结构之后，观察蜂便对采蜜蜂所在的蜜源进行考察，根据两种打分函数计算对应网络结构的适应度值，然后利用轮盘赌的方式随机选择一个蜜源网络结构，适应度越大的网络结构被认为邻域拥有更优网络结构的可能性越大，被选到的几率也就越大，采蜜蜂对选中网络结构的邻域进行搜索，重复采蜜蜂阶段的操作。
进一步，在所述S4步骤中，侦察蜂阶段具体为：侦察蜂在采蜜蜂和观察蜂阶段中扮演监督员的角色，记录每个蜜源网络结构的采蜜次数以及最优蜜源网络结构，当某一蜜源网络结构达到最大采蜜次数，则放弃该蜜源，并通过步骤1重新生成新的蜜源网络结构，假如每只采蜜蜂均将对应蜜源开采至最大采蜜次数，则找到最优蜜源，即最优网络结构，此时算法结束。
进一步，适应度函数是用于评判网络结构好坏的标准，因此采用BIC与MIT联合打分的方法对网络的优劣进行判断，同时，由于BIC打分函数可以处理的父节点个数有限，因此采用可分解BIC打分方法处理较多个数的父节点，进而可以处理规模较大的网络结构。
与现有技术相比，本发明的有益效果是：一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，首先利用扩张、收缩、对称性检测三个阶段，通过条件互信息计算得到节点的马尔可夫毯，构建初始蜜源网络结构，然后，基于初始蜜源通过随机加边、减边、逆转边3种操作生成新的蜜源，直到达到最大初始蜜源数目，利用人工蜂群算法的三种操作(采蜜蜂、观察蜂和侦察蜂)与贝叶斯网络的BIC与MIT打分方法，对贝叶斯网络结构进行演化，找到最优的网络结构，快速准确的获取到影响表型性状的上位性基因位点，辅助基因功能挖掘。
附图说明
图1为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的贝叶斯网络矩阵编码表示图；
图2为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的贝叶斯网络三种操作示意图；
图3为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的具体实施的流程图；
图4为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的不同方法2位点上位性检测的准确率比较；
图5为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的不同方法2位点上位性检测的F1-score比较；
图6为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的不同方法2位点上位性检测的FPR比较；
图7为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的不同方法3位点上位性检测的准确率比较；
图8为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的不同方法3位点上位性检测的F1-score比较；
图9为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的不同方法3位点上位性检测的FPR比较；
图10为本发明实施例提供的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法的可分解BIC打分方法表示图；
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
本发明实施例提供一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，具体为：
1、用0，1，2形式的数据表示基因型数据，如对SNP基因型为AT的数据表示如下：AA用0表示，TT用2表示，AT/TA用1表示。0表示纯合子常见基因型，1表示杂合子，2表示纯合子少见基因型。图10所示为4个SNP(SNP1～SNP4)对应10个样本的基因型数据，最后一列Class表示表型性状，其中Class＝1表示case(患病)，Class＝0表示control(对照)。
2、将目标节点分为表型Class节点和非Class节点两种情况，利用扩张、收缩、对称性检测三个阶段，通过条件互信息计算获取不同节点的马尔可夫毯。
(1)扩张阶段：分为目标节点为非Class节点与目标节点为Class节点两种情况进行处理。构建表型Class节点的马尔可夫毯时，在初始阶段假设任意一个节点SNP1已经加入Class节点的马尔可夫毯中，通过Eq.(1)计算另外任意一个节点SNP2与已加入马尔可夫毯的SNP1节点在Class条件下的条件互信息I(SNP1,SNP2|Class)。构建非Class节点的马尔可夫毯时，在初始阶段假设表型Class节点已经加入到目标SNP1节点的马尔可夫毯中，通过Eq.(1)计算任意一个节点SNP2与Class在目标节点SNP1条件下的条件互信息I(SNP2,Class|SNP1)。由于G-test服从卡方分布，而G(X,Y|Z)＝I(X,Y|Z)*2*m，其中m为样本数，通过换算可以将互信息转换成G-test检验，再通过特定阈值筛选节点，从而得到不同节点的马尔可夫毯。
(2)收缩阶段：进行去除多余节点操作。对于任意目标节点SNP1(或者Class节点)，计算在去除节点SNP2的马尔可夫毯的条件下目标节点SNP1与SNP2的条件互信息I(SNP1,SNP2|mb(SNP1))。若条件互信息大于阈值，则删除节点SNP2并更新节点SNP1的马尔可夫毯。
(3)对称性检测：确保任意目标节点马尔可夫毯中的节点之间相互依赖。假设SNP1存在于SNP2的马尔可夫毯中，SNP1与SNP2相互依赖。与之对应，SNP2应该存在于SNP1的马尔可夫毯中。倘若两个节点的马尔可夫毯并不对称，那么分别在马尔可夫毯中删除另一节点。
3、在贝叶斯网络结构的空间中进行搜索时，本发明采用的人工蜂群优化算法中的每个蜜源对应一个贝叶斯网络结构。用邻接矩阵对贝叶斯网络个体进行表示，设网络中SNP位点和表型性状的个数为n，每个个体表示为n×n的邻接矩阵C。采用0/1的编码方案，如果节点i是节点j的父节点，Cij＝1，否则，Cij＝0。如图1所示。
人工蜂群算法中的初始蜜源(贝叶斯网络)由SNP节点、Class节点和节点间的边组成。基于各SNP节点以及Class节点的马尔可夫毯构建的初始网络(初始蜜源)，通过随机增加边、删除边、逆转边操作生成下一个网络个体，直到网络个体数量达到最大初始蜜源数。人工蜂群算法以初始蜜源作为起始点开始迭代。
4、在初始蜜源的基础上，通过人工蜂群优化算法的三种蜜蜂行为(采蜜蜂、观察蜂和侦察蜂)与贝叶斯网络的打分函数，对包括SNP位点和表型性状的贝叶斯网络进行演化，得到网络结构的最优解，进而快速准确的获取到影响表型性状的上位性基因位点(图3)。
(1)采蜜蜂阶段
首先采蜜蜂对初始蜜源的邻域进行搜索：通过加边、减边、逆转边3种操作(图2)对网络结构进行更新。然后利用BIC和MIT两种评分函数同时对新的网络结构(蜜源)进行打分，同时与初始网络结构(蜜源)的分数进行比较。
如果新蜜源网络结构优于初始蜜源，则将该网络结构进行保存，并记录网络结构提高的分数。最后，根据提高的分数使用轮盘赌的方式随机选择新蜜源网络结构。设第i个网络的适应度值为fi，则i被选中的概率Pi如Eq.(2)所示，其中N表示种群的大小。在该方法中，使网络分数提高越大的操作被选中的概率越大，这不仅可以保证网络演化的正确性，而且可以避免陷入局部最优。假如轮盘赌方法选择放弃本次加边、减边或者逆转边操作，则说明该蜜源邻域已无更优的网络结构。表示该网络结构已是最优的网络结构，采蜜蜂对其进行采蜜。
BIC打分函数的计算公式如Eq.(3)所示，其中m表示样本的总数量，n表示变量的个数，ri表示第i个变量的取值个数，qi表示第i个变量的父变量的组合个数，mijk表示第i个变量取第k个值，其父变量取第j个组合的样本个数。
MIT打分函数的计算公式如Eq.(4)所示，其中D表示特定数据集，G表示网络结构，n表示节点数目，N表示样本个数。I(Xi,∏Xi)表示目标节点Xi与其父节点集合∏Xi的互信息值，表示惩罚项，其中Si表示节点Xi的父节点总数，表示在自由度为置信区间为α时卡方分布的值。Eq.(4)中自由度的定义如Eq.(5)所示。
其中表示节点Xi的父节点排列顺序，其父节点排列顺序由父节点取值个数多少决定。取值个数越多，排列顺序越靠前。表示节点Xi某个父节点的取值个数，表示节点Xi最后一个父节点的取值个数。在Eq.(5)中，分为如下两种情况求自由度(1).当父节点个数为1时，通过(目标节点取值个数-1)*(父节点取值个数-1)进行计算。(2).当父节点个数大于1时，根据之前的父节点排列顺序，依次相乘父节点取值个数。其中目标节点取值个数和最后一个父节点取值个数需要减1后进行计算。
给定有N个样本的数据集D，假设X和Y在给定Z的情况下条件独立成立，可得2NI(X,Y|Z)近似于自由度为L＝(rX-1)*(rY-1)*rZ的卡方分布X2(L)，其中rX,rY和rZ分别表示变量X,Y和Z的取值个数。由此可得2NI(Xi,Xij|{Xi1,…,Xi(j-1)})近似于自由度为lij的卡方分布X2(lij)。给定置信区间α，当满足Eq.(6)时，可得进而得到Eq.(4)中惩罚项
(2)观察蜂阶段
观察蜂阶段和采蜜蜂阶段交替进行，每当采蜜蜂搜索一次邻域得到较优的网络结构之后，观察蜂便对采蜜蜂所在的蜜源进行考察(图3中阶段2)，根据两种打分函数计算对应网络结构的适应度值。然后利用轮盘赌的方式随机选择一个蜜源网络结构，适应度越大的网络结构被认为邻域拥有更优网络结构的可能性越大，被选到的几率也就越大。采蜜蜂对选中网络结构的邻域进行搜索(图3中阶段2-阶段3)，重复采蜜蜂阶段的操作。
(3)侦察蜂阶段
侦察蜂在采蜜蜂和观察蜂阶段中扮演监督员的角色，记录每个蜜源网络结构的采蜜次数以及最优蜜源网络结构。当某一蜜源网络结构达到最大采蜜次数，则放弃该蜜源，并通过步骤2中的方法重新生成新蜜源(图3中阶段4)。
5、可分解BIC评分函数
BIC评分函数的时间复杂度会随着父节点个数的增长而成指数增长，该方法无法对父节点个数较多的节点进行评分。为了解决此问题，本发明利用可分解BIC评分函数进行处理，其计算方式如Eq.(7)-Eq.(9)所示。
BIC*(X,Π1,Π2)＝BIC(X,Π1)+BIC(X,Π2)+inter(X,Π1,Π2) (7)
BIC(X,Π1∪Π2)＝BIC*(X,Π1,Π2)+N·(I(Π1,Π2|X)-I(Π1,Π2)) (9)
在Eq.(7)中，将X节点的父节点集合∏分解为∏1和∏2，Π1∪Π2＝Π，通过Eq.(8)和Eq.(9)计算得到近似BIC分数。在Eq.(9)计算互信息I(∏1,∏2)和条件互信息I(∏1,∏2|X)的过程中，∏1,∏2以及X的组合个数随父节点的增多呈爆炸式增长。计算机内存往往无法存储这么多组合数对应取值个数的频数，因此可分解BIC评分方法仍具有局限性，且实验结果表明该方法无法处理18个以上父节点。
考虑到显性上位性数据属于典型的高维少样本型数据，往往没有涵括所有的取值组合，因此本发明采用如下方法进行简化处理：仅仅统计所有样本包括的组合的频数并进行存储，未在样本内出现的组合的频数计为0，从而不对所有组合的频数进行统计。如图10所示，SNP的取值有以下三种：0/1/2，理论上两个SNP位点的取值组合类型有9种。而实际中，此两个SNP位点的取值类型往往少于9种。如图10中所示，SNP2和SNP3的组合类型数就只有6种。仅仅对这6种组合类型进行频数统计，而不是遍历所有的取值组合，从而大大减少计算的时间和空间复杂度。
6、算法终止
当达到最大迭代次数或者每只采蜜蜂均将对应蜜源开采至最大开采次数，则找到最优蜜源(最优网络结构)，算法结束。
7、通过实验说明人工蜂群优化贝叶斯网络的上位性位点挖掘方法的高效性，分别比较二节点和三节点上位性检测的准确率、F1-score、FPR和运行时间。
下面是应用本发明中的方法在GAMETES软件生成数据集上进行上位性位点挖掘的实例，通过相关的实验来详细说明本发明方法挖掘上位性位点的高效性。GAMETES软件是一款业界常用的用于生成Epistasis模拟数据的软件[9]，该软件可以快速准确地生成Epistasis模拟数据，通过改变不同的参数生成特定的两位点甚至多位点Epistasis模型。可以设置的参数包括：SNP位点的个数、遗传率(heritability,h2)、最小等位基因频率(MAF)以及患病率(prevalence)等。生成模拟数据的文件中第1行为位点名称，最后1列为Class标签，1表示患病，0表示对照。基因型数据用0，1，2表示，0表示纯合子常见基因型，1表示杂合子，2表示纯合子少见基因型。
将本发明中人工蜂群优化贝叶斯网络的上位性位点挖掘方法记为BnBeeEpi，实验比较的上位性检测方法包括以下几种：BEAM,AntEpiSeeker,BOOST，ESMO和未加入人工蜂群算法仅利用马尔可夫毯获取网络的方法omb-Fast。通过设置不同的遗传率h2(0.05,0.15,0.25,0.35,0.4)和最小等位基因频率MAF(0.1,0.2,0.3,0.4)，采用GAMETES软件生成了不同的数据集，每个参数设置下的数据集包括100个文件。利用了下述几个评价指标对上位性检测的效果进行比较：准确率(Eq.(10)),精确率(Eq.(11)),召回率(Eq.(12)),F-score(Eq.(13)),假阳性率(FalsePositiveRate,FPR)(Eq.(14))。其中准确率用于评价输出结果中能否包含所有正确结果。召回率是指输出结果中真阳性的数目和数据中总共真阳性数目的比率。在上位性挖掘研究中，较高的召回率意味着不管SNP组合数有多少，该算法也能检测到较全的SNP组合。精确率是指输出结果中真阳性的数目占输出总数的多少，较高的精确率意味着该算法的输出结果与标准结果吻合度较高。F1-score综合考虑了精确率与召回率这两个指标。FPR表示输出中假阳性结果的数目占数据中假阳性数目的多少。
实验1.2位点上位性检测的准确率和效率对比
在本实验中，比较了设置不同遗传率h2和最小等位基因频率MAF情况下上位性位点挖掘的准确率，F1-score以及FPR。图4，图5和图6分别给出了不同方法2位点上位性挖掘的准确率、F1-score和FPR的比较结果。表1给出了不同方法的2位点上位性挖掘的时间对比。
通过图4可见，BEAM方法的上位性检测准确率要明显低于其他几种方法，其他方法的准确率基本为100％。由图5可得，BOOST和BnBeeEpi的F1-score相差不大，两者的F1-score明显优于AntEpiSeeker,BEAM以及ESMO。omb-Fast的F1-score略低于BOOST和BnBeeEpi方法。由图6可得，BnBeeEpi的假阳性率(FPR)明显低于AntEpiSeeker以及ESMO方法，稍微高于BEAM和BOOST方法。但BEAM方法的准确性比较低，而BOOST只能用于两位点的上位性检测，并不能处理多位点上位性检测问题。另外，可见omb-Fast的假阳性率(FPR)稍微高于BEAM,BOOST与BnBeeEpi方法。并且由表1可以看出，omb-Fast的运行速度非常快，其速度优于BOOST。BnBeeEpi因为使用了人工蜂群算法进行优化处理增加了迭代次数，减少了假阳性，但牺牲了运行时间。
表1不同方法2-locus上位性检测的时间对比(单位:分)
综上可得，BnBeeEpi方法在保证检测准确率的情况下，有较好的F1-score以及较低的假阳性率(FPR)，整体上优于其他方法。另外，可见omb-Fast和BnBeeEpi的检测准确率均为100％，且omb-Fast方法的运行效率很高，该方法能快速的学习得到基本正确的网络结构。但BnBeeEpi方法的F1-score和FPR明显优于omb-Fast，这说明在学习得到基本网络结构的基础上，利用人工蜂群优化方法可以大大的减少上位性检测的假阳性率。但BnBeeEpi方法的上位性检测的效率远远低于omb-Fast方法，可见利用人工蜂群优化贝叶斯网络花费了一定的时间。
实验2.3位点上位性检测的准确率和效率对比
在本实验中，比较了不同方法在设置不同遗传率h2和最小等位基因频率MAF情况下的3位点上位性挖掘的准确率、F1-score以及FPR，如图7-图9所示。表2给出了不同方法的3位点上位性挖掘的时间对比。
由图7可得，BEAM方法的3-locus上位性检测的准确性最低，ESMO,omb-Fast与BnBeeEpi三种方法的准确性相差不大。由图8，图9以及表2可得，BnBeeEpi方法的F1-score和假阳性率(FPR)明显优于其他几种方法。虽然omb-Fast方法的准确性略高于BnBeeEpi方法，而且速度也非常快，但后者的假阳性率大大减少，可见将人工蜂群用于对贝叶斯网络进行演化，进而用于上位性位点检测方法的有效性。
表2不同方法3-locus上位性检测的时间对比(单位:分)
8、BIC可分解打分函数的实验对比
实验3.BIC与可分解BIC打分函数可处理父节点个数对比
分别在每组参数生成的模拟数据中随机抽取2个数据集进行实验，将原始的BIC与可分解BIC打分函数进行对比。由表3可得，BIC打分函数能够处理的最大父节点数为18个，可分解BIC打分函数能够处理的最大父节点数为30个。
表3采用可分解BIC打分方法的结果对比
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (10)

1.一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于，包括如下步骤：
S1，利用扩张、收缩、对称性检测三个阶段，将目标节点分为表型Class节点和非Class节点两种情况，通过条件互信息计算获取不同节点的马尔可夫毯；
S2，根据得到不同节点的马尔可夫毯，构建包括SNP位点和表型Class的网络，得到初始蜜源，即初始网络结构；
S3，在不生成环的前提下，基于初始蜜源网络个体，通过随机增加边、删除边、逆转边的操作生成新的网络个体，直到网络个体数量达到最大初始蜜源数目；
S4，在步骤3所得初始蜜源的基础上，通过人工蜂群优化算法的三种蜜蜂行为与贝叶斯网络的BIC和MIT两种打分函数，对包括SNP位点和表型Class的贝叶斯网络进行演化，得到较优的贝叶斯网络结构，从而准确快速地获取到影响表型性状的上位性基因位点，其中三种蜜蜂分别为采蜜蜂、观察蜂和侦察蜂。
2.如权利要求1所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于：在所述S1步骤前，将基因型数据表示为0，1，2形式的数据，对于特定SNP基因型数据，例如AT，进行如下表示：AA用0表示，TT用2表示，AT/TA用1表示，其中0表示纯合子常见基因型，1表示杂合子，2表示纯合子少见基因型，表型Class取值为0或1，0表示为未患病，1表示为患病。
3.如权利要求1所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于，在所述S1步骤中，扩张阶段具体为：在构建非Class节点的马尔可夫毯时，假设表型Class节点已经加入到目标SNP节点的马尔可夫毯中，计算任意一个SNP节点与Class在目标SNP节点条件下的条件互信息，并将所述条件互信息小于阈值的SNP位点加入目标节点的马尔可夫毯中；在构建Class节点的马尔可夫毯时，假设任意一个SNP节点已经加入表型Class节点的马尔可夫毯中，计算另外任意一个SNP节点与已加入马尔可夫毯的SNP节点在Class条件下的条件互信息，并将小于阈值的SNP位点加入Class节点的马尔可夫毯中。
4.如权利要求1所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于，在所述S1步骤中，收缩阶段具体为：去除多余节点；对于任意目标节点SNP1或者Class节点，计算在去除SNP2的马尔可夫毯的条件下目标节点SNP1与SNP2的条件互信息，若所述条件互信息大于阈值，则删除SNP2并更新节点SNP1的马尔可夫毯。
5.如权利要求1所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于，在所述S1步骤中，对称性检测具体为：由于本发明中的方法是基于条件独立性检测，因此任意目标节点马尔可夫毯中的节点之间应该相互依赖，假设SNP1存在于SNP2的马尔可夫毯中，SNP1与SNP2应该相互依赖，与之对应，SNP2应该存在于SNP1的马尔可夫毯中，假如两个节点的马尔可夫毯并不对称，那么分别在两个节点的马尔可夫毯中删除另外一节点。
6.如权利要求1所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于，在所述S4步骤中，采蜜蜂阶段具体为：
S40，首先采蜜蜂对初始蜜源的邻域进行搜索，即通过加边、减边、逆转边三种操作对网络结构进行更新；
S41，然后利用BIC和MIT两种评分函数同时对新的网络结构进行打分，并与初始网络结构的分数进行比较，若新的网络结构优于初始网络结构，则将该网络结构进行保存，并记录该网络结构提高的分数；
S42，最后，根据提高的分数使用轮盘赌的方式随机选择新的网络结构，此时采蜜蜂对其进行采蜜；
以上的网络结构即蜜源。
7.如权利要求6所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于：在所述采蜜蜂阶段的方法中，使网络分数提高越大的操作被选择的概率越大，这不仅可以保证网络演变的正确性，又不至于陷入局部最优，假如轮盘赌方法选择放弃本次加边、减边或者逆转边操作，则说明该蜜源邻域已无更优网络结构，表示该网络结构已是最优的网络结构，此时采蜜蜂对其进行采蜜。
8.如权利要求6所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于，在所述S4步骤中，观察蜂阶段具体为：观察蜂阶段和采蜜蜂阶段交替进行，每当采蜜蜂搜索一次邻域得到较优的网络结构之后，观察蜂便对采蜜蜂所在的蜜源进行考察，根据两种打分函数计算对应网络结构的适应度值，然后利用轮盘赌的方式随机选择一个蜜源网络结构，适应度越大的网络结构被认为邻域拥有更优网络结构的可能性越大，被选到的几率也就越大，采蜜蜂对选中网络结构的邻域进行搜索，重复采蜜蜂阶段的操作。
9.如权利要求8所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于，在所述S4步骤中，侦察蜂阶段具体为：侦察蜂在采蜜蜂和观察蜂阶段中扮演监督员的角色，记录每个蜜源网络结构的采蜜次数以及最优蜜源网络结构，当某一蜜源网络结构达到最大采蜜次数，则放弃该蜜源，并通过步骤1重新生成新的蜜源网络结构，假如每只采蜜蜂均将对应蜜源开采至最大采蜜次数，则找到最优蜜源，即最优网络结构，此时算法结束。
10.如权利要求1所述的一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法，其特征在于：适应度函数是用于评判网络结构好坏的标准，因此采用BIC与MIT联合打分的方法对网络的优劣进行判断，同时，由于BIC打分函数可以处理的父节点个数有限，因此采用可分解BIC打分方法处理较多个数的父节点，进而可以处理规模较大的网络结构。
CN201910859461.7A 2019-09-11 2019-09-11 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法 Active CN110570909B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910859461.7A CN110570909B (zh) 2019-09-11 2019-09-11 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910859461.7A CN110570909B (zh) 2019-09-11 2019-09-11 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法

Publications (2)

Publication Number Publication Date
CN110570909A true CN110570909A (zh) 2019-12-13
CN110570909B CN110570909B (zh) 2023-03-03

Family

ID=68779336

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910859461.7A Active CN110570909B (zh) 2019-09-11 2019-09-11 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法

Country Status (1)

Country Link
CN (1) CN110570909B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111833967A (zh) * 2020-07-10 2020-10-27 华中农业大学 基于k-tree优化贝叶斯网络的上位性位点挖掘方法
CN114550827A (zh) * 2022-01-14 2022-05-27 山东师范大学 一种基因序列比对方法及系统
CN115577313A (zh) * 2022-12-09 2023-01-06 中国南方电网有限责任公司超高压输电公司广州局 电力电网故障诊断方法、装置、计算机设备及存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120251502A1 (en) * 2008-10-24 2012-10-04 The Government of the US as Represented by the Secretary of the Dept. of health Human Ebola Virus Species and Compositions and Methods Thereof
CN105205344A (zh) * 2015-05-18 2015-12-30 上海交通大学 基于多目标蚁群优化算法的基因位点挖掘方法
CN109448794A (zh) * 2018-10-31 2019-03-08 华中农业大学 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法
CN109472320A (zh) * 2018-11-30 2019-03-15 南京农业大学 不确定条件下作物生育期模型品种参数自动校正框架

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120251502A1 (en) * 2008-10-24 2012-10-04 The Government of the US as Represented by the Secretary of the Dept. of health Human Ebola Virus Species and Compositions and Methods Thereof
CN105205344A (zh) * 2015-05-18 2015-12-30 上海交通大学 基于多目标蚁群优化算法的基因位点挖掘方法
CN109448794A (zh) * 2018-10-31 2019-03-08 华中农业大学 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法
CN109472320A (zh) * 2018-11-30 2019-03-15 南京农业大学 不确定条件下作物生育期模型品种参数自动校正框架

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
NATHAN FORTIER等: "《Abductive inference in Bayesian networks using distributed overlapping swarm intelligence》", 《SOFT COMPUTING》 *
张平 等: "《基于人工蜂群算法的贝叶斯网络结构学习》", 《智能系统学报》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111833967A (zh) * 2020-07-10 2020-10-27 华中农业大学 基于k-tree优化贝叶斯网络的上位性位点挖掘方法
CN111833967B (zh) * 2020-07-10 2022-05-20 华中农业大学 基于k-tree优化贝叶斯网络的上位性位点挖掘方法
CN114550827A (zh) * 2022-01-14 2022-05-27 山东师范大学 一种基因序列比对方法及系统
CN114550827B (zh) * 2022-01-14 2022-11-22 山东师范大学 一种基因序列比对方法及系统
CN115577313A (zh) * 2022-12-09 2023-01-06 中国南方电网有限责任公司超高压输电公司广州局 电力电网故障诊断方法、装置、计算机设备及存储介质

Also Published As

Publication number Publication date
CN110570909B (zh) 2023-03-03

Similar Documents

Publication Publication Date Title
CN110570909B (zh) 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法
JP6888123B2 (ja) 深層畳み込みニューラルネットワークを事前訓練するための深層学習ベースの技術
JP7247253B2 (ja) 経験的バリアントスコア（ｅｖｓ）ベースの深層学習バリアントコーラ
KR102223129B1 (ko) 심층 학습 기반 비정상 스플라이싱 검출
Mejía-Guerra et al. A k-mer grammar analysis to uncover maize regulatory architecture
CN109448794B (zh) 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法
Boekhout et al. The evolving species concepts used for yeasts: from phenotypes and genomes to speciation networks
KR20230084319A (ko) 심층 컨볼루션 신경망을 트레이닝하기 위한 심층 학습-기반 기술
WO2020014280A1 (en) DEEP LEARNING-BASED FRAMEWORK FOR IDENTIFYING SEQUENCE PATTERNS THAT CAUSE SEQUENCE-SPECIFIC ERRORS (SSEs)
AU2019272065B2 (en) Deep learning-based framework for identifying sequence patterns that cause sequence-specific errors (SSEs)
CN107025386B (zh) 一种基于深度学习算法进行基因关联分析的方法
Lee et al. BNTagger: improved tagging SNP selection using Bayesian networks
Whitehouse et al. Timesweeper: accurately identifying selective sweeps using population genomic time series
CN117637033A (zh) 一种基于遗传知识库的新致病基因位点机器学习挖掘方法
Yang et al. BnBeeEpi: An approach of epistasis mining based on artificial bee colony algorithm optimizing Bayesian network
CN116153396A (zh) 一种基于迁移学习的非编码变异预测方法
Lopez et al. Extracting biological knowledge by fuzzy association rule mining
Zhang et al. Tree‐based linkage and association analyses of asthma
CN116994645B (zh) 基于交互式推理网络的piRNA与mRNA靶标对的预测方法
CN115995262B (zh) 基于随机森林及lasso回归解析玉米遗传机理的方法
Huang Sparse model learning for inferring genotype and phenotype associations
Abeysinghe et al. A Note on Plant Virus Images for Use in Machine Learning.
Sitarcık et al. Flower Pollination Algorithm for Detection of Epistasis Associated with a Phenotype
Mejía-Guerra et al. k-mer grammar uncovers maize regulatory architecture
Chan EVALUATING AND CREATING GENOMIC TOOLS FOR CASSAVA BREEDING

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