CN111833967A - 基于k-tree优化贝叶斯网络的上位性位点挖掘方法 - Google Patents
基于k-tree优化贝叶斯网络的上位性位点挖掘方法 Download PDFInfo
- Publication number
- CN111833967A CN111833967A CN202010683358.4A CN202010683358A CN111833967A CN 111833967 A CN111833967 A CN 111833967A CN 202010683358 A CN202010683358 A CN 202010683358A CN 111833967 A CN111833967 A CN 111833967A
- Authority
- CN
- China
- Prior art keywords
- node
- tree
- nodes
- snp
- class
- 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
Images
Classifications
-
- 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
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Biophysics (AREA)
- Software Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Data Mining & Analysis (AREA)
- Databases & Information Systems (AREA)
- Epidemiology (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- Public Health (AREA)
- Bioethics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- General Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
Abstract
本发明提供一种基于k‑tree优化贝叶斯网络的上位性位点挖掘方法。均匀抽样蒲公英代码,得到每个节点的邻居节点集合,构建对应的特征树结构;利用与k‑clique中节点进行合并组成新节点的方法,构建得到包括大规模SNP位点和表型性状节点的Rényi k‑tree,得到相应的k‑tree;利用基于度选择的树分解算法将k‑tree对应的图分解为不同的k‑clique;利用优化马尔科夫毯的Fast‑IAMB算法获取不同节点的马尔可夫毯,学习得到不同k‑clique对应的子贝叶斯网络结构;合并子网络得到整体网络结构。重复以上步骤,得到包括SNP位点和表型性状节点的网络结构,进而得到影响表型性状的上位性位点。本发明可以帮助生物学研究者得到影响特定表型性状的上位性基因位点,为不同物种复杂数量性状的遗传基础解析提供借鉴。
Description
技术领域
本发明属于生物信息技术领域,特别是涉及一种基于k-tree优化贝叶斯网络的上位性位点挖掘方法。
背景技术
随着科学技术的发展和医疗水平的提升,过去长期困扰人们的一些传染病基本得到了控制,而复杂疾病成为了目前影响人类健康的主要疾病。复杂疾病占人类疾病的大约80%以上,对人类健康造成了极大的伤害。哮喘、癌症、糖尿病、高血压、老年痴呆症、类风湿性关节炎、精神分裂症、心脏病、心血管疾病、肥胖、肿瘤等常见慢性疾病,统称为复杂疾病。目前传统的遗传流行病和生物医学模式向可预防的精准化医学模式转变,能够为预防和治疗困扰人类多年的复杂疾病提出新的解决方案。
复杂疾病不符合孟德尔遗传定律,潜在的遗传机理十分复杂。除了基因主效应的影响,还包括基因-基因相互作用,基因-环境相互作用等不符合孟德尔遗传定律的因素,给研究者带来了困难和挑战。因此,急需阐明复杂疾病的致病原因及遗传机制,为复杂疾病的诊断和治疗提供科学依据,进而为人类健康提供保障。通过生物学大量实验研究发现,控制生物复杂性状的主要原因是基因与基因之间的相互作用。孟德尔经典实验只是在简单性状上效果良好,仅仅能解释很小部分的遗传变异。上位性主要指SNP之间的作用,上位性效应检测可以识别出复杂性状位点间相关联的遗传信息,从而解释复杂性状背后的遗传变异,揭露出遗传机理,进而在一定程度上解决“遗传丢失性”的现象。另外,目前全基因组关联分析方法被成功应用于挖掘检测影响表型性状的基因位点,但是该方法只能解释小部分的遗传变异。其中主要原因是该方法只能检测主要基因,而忽略了基因之间的相互作用,即上位性。因此,在大规模的全基因组数据范围内提出更有效,更准确的上位性检测方法具有重要意义,也对复杂疾病致病机理的发现、诊断、治疗和预防有着非常重要的作用。
发明内容
本发明目的在于针对现有技术的不足,提出一种基于k-tree优化贝叶斯网络的上位性位点挖掘方法。具体技术方案如下:
一种基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:包含以下步骤:
步骤1,将基因型SNP和表型Class看成节点,并将SNP基因型数据表示为0/1/2形式的数据,表型Class表示为0/1型的数据;
步骤2,基于节点个数n和树宽k生成表示网络图结构的蒲公英代码,得到节点集合和边标记集合;然后识别环结构并进行去环操作,得到每个节点的邻居节点集合,进而构建得到蒲公英代码对应的特征树结构;
步骤2.1,基于节点个数n和树宽k,通过随机生成蒲公英代码,求得节点集合p和边标记集合l;
步骤2.2,识别并去掉环结构;当节点本身成环时,直接去掉该环;当多个节点形成环时,通过交换最大编号与最小编号节点位置的方法去掉环结构,从而更新节点集合p和边标记集合l;
步骤2.3,找到每个节点的邻居节点,构建特征树;首先根据节点集合p,得到不同节点i的邻居节点用around[i]表示,计算得到neighbor[p[i]]的值;然后判断around与neighbor中对应位置的元素个数是否相等;如果相等,则计算around[i]-neighbor[i]得到节点i的邻居节点;如果不相等,交换neighbor中第一个位置以及第一个偶数位置的值,然后通过计算around[i]-neighbor[i]得到节点i的邻居节点;
步骤2.4,基于更新后的边标记集合l,根据从右侧分支到左侧分支的顺序,对节点间的边进行标记,包括将边标记为a,b,c,其中与节点0相连的边标记为ε,进而生成特征树结构;
步骤3,根据生成的特征树,生成特征树的骨架,进而将特定节点v与对应k-clique中的节点相连生成Rényi k-tree;
步骤3.1,对于节点数n和树宽k,以k=3参考,用{n,n+1,n+2}表示根节点0,k-clique表示朝着根方向的特定节点的邻居节点集;将节点v与对应k-clique中节点进行合并组成新的节点v∪C,C表示k-clique,得到特征树的骨架;
步骤3.2,将{n,n+1,n+2}作为根节点,对特征树骨架中的形式为v∪C的节点进行遍历,依次连接v与C中各个节点,从而得到Rényi k-tree;
步骤4,生成k-tree,在Rényi k-tree的基础上,得到各个节点的邻居节点,进的得到包括SNP和Class节点的k-tree结构图;
步骤5,采用基于度选择的算法将k-tree对应的图进行树分解,得到不同的团;
步骤6,利用马尔科夫优化的贝叶斯网络方法构建不同k-clique内节点对应的子网络;将目标节点分为表型Class节点和非Class节点两种情况,得到不同团对应的子网络结构;
步骤6.1扩张阶段:对于表型Class节点,在初始阶段假设任意一个节点SNP1已经加入Class节点的马尔可夫毯中,计算另外任意一个节点SNP2与已加入马尔可夫毯的SNP1节点在Class条件下的条件互信息I(SNP1,SNP2|Class),将条件互信息小于阈值的SNP2位点加入Class节点的马尔可夫毯中;对于非Class节点,在初始阶段假设表型Class节点已经加入到目标SNP1节点的马尔可夫毯中,计算任意一个节点SNP2与Class在目标节点SNP1条件下的条件互信息I(SNP2,Class|SNP1);将条件互信息小于阈值的SNP2位点加入目标节点SNP1的马尔可夫毯中;
步骤6.2收缩阶段:去除马尔科夫毯中多余节点;对于任意目标节点SNP1或者Class节点,计算在去除节点SNP2的马尔可夫毯的条件下目标节点SNP1与SNP2的条件互信息;若条件互信息大于阈值,则删除节点SNP2并更新节点SNP1的马尔可夫毯;
步骤6.3对称性检测:确保任意目标节点马尔可夫毯中的节点之间相互依赖;假设SNP1存在于SNP2的马尔可夫毯中,SNP1与SNP2相互依赖;与之对应,SNP2应该存在于SNP1的马尔可夫毯中;倘若两个节点的马尔可夫毯并不对称,则分别在每个节点的马尔可夫毯中删除另一节点;
步骤7,在构建得到不同k-clique内节点对应子网络的基础上,对各个子网络进行合并,得到包括SNP位点和表型性状节点网络结构;
步骤8,重复执行步骤1-步骤6,利用均匀抽样蒲公英代码的方法将网络图的构建转换为特定编码的生成问题,构建得到包括SNP位点和表型性状节点的网络结构;然后根据网络中边出现的次数,得到最终的包括SNP位点和表型性状节点网络结构,进而根据节点间的联系,得到影响表型性状的上位性位点。
具体地,在步骤1中,将SNP基因型数据用0/1/2表示,假设AT为参考对象,进行如下表示:AA用0表示,TT用2表示,AT/TA用1表示,其中0表示纯合子常见基因型,1表示杂合子,2表示纯合子少见基因型;表型Class用0/1表示,0表示为未患病,1表示为患病。
具体地,在步骤3.1中,根据除了与根节点相连的k-clique,与其它节点相连的k-clique不能重复的规则,将节点v与对应k-clique中节点进行合并组成新的节点,得到特征树的骨架。
具体地,在步骤5中,对于图的每一条边,至少存在一个团包含边的两个顶点;假设t1,t2,t3是图中节点,t2在t1到t3的路径上,用Xt1,Xt2,Xt3表示不同的团;首先找到度最小的节点,根据树分解规则(节点v∈Xt1,v∈Xt3,可得v∈Xt2)对k-tree对应的图进行树分解,得到不同的团。
具体地,在步骤6中,考虑到上位性位点挖掘问题的特殊性,将表型性状加入到每个团,进行子网络的构建;利用扩张、收缩、对称性检测三个阶段,通过条件互信息计算获取不同节点的马尔可夫毯,进而构建得到不同团对应的子网络结构。
具体地,在步骤8中,通过均匀抽样蒲公英代码重复执行以上特征树构建,Rényik-tree生成,k-tree生成,k-tree分解生成不同的团,子网络构建,子网络合并步骤,得到最终包括SNP位点和表型性状节点的网络结构,进而得到影响表型性状的上位性位点。
本发明提供的解决方案给出了一种基于k-tree优化贝叶斯网络构建包括SNP位点和表型性状Class的网络结构,进而挖掘上位性位点的方法。本发明所提供的技术流程如图1所示。本发明根据节点个数n和树宽k,得到表示网络图结构的蒲公英代码。然后进行去环结构操作,得到每个节点的邻居节点集合,进而构建得到蒲公英代码对应的特征树结构。利用与k-clique中节点进行合并组成新节点的方法,并依据与非根节点相连的k-clique不能重复的原则,构建得到包括大规模SNP位点和表型性状节点的Rényi k-tree,进而构建得到相应的k-tree。然后利用基于度选择的树分解算法将k-tree对应的图分解为不同的k-clique。通过扩张、收缩、对称性检测三个阶段,利用优化马尔科夫毯的Fast-IAMB算法获取不同节点的马尔可夫毯,快速准确的学习得到不同k-clique对应的子贝叶斯网络结构。最后进行子网络合并,得到包括SNP位点和表型性状节点的网络结构。重复执行以上步骤,利用均匀抽样蒲公英代码的方法,通过生成特征树,生成Rényi k-tree,生成k-tree,k-tree分解生成不同的团,子网络构建,子网络合并等方法得到包括SNP位点和表型性状节点网络结构,进而根据节点间的联系,得到影响表型性状的上位性位点。本发明可以帮助生物学研究者得到影响特定表型性状的上位性基因位点,进而辅助基因功能挖掘,以及为不同物种复杂数量性状的遗传基础解析提供借鉴。
附图说明
图1 为本发明具体实施的流程示意图;
图2 Rényi k-tree的生成图;
图3 k-tree图结构;
图4 k-tree图分解过程图;
具体实施方式
本发明为了解决其技术问题采用了如下的技术方案:
一种基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:包含以下步骤:
步骤1,将基因型SNP和表型Class看成节点,并将SNP基因型数据表示为0/1/2形式的数据,表型Class表示为0/1型的数据;
步骤2,基于节点个数n和树宽k生成表示网络图结构的蒲公英代码,得到节点集合和边标记集合;然后识别环结构并进行去环操作,得到每个节点的邻居节点集合,进而构建得到蒲公英代码对应的特征树结构;
步骤2.1,基于节点个数n和树宽k,通过随机生成蒲公英代码,求得节点集合p和边标记集合l;
步骤2.2,识别并去掉环结构;当节点本身成环时,直接去掉该环;当多个节点形成环时,通过交换最大编号与最小编号节点位置的方法去掉环结构,从而更新节点集合p和边标记集合l;
步骤2.3,找到每个节点的邻居节点,构建特征树;首先根据节点集合p,得到不同节点i的邻居节点用around[i]表示,计算得到neighbor[p[i]]的值;然后判断around与neighbor中对应位置的元素个数是否相等;如果相等,则计算around[i]-neighbor[i]得到节点i的邻居节点;如果不相等,交换neighbor中第一个位置以及第一个偶数位置的值,然后通过计算around[i]-neighbor[i]得到节点i的邻居节点;
步骤2.4,基于更新后的边标记集合l,根据从右侧分支到左侧分支的顺序,对节点间的边进行标记,包括将边标记为a,b,c,其中与节点0相连的边标记为ε,进而生成特征树结构;
步骤3,根据生成的特征树,生成特征树的骨架,进而将特定节点v与对应k-clique中的节点相连生成Rényi k-tree;
步骤3.1,对于节点数n和树宽k,以k=3参考,用{n,n+1,n+2}表示根节点0,k-clique表示朝着根方向的特定节点的邻居节点集;将节点v与对应k-clique中节点进行合并组成新的节点v∪C,C表示k-clique,得到特征树的骨架;
步骤3.2,将{n,n+1,n+2}作为根节点,对特征树骨架中的形式为v∪C的节点进行遍历,依次连接v与C中各个节点,从而得到Rényi k-tree;
步骤4,生成k-tree,在Rényi k-tree的基础上,得到各个节点的邻居节点,进的得到包括SNP和Class节点的k-tree结构图;
步骤5,采用基于度选择的算法将k-tree对应的图进行树分解,得到不同的团;
步骤6,利用马尔科夫优化的贝叶斯网络方法构建不同k-clique内节点对应的子网络;将目标节点分为表型Class节点和非Class节点两种情况,得到不同团对应的子网络结构;
步骤6.1扩张阶段:对于表型Class节点,在初始阶段假设任意一个节点SNP1已经加入Class节点的马尔可夫毯中,计算另外任意一个节点SNP2与已加入马尔可夫毯的SNP1节点在Class条件下的条件互信息I(SNP1,SNP2|Class),将条件互信息小于阈值的SNP2位点加入Class节点的马尔可夫毯中;对于非Class节点,在初始阶段假设表型Class节点已经加入到目标SNP1节点的马尔可夫毯中,计算任意一个节点SNP2与Class在目标节点SNP1条件下的条件互信息I(SNP2,Class|SNP1);将条件互信息小于阈值的SNP2位点加入目标节点SNP1的马尔可夫毯中;
步骤6.2收缩阶段:去除马尔科夫毯中多余节点;对于任意目标节点SNP1或者Class节点,计算在去除节点SNP2的马尔可夫毯的条件下目标节点SNP1与SNP2的条件互信息;若条件互信息大于阈值,则删除节点SNP2并更新节点SNP1的马尔可夫毯;
步骤6.3对称性检测:确保任意目标节点马尔可夫毯中的节点之间相互依赖;假设SNP1存在于SNP2的马尔可夫毯中,SNP1与SNP2相互依赖;与之对应,SNP2应该存在于SNP1的马尔可夫毯中;倘若两个节点的马尔可夫毯并不对称,则分别在每个节点的马尔可夫毯中删除另一节点;
步骤7,在构建得到不同k-clique内节点对应子网络的基础上,对各个子网络进行合并,得到包括SNP位点和表型性状节点网络结构;
步骤8,重复执行步骤1-步骤6,利用均匀抽样蒲公英代码的方法将网络图的构建转换为特定编码的生成问题,构建得到包括SNP位点和表型性状节点的网络结构;然后根据网络中边出现的次数,得到最终的包括SNP位点和表型性状节点网络结构,进而根据节点间的联系,得到影响表型性状的上位性位点。
具体地,在步骤1中,将SNP基因型数据用0/1/2表示,假设AT为参考对象,进行如下表示:AA用0表示,TT用2表示,AT/TA用1表示,其中0表示纯合子常见基因型,1表示杂合子,2表示纯合子少见基因型;表型Class用0/1表示,0表示为未患病,1表示为患病。
具体地,在步骤3.1中,根据除了与根节点相连的k-clique,与其它节点相连的k-clique不能重复的规则,将节点v与对应k-clique中节点进行合并组成新的节点,得到特征树的骨架。
具体地,在步骤5中,对于图的每一条边,至少存在一个团包含边的两个顶点;假设t1,t2,t3是图中节点,t2在t1到t3的路径上,用Xt1,Xt2,Xt3表示不同的团;首先找到度最小的节点,根据树分解规则(节点v∈Xt1,v∈Xt3,可得v∈Xt2)对k-tree对应的图进行树分解,得到不同的团。
具体地,在步骤6中,考虑到上位性位点挖掘问题的特殊性,将表型性状加入到每个团,进行子网络的构建;利用扩张、收缩、对称性检测三个阶段,通过条件互信息计算获取不同节点的马尔可夫毯,进而构建得到不同团对应的子网络结构。
具体地,在步骤8中,通过均匀抽样蒲公英代码重复执行以上特征树构建,Rényik-tree生成,k-tree生成,k-tree分解生成不同的团,子网络构建,子网络合并步骤,得到最终包括SNP位点和表型性状节点的网络结构,进而得到影响表型性状的上位性位点。
本发明的具体实现过程,还可以采用如下步骤:
1、用0/1/2形式的数据表示基因型数据,如对SNP基因型为AT的数据表示如下:AA用0表示,TT用2表示,AT/TA用1表示。0表示纯合子常见基因型,1表示杂合子,2表示纯合子少见基因型。Class表示表型性状,其中Class=1表示case(患病),Class=0表示control(未患病,对照)。将SNP和Class看作贝叶斯网络中的节点。
2、将蒲公英代码表示为(Q,S),n表示变量个数,k表示树宽,S为(n-k-2)*2的整数矩阵,其中元素(i,j),1≤i≤n-k and 1≤j≤k,或者为(0,ò),ò为不在{0,1,..,n-1}中任意数。对于n=9,k=3,Q=[0,1,8],(9-3-2)*2的S如Eq.(1)所示。
计算得到以下参数:p,ε,m,s。
(1)p为不在Q中的最小编号节点,得到p=2。
(2)ε为从蒲公英代码(Q,S)中返回的来自Q的向量值组合。计算的规则是:首先将[Q[i]]赋值为n-k+i。然后进行判断,如果qi是Q中最小的节点,则对于每一个不属于Q的节点则对于节点t∈{n-k+1,...,n}-Q,
(3)m表示向量ε反转后的向量。例如,对于ε=[6,7,2,3,4,5,0,1,8],得到m=[8,1,0,5,4,3,2,7,6]。
基于得到的ε,m和s,设节点集合H=[0,1,2,3,4,5,6,7,8],边标记集合L=[ε,a,b,c]。根据s=1,将(0,-1)插入矩阵S第一行,得到新的矩阵如Eq.(2)所示。其中0表示初始虚拟根节点,-1为设置的初始值。
然后利用以下步骤解码矩阵S,进而构建得到特征树。
(1)求节点的相邻节点集合和边标记集合。
用p表示节点集合,l表示边的标记集合。设置jump=v=0,对节点v从0到n进行遍历,将标记为v-jump的节点加入节点集合,v-jump对应的标记加入标记集合。为了表示方便,将用数字形式表示的边的标记用对应的字母表示。
对于上述例子,得到v=0,jump=0+1=1。v=1,jump=1,v-jump=0,将节点0加入到p中,将0对应的字母ε加入到l中。v=2,jump=1,v-jump=1,节点1加入到p中,将1对应的字母a加入到l中。按照此计算过程,得到p=[0,1,3,4,2,5,6,8,7],l=[ε,a,b,c]。
(2)判断并去掉环结构。当节点本身成环时,直接去掉该环。当多个节点形成环时,通过交换最大编号与最小编号节点位置的方法去掉环结构。
(3)找到每个节点的邻居节点,构建特征树。
首先根据节点集合p,得到不同节点i的邻居节点around[i],计算得到neighbor[p[i]]的值。然后判断around与neighbor中对应位置的元素个数是否相等。如果相等,则计算around[i]-neighbor[i]得到节点i的邻居节点。如果不相等,将neighbor中相应位置的元素加入的值,交换neighbor中第一个位置以及第一个偶数位置的值,通过计算around[i]-neighbor[i]得到节点i的邻居节点。
例如,假设去环操作后得到p=[0,0,2,4,3,8,6,5,7],l=[ε,a,b,c]。当i=6时,
j=6,around[6]={5,8}
j=7,around[6]=around[6]∪7=[5,8,7]
j=8,around[6]=around[6]∪8=[5,8]
依次根据i,j的值计算可以得到其余around集合,然后从0到len(p[v])=9进行遍历,得到如下。
neighbor[p[0]]=neighbor[0]=neighbor[p[0]]∪0=neighbor[[0,0]]
neighbor[p[1]]=neighbor[0]=neighbor[[0,1]]
neighbor[p[2]]=neighbor[2]=neighbor[[2,2]]
依次计算得到neighbor值,around第一个位置元素为[5,8,7],与neighbor第一个位置[0,0]个数不相等,则使neighbor中的[0,0]元组加入得到[0,0,5]。然后交换neighbor中第一个位置与第一个偶数位置的值2,得到集合[2,0,5]。得到around[5,8,7]-neighbor[2,0,5]=[5-2,8-0,7-5],得到与节点0相连的节点集合为[3,8,2]。然后对第二个位置的元组个数进行判断,[5,8]与[0,1]个数相等,则计算around[5,8]-neighbor[0,1]得到与1相连的节点集合为[5,7]。依次求出其它节点相连的其他节点,按照从右侧分支到左侧分支的顺序,将节点间的边标记为a,b,c,其中与节点0相连的边标记为ε,得到如图2(a)所示的特征树。
3、生成Rényi k-tree。特征树和它的Rényi k-tree具有一对一关系,由特征树生成Rényi k-tree主要包括以下步骤:
(1)对于节点数n和树宽k,例如k=3,用{n,n+1,n+2}表示根节点0,k-clique表示朝着根方向的特定节点的邻居节点集。根据除了与根节点相连的k-clique,与其它节点相连的k-clique不能重复的规则。从根节点到叶子节点的顺序,对特征树中的节点v开始遍历。将节点v与对应k-clique中节点进行合并组成新的节点,表示为v∪C(C表示k-clique)。然后,将节点v与新构建的节点进行相连,得到特征树的骨架,其中与根节点相连的边标记为ε。
例如,对于图2(a)中的特征树,给定根为(9,10,11),得到{3}∪{9,11,10},{8}∪{9,11,10},{2}∪{9,11,10}。对于图2(a)中的节点5和6,得到{5}∪{8,9,10}。对于节点6,若合并节点得到{6}∪{8,9,10},则不满足与非根节点相连的k-clique不能重复的规则,因此得到{6}∪{8,9,11}。同理,得到{4}∪{2,10,11},{1}∪{5,8,9},{7}∪{1,5,8},从而得到特征树的骨架,如图2(b)所示。
(2)将{n,n+1,n+2}作为根节点,对特征树骨架中的形式为v∪C的节点进行遍历,依次连接v与C中各个节点,从而得到Rényi k-tree。
例如,对于图2(b)中的树结构,与节点3相连的节点包括9,10,11,与5相连的节点为8,9,10,与7相连的节点为1,5,8等,从而得到如图2(c)所示的Rényi k-tree。
4、生成k-tree。其根节点看作Q,根据节点个数n和树宽k,用二维数组t表示Rényik-tree,对节点u∈[0,n]进行遍历,变量i从0到t[u]进行循环遍历,得到v=t[u,i],利用步骤1中的方法计算得到与进而得到各节点的邻居节点
例如,对于Q={9,10,11},n=11,k=3,根据图2(c)得到二维数组t=[[1,6,8][3,9,5,7,8][4,6,9,10,11],...]。
根据上述计算过程,得到各个节点的邻居节点,从而得到k-tree的部分图如图3所示。
5、k-tree分解生成不同的团。对于图的每一条边,至少存在一个团包含边的两个顶点。假设t1,t2,t3是图中节点,t2在t1到t3的路径上,用Xt1,Xt2,Xt3表示不同的团。首先找到度最小的节点,根据树分解规则(节点v∈Xt1,v∈Xt3,可得v∈Xt2)对k-tree对应的图进行树分解,得到不同的团。
例如,如图4(a)所示的k-tree对应的图结构,节点8为度最小的节点,与节点7相连组成团,如图4(b)所示。由于树分解的结果不唯一,在此例子中,我们选择这两个节点组成团。由于不同节点有相同度等情况的存在,树分解的结果不唯一。然后选择节点1作为下一个度最小的节点,根据树分解规则,对于图的每一条边,至少存在一个团包含边的两个顶点。选择与1相连的节点2,4组成团,如图4(c)所示。根据树分解规则,节点5在路径3->5->7上,选择节点4,5,7组成团,节点3,5,6组成团,节点5,6,7组成团,从而满足节点5属于团{3,5,6},同时属于团{5,6,7},那么节点5一定属于团{4,5,7}。根据树分解规则所有团中的节点需要涵盖图中的所有节点,选择节点6,7,9组成团,得到如图4(e)所示的团。
6、构建子网络。利用扩张、收缩、对称性检测三个阶段,将目标节点分为表型Class节点和非Class节点两种情况,通过条件互信息计算,从优化马尔科夫毯的角度获取不同节点的马尔可夫毯,从而构建不同团内节点对应的子网络结构。
(1)扩张阶段:对于表型Class节点,在初始阶段假设任意一个节点SNP1已经加入Class节点的马尔可夫毯中,通过Eq.(1)计算另外任意一个节点SNP2与已加入马尔可夫毯的SNP1节点在Class条件下的条件互信息I(SNP1,SNP2|Class)。对非Class节点,在初始阶段假设表型Class节点已经加入到目标SNP1节点的马尔可夫毯中,通过Eq.(3)计算任意一个节点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的马尔可夫毯中。假如两个节点的马尔可夫毯并不对称,那么分别在两个节点的马尔可夫毯中删除另外一节点。
7、子网络合并构建得到网络结构。在利用上述步骤构建得到不同k-clique内节点对应的子网络的基础上,对各个子网络进行合并,得到最终的包括SNP位点和表型性状节点的网络结构。
8、重复多次执行以上步骤1-步骤6,构建得到多个包括大规模SNP位点和表型性状节点的网络结构。根据网络中边出现的次数,得到最终包括SNP位点和表型性状节点的网络结构。根据网络中节点间的联系,得到影响表型性状的上位性位点,算法结束。
本文中所做的步骤围绕本发明方法进行说明,本发明所属技术领域的技术人员可以对其做各种各样的修改或补充或采用类似的方式替代,但并不会偏离本发明的精神或者超越本方法所定义的范围。
Claims (6)
1.一种基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:包含以下步骤:
步骤1,将基因型SNP和表型Class看成节点,并将SNP基因型数据表示为0/1/2形式的数据,表型Class表示为0/1型的数据;
步骤2,基于节点个数n和树宽k生成表示网络图结构的蒲公英代码,得到节点集合和边标记集合;然后识别环结构并进行去环操作,得到每个节点的邻居节点集合,进而构建得到蒲公英代码对应的特征树结构;
步骤2.1,基于节点个数n和树宽k,通过随机生成蒲公英代码,求得节点集合p和边标记集合l;
步骤2.2,识别并去掉环结构;当节点本身成环时,直接去掉该环;当多个节点形成环时,通过交换最大编号与最小编号节点位置的方法去掉环结构,从而更新节点集合p和边标记集合l;
步骤2.3,找到每个节点的邻居节点,构建特征树;首先根据节点集合p,得到不同节点i的邻居节点用around[i]表示,计算得到neighbor[p[i]]的值;然后判断around与neighbor中对应位置的元素个数是否相等;如果相等,则计算around[i]-neighbor[i]得到节点i的邻居节点;如果不相等,交换neighbor中第一个位置以及第一个偶数位置的值,然后通过计算around[i]-neighbor[i]得到节点i的邻居节点;
步骤2.4,基于更新后的边标记集合l,根据从右侧分支到左侧分支的顺序,对节点间的边进行标记,包括将边标记为a,b,c,其中与节点0相连的边标记为ε,进而生成特征树结构;
步骤3,根据生成的特征树,生成特征树的骨架,进而将特定节点v与对应k-clique中的节点相连生成Rényi k-tree;
步骤3.1,对于节点数n和树宽k,以k=3参考,用{n,n+1,n+2}表示根节点0,k-clique表示朝着根方向的特定节点的邻居节点集;将节点v与对应k-clique中节点进行合并组成新的节点v∪C,C表示k-clique,得到特征树的骨架;
步骤3.2,将{n,n+1,n+2}作为根节点,对特征树骨架中的形式为v∪C的节点进行遍历,依次连接v与C中各个节点,从而得到Rényi k-tree;
步骤4,生成k-tree,在Rényi k-tree的基础上,得到各个节点的邻居节点,进的得到包括SNP和Class节点的k-tree结构图;
步骤5,采用基于度选择的算法将k-tree对应的图进行树分解,得到不同的团;
步骤6,利用马尔科夫优化的贝叶斯网络方法构建不同k-clique内节点对应的子网络;将目标节点分为表型Class节点和非Class节点两种情况,得到不同团对应的子网络结构;
步骤6.1扩张阶段:对于表型Class节点,在初始阶段假设任意一个节点SNP1已经加入Class节点的马尔可夫毯中,计算另外任意一个节点SNP2与已加入马尔可夫毯的SNP1节点在Class条件下的条件互信息I(SNP1,SNP2|Class),将条件互信息小于阈值的SNP2位点加入Class节点的马尔可夫毯中;对于非Class节点,在初始阶段假设表型Class节点已经加入到目标SNP1节点的马尔可夫毯中,计算任意一个节点SNP2与Class在目标节点SNP1条件下的条件互信息I(SNP2,Class|SNP1);将条件互信息小于阈值的SNP2位点加入目标节点SNP1的马尔可夫毯中;
步骤6.2收缩阶段:去除马尔科夫毯中多余节点;对于任意目标节点SNP1或者Class节点,计算在去除节点SNP2的马尔可夫毯的条件下目标节点SNP1与SNP2的条件互信息;若条件互信息大于阈值,则删除节点SNP2并更新节点SNP1的马尔可夫毯;
步骤6.3对称性检测:确保任意目标节点马尔可夫毯中的节点之间相互依赖;假设SNP1存在于SNP2的马尔可夫毯中,SNP1与SNP2相互依赖;与之对应,SNP2应该存在于SNP1的马尔可夫毯中;倘若两个节点的马尔可夫毯并不对称,则分别在每个节点的马尔可夫毯中删除另一节点;
步骤7,在构建得到不同k-clique内节点对应子网络的基础上,对各个子网络进行合并,得到包括SNP位点和表型性状节点网络结构;
步骤8,重复执行步骤1-步骤6,利用均匀抽样蒲公英代码的方法将网络图的构建转换为特定编码的生成问题,构建得到包括SNP位点和表型性状节点的网络结构;然后根据网络中边出现的次数,得到最终的包括SNP位点和表型性状节点网络结构,进而根据节点间的联系,得到影响表型性状的上位性位点。
2.根据权利要求1所述的基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:在步骤1中,将SNP基因型数据用0/1/2表示,假设AT为参考对象,进行如下表示:AA用0表示,TT用2表示,AT/TA用1表示,其中0表示纯合子常见基因型,1表示杂合子,2表示纯合子少见基因型;表型Class用0/1表示,0表示为未患病,1表示为患病。
3.根据权利要求1所述的基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:在步骤3.1中,根据除了与根节点相连的k-clique,与其它节点相连的k-clique不能重复的规则,将节点v与对应k-clique中节点进行合并组成新的节点,得到特征树的骨架。
4.根据权利要求1所述的基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:在步骤5中,对于图的每一条边,至少存在一个团包含边的两个顶点;假设t1,t2,t3是图中节点,t2在t1到t3的路径上,用Xt1,Xt2,Xt3表示不同的团;首先找到度最小的节点,根据树分解规则(节点v∈Xt1,v∈Xt3,可得v∈Xt2)对k-tree对应的图进行树分解,得到不同的团。
5.根据权利要求1所述的基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:在步骤6中,考虑到上位性位点挖掘问题的特殊性,将表型性状加入到每个团,进行子网络的构建;利用扩张、收缩、对称性检测三个阶段,通过条件互信息计算获取不同节点的马尔可夫毯,进而构建得到不同团对应的子网络结构。
6.根据权利要求1所述的基于k-tree优化贝叶斯网络的上位性位点挖掘方法,其特征在于:在步骤8中,通过均匀抽样蒲公英代码重复执行以上特征树构建,Rényi k-tree生成,k-tree生成,k-tree分解生成不同的团,子网络构建,子网络合并步骤,得到最终包括SNP位点和表型性状节点的网络结构,进而得到影响表型性状的上位性位点。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010683358.4A CN111833967B (zh) | 2020-07-10 | 2020-07-10 | 基于k-tree优化贝叶斯网络的上位性位点挖掘方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010683358.4A CN111833967B (zh) | 2020-07-10 | 2020-07-10 | 基于k-tree优化贝叶斯网络的上位性位点挖掘方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111833967A true CN111833967A (zh) | 2020-10-27 |
CN111833967B CN111833967B (zh) | 2022-05-20 |
Family
ID=72923360
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010683358.4A Active CN111833967B (zh) | 2020-07-10 | 2020-07-10 | 基于k-tree优化贝叶斯网络的上位性位点挖掘方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111833967B (zh) |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20070250463A1 (en) * | 2005-10-31 | 2007-10-25 | Tsutomu Sawada | Learning apparatus and method |
US20130024417A1 (en) * | 2011-07-18 | 2013-01-24 | Andreas Joanni | Method, system and computer program product for automatic generation of bayesian networks from system reliability models |
CN104946765A (zh) * | 2015-06-25 | 2015-09-30 | 华中农业大学 | 基于基因组测序的体细胞突变位点挖掘方法 |
CN105205344A (zh) * | 2015-05-18 | 2015-12-30 | 上海交通大学 | 基于多目标蚁群优化算法的基因位点挖掘方法 |
CN109411023A (zh) * | 2018-09-30 | 2019-03-01 | 华中农业大学 | 一种基于贝叶斯网络推理的基因间交互关系挖掘方法 |
CN109448794A (zh) * | 2018-10-31 | 2019-03-08 | 华中农业大学 | 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法 |
US20190341127A1 (en) * | 2018-05-03 | 2019-11-07 | The Chinese University Of Hong Kong | Size-tagged preferred ends and orientation-aware analysis for measuring properties of cell-free mixtures |
CN110570909A (zh) * | 2019-09-11 | 2019-12-13 | 华中农业大学 | 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法 |
-
2020
- 2020-07-10 CN CN202010683358.4A patent/CN111833967B/zh active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20070250463A1 (en) * | 2005-10-31 | 2007-10-25 | Tsutomu Sawada | Learning apparatus and method |
US20130024417A1 (en) * | 2011-07-18 | 2013-01-24 | Andreas Joanni | Method, system and computer program product for automatic generation of bayesian networks from system reliability models |
CN105205344A (zh) * | 2015-05-18 | 2015-12-30 | 上海交通大学 | 基于多目标蚁群优化算法的基因位点挖掘方法 |
CN104946765A (zh) * | 2015-06-25 | 2015-09-30 | 华中农业大学 | 基于基因组测序的体细胞突变位点挖掘方法 |
US20190341127A1 (en) * | 2018-05-03 | 2019-11-07 | The Chinese University Of Hong Kong | Size-tagged preferred ends and orientation-aware analysis for measuring properties of cell-free mixtures |
CN109411023A (zh) * | 2018-09-30 | 2019-03-01 | 华中农业大学 | 一种基于贝叶斯网络推理的基因间交互关系挖掘方法 |
CN109448794A (zh) * | 2018-10-31 | 2019-03-08 | 华中农业大学 | 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法 |
CN110570909A (zh) * | 2019-09-11 | 2019-12-13 | 华中农业大学 | 一种人工蜂群优化贝叶斯网络的上位性位点挖掘方法 |
Non-Patent Citations (2)
Title |
---|
FRANCESCO MERCATI等: "High-throughput 18k SNP array to assess genetic variability of the main grapevine cultivars from sicily", 《TREE GENETIC & GENOMES》 * |
张琪: "基于HBV病毒序列的突变位点挖掘与系统进化研究", 《中国优秀博硕士学位论文全文数据库(硕士)医药卫生科技辑》 * |
Also Published As
Publication number | Publication date |
---|---|
CN111833967B (zh) | 2022-05-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
KR102314219B1 (ko) | 심층 컨볼루션 신경망의 앙상블을 트레이닝하기 위한 반감독 학습 | |
Tang et al. | Mitochondrial phylogenomics of the Hymenoptera | |
Kautt et al. | Multispecies outcomes of sympatric speciation after admixture with the source population in two radiations of Nicaraguan crater lake cichlids | |
Richardson et al. | Statistical methods in integrative genomics | |
Yang et al. | CMDR based differential evolution identifies the epistatic interaction in genome-wide association studies | |
KR102628141B1 (ko) | 서열-특정 오류(sse)를 유발시키는 서열 패턴을 식별하기 위한 심층 학습-기반 프레임워크 | |
CN109448794B (zh) | 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法 | |
CN111833967B (zh) | 基于k-tree优化贝叶斯网络的上位性位点挖掘方法 | |
Masutani et al. | Investigating the mitochondrial genomic landscape of Arabidopsis thaliana by long-read sequencing | |
Gu et al. | Phylogeny and species delimitation of the genus Longgenacris and Fruhstorferiola viridifemorata species group (Orthoptera: Acrididae: Melanoplinae) based on molecular evidence | |
CN109493919B (zh) | 基于条件概率的基因型指派方法 | |
CN107058298A (zh) | 一种基于人工减数分裂的辅助基因组组装方法 | |
Sell | Addressing challenges of ancient DNA sequence data obtained with next generation methods | |
Kwarciak et al. | Tabu search algorithm for DNA sequencing by hybridization with multiplicity information available | |
Rastas et al. | Haplotype inference via hierarchical genotype parsing | |
Garg | Computational haplotyping: theory and practice | |
Reddy et al. | Untangling taxonomic confusion and diversification patterns of the Streak-breasted Scimitar Babblers (Timaliidae: Pomatorhinus ruficollis complex) in southern Asia | |
Delgado et al. | Viral Fitness Landscapes Based on Self-organizing Maps | |
Bucur | A stochastic de novo assembly algorithm for viral-sized genomes obtains correct genomes and builds consensus | |
Trujillo et al. | Getting higher on rugged landscapes: Inversion mutations open access to fitter adaptive peaks in NK fitness landscapes | |
Wu et al. | A practical algorithm based on particle swarm optimization for haplotype reconstruction | |
Płoński et al. | Quick path finding—Quick algorithmic solution for unambiguous labeling of phylogenetic tree nodes | |
Hossain et al. | An extension of heuristic algorithm for reconstructing multiple haplotypes with minimum error correction | |
Ebler | Design and application of methods for genome inference | |
Saha | Computational methods to study gene regulation in humans using DNA and RNA sequencing data |
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 |