CN114925837A - 基于混合熵优化互信息的基因调控网络构建方法 - Google Patents
基于混合熵优化互信息的基因调控网络构建方法 Download PDFInfo
- Publication number
- CN114925837A CN114925837A CN202210294214.9A CN202210294214A CN114925837A CN 114925837 A CN114925837 A CN 114925837A CN 202210294214 A CN202210294214 A CN 202210294214A CN 114925837 A CN114925837 A CN 114925837A
- Authority
- CN
- China
- Prior art keywords
- gene
- mutual information
- genes
- network
- entropy
- 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
- 108090000623 proteins and genes Proteins 0.000 title claims abstract description 233
- 230000033228 biological regulation Effects 0.000 title claims abstract description 83
- 238000005457 optimization Methods 0.000 title claims abstract description 35
- 238000010276 construction Methods 0.000 title claims abstract description 30
- 238000009826 distribution Methods 0.000 claims abstract description 62
- 239000011159 matrix material Substances 0.000 claims abstract description 40
- 230000008602 contraction Effects 0.000 claims abstract description 34
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 19
- 238000012216 screening Methods 0.000 claims abstract description 15
- 230000014509 gene expression Effects 0.000 claims abstract description 12
- 238000006243 chemical reaction Methods 0.000 claims abstract description 4
- 238000009795 derivation Methods 0.000 claims abstract description 4
- 238000000034 method Methods 0.000 claims description 96
- 230000001105 regulatory effect Effects 0.000 claims description 56
- 238000004364 calculation method Methods 0.000 claims description 24
- 230000008569 process Effects 0.000 claims description 12
- 239000013598 vector Substances 0.000 claims description 9
- 238000007476 Maximum Likelihood Methods 0.000 claims description 8
- 230000001186 cumulative effect Effects 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 6
- 230000003993 interaction Effects 0.000 claims description 6
- -1 Z)) is from P (X Inorganic materials 0.000 claims description 3
- 230000001364 causal effect Effects 0.000 claims description 2
- 230000002068 genetic effect Effects 0.000 claims description 2
- 238000009827 uniform distribution Methods 0.000 claims description 2
- 239000006185 dispersion Substances 0.000 abstract description 4
- 241000588724 Escherichia coli Species 0.000 description 28
- 238000002474 experimental method Methods 0.000 description 20
- 230000006870 function Effects 0.000 description 11
- 230000033616 DNA repair Effects 0.000 description 7
- 230000037361 pathway Effects 0.000 description 5
- 108091023040 Transcription factor Proteins 0.000 description 4
- 102000040945 Transcription factor Human genes 0.000 description 4
- 210000004027 cell Anatomy 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 108700026220 vif Genes Proteins 0.000 description 4
- 101100278439 Archaeoglobus fulgidus (strain ATCC 49558 / DSM 4304 / JCM 9628 / NBRC 100126 / VC-16) pol gene Proteins 0.000 description 3
- 210000003494 hepatocyte Anatomy 0.000 description 3
- 101150005648 polB gene Proteins 0.000 description 3
- 101150061752 ruvA gene Proteins 0.000 description 3
- 101150062749 uvrY gene Proteins 0.000 description 3
- 101100235354 Pseudomonas putida (strain ATCC 47054 / DSM 6125 / CFBP 8728 / NCIMB 11950 / KT2440) lexA1 gene Proteins 0.000 description 2
- 240000004808 Saccharomyces cerevisiae Species 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 2
- 101150047523 lexA gene Proteins 0.000 description 2
- 238000010801 machine learning Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 101150079601 recA gene Proteins 0.000 description 2
- 230000020509 sex determination Effects 0.000 description 2
- 101150005573 uvrA gene Proteins 0.000 description 2
- 101150073340 uvrD gene Proteins 0.000 description 2
- 101150106774 9 gene Proteins 0.000 description 1
- 238000001353 Chip-sequencing Methods 0.000 description 1
- 108020004414 DNA Proteins 0.000 description 1
- 241000282412 Homo Species 0.000 description 1
- 238000003559 RNA-seq method Methods 0.000 description 1
- 108700005075 Regulator Genes Proteins 0.000 description 1
- 230000027151 SOS response Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000004071 biological effect Effects 0.000 description 1
- 238000000205 computational method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000004034 genetic regulation Effects 0.000 description 1
- 210000002149 gonad Anatomy 0.000 description 1
- 230000002710 gonadal effect Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000002493 microarray Methods 0.000 description 1
- 230000000813 microbial effect Effects 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000002904 solvent Substances 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 101150046028 umuD gene Proteins 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N7/00—Computing arrangements based on specific mathematical models
- G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N5/00—Computing arrangements using knowledge-based models
- G06N5/04—Inference or reasoning models
- G06N5/041—Abduction
-
- 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
- Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
- Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
- Y04S10/00—Systems supporting electrical power generation, transmission or distribution
- Y04S10/50—Systems or methods supporting the power network operation or management, involving a certain degree of interaction with the load-side end user applications
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Artificial Intelligence (AREA)
- Probability & Statistics with Applications (AREA)
- Evolutionary Computation (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computational Linguistics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Operations Research (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于混合熵优化互信息的基因调控网络构建方法，该方法包括：根据设置的离散度对基因表达连续型数据进行离散化处理，根据真实概率与James‑Stein估计概率产生的均方误差MSE计算得到收缩强度λ；根据概率与熵值的转换公式，得到James‑Stein估计熵值；通过概率分布的β矩求导简化Dirichlet先验分布下的贝叶斯熵值估计，将两种熵值估计器得到的值转换为互信息矩阵；计算与互信息矩阵类似的Z‑score矩阵，将两个矩阵组合得到初始基因调控网络；根据路径一致算法进行遍历，通过动态阈值对基因调控网络中基因间关系进一步进行筛选，得到最终的基因调控网络。本发明解决了现有技术中构建网络存在大量错误调控关系的问题，得到了更准确的基因调控网络。
Description
技术领域
本发明涉及生物信息技术领域，尤其涉及一种基于混合熵优化互信息的基因调控网络构建方法。
背景技术
推断基因调控网络(GRN)，也称为逆向工程，是计算生物学中的一个关键问题。基因调控网络主要是指细胞内或特定基因组内基因之间相互作用形成的网络，尤其是由基因调控引起的相互作用。准确推断GRN对于理解生物活性和发现新的Pathway途径至关重要。随着ChIP-seq,RNA-seq等新一代测序技术(NGS)的快速发展，产生了大量的基因表达数据，这为基因调控关系的挖掘提供了前所未有机会。因此，如何在生物分子水平上构建基因调控网络是生物学研究中的一个复杂问题。
基于基因表达数据，通过逆向工程构建适当的数学或机器学习模型是推断基因调控网络的常用方法。迄今为止，研究者们已经开发了多种推断基因调控网络的计算方法，这些方法大致分为基于模型的方法和基于信息论的方法。其中基于模型的方法可以大致分为三类：统计方法、机器学习方法和概率图模型方法。互信息是构建基因调控网络最广泛使用的信息论方法。该方法不需要先验知识，主要通过计算(条件)互信息确定基因之间的调控关系。另外，互信息方法可以处理大量的基因，对基因间的非线性关系有较好的学习效果。目前常用的基于互信息的基因网络构建方法包括ARACNE，CLR， MRNET，NARROMI和CMI2等方法。然而，互信息方法在一定程度上高估了基因之间的直接调控关系。另外，互信息无法区分间接和直接的基因调控关系，导致较高的假阳性率。尽管条件互信息(CMI)可以区分直接和间接的调控关系，但该方法往往低估基因调控关系的强度，进而影响基因调控网络构建的准确性。
发明内容
本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种基于混合熵优化互信息的基因调控网络构建方法。
本发明解决其技术问题所采用的技术方案是：
本发明提供一种基于混合熵优化互信息的基因调控网络构建方法，该方法包括以下步骤：
S1，首先根据设置的离散度对基因表达连续型数据进行离散化处理，计算并得到每个基因对应的计数向量。根据真实概率与James-Stein估计概率产生的均方误差MSE计算得到收缩强度λ，以保证风险函数最小化。
S2，根据得到的收缩强度λ，计算无先验分布James-Stein估计的概率，根据概率与熵值的转换公式，得到James-Stein估计熵值。通过概率分布的β矩求导简化Dirichlet先验分布下的贝叶斯熵值估计，其中先验参数a通过收缩强度λ来计算，具体为当λ等于某一定值时可以使基于Dirichlet先验计算的概率等于贝叶斯估计下的概率，由此得到先验参数a与收缩强度λ的对应关系。根据熵值和互信息计算公式，将两种熵值估计器得到的值转换为互信息矩阵。
S3，在对两种熵值估计的互信息矩阵进行优化的基础上，计算每个基因的Z-score，进而计算每对基因的Z-score，得到与互信息矩阵类似的Z-score 矩阵，将两个矩阵组合得到初始基因调控网络。
S4，根据初始阈值筛选初始的基因调控网络中的调控关系，然后根据路径一致算法进行遍历，通过动态阈值对基因调控网络中基因间关系进一步进行筛选，得到最终的基因调控网络。
进一步，在所述S2步骤中，James-Stein收缩估计适用于高维数据集的基因网络推理计算。也就是说，James-Stein收缩估计对具有少量样本的数据集具有更好的计算效果。James-Stein收缩估计通过增加两个不同模型的权重以确保均方误差最小化。这两种模型分别为低偏差、高方差的高维模型和低偏差、高方差的低维模型。
进一步，在所述S2步骤中，基于Dirichlet先验分布的贝叶斯估计具体为 Beta分布的高维推广，在确定先验分布的参数a时，相当于向所有维度对应的单元格添加一个伪计数。为了得到先验参数a的值，本发明考虑了先验参数与收缩强度λ之间的关系。
进一步，在所述S3步骤中，标准化互信息矩阵具体为通过考虑基因间互信息的累积分布，进而消除网络中基因之间的一些间接调控关系。然后，根据每对基因的互信息，计算得到与互信息的经验分布相关的Z-score。Z-score 矩阵将替代互信息矩阵作为基因调控网络的初始权重矩阵，便于下一步去除网络中冗余边的操作。
进一步，在所述S4步骤中，在通过路径一致性算法遍历整个基因调控网络的过程中，计算不同基因对在多个基因影响下的条件互包含信息(CMI2)。根据计算得到的动态阈值更新当前阈值，并删除低于该阈值的基因间关系，判断并更新整个邻接矩阵，直到不删除任何关系为止。通过上述方式逐渐删除冗余边，进而形成一个更准确的基因调控网络。
进一步，在所述S4步骤中，动态阈值筛选过程具体为使用动态设置阈值方法来筛选基因之间的冗余调控关系。阈值越大，对基因调控关系的控制就越严格。首先根据设置的初始阈值初步筛选并删除冗余关系，同时设置order 值，即可以从与基因对相关联的所有基因中抽取的最大值，从1到order随机选择当前基因对下相关联的基因，作为间接调控影响来计算条件互包含信息 (CMI2)。随机选择产生的组合数越多，该基因对受到间接调控的可能性就越大，相应阈值也应该增加，因此通过设置动态阈值实现实时删除图中的冗余边。
与现有技术相比，本发明产生的有益效果是：本发明提供一种基于混合熵优化互信息的基因调控网络构建方法，通过基因表达数据构建互信息矩阵，对互信息矩阵和初始基因调控网络优化，进而构建准确性更高的基因调控网络的方法。首先将无先验分布的James-Stein熵估计和贝叶斯熵估计与 Dirichlet先验分布相结合，计算基因之间的互信息以获得初始基因调控网络。这有利于获得更多正确的基因间调控关系。为了考虑上下文信息对基因间调控关系的影响，进一步使用相关性算法的上下文似然优化互信息矩阵分布。通过以上两步，可以得到一个基因调控关系较准确的初始基因调控网络。由于互信息不能有效处理基因间的间接调控关系，构建的基因网络中仍包含大量错误的调控关系，假阳性率较高。为了解决此问题，主要通过考虑互信息的累积分布来消除基因之间的一些间接调控关系。为了删除冗余和不正确的基因间调控关系，通过计算多个基因影响下基因之间的条件互包含信息进一步优化整个网络。经过上述多次迭代处理后，最初的基因调控网络逐渐优化，进而得到更准确的基因调控网络。
附图说明
下面将结合附图及实施例对本发明作进一步说明，附图中：
图1为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的两个真实大肠杆菌数据集的标准网络。
图2为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法MEOMI的框架。
图3为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的去除冗余关系的过程。
图4为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法DREAM3 Size50数据集上的AUPR和AUROC比较结果。
图5为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的CMI2NI和MEOMI在DREAM3 Size50数据集的TPR比较结果。
图6为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的DREAM3 Size100数据集的AUPR和AUROC比较结果。
图7为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法CMI2NI,NARROMI和MEOMI在DREAM3 Size100数据集的TPR 比较结果。
图8为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法在DREAM5数据集中9种方法的F1-score,AUPR和AUC比较。
图9为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的SOS数据集的AUPR和AUROC比较结果。
图10为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的大肠杆菌SOS通路网络(Data1)的预测网络。
图11为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的大肠杆菌SOS DNA修复网络的预测网络。
图12为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的大肠杆菌SOS DNA修复网络在不同数据集中8个基因的分布。
图13为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的9种不同方法在GSD(drop rate＝0)和GSD(drop rate＝50)数据集的 TPR,PPV,F1-score,FPR,ACC,MCC,AUPR和AUC比较结果。
图14为本发明实施例提供的一种基于混合熵优化互信息的基因调控网络构建方法的hHEP(Top-100)、hHEP(Top-200)和hHEP(Top-500)数据集的TPR, PPV,F1-score,FPR,ACC,MCC,AUPR和AUC的比较结果。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
请参阅图1和图2，本发明实施例提供一种基于混合熵优化互信息的基因调控网络构建方法，包括如下步骤：S1，首先根据设置的离散度对基因表达连续型数据进行离散化处理，计算得到每个基因对应的计数向量。根据真实概率与James-Stein估计概率产生的均方误差MSE计算得到收缩强度λ，以保证风险函数最小化。S2，根据得到的收缩强度λ，计算无先验分布James-Stein 估计的概率，根据概率与熵值的转换公式，得到James-Stein估计熵值。通过概率分布的β矩求导简化Dirichlet先验分布下的贝叶斯熵值估计，其中先验参数a取决于收缩强度λ的值，其中先验参数a通过收缩强度λ来计算，具体为当λ等于某一定值时可以使基于Dirichlet先验计算的概率等于贝叶斯估计下的概率，由此得到先验参数a与收缩强度λ的对应关系。根据熵值和互信息计算公式，将两种熵值估计器得到的值转换为互信息矩阵。S3，对两种熵值估计得到的互信息矩阵进行优化，计算每个基因的Z-score，然后计算每对基因的Z-score，得到与互信息矩阵类似的Z-score矩阵，将两个矩阵组合得到初始基因调控网络。S4，根据初始阈值筛选初始基因调控网络的调控关系，然后根据路径一致算法进行遍历，通过动态调整阈值的方法对基因调控网络中的关系进一步进行筛选，得到最终的基因调控网络。
在本实施例中，为了计算基因之间的互信息，首先对连续基因表达数据集进行离散化。如果随机变量X的值分布在区间[a,b]中，则根据区间的大小将该区间划分为等距的子区间。子区间的数量表示为bin，离散化后的随机变量X如Eq.(1)所示。
X＝[X1,X2,X3,......,Xn] (1)
经过离散化操作后，随机变量X的n个变量Xi分布在K个bin中，其中 K表示分布概率大于0的bin数量。每个bin对应的索引向量为xi，随机变量 X对应的索引向量如Eq.(2)所示。
χ＝[x1,x2,x3,......,xK] (2)
此外，收缩目标t没有方差，但有更高的偏差。确定最佳收缩强度λ的第一步是选择合适的损失函数，在本实施例中，使用平方误差作为损失函数。第二步是最小化风险函数R(λ)，在本实施例中，使用均方误差进行计算，如 Eq.(3)所示。
然后得到最小化MSE的收缩强度λ，如Eq.(4)所示。
给定
和
可以得到Eq.(5)。
为了避免过度收缩，当
令
当出现负收缩时，令
因此可以保证
作为本发明实施例的优化方案，请参阅图1和图2，在所述S2步骤中， James-Stein收缩估计适用于高维数据集的基因调控网络推理计算。也就是说， James-Stein收缩估计对具有少量样本的数据集具有较好的计算效果。 James-Stein收缩估计通过增加两个不同模型的权重来确保均方误差最小化。这两种模型分别为低偏差、高方差的高维模型和低偏差、高方差的低维模型。在本实施例中，基于James-Stein估计将ML估计值降低到tk，并通过最小化均方误差(MSE)计算收缩强度λ。这有助于提高计算精度并获得比其他单一估计器更好的性能，如Eq.(6)所示。
在Eq.(6)中，λ为收缩强度，其值在0～1之间。tk为对应的收缩目标。一般，
对应于最大熵目标的均匀分布。此外，收缩的目标tk没有方差，但有较高的偏差。
最大似然估计(ML)是一个典型的熵估计方法，如Eq.(7)所示。
每个bin的计数n(xi)和频率θ(xi)之间关联的多项式分布如Eq.(8)所示。
通过取Eq.(8)右侧似然函数的最大值，对参数θ(xi)进行ML估计。得到θ(xi)的最大似然估计值，记为
进而可以证明方差
偏差
和
ML估计是一个具有高方差的无偏模型。
最大似然估计的熵计算是基于渐近理论，当样本量较小时，bin中的计数更可能为0。进而最大似然估计的概率也会为0，这将影响计算的精度。
作为本发明实施例的优化方案，请参阅图1和图2，在所述S2步骤中，基于Dirichlet先验分布的贝叶斯估计为Beta分布的高维推广，在确定先验分布的参数a时，相当于向所有维度对应的单元格添加一个伪计数。为了得到先验参数a的值，本实施例考虑了先验参数与收缩强度λ之间的关系。在本实施例中，Dirichlet分布也被称为多元Beta分布，可以用来估计离散数据的熵，它是Beta分布的高维推广。参数为a1～aK，维数为K的狄利克雷先验分布的密度形式用Eq.(9)表示。
其中ai是事件θi的先验参数，Γ(·)是伽玛函数，δ的定义如Eq.(10)所示。
如果没有先验知识，先验参数ai的取值默认是相等的，表示每件事情在未知情况下发生的概率是相等的，故得到a1＝a2＝...＝aK＝a，且
基于狄利克雷先验分布的贝叶斯估计如Eq.(11)所示。
一般都是根据Eq.(11)中估计的概率θ直接估计熵值。然而，用该方法得到的熵估计的精度是未知的。因此，本实施例直接从现有的数据中估计熵。根据估计概率θ计算的熵
不等于直接基于索引向量计算的估计熵
为了快速有效地得到估计的熵，本实施例首先得到概率分布的第β个矩的方程，如Eq.(12)所示。
其中β对应矩，xj表示每个bin的计数，K为大于0的bin数。根据
得到Eq.(13)。
其中
是普西函数。可以得到，该计算范式取决于离散数据计数向量以及设置的先验参数，计算过程简单且易实现。其中最为关键的是确定先验分布的参数a，相当于给K个单元格增加一个伪计数。为了得到先验参数a的值，本实施例考虑了先验参数与收缩强度λ之间的关系。
通过设置
和
可以得到由James-Stein收缩估计的收缩强度λ与Dirichlet先验参数a之间的相应关系，如Eq.(14)所示。
应用于熵值估计的ai有多种选择情况。当ai＝0时，即无先验知识，得到的熵值估计为经典的ML估计解。当
得到的估计值为最大最小似然估计；当
这些值都是在未知先验下的合理估计。
作为本发明实施例的优化方案，请参阅图1和图2，在所述S3步骤中，标准化互信息矩阵具体为通过考虑基因间互信息的累积分布，进而消除网络中基因之间的一些间接调控关系。然后，根据每对基因的互信息，计算得到一个与互信息的经验分布相关的Z-score。Z-score矩阵将替代互信息矩阵作为基因调控网络的初始权重矩阵，以便下一步去除冗余关系处理。在本实施例中，为了获得更准确的初始基因调控网络，对得到的互信息矩阵MIM进一步进行处理。主要通过考虑基因间互信息的累积分布，消除基因调控网络中基因之间的一些间接调控关系。Eq.(15)定义了单个基因Gi的Z-score Zi。
在Eq.(15)中，I(Xi；Xj)表示基因对Gi和Gj间的互信息。任意两个基因Gi和Gj对应的随机变量分别记为Xi和Xj。将MIM中所有与Gi相关的Gj的互信息形成经验分布，标准化处理该分布后得到Zi。在Eq.(15)中，μi和σi分别表示均值和标准差。为了获得单个基因的Zi，计算基因Gi和Gj的似然互信息得分，如Eq.(16)所示。
最后，根据每对基因的互信息得到一个与互信息的经验分布相关的Zij分数矩阵。Zij分数矩阵将MIM作为基因调控网络的初始权重矩阵，以便在下一步中进行去除冗余关系处理。
作为本发明实施例的优化方案，请参阅图1和图2，在所述S4步骤中，根据路径一致算法进行遍历具体为通过路径一致性算法遍历整个基因调控网络。该算法计算了基因对在多个基因影响下的条件互包含信息(CMI2)。根据计算得到的动态阈值更新当前阈值，并删除低于该阈值的关系。判断并更新整个邻接矩阵，直到不再删除任何边为止。通过上述方法，逐渐去除冗余边，进而形成一个更准确的基因调控网络。在本实施例中，由于互信息方法不能有效地处理基因间的间接调控关系，导致所构建的基因调控网络中的假阳性率较高。为了解决这一问题，本实施例进一步使用条件互包含信息计算方法逐步去除冗余关系。具体来说，通过路径一致性算法遍历整个基因调控网络。通过计算任一基因对在多个基因影响下的条件互包含信息(CMI2)。逐渐去除冗余关系，进而形成一个更准确的基因调控网络。具体的实现过程如图2所示。
如图2所示，对于基因G1，在初始基因调控网络中，G1和G2之间存在一条边。然后搜索与G1和G2具有连接边的基因Gi。假设满足上述条件的基因Gi的数目为Order_Max，即与G1和G2相互作用的最大基因数为Order_Max。如果Order＝0，则根据初始阈值λ0获得初始基因调控网络，然后逐渐增加Order 以增量更新基因调控网络。
具体过程为，从Order＝1开始，计算所有Gi的条件交互信息CMI2(G1,G2 |Gi)，返回最大值为MAX(CMI2(G1,G2|Gi))。如果MAX(CMI2(G1,G2|Gi))小于阈值λ1，则删除G1和G2之间的边，并将基因邻接矩阵中G1(1,2)的值设置为0。如果MAX(CMI2(G1,G2|Gi))大于阈值λ1，增加Order的值。计算Order 基因的CMI2(G1,G2|Gi)，并按照上述步骤依次将计算结果与阈值λi进行比较，直到Order＝Order_Max或Order达到Order_Set的最大值。按照上述步骤遍历所有基因，逐渐删除不正确的基因间关系，从而形成较准确的基因调控网络。
对于给定Z的随机变量X和Y，CMI2的计算方式如Eq.(17)所示。
在Eq.(17)中，DKL(P(X,Y,Z)||PX→Y(X,Y,Z))是从P(X,Y,Z)到PX→Y(x,y,z)的 KL散度，用于测量两个分布之间的距离。P(X,Y,Z)表示X,Y和Z之间的联合概率。PX→Y(x,y,z)为X到Y的干扰概率。X和Y之间的因果变量强度见Eq. (18)。
CX→Y(X；Y|Z)＝DKL(P(X,Y,Z)||PX→Y(X,Y,Z)) (18)
CX→Y(X；Y|Z)是不对称的，CMI2(X；Y|Z)是CX→Y(X；Y|Z)和CY→X(Y；X|Z)的平均值。在扩展Eq.(18)后，可以得到Eq.(19)。
DKL(P(X,Y,Z)||PY→X(X,Y,Z))＝I(X；Y|Z)+DKL(P(X|Z)||PY→X(X|Z)) (19)
通过结合Eq.(20)和Eq.(17)，可以得到CMI2(X；Y|Z)，如Eq.(20)所示。
其中PX→Y(Y|Z)的计算方法如Eq.(21)所示。
为了估计Eq.(21)中的密度函数，本实施例使用高斯核密度估计方法计算相应的密度函数。当使用高斯核密度法估计概率密度函数时，主要应用基因表达数据受多元高斯分布影响的广泛假设进行计算，这有助于更有效的估计 CMI2(X；Y|Z)。
作为本发明实施例的优化方案，请参阅图1和图2，在所述S4步骤中，动态阈值筛选具体为使用动态设置阈值方法来筛选基因之间的冗余调控关系。阈值越大，对基因调控关系的控制就越严格。首先根据设置的初始阈值初步筛选并删除冗余边，同时设置order值，即可以从与基因对相关联的所有基因中抽取的最大值。从1到order随机选择当前基因对下相关联的基因，作为间接调控影响来计算条件互包含信息(CMI2)。随机选择产生的组合数越多，该基因对受到间接调控的可能性就越大，同时增加相应阈值，因此通过设置动态阈值实现实时删除冗余边。在本该阈值对实验结果有很大的影响，本实施例使用动态设置阈值的方法来筛选基因之间的冗余调控关系。阈值越大，对基因调控关系的控制就越严格。考虑到多种间接调控关系对基因间直接调控关系的影响，在计算各阶CMI2(Gi,Gj|Gk)时，同时考虑了从Order_Max中提取的阶基因的初始阈值α和组合数m。显然，m依赖于Order的值和所有相关基因Order_Max的数量。m的值越大，受到的间接调控就越大，从而增加相应的阈值。
在另外的实施例中，通过实验来说明基于混合熵优化互信息的基因调控网络构建方法的有效性，分别比较了DREAM challenge模拟数据集(DREAM3 和DREAM5)、3个真实的大肠杆菌数据集(大肠杆菌SOS通路网络、大肠杆菌SOS DNA修复网络和大肠杆菌Community网络)和2个人类数据集(GSD 和hHEP)的多种指标。
为了对不同方法的性能进行比较，使用真阳性率(TPR)，假阳性率(FPR)，阳性预测值(PPV)，整体准确率(ACC)，F1-score和Matthews相关系数(MCC) 进行评价，如Eq.(22)～Eq.(27)所示。
TPR＝TP/(TP+FN) (22)
FPR＝FP/(TN+FP) (23)
PPV＝TP/(TP+FP) (24)
ACC＝(TP+TN)/(TP+FP+TN+FN) (25)
F1-score＝2PPV×TPR/(PPV+TPR) (26)
其中TP是正确识别的边数，TN是正确识别的非标准网络边数，FP是错误识别的边数，FN是错误识别的非标准网络边数。通过从大到小设置不同的阈值，计算得到一系列TP,FP,TN和FN，进而计算得到绘制接收机工作特性 (ROC)曲线的TPR和FPR值和用于绘制PR曲线的TPR和PPV。计算ROC曲线下的面积(AUC)作为比较不同算法的另一个度量指标。PR(Precision-Recall) 曲线可以综合考虑TP和FP，对于极度不平衡的数据PR比ROC更实用，因此将其作为衡量不同算法性能的一个度量指标。
将本发明中基于混合熵优化互信息的基因调控网络构建方法记为 MEOMI，比较了MEOMI与目前一些有代表性的基因网络构建方法，包括 GENIE3,CLR,ARECNE,MRNET,CMI2NI,NARROMI和BiXGBoost。
实验1.DREAM3 Challenge Size50数据集的实验结果
DREAM challenge包括了一系列从真实物种源网络中选取的带有噪声和基准网络的基因表达数据。在本实验中，从DREAM3 challenge数据集中选择了10个Size50的数据集进行比较。这些数据集中有50个基因和50个样本，它们包括5组数据：Ecoli1,Ecoli2,Ecoli3,Yeast1和Yeast2。每组包括两种数据类型：heterozygous和null-mutants。对于每种数据类型，使用5个数据集进行实验比较。Size50数据集的详细信息如表1所示。
表1 DREAM3中Size50数据集信息
MEOMI方法具有三个参数。第一个参数是bins，表示数据集的分散度。根据初步实验结果，参数bins一般设置为5。第二个参数order是控制基因对直接调控的多基因数量。当order达到某个阈值时，该参数终止算法。order 的取值对算法的运行时间影响很大，设置较大的order会导致MEOMI的效率低下。一般来说，设置order＝2或order＝3。第三个参数λ是用于筛选基因调控网络中边的初始阈值λ0。根据初步实验结果，λ的范围一般设置为0.000001～0.01。此外，实验结果显示λ的值主要与数据集中的基因数量有关。基因数量越多，λ值越大。对于Size50数据集，在MEOMI算法中设置λ＝0.000001。图3显示了9种不同方法的AUPR和AUROC比较结果。
在图3(B)中，MEOMI方法的AUROC在10个数据集中明显高于其他8 种方法。在Data1和Data10中MEOMI方法的AUC略低于GENIE3方法的 AUC，但这两种方法之间的差异非常小。图3(A)显示CMI2NI和MEOMI方法的AUPR值比较高，总体上高于其他7种方法。除了少数数据集外，MEOMI 方法的AUPR与CMI2NI方法相似，均优于其他方法的AUPR。尽管CMI2NI 在AUPR上的表现优于MEOMI，但其在10个数据集上的TPR均低于MEOMI 的TPR。如图4所示。尽管CMI2NI的AUPR比MEOMI高，但CMI2NI丢失了更多正确边。总体而言，Size50数据集中MEOMI的AUPR和AUROC显着高于其他方法。
实验2.DREAM3 Challenge Size100数据集的实验结果
在本实验中，从DREAM3 challenge数据集中选择了10个Size100的数据集进行比较。这些数据集中有100个基因和100个样本，包括5组数据：Ecoli1,Ecoli2,Ecoli3,Yeast1和Yeast2。每组包括两种数据类型：heterozygous 和null-mutants。对于每种数据类型，使用5个数据集进行实验比较，采用 Size100数据集的详细信息如表2所示。对于Size100数据集，在MEOMI中，λ同样设置为0.000001。图5显示了9种不同方法的AUPR和AUROC比较结果。
表2 DREAM3中Size100数据集信息
图5(B)显示MEOMI方法的AUROC在不同数据集上优于其他方法，整体结果趋势与GENIE3-RF-sqrt大致相同，在部分heterozygous数据上的表现略低于GENIE3-RF-sqrt。在5个heterozygous数据集中，所有方法的AUC均小于0.6，表明该实验数据的特征不是很明显。此外，可以看到GENIE3-RF-sqrt 方法的结果优于GENIE3-RF-all方法，这表示随机选择特征对heterozygous数据集的实验结果有更大的影响。在图5(A)中，MEOMI的AUPR在某些数据集上略低于CMI2NI和NARROMI。但MEOMI的AUPR优于其他6种方法。虽然CMI2NI和NARROMI在AUPR方面的表现优于MEOMI，但两种方法在 10个数据集上对应的TPR低于MEOMI。如图6所示，尽管CMI2NI和 NARROMI的AUPR比MEOMI高，但这两种方法丢失了更多正确的边。总体而言，MEOMI的AUPR和AUROC在Size100数据集上优于其他8种方法。
实验3.DREAM5 Challenge Size1643数据集的实验结果
DREAM5数据集包括1643个基因、805个样本和195个转录因子。考虑到标准网络中仅涉及转录因子，仅使用与195个转录因子相关的调控关系进行实验比较。DREAM5数据集的实验结果如图7所示，如上所述，参数λ主要与数据集中的基因数量有关。基因数量越多，λ值越大。对于包含比DREAM3 数据集更多数目基因的DREAM5数据集，将λ设置为0.001。
图7显示MEOMI的F1-score在DREAM5数据集上低于NARROMI, CMI2NI和ARACNE，但高于其他方法。NARROMI,CMI2NI和ARACNE都使用互信息进行计算以构建基因调控网络，但这些方法获得的正确边太少。这三种方法的TPR分别为0.2459,0.2985和0.1579，远低于MEOMI方法的 TPR(0.9424)。此外，MEOMI的AUC和AUPR均高于其他8种方法。可见 MEOMI在处理大量基因时具有较好的学习准确度。
实验4.大肠杆菌SOS网络的实验结果
大肠杆菌(E.coli)SOS网络被广泛用于评估基因调控网络构建算法的优劣。本发明使用三个真实的大肠杆菌(E.coli)数据集对9种基因网络构建方法进行了比较和分析。第一个数据集包括大肠杆菌SOS通路中的9个基因、24 条边和9个样本。标准网络中有9个基因(G1～G9)，包括参与SOS反应的主要介质、已知调节基因和sigma因子基因，标准网络如图8(A)所示。第二个大肠杆菌SOS数据集包括参与DNA修复过程的8个基因(uvrD,lexA,umuD,recA,uvrA,uvrY,ruvA和polB)和50个样本，标准网络如图8(B)所示。本发明主要使用不考虑方向和自环的标准网络进行实验比较，如图8(C)所示。该数据集是通过在实验中设置一定的时间步长(步长为50)，在不同紫外光强度下对基因间隔进行采样而生成的，总共进行了4次实验，数据集的预处理保留了第一个时间点(0)。三个真实大肠杆菌数据集的信息如表3所示。
表3三种大肠杆菌真实数据集信息
初步实验结果表明，对于基因和样本数量较少的数据集，往往设置较大的离散度(bins)。将大肠杆菌SOS DNA通路网络和大肠杆菌SOS DNA修复网络数据集分别设置为8和10。此外，两个数据集的参数λ分别设置为0.016 和0.003。本实验比较了不同方法在大肠杆菌SOS数据集的TP,FP,TN,FN, TPR,PPV,F1-score,FPR,ACC,MCC,AUPR和AUC。图9显示了9种方法 (GENIE3-RF-sqrt,GENIE3-RF-all,CLR,ARACNE,MRNET,CMI2NI,NARROMI, BiXGBoost和MEOMI)的AUPR和AUROC比较结果。
如图9所示，对于大肠杆菌SOS通路网络数据集Data1，MEOMI的AUC 和AUPR优于其他方法。对于大肠杆菌SOS DNA修复网络数据集Data2，除 Data2-2数据集外，MEOMI的AUC和AUPR优于其他方法。具体来讲，MEOMI 在Data2-2的准确率低于ARACNE,CLR和MRNET，MEOMI与CMI2NI的准确率相差不大。利用MEOMI预测的大肠杆菌SOS网络(Data1和Data2)如图 10和图11所示。实线表示FP(错误预测的边数)，虚线表示TP(正确预测的边数)。可见MEOMI预测的网络具有比较高的准确率。
为了进一步分析MEOMI在Data2-2数据集中学习效果不好的原因，绘制了SOS DNA修复网络的4个数据集中每个基因(polB,ruvA,uvrY,uvrA,recA, umuDC,lexA和uvrD)对应的数据分布，结果如图12所示。
在图12中，4个数据集分布差异最大的三个基因分别是uvrY,ruvA和polB。进一步分析MEOMI构建的基因网络后，发现与这三个基因相关的边的得分在标准网络所有边中是最低的，如表4所示。因此，可见Data2-2中MEOMI方法的AUC和AUPR性能可能与数据分布有关。
表4 Data2-2的标准网络中基因对的预测得分
实验5.大肠杆菌Community网络的实验结果
大肠杆菌Community网络数据集是从许多微生物微阵列(M3D)数据库 (version4,build 6)中得到的。大肠杆菌Community网络数据集包含4297个基因的907个芯片的测量值，包括175个转录因子和4564个调控关系。参数设置为λ＝0.0001，order＝4。比较了大肠杆菌Community网络中不同方法的TP,FP, TN,FN,TPR,PPV,F1-score,FPR,ACC,MCC,AUPR和AUC。部分实验结果如表5所示。
如表5所示，MEOMI的AUC和AUPR仅略低于CLR和MRNET，但高于其他6种方法。
表5大肠杆菌Community网络的实验结果
MEOMI,CLR和MRNET方法主要基于互信息理论构建基因调控网络。虽然MEOMI在某些参数上略逊于CLR和MRNET，但MEOMI在保证TP变化不大的同时可以大大降低FP。表明MEOMI可以删除更多的网络中假阳性边，有效地减少了错误的基因间接调控关系的数量。此外，可见MEOMI的 F1-score低于ARACNE和NARROMI，并且与其他方法相差不大。然而ARACNE和NARROMI的TP远小于MEOMI，表明这两种方法无法获得较多数目正确的边。总体而言，MEOMI在处理具有大量基因的真实大肠杆菌数据集中仍然具有优势。
实验6.性腺性别测定(GSD)数据集实验结果
本发明使用两个人类数据集进行了实验验证：性腺性别测定(GSD)和人类成熟肝细胞(hHEP)。GSD数据集包括19个基因、79边和2000个样本。采用 dropout率q＝0和q＝50的GSD数据集进行实验比较。图13显示了9种不同方法在GSD(drop rate＝0)和GSD(drop rate＝50)数据集的TPR,PPV,F1-score, FPR,ACC,MCC,AUPR和AUC比较结果。
需要说明的是，ARACNE和NARROMI的TPR明显远低于其他7种方法，表明ARACNE和NARROMI无法获得更多的正确边。因此，在本实验中不对这两种方法进行比较。如图13所示，MEOMI的性能总体上优于 GENIE3-RF-sqrt,GENIE3-RF-all,CMI2NI和BiXGBoost。虽然MEOMI的TPR 低于CLR和MRNET，但两种方法的FPR明显高于MEOMI。
实验7.人成熟肝细胞(hHEP)数据集实验结果
人类成熟肝细胞(hHEP)数据集包括11515个基因和425个样本。分别选取与9,17,33个TF基因相关的前100,200,500个变异最大的基因进行实验比较。主要使用从STRING数据库中获得的标准网络进行验证和分析。图14 显示了hHEP(Top-100)、hHEP(Top-200)和hHEP(Top-500)数据集的TPR,PPV, F1-score,FPR,ACC,MCC,AUPR和AUC的比较结果。同理，ARACNE和 NARROMI的TPR明显低于其他7种方法，本实验不对这两种方法的结果进行比较。
通过实验结果可得MEOMI和CMI2NI的性能总体上优于其他5种方法。虽然MEOMI在hHEP(前500)数据集上的TPR低于MRNET，但MRNET的FPR 明显低于MEOMI。在三个数据集中，CLR和MRNET的FPR均高于其他5种算法。总之，可见MEOMI在处理人类较多数目的基因时比其他算法具有更好的学习效果。
应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。
Claims (10)
1.一种基于混合熵优化互信息的基因调控网络构建方法，其特征在于，该方法包括以下步骤：
S1、获取基因表达连续型数据，根据设置的离散度对基因表达连续型数据进行离散化处理，计算并得到每个基因对应的计数向量；根据真实概率与James-Stein估计概率产生的均方误差MSE计算得到收缩强度λ；
S2、根据得到的收缩强度λ，计算无先验分布James-Stein估计的概率，根据概率与熵值的转换公式，得到James-Stein估计熵值；通过概率分布的β矩求导简化Dirichlet先验分布下的贝叶斯熵值估计，其中先验参数a通过收缩强度λ来计算，具体为当λ等于某一定值时使基于Dirichlet先验计算的概率等于贝叶斯估计下的概率，由此得到先验参数a与收缩强度λ的对应关系；根据熵值和互信息计算公式，将两种熵值估计器得到的值转换为互信息矩阵；
S3、在对两种熵值估计的互信息矩阵进行优化的基础上，计算每个基因的Z-score，进而计算每对基因的Z-score，得到与互信息矩阵类似的Z-score矩阵，将两个矩阵组合得到初始基因调控网络；
S4、根据初始阈值筛选初始的基因调控网络中的调控关系，然后根据路径一致算法进行遍历，通过动态阈值对基因调控网络中基因间关系进一步进行筛选，得到最终的基因调控网络。
2.根据权利要求1所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S1步骤中：
首先对连续基因表达数据集进行离散化；如果随机变量X的值分布在区间[a,b]中，则根据区间的大小将该区间划分为等距的子区间；子区间的数量表示为bin，离散化后的随机变量X如Eq.(1)所示：
X＝[X1,X2,X3,......,Xn] (1)
经过离散化操作后，随机变量X的n个变量Xi分布在K个bin中，其中K表示分布概率大于0的bin数量；每个bin对应的索引向量为xi，随机变量X对应的索引向量如Eq.(2)所示：
χ＝[x1,x2,x3,......,xK] (2)
此外，收缩目标t没有方差，但有更高的偏差；确定最佳收缩强度λ的第一步是选择合适的损失函数，使用平方误差作为损失函数；第二步是最小化风险函数R(λ)，使用均方误差进行计算，如Eq.(3)所示：
然后得到最小化MSE的收缩强度λ，如Eq.(4)所示：
给定
和
得到Eq.(5)：
为了避免过度收缩，当
令
当出现负收缩时，令
保证
3.根据权利要求1所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S2步骤中：
James-Stein收缩估计适用于高维数据集的基因网络推理计算，也就是说，James-Stein收缩估计对具有少量样本的数据集具有更好的计算效果，James-Stein收缩估计通过增加两个不同模型的权重以确保均方误差最小化，这两种模型分别为低偏差、高方差的高维模型和低偏差、高方差的低维模型；
基于Dirichlet先验分布的贝叶斯估计具体为Beta分布的高维推广，在确定先验分布的参数a时，相当于向所有维度对应的单元格添加一个伪计数，为了得到先验参数a的值，考虑了先验参数与收缩强度λ之间的关系。
4.根据权利要求3所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S2步骤中James-Stein收缩估计的具体方法为：
基于James-Stein估计将ML估计值降低到tk，并通过最小化均方误差MSE计算收缩强度λ；如Eq.(6)所示：
在Eq.(6)中，λ为收缩强度，其值在0～1之间；tk为对应的收缩目标；
对应于最大熵目标的均匀分布；
最大似然估计ML如Eq.(7)所示：
每个bin的计数n(xi)和频率θ(xi)之间关联的多项式分布如Eq.(8)所示：
通过取Eq.(8)右侧似然函数的最大值，对参数θ(xi)进行ML估计；得到θ(xi)的最大似然估计值，记为
进而证明方差
偏差
和
ML估计是一个具有高方差的无偏模型。
5.根据权利要求4所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S2步骤中基于Dirichlet先验分布的贝叶斯估计具体为Beta分布的高维推广的具体方法为：
Dirichlet分布用来估计离散数据的熵，它是Beta分布的高维推广；参数为a1～aK，维数为K的狄利克雷先验分布的密度形式用Eq.(9)表示：
其中ai是事件θi的先验参数，Γ(·)是伽玛函数，δ的定义如Eq.(10)所示：
如果没有先验知识，先验参数ai的取值默认是相等的，表示每件事情在未知情况下发生的概率是相等的，故得到a1＝a2＝...＝aK＝a，且
基于狄利克雷先验分布的贝叶斯估计如Eq.(11)所示：
首先得到概率分布的第β个矩的方程，如Eq.(12)所示：
其中β对应矩，xj表示每个bin的计数，K为大于0的bin数；根据
得到Eq.(13)：
其中
是普西函数；
通过设置
和
得到由James-Stein收缩估计的收缩强度λ与Dirichlet先验参数a之间的相应关系，如Eq.(14)所示：
应用于熵值估计的ai有多种选择情况；当ai＝0时，即无先验知识，得到的熵值估计为经典的ML估计解；当
得到的估计值为最大最小似然估计；当
这些值都是在未知先验下的合理估计。
6.根据权利要求1所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S3步骤中：
标准化互信息矩阵具体为通过考虑基因间互信息的累积分布，进而消除网络中基因之间的一些间接调控关系；然后，根据每对基因的互信息，计算得到与互信息的经验分布相关的Z-score；Z-score矩阵将替代互信息矩阵作为基因调控网络的初始权重矩阵，便于下一步去除网络中冗余边的操作。
7.根据权利要求6所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S3步骤中：
通过考虑基因间互信息的累积分布，消除基因调控网络中基因之间的一些间接调控关系；Eq.(15)定义了单个基因Gi的Z-score Zi：
在Eq.(15)中，I(Xi；Xj)表示基因对Gi和Gj间的互信息；任意两个基因Gi和Gj对应的随机变量分别记为Xi和Xj；将MIM中所有与Gi相关的Gj的互信息形成经验分布，标准化处理该分布后得到Zi；μi和σi分别表示均值和标准差；为了获得单个基因的Zi，计算基因Gi和Gj的似然互信息得分，如Eq.(16)所示：
最后，根据每对基因的互信息得到一个与互信息的经验分布相关的Zij分数矩阵；Zij分数矩阵将MIM作为基因调控网络的初始权重矩阵。
8.根据权利要求1所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S4步骤中：
在通过路径一致性算法遍历整个基因调控网络的过程中，计算不同基因对在多个基因影响下的条件互包含信息CMI2；根据计算得到的动态阈值更新当前阈值，并删除低于该阈值的基因间关系，判断并更新整个邻接矩阵，直到不删除任何关系为止；通过上述方式逐渐删除冗余边，进而形成一个更准确的基因调控网络；
动态阈值筛选过程具体为：使用动态设置阈值方法来筛选基因之间的冗余调控关系，阈值越大，对基因调控关系的控制就越严格；首先根据设置的初始阈值初步筛选并删除冗余关系，同时设置order值，即从与基因对相关联的所有基因中抽取的最大值，从1到order随机选择当前基因对下相关联的基因，作为间接调控影响来计算条件互包含信息CMI2；随机选择产生的组合数越多，该基因对受到间接调控的可能性就越大，相应阈值也应该增加，因此通过设置动态阈值实现实时删除图中的冗余边。
9.根据权利要求8所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S4步骤中，计算不同基因对在多个基因影响下的条件互包含信息CMI2的具体方法为：
对于基因G1，在初始基因调控网络中，G1和G2之间存在一条边；然后搜索与G1和G2具有连接边的基因Gi；假设满足上述条件的基因Gi的数目为Order_Max，即与G1和G2相互作用的最大基因数为Order_Max；如果Order＝0，则根据初始阈值λ0获得初始基因调控网络，然后逐渐增加Order以增量更新基因调控网络；
具体过程为，从Order＝1开始，计算所有Gi的条件交互信息CMI2(G1,G2|Gi)，返回最大值为MAX(CMI2(G1,G2|Gi))；如果MAX(CMI2(G1,G2|Gi))小于阈值λ1，则删除G1和G2之间的边，并将基因邻接矩阵中G1(1,2)的值设置为0；如果MAX(CMI2(G1,G2|Gi))大于阈值λ1，增加Order的值；计算Order基因的CMI2(G1,G2|Gi)，并依次将计算结果与阈值λi进行比较，直到Order＝Order_Max或Order达到Order_Set的最大值；遍历所有基因，逐渐删除不正确的基因间关系，从而形成较准确的基因调控网络；
对于给定Z的随机变量X和Y，CMI2的计算方式如Eq.(17)所示：
在Eq.(17)中，DKL(P(X,Y,Z)||PX→Y(X,Y,Z))是从P(X,Y,Z)到PX→Y(x,y,z)的KL散度，用于测量两个分布之间的距离；P(X,Y,Z)表示X,Y和Z之间的联合概率；PX→Y(x,y,z)为X到Y的干扰概率；X和Y之间的因果变量强度见Eq.(18)：
CX→Y(X；Y|Z)＝DKL(P(X,Y,Z)||PX→Y(X,Y,Z)) (18)
CX→Y(X；Y|Z)是不对称的，CMI2(X；Y|Z)是CX→Y(X；Y|Z)和CY→X(Y；X|Z)的平均值；在扩展Eq.(18)后，得到Eq.(19)：
DKL(P(X,Y,Z)||PY→X(X,Y,Z))＝I(X；Y|Z)+DKL(P(X|Z)||PY→X(X|Z)) (19)
通过结合Eq.(20)和Eq.(17)，得到CMI2(X；Y|Z)，如Eq.(20)所示：
其中PX→Y(Y|Z)的计算方法如Eq.(21)所示：
为了估计Eq.(21)中的密度函数，使用高斯核密度估计方法计算相应的密度函数；应用基因表达数据受多元高斯分布影响的广泛假设进行计算，这有助于更有效的估计CMI2(X；Y|Z)。
10.根据权利要求9所述的基于混合熵优化互信息的基因调控网络构建方法，其特征在于，在所述S4步骤中，动态阈值筛选的具体方法为：
使用动态设置阈值的方法来筛选基因之间的冗余调控关系；阈值越大，对基因调控关系的控制就越严格；考虑到多种间接调控关系对基因间直接调控关系的影响，在计算各阶CMI2(Gi,Gj|Gk)时，同时考虑了从Order_Max中提取的阶基因的初始阈值α和组合数m；m依赖于Order的值和所有相关基因Order_Max的数量；m的值越大，受到的间接调控就越大，从而增加相应的阈值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210294214.9A CN114925837B (zh) | 2022-03-23 | 2022-03-23 | 基于混合熵优化互信息的基因调控网络构建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210294214.9A CN114925837B (zh) | 2022-03-23 | 2022-03-23 | 基于混合熵优化互信息的基因调控网络构建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114925837A true CN114925837A (zh) | 2022-08-19 |
| CN114925837B CN114925837B (zh) | 2024-04-16 |
Family
ID=82805740
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210294214.9A Active CN114925837B (zh) | 2022-03-23 | 2022-03-23 | 基于混合熵优化互信息的基因调控网络构建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114925837B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116129992A (zh) * | 2023-04-17 | 2023-05-16 | 之江实验室 | 基于图神经网络的基因调控网络构建方法及系统 |
| CN117409962A (zh) * | 2023-12-14 | 2024-01-16 | 北京科技大学 | 一种基于基因调控网络的微生物标记物的筛选方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20020019705A1 (en) * | 1998-10-21 | 2002-02-14 | Kauffman Stuart A. | System and method for analysis of genetic networks |
| CN110555530A (zh) * | 2019-09-02 | 2019-12-10 | 东北大学 | 一种基于分布式的大规模基因调控网络构建方法 |
| CN111223523A (zh) * | 2020-01-06 | 2020-06-02 | 中南大学 | 基于多时滞因果熵的基因调控网络构建方法及系统 |
| US20220068478A1 (en) * | 2020-09-01 | 2022-03-03 | NEC Laboratories Europe GmbH | Scalable, accurate and reliable measure of variable dependence and independence, and utilization of the measure to train a neural network |
-
2022
- 2022-03-23 CN CN202210294214.9A patent/CN114925837B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20020019705A1 (en) * | 1998-10-21 | 2002-02-14 | Kauffman Stuart A. | System and method for analysis of genetic networks |
| CN110555530A (zh) * | 2019-09-02 | 2019-12-10 | 东北大学 | 一种基于分布式的大规模基因调控网络构建方法 |
| CN111223523A (zh) * | 2020-01-06 | 2020-06-02 | 中南大学 | 基于多时滞因果熵的基因调控网络构建方法及系统 |
| US20220068478A1 (en) * | 2020-09-01 | 2022-03-03 | NEC Laboratories Europe GmbH | Scalable, accurate and reliable measure of variable dependence and independence, and utilization of the measure to train a neural network |
Non-Patent Citations (1)
| Title |
|---|
| 缑葵香;宫秀军;汤莉;: "基于时序互信息构建基因调控网络", 天津大学学报, no. 07, 15 July 2010 (2010-07-15) * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116129992A (zh) * | 2023-04-17 | 2023-05-16 | 之江实验室 | 基于图神经网络的基因调控网络构建方法及系统 |
| CN117409962A (zh) * | 2023-12-14 | 2024-01-16 | 北京科技大学 | 一种基于基因调控网络的微生物标记物的筛选方法 |
| CN117409962B (zh) * | 2023-12-14 | 2024-03-29 | 北京科技大学 | 一种基于基因调控网络的微生物标记物的筛选方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114925837B (zh) | 2024-04-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11568957B2 (en) | Methods and systems for copy number variant detection | |
| CN114925837A (zh) | 基于混合熵优化互信息的基因调控网络构建方法 | |
| CN110675912B (zh) | 一种基于结构预测的基因调控网络构建方法 | |
| CN113488104B (zh) | 基于局部和全局的网络中心性分析的癌症驱动基因预测方法及系统 | |
| Siddiqi et al. | A new heuristic for the data clustering problem | |
| CN115019883A (zh) | 一种基于多网络图卷积的癌症驱动基因识别方法 | |
| Ruffieux et al. | A global-local approach for detecting hotspots in multiple-response regression | |
| CN111832817A (zh) | 基于mcp罚函数的小世界回声状态网络时间序列预测方法 | |
| Lei et al. | An approach of gene regulatory network construction using mixed entropy optimizing context-related likelihood mutual information | |
| Yelmen et al. | Deep convolutional and conditional neural networks for large-scale genomic data generation | |
| Dudoit et al. | Loss-based estimation with cross-validation: Applications to microarray data analysis | |
| Das et al. | Applying restrained genetic algorithm for attribute reduction using attribute dependency and discernibility matrix | |
| Chen et al. | Genome-scale protein function prediction in yeast Saccharomyces cerevisiae through integrating multiple sources of high-throughput data | |
| Sitarčík et al. | epiBAT: Multi-objective bat algorithm for detection of epistatic interactions | |
| Tareq et al. | A new density-based method for clustering data stream using genetic algorithm | |
| Li et al. | Functional dissection of regulatory models using gene expression data of deletion mutants | |
| CN113450872B (zh) | 磷酸化位点特异激酶的预测方法 | |
| Lewin et al. | Bayesian methods for gene expression analysis | |
| CN109360602B (zh) | 基于模糊优先性的dna编码序列设计方法及装置 | |
| CN117727373B (zh) | 基于样本和特征双加权的特征约简中智c-均值聚类方法 | |
| Mohammadi et al. | Dealing with missing values in microarray data | |
| Hequet | Biologically-informed interpretable deep learning techniques for BMI prediction and gene interaction detection | |
| Yonar et al. | Evaluation and comparison of metaheuristic methods to estimate the parameters of gamma distribution | |
| Vilsen | Statistical modelling of Massively Parallel Sequencing data in forensic genetics | |
| Robin et al. | Applications in genomics |
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 |