CN113066522A - 一种基于模块化识别的基因网络推理方法 - Google Patents
一种基于模块化识别的基因网络推理方法 Download PDFInfo
- Publication number
- CN113066522A CN113066522A CN202110309281.9A CN202110309281A CN113066522A CN 113066522 A CN113066522 A CN 113066522A CN 202110309281 A CN202110309281 A CN 202110309281A CN 113066522 A CN113066522 A CN 113066522A
- Authority
- CN
- China
- Prior art keywords
- gene
- module
- network
- genes
- modules
- 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 279
- 238000000034 method Methods 0.000 title claims abstract description 64
- 230000033228 biological regulation Effects 0.000 claims abstract description 59
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 52
- 230000014509 gene expression Effects 0.000 claims abstract description 49
- 230000006870 function Effects 0.000 claims abstract description 18
- 238000010606 normalization Methods 0.000 claims abstract description 12
- 238000012545 processing Methods 0.000 claims abstract description 7
- 238000012549 training Methods 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 59
- 230000001105 regulatory effect Effects 0.000 claims description 34
- 238000012417 linear regression Methods 0.000 claims description 21
- 238000004364 calculation method Methods 0.000 claims description 20
- 230000008569 process Effects 0.000 claims description 17
- 239000013598 vector Substances 0.000 claims description 15
- 230000002087 whitening effect Effects 0.000 claims description 11
- 230000008859 change Effects 0.000 claims description 9
- 230000001174 ascending effect Effects 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 108091008025 regulatory factors Proteins 0.000 claims description 6
- 102000037983 regulatory factors Human genes 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 5
- 239000013604 expression vector Substances 0.000 claims description 5
- 238000005457 optimization Methods 0.000 claims description 5
- 238000011160 research Methods 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 4
- 238000012174 single-cell RNA sequencing Methods 0.000 claims description 4
- 108700005081 Overlapping Genes Proteins 0.000 claims description 3
- 238000002493 microarray Methods 0.000 claims description 3
- 238000012887 quadratic function Methods 0.000 claims description 3
- 238000011084 recovery Methods 0.000 claims description 3
- 238000000551 statistical hypothesis test Methods 0.000 claims description 3
- 239000000203 mixture Substances 0.000 claims description 2
- 238000004088 simulation Methods 0.000 claims description 2
- 230000004927 fusion Effects 0.000 abstract description 2
- 238000012163 sequencing technique Methods 0.000 abstract 1
- 241000588724 Escherichia coli Species 0.000 description 8
- 230000000694 effects Effects 0.000 description 7
- 238000012360 testing method Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 4
- 230000002068 genetic effect Effects 0.000 description 4
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 210000004027 cell Anatomy 0.000 description 3
- 240000004808 Saccharomyces cerevisiae Species 0.000 description 2
- 230000024245 cell differentiation Effects 0.000 description 2
- 238000003064 k means clustering Methods 0.000 description 2
- 230000001717 pathogenic effect Effects 0.000 description 2
- 101150044182 8 gene Proteins 0.000 description 1
- 101150040974 Set gene Proteins 0.000 description 1
- 230000004913 activation Effects 0.000 description 1
- 230000006399 behavior Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 125000004432 carbon atom Chemical group C* 0.000 description 1
- 230000011712 cell development Effects 0.000 description 1
- 230000032823 cell division Effects 0.000 description 1
- 230000010261 cell growth Effects 0.000 description 1
- 238000003066 decision tree Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 238000001415 gene therapy Methods 0.000 description 1
- 230000002710 gonadal effect Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035790 physiological processes and functions Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000008844 regulatory mechanism Effects 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 230000020509 sex determination Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 210000000130 stem cell Anatomy 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000002626 targeted therapy Methods 0.000 description 1
- 238000013518 transcription Methods 0.000 description 1
- 230000035897 transcription Effects 0.000 description 1
- 230000002103 transcriptional effect Effects 0.000 description 1
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
- G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
-
- 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)
- Bioinformatics & Cheminformatics (AREA)
- Medical Informatics (AREA)
- Biophysics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- General Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Physiology (AREA)
- Artificial Intelligence (AREA)
- Bioethics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Data Mining & Analysis (AREA)
- Databases & Information Systems (AREA)
- Epidemiology (AREA)
- Evolutionary Computation (AREA)
- Public Health (AREA)
- Software Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
Abstract
本发明公开了一种基于模块化识别的基因网络推理方法。该方法将n个基因的表达信息作为训练集，其中每个基因有m个样本；采用ICA‑FDR算法对基因进行模块识别并将n个基因划分至不同基因模块中；基因模块内部的调控关系采用以梯度提升树为基础的算法进行推理，基因模块间的调控关系采用以稀疏回归为基础的算法进行推理，得出每个基因对的相关性得分；针对每个基因模块内与模块间推理得到的相关性得分分别进行归一化处理后合并，并按照得分降序排序得到最终基因调控网络。本发明为基因模块识别与基因网络推理提供了一个无缝衔接的融合框架，提高了基因模块识别的准确性，增加了基因调控网络功能的可解释性。
Description
技术领域
本发明属于生物信息学中的基因调控网络推理领域，尤其涉及一种基于模块化识别的基因网络推理方法。
背景技术
如何在转录水平上准确阐明调控因子与目标基因之间的调控关系，是近年来计算生物学和生物信息学的核心挑战之一。更准确地识别转录调控因子与目标基因之间的调控关系对于探究细胞生长和分裂、细胞分化和发育等活动规律至关重要。此外，基因调控网络对于现代医学研究提供了有力帮助，能够从生命活动最底层—基因控制的角度模拟和预测致病基因，这将有助于医护人员对病人精准地实施靶向治疗，并有助于相应药物的研发。然而，基因调控网络推理方法大部分针对整个基因网络进行推理，但是这些方法仅推断了基因网络的拓扑结构却缺乏了明确的生物学解释，从而限制了基因网络在致病基因预测和基因治疗等领域的应用。研究表明基因调控网络结构具有模块性。在一定程度上，基因调控网络的模块化分析可以辅助我们理解整体网络的动力学性质，探索某些未知的基因功能模块，并指导基因网络的重构。然而，目前大多数经典基因模块识别的算法应用的先决条件是需要知道基因网络的拓扑结构，无法由数据驱动角度切入研究基因模块。
发明内容
基于背景技术存在的问题，本发明提出了一种基于模块化识别的高效基因网络推理方法。对基因组学和转录组学数据的深入挖掘，结合基因网络固有的模块特点，从数据驱动推理出具有社区结构的基因调控网络，为基因模块识别与基因网络推理提供了一个无缝衔接的融合框架。
为了解决上述问题，本发明采用如下技术方案：一种基于模块化识别的基因网络推理方法，该方法包括以下步骤：
S1:针对基因调控网络的推理问题，该方法从数据驱动角度切入进行基因调控网络的推理过程。假设研究对象为包含n个基因的基因调控网络，每个基因有m个表达数据样本。该方法将大小为m×n的基因表达矩阵X作为训练数据集。
S2:采用ICA-FDR算法对基因的表达矩阵进行分析，识别表达功能相近的基因。采用FastICA算法将表达矩阵X分解为混合矩阵A和源矩阵S。源矩阵S蕴含了基因数据之间的功能相似性信息，对源矩阵S进行FDR计算分析，将n个基因划分至不同的基因模块中，完成基因模块识别。
S3:针对每个基因模块内基因之间的较为紧密的调控关系，该方法采用以梯度提升树为基础的算法进行子网络推理。梯度提升树算法计算每一个基因模块中调控因子与目标基因之间的相关性得分，并根据相关性得分进行降序排列，存在调控关系的可能性较大的调控边将集中在得分列表前端。
S4:针对两个基因分别属于不同基因模块之间的基因对，该方法采用以稀疏线性回归为基础的算法进行模块间网络推理。采用稀疏线性回归算法计算调控因子与目标基因之间的相关性得分，并根据相关性得分进行降序排列，存在调控关系的可能性较大的调控边将集中在得分列表前端。
S5:对模块内基因子网络的调控边得分与模块间网络的调控边得分进行归一化处理后合并，并按降序重新排列合并后的得分列表，得到总相关性得分列表。
S6:根据基因调控网络的复杂度及研究需求，选取经验阈值，将相关性得分大于等于此阈值的调控边保留，小于此阈值的调控边舍弃，完成基因调控网络推理过程。
进一步地，所述步骤S1中的基因表达矩阵的数据种类具体包括但不限于：
S11:时序微阵列表达数据；
S12:单细胞RNA测序(scRNA-seq)数据；
S13:基因表达仿真数据。
进一步地，所述步骤S2中ICA-FDR算法的基本步骤包括：
S21:确定基因模块的数量n_comps，对基因表达矩阵进行白化预处理，采用FastICA算法寻找一个迭代优化方向ω，使得此方向迭代的非高斯性最大。迭代收敛后，将表达矩阵X分解为混合矩阵A和源矩阵S，X＝A×S，源矩阵中即包含基因与基因模块之间的关联信息。
S22:采用错误发现率(FDR)分析源矩阵S，将表达功能相近的基因划分至同一基因模块中。
进一步地，所述步骤S21中的具体计算方法的步骤为：
S211:对表达矩阵进行白化预处理的计算公式为：
其中，x代表基因的表达向量，若基因向量的协方差矩阵为E{xxT}，那么其特征分解式为E{xxT}＝EDET。其中，E是协方差矩阵E{xxT}的特征向量组成的矩阵，D＝diag(d1,...,dn)是其特征值的对角矩阵。白化处理使学习过程的参数减少至二分之一，使原始矩阵的复杂度由n2变为
S212:经过白化处理后，采用FastICA找到一个最优优化方向ω使得非高斯性JG(ω)最大化，非高斯性由负熵衡量：
JG(ω)＝H(ωgauss)-H(ω)
上式中ωgauss是高斯分布，H(·)代表熵值，当JG(ω)越大非高斯性将越大。若设ν是单位方差并且零均值的高斯变量，G(·)是用于提高参数估计鲁棒性的非二次函数，例如：
其中，a1是任意常数；那么非高斯性的具体计算式为：
进一步地，所述步骤S22中的具体计算方法为：
S221:对于源矩阵S当中的基因数据进行统计学假设检验，分析每个基因模块中基因的集中程度并得到P-value，对所有候选基因的P-value进行升序排列，P-value的序号记为i，并计算相应的Q值：
上式中pik代表在第k个基因模块中，升序排列后的第i个P-value，n是候选基因个数。
S222:根据基因数量规模及基因模块的个数选取阈值q_cuoff，将每个基因模块中Q值小于此阈值q_cuoff的基因分配给此基因模块。一个基因可能属于多个基因模块，即允许存在基因模块重叠现象。
进一步地，所述步骤S4中稀疏线性回归建模方法为：
S41:针对分别属于不同模块的两个基因之间的调控关系，建立如下稀疏线性回归模型：
Et＝αr1,tEr1+αr2,tEr2+…+αrt,tErt+βt
定义基因集合为G＝{g1,g2,…,gn}，其中Et代表被调控目标基因gt的数据表达向量，并且gt∈G。与基因gt不属于同一模块的潜在调控因子集合表示为G-t，其表达式为G-t＝{gr1,gr2,…,grt}，那么Er1,Er2,…,Ert为潜在调控因子G-t的表达数据向量。αrt,t代表潜在调控因子grt对目标基因gt的回归系数，βt代表回归过程中的噪声向量。
S42:采用稀疏线性回归的相关目标函数来进行α,β参数的拟合，收敛后α即为调控因子与目标基因之间的相关性得分。稀疏线性回归的相关目标函数包括但不仅限于L1、L2等。
进一步地，所述步骤S5所采用的相关性得分归一化方法为最大最小归一化，其计算公式为：
式中si代表调控边i所对应的相关性得分，s为相关性得分向量。采用此归一化方法分别处理所有模块内子网络和模块间网络的相关性得分列表。
进一步地，根据步骤S2识别出的基因模块，采用Frr指标评价基因模块识别的准确性；根据步骤S5得到的总相关性得分列表，采用AUROC与AUPR指标评价基因网络推理的准确性。
进一步地，评估模块识别的准确性，采用Frr指标将识别结果与真实基因模块进行对比，计算方法如下：
其中Recovery代表真实基因模块拟合模块识别结果的程度，Relevance代表模块识别结果拟合真实基因模块的程度，具体公式如下：
其中，
和M分别代表真实基因模块集合与识别模块集合，Jaccard(·)则用来度量两个基因模块集合之间的相似度，Jaccard距离越大，真实模块与识别模块之间的相似度越低。
进一步地，评估基因调控网络推理的准确性，根据总相关性评分表绘制ROC曲线与PR曲线。给定一个分数阈值，评分表中的推理结果与真实基因调控网络对比，可以分为四类：真阳性(True Positive,TP)，真阴性(True Negative,TN)，假阳性(False Positive,FP)，假阴性(False Negative,FN)。其中：
随着阈值变化，ROC曲线表现了FPR和TPR之间的变化趋势，PR曲线表示了recall和precision的变化趋势，分别计算ROC与PR的曲线下面积得到AUROC与AUPR，参数越大则表明推理结果越准确。
与现有技术相比，本发明的有益效果是：
(1)提供了一种具有基因模块识别功能的基因调控网络推理算法。采用以ICA为基础的统计学分析将功能相近的基因划分至同一模块，为研究人员进一步研究基因调控网络背后更深层次的生物物理意义提供了有效参考，增加了基因调控网络的生物学可解释性。
(2)集成梯度提升树准确性与稀疏线性回归高效性的优势，对基因模块内信息交流较为频繁的模块内网络采用梯度提升树算法进行推理，对基因模块间调控关系较为稀疏的网络采用稀疏线性回归的高效算法进行推理。本发明进一步提升了基因调控网络推理的速度与准确性。
附图说明
图1是本发明具体实施例中基于模块化识别的基因网络推理流程图；
图2是本发明具体实施例中基于ICA-FDR算法的基因模块识别原理示意图；
图3是本发明具体实施例中SCODE数据集基因网络推理的ROC曲线与PR曲线。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于任何熟悉本技术领域的技术人员理解本发明，但本发明的保护范围不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。
下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。
本实施例运行的硬件环境：笔记本电脑一台，CPU：2.6GHz，RAM：8.0GB；软件环境：Python 3.6，R 3.3.3；操作平台：Win10。
本实施例在时序数据集与单细胞数据集上测试了所提出方法的效果，数据集包括：BEELINE项目中的gonadal sex determination(GSD)数据集、SCODE项目中的PrE、MEF、DE三个真实单细胞数据集、DREAM5挑战中的E.coli和Yeast数据集。实验数据集的详细信息如表1所示。
表1实验数据集详细信息
本实施例提出的一种基于模块化识别的基因网络推理方法，流程图如附图1所示，以时序数据的DREAM5的E.coli数据集为例说明本方法的实现流程，具体包含以下步骤：
S1:对DREAM5的E.coli公共数据集进行数据预处理，E.coli的基因调控网络为包含4511个基因的基因调控网络，每个基因有805个基因表达数据样本。因此将大小为m×n＝805×4511的基因表达矩阵X作为训练数据集。此外，E.coli数据属于时序微阵列表达数据，表达数据如图1中流程A所示。
S2:采用ICA-FDR算法对基因的表达矩阵进行分析，识别表达功能相近的基因。采用FastICA算法将表达矩阵X分解为混合矩阵A和源矩阵S。源矩阵S蕴含了基因数据之间的功能相似性信息，对源矩阵S进行FDR计算分析，将n＝4511个基因划分至不同的基因子集(基因模块)中。其具体技术方案为：
S21:根据E.coli网络数据集中存在4511个基因，此基因调控网络规模较大，因此确定基因模块的数量n_comps＝80，并对基因表达矩阵进行白化预处理，采用FastICA算法寻找一个迭代优化方向ω，使得此方向迭代的非高斯性最大。迭代收敛后，将表达矩阵分解为混合矩阵和源矩阵X＝A×S，源矩阵S中即包含基因与基因模块之间的关联信息。矩阵中基因表达值可由公式表示为：
上式中xij代表基因j在样本i的表达值，aik代表样本i对模块k的贡献度，skj则代表基因j在基因模块k中的激活水平。
S211:对表达矩阵进行白化预处理的计算公式为：
其中，x代表基因的表达向量，若基因向量的协方差矩阵为E{xxT}，那么其特征分解式为E{xxT}＝EDET。其中，E是协方差矩阵E{xxT}的特征向量组成的矩阵，D＝diag(d1,...,dn)是其特征值的对角矩阵。白化处理使学习过程的参数减少至二分之一，使原始矩阵的复杂度由45112变为
S212:经过白化处理后，采用FastICA找到一个最优优化方向ω使得非高斯性JG(ω)最大化，非高斯性由负熵衡量：
JG(ω)＝H(ωgauss)-H(ω)
上式中ωgauss是高斯分布，H(·)代表熵值，当JG(ω)越大非高斯性将越大。本方法具体表达式为：
上式中，ν是单位方差并且零均值的高斯变量，G(·)是用于提高参数估计鲁棒性的非二次函数，例如：
其中，a1是任意常数。
S22:采用错误发现率(FDR)分析源矩阵S，将表达功能相近的基因划分至同一基因模块中，其具体技术方案为：
S221:对于源矩阵S当中的基因数据进行统计学假设检验，分析每个基因模块中基因的集中程度并得到P-value，对所有候选基因的P-value进行升序排列，P-value的序号记为i，并计算相应的Q值：
上式中pik代表在第k个基因模块中，升序排列后的第i个P-value，n是候选基因个数，在此实例中即为n＝4511。
S222:根据E.coli基因数量规模及基因模块的个数选取阈值q_cuoff，将每个基因模块中Q值小于此阈值q_cuoff的基因分配给此基因模块。一个基因可能属于多个基因模块，即允许存在基因模块重叠现象。此实例中对E.coli网络选取阈值为：q_cutoff＝1×10-3。
综上所述，ICA-FDR算法进行基因模块识别的原理如图2所示。
S3:采用ICA-FDR2算法对基因表达矩阵进行分析，作为ICA-FDR算法的对比实验。
ICA-FDR2算法具体为：采用ICA算法生成每个基因模块的过程中，将根据基因权重的正负，分别将其划入不同的基因模块，后续FDR分析同步骤S22。即ICA-FDR2生成的基因模块数量为ICA-FDR的2倍。
S4:针对每个基因模块内基因之间的较为紧密的调控关系，该方法采用以梯度提升树为基础的算法进行子网络推理。梯度提升树算法计算每一个基因模块中调控因子与目标基因之间的相关性得分，并根据相关性得分进行降序排列，存在调控关系的可能性较大的调控边将集中在得分列表前端，如图1中流程C所示。
在基于梯度提升树回归的基因调控网络推理方法中，目标基因j的行为受其他调控因子控制的数学模型描述如下式：
上式中，
代表在样本s中的基因表达向量，因此基因j即为目标基因，m代表样本个数；
为除目标基因j之外的所有此模块内的基因在样本s中的表达向量，这些基因均为j的潜在调控因子；kn代表目标基因j所在基因模块的基因数量，即
中共有(kn-1)个基因；εs代表均值为0的随机噪声；f(·)函数由决策树集合的建立方法决定。对于每个树节点φ，分支后的方差的减少值计算如下：
I(φ)＝SVar(S)-SlVar(Sl)-SrVar(Sr)
上式中S是在树节点φ中包含的样本集合，Sl和Sr分别代表左子树与右子树样本集合，Var(·)代表树分枝方差。
S5:针对两个基因分别属于不同基因模块之间的基因对，该方法采用以稀疏线性回归为基础的算法进行模块间网络推理。采用稀疏线性回归算法计算调控因子与目标基因之间的相关性得分，并根据相关性得分进行降序排列，存在调控关系的可能性较大的调控边将集中在得分列表前端，如图1中流程D所示。稀疏线性回归的模块间推理具体技术方案为：
S51:针对分别属于不同模块的两个基因之间的调控关系，建立如下稀疏线性回归模型：
Et＝αr1,tEr1+αr2,tEr2+…+αrt,tErt+βt
定义基因集合为G＝{g1,g2,…,gn}，其中Et代表被调控目标基因gt的数据表达向量，并且gt∈G。与基因gt不属于同一模块的潜在调控因子集合表示为G-t，其表达式为G-t＝{gr1,gr2,…,grt}，那么Er1,Er2,…,Ert为潜在调控因子G-t的表达数据向量。αrt,t代表潜在调控因子grt对目标基因gt的回归系数，βt代表回归过程中的噪声向量。
S52:采用稀疏线性回归的相关目标函数进行α,β参数的拟合，收敛后α即为调控因子与目标基因之间的相关性得分。稀疏线性回归的相关目标函数包括但不仅限于L1、L2等，在此实例中目标函数的表达式为：
S6:对模块内基因子网络的调控边得分与模块间网络的调控边得分进行归一化处理后合并，并按降序重新排列合并后的得分列表，得到总相关性得分列表。所采用的相关性得分归一化方法为最大最小归一化，如图1中流程E所示，其计算公式为：
上式中si代表调控边i所对应的相关性得分，s为相关性得分向量。采用此归一化方法分别处理所有模块内子网络和模块间网络的相关性得分列表。
S7:根据基因调控网络的复杂度及研究需求，选取合适的阈值分数，将相关性得分大于等于此阈值的调控边保留，小于此阈值的调控边舍弃，完成基因调控网络推理过程。
S8:根据步骤S2识别出的基因模块，采用Frr指标评价基因模块识别的准确性；根据步骤S6得到的总相关性得分列表，采用AUROC与AUPR指标评价基因网络推理的准确性。具体技术方案为：
S81:评估模块识别的准确性，采用Frr指标将识别结果与真实基因模块进行对比，计算方法如下：
其中Recovery代表真实基因模块拟合模块识别结果的程度，Relevance代表模块识别结果拟合真实基因模块的程度，具体公式如下：
其中，
和M分别代表真实基因模块集合与识别模块集合，Jaccard(·)则用来度量两个基因模块集合之间的相似度，Jaccard距离越大，真实模块与识别模块之间的相似度越低。
S82:评估基因调控网络推理的准确性，根据总相关性评分表绘制ROC曲线与PR曲线。给定一个分数阈值，评分表中的推理结果与真实基因调控网络对比，可以分为四类：真阳性(True Positive,TP)，真阴性(True Negative,TN)，假阳性(False Positive,FP)，假阴性(False Negative,FN)。其中：
随着阈值变化，ROC曲线表现了FPR和TPR之间的变化趋势，PR曲线表示了recall和precision的变化趋势，分别计算ROC与PR的曲线下面积得到AUROC与AUPR，数值越大则表明推理结果越准确。
本发明的主要贡献是为基因模块识别与基因网络推理提出了一个无缝的框架，实现了基于已识别的基因模块执行基因网络推断。
为了从基因表达数据中检测出功能相近的基因并划分至同一模块，在本发明方法中应用了基于ICA的分解算法。为了检验本发明的ICA-FDR算法模块识别的准确度，选取了若干模块识别算法包括：ICA-FDR2、ICA-zscore、PCA、K-means，在不同数据集上进行实验并用Frr指标将识别结果与数据集提供的真实基因模块(模块金标准)进行对比，实验结果如表2-4所示：
表2 DREAM5网络基因模块识别Frr评估结果
如上表中，最优Frr值由加粗字体表示，Minimal，Strict和Interconnected是基因模块金标准的三种不同的模块定义。从表2中可以看出，在三种基因模块定义下，基于ICA或PCA分解的方法均明显优于K-means聚类。从三种模块定义的角度进行分析，Minimal和Strict下获得的Frr指标又明显高于Interconnected组。
为了直观展示模块识别的效果，以E.coli网络为例对ICA-FDR识别出的不同基因调控模块进行不同形状的标注。如图1所示，节点代表基因，而边则代表两个基因之间的调控关系。
表3 BEELINE网络基因模块识别Frr评估结果
表4 SCODE网络基因模块识别Frr评估结果
表3和表4分别是在BEELINE合成数据集和SCODE真实单细胞数据集上测试的结果，Frr指标表明ICA-FDR算法与PCA分解和k均值聚类相比，在检测基因模块准确性方面具有明显优势，与DREAM5数据集的结果基本保持一致。由此可见，本发明为基因模块识别提供了一种数据驱动的解决方案，即使没有精确的生物学注释，也可以直接从基因表达数据中挖掘功能模块，对于加深网络调控机制的理解至关重要。
为了检验本发明提出的方法的准确性，需要将本方法推理得出的基因调控网络与金标准网络进行对比，金标准网络即为数据集所对应的真实基因调控网络，其中标注了哪些基因之间存在调控关系。绘制推理结果ROC与PR曲线并计算曲线下面积，面积取值范围为[0,1]，面积越大表明推理效果越好。
对比方法选用了以线性回归为基础的Ridge、Linear Regression、TIGRESS，在DREAM5Challenge中被提出的以GBDT框架为基础的GRNBoost2算法、基于互信息的CLR算法。将表1中的8个基因网络进行基因网络推理，以ROC与PR曲线下面积AUROC与AUPR两个指标衡量网络推理准确性及有效性，实验结果如下表所示：
表5DREAM5基因网络推理结果对比
目前在基因网络推理领域，基于GBDT框架的GRNBoost2算法的准确度在众多测试集上表现良好，可以视为顶级基因网络推理算法。在计算时间上，DREAM5 E.coli网络花费了41h 46min来完成GRNBoost2推理过程，而本方法只花费了1h 8min；使用GRNBoost2推理，DREAM5 Yeast网络花费了33h 36min完成GRNBoost2的推理过程，而本方法仅花费了1h9min。从表5中的准确性测试表明，本方法与GRNBoost2算法相比AUROC与AUPR几乎处于同一水平。因此，与GRNBoost2相比，本方法显著降低了算法复杂度，提高了网络推理的速度，并且未引起网络推理准确性的明显下降。
表6 BEELINE基因网络推理结果对比
表7 SCODE基因网络推理结果对比
表6将本方法与其他经典网络推理算法进行对比，结果表明在BEELINE网络中，本方法可以在这4种方法中达到较高的准确度，并基本与GRNBoost2算法推理效果持平。针对SCODE项目中真实单细胞数据集的基因网络推理，AUPR和AUROC计算结果表明，本算法在PrE和MEF网络上可以达到最优效果。图3详细绘制了几种对比算法的PR和ROC曲线，PR曲线的形状表明，与其他算法相比由本方法推断的降序得分列表的顶部具有更多符合真实网络结构的调控关系。对于涉及人类干细胞分化的DE数据集，本方法的准确性指数低于现有的GRNBoost2，其原因可能因为人类基因调控是较为复杂生理过程，导致推理准确度普遍较低。
上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。
Claims (10)
1.一种基于模块化识别的基因网络推理方法，其特征在于，包括以下步骤：
S1:针对基因调控网络的推理问题，从数据驱动角度切入进行基因调控网络的推理过程。假设研究对象为包含n个基因的基因调控网络，每个基因有m个表达数据样本，将大小为m×n的基因表达矩阵X作为训练数据集。
S2:采用ICA-FDR算法对基因的表达矩阵进行分析，识别表达功能相近的基因。采用FastICA算法将表达矩阵X分解为混合矩阵A和源矩阵S。源矩阵S蕴含了基因数据之间的功能相似性信息，对源矩阵S进行FDR计算分析，将n个基因划分至不同的基因模块中，完成基因模块识别。
S3:针对每个基因模块内基因之间的较为紧密的调控关系，采用以梯度提升树为基础的算法进行子网络推理。梯度提升树算法计算每一个基因模块中调控因子与目标基因之间的相关性得分，并根据相关性得分进行降序排列，存在调控关系的可能性较大的调控边将集中在得分列表前端。
S4:针对两个基因分别属于不同基因模块之间的基因对，采用以稀疏线性回归为基础的算法进行模块间网络推理。采用稀疏线性回归算法计算调控因子与目标基因之间的相关性得分，并根据相关性得分进行降序排列，存在调控关系的可能性较大的调控边将集中在得分列表前端。
S5:对模块内基因子网络的调控边得分与模块间网络的调控边得分进行归一化处理后合并，并按降序重新排列合并后的得分列表，得到总相关性得分列表。
S6:根据基因调控网络的复杂度及研究需求，选取经验阈值，将相关性得分大于等于此阈值的调控边保留，小于此阈值的调控边舍弃，完成基因调控网络推理过程。
2.根据权利要求1所述的一种基于模块化识别的基因网络推理方法，其特征在于，所述步骤S1中的基因表达矩阵的数据种类具体包括但不限于：
S11:时序微阵列表达数据；
S12:单细胞RNA测序(scRNA-seq)数据；
S13:基因表达仿真数据。
3.根据权利要求1所述的一种基于模块化识别的基因网络推理方法，其特征在于，所述步骤S2中ICA-FDR算法的基本步骤包括：
S21:确定基因模块的数量n_comps，对基因表达矩阵进行白化预处理，采用FastICA算法寻找一个迭代优化方向ω，使得此方向迭代的非高斯性最大。迭代收敛后，将表达矩阵X分解为混合矩阵A和源矩阵S，X＝A×S，源矩阵中即包含基因与基因模块之间的关联信息。
S22:采用错误发现率(FDR)分析源矩阵S，将表达功能相近的基因划分至同一基因模块中。
4.根据权利要求3所述的一种基于模块化识别的基因网络推理方法，其特征在于，所述步骤S21中的具体计算方法的步骤为：
S211:对表达矩阵进行白化预处理的计算公式为：
其中，x代表基因的表达向量，若基因向量的协方差矩阵为E{xxT}，那么其特征分解式为E{xxT}＝EDET。其中，E是协方差矩阵E{xxT}的特征向量组成的矩阵，D＝diag(d1,...,dn)是其特征值的对角矩阵。
S212:经过白化处理后，采用FastICA找到一个最优优化方向ω使得非高斯性JG(ω)最大化，非高斯性由负熵衡量：
JG(ω)＝H(ωgauss)-H(ω)
上式中ωgauss是高斯分布，H(·)代表熵值，当JG(ω)越大非高斯性将越大。
若设ν是单位方差并且零均值的高斯变量，G(·)是用于提高参数估计鲁棒性的非二次函数，非高斯性的具体计算式为：
5.根据权利要求3所述的一种基于模块化识别的基因网络推理方法，其特征在于，所述步骤S22中的具体计算方法为：
S221:对于源矩阵S当中的基因数据进行统计学假设检验，分析每个基因模块中基因的集中程度并得到P-value，对所有候选基因的P-value进行升序排列，P-value的序号记为i，并计算相应的Q值：
上式中pik代表在第k个基因模块中，升序排列后的第i个P-value，n是候选基因个数。
S222:根据基因数量规模及基因模块的个数选取阈值q_cuoff，将每个基因模块中Q值小于此阈值q_cuoff的基因分配给此基因模块。一个基因可能属于多个基因模块，即允许存在基因模块重叠现象。
6.根据权利要求1所述的一种基于模块化识别的基因网络推理方法，其特征在于，所述步骤S4中稀疏线性回归建模方法为：
S41:针对分别属于不同模块的两个基因之间的调控关系，建立如下稀疏线性回归模型：
Et＝αr1,tEr1+αr2,tEr2+…+αrt,tErt+βt
定义基因集合为G＝{g1,g2,…,gn}，其中Et代表被调控目标基因gt的数据表达向量，并且gt∈G。与基因gt不属于同一模块的潜在调控因子集合表示为G-t，其表达式为G-t＝{gr1,gr2,…,grt}，那么Er1,Er2,…,Ert为潜在调控因子G-t的表达数据向量。αrt,t代表潜在调控因子grt对目标基因gt的回归系数，βt代表回归过程中的噪声向量。
S42:采用稀疏线性回归的相关目标函数来进行α,β参数的拟合，收敛后α即为调控因子与目标基因之间的相关性得分。
7.根据权利要求1所述的一种基于模块化识别的基因网络推理方法，其特征在于，所述步骤S5所采用的相关性得分归一化方法为最大最小归一化，其计算公式为：
式中si代表调控边i所对应的相关性得分，s为相关性得分向量。
8.根据权利要求1所述的一种基于模块化识别的基因网络推理方法，其特征在于，根据步骤S2识别出的基因模块，采用Frr指标评价基因模块识别的准确性；根据步骤S5得到的总相关性得分列表，采用AUROC与AUPR指标评价基因网络推理的准确性。
9.根据权利要求8所述的一种基于模块化识别的基因网络推理方法，其特征在于，评估模块识别的准确性，采用Frr指标将识别结果与真实基因模块进行对比，计算方法如下：
其中Recovery代表真实基因模块拟合模块识别结果的程度，Relevance代表模块识别结果拟合真实基因模块的程度，具体公式如下：
其中，
和M分别代表真实基因模块集合与识别模块集合，Jaccard(·)则用来度量两个基因模块集合之间的相似度，Jaccard距离越大，真实模块与识别模块之间的相似度越低。
10.根据权利要求8所述的一种基于模块化识别的基因网络推理方法，其特征在于，评估基因调控网络推理的准确性，根据总相关性评分表绘制ROC曲线与PR曲线。给定一个分数阈值，评分表中的推理结果与真实基因调控网络对比，可以分为四类：真阳性(TruePositive,TP)，真阴性(True Negative,TN)，假阳性(False Positive,FP)，假阴性(FalseNegative,FN)。其中：
随着阈值变化，ROC曲线表现了FPR和TPR之间的变化趋势，PR曲线表示了recall和precision的变化趋势，分别计算ROC与PR的曲线下面积得到AUROC与AUPR，参数越大则表明推理结果越准确。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110309281.9A CN113066522B (zh) | 2021-03-23 | 2021-03-23 | 一种基于模块化识别的基因网络推理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110309281.9A CN113066522B (zh) | 2021-03-23 | 2021-03-23 | 一种基于模块化识别的基因网络推理方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113066522A true CN113066522A (zh) | 2021-07-02 |
| CN113066522B CN113066522B (zh) | 2022-07-12 |
Family
ID=76563089
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110309281.9A Active CN113066522B (zh) | 2021-03-23 | 2021-03-23 | 一种基于模块化识别的基因网络推理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113066522B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113947149A (zh) * | 2021-10-19 | 2022-01-18 | 大理大学 | 基因模块群的相似性度量方法、装置、电子设备及存储介质 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1452993A1 (en) * | 2002-12-23 | 2004-09-01 | STMicroelectronics S.r.l. | Method of analysis of a table of data relating to expressions of genes and relative identification system of co-expressed and co-regulated groups of genes |
| CN109409522A (zh) * | 2018-08-29 | 2019-03-01 | 浙江大学 | 一种基于集成学习的生物网络推理算法 |
| CN110111848A (zh) * | 2019-05-08 | 2019-08-09 | 南京鼓楼医院 | 一种基于rnn-cnn神经网络融合算法的人体周期表达基因识别方法 |
-
2021
- 2021-03-23 CN CN202110309281.9A patent/CN113066522B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1452993A1 (en) * | 2002-12-23 | 2004-09-01 | STMicroelectronics S.r.l. | Method of analysis of a table of data relating to expressions of genes and relative identification system of co-expressed and co-regulated groups of genes |
| CN109409522A (zh) * | 2018-08-29 | 2019-03-01 | 浙江大学 | 一种基于集成学习的生物网络推理算法 |
| CN110111848A (zh) * | 2019-05-08 | 2019-08-09 | 南京鼓楼医院 | 一种基于rnn-cnn神经网络融合算法的人体周期表达基因识别方法 |
Non-Patent Citations (1)
| Title |
|---|
| 孔薇等: "基于转录调控网络研究的乳腺癌与系统性红斑狼疮免疫系统发病机理探寻", 《青岛科技大学学报(自然科学版)》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113947149A (zh) * | 2021-10-19 | 2022-01-18 | 大理大学 | 基因模块群的相似性度量方法、装置、电子设备及存储介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113066522B (zh) | 2022-07-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| You et al. | Feature selection for high-dimensional multi-category data using PLS-based local recursive feature elimination | |
| CN109411023B (zh) | 一种基于贝叶斯网络推理的基因间交互关系挖掘方法 | |
| CN107784196B (zh) | 基于人工鱼群优化算法识别关键蛋白质的方法 | |
| CN112215259B (zh) | 基因选择方法和装置 | |
| Alagukumar et al. | Classification of microarray gene expression data using associative classification | |
| CN113066522B (zh) | 一种基于模块化识别的基因网络推理方法 | |
| CN111429970B (zh) | 基于极端梯度提升方法进行特征选择来获取多基因风险评分的方法及系统 | |
| Chowdhury et al. | UICPC: centrality-based clustering for scRNA-seq data analysis without user input | |
| Akalın et al. | Classification of exon and intron regions on dna sequences with hybrid use of sbert and anfis approaches | |
| Sato et al. | A non-parametric Bayesian approach for predicting RNA secondary structures | |
| CN113223655A (zh) | 基于变分自编码器的药物-疾病关联预测方法 | |
| Bagyamani et al. | Biological significance of gene expression data using similarity based biclustering algorithm | |
| Koumadorakis et al. | Gene Regulatory Network Reconstruction Using Single-Cell RNA-Sequencing | |
| Ghai et al. | Proximity measurement technique for gene expression data | |
| CN114238561B (zh) | 基于三元损失训练策略的生物医学实体关系抽取方法 | |
| Rehman et al. | An Optimized Neural Network with Bat Algorithm for DNA Sequence Classification | |
| Liu et al. | Function Clustering Self-Organization Maps (FCSOMs) for mining differentially expressed genes in Drosophila and its correlation with the growth medium | |
| Sun et al. | An enhanced LRMC method for drug repositioning via gcn-based HIN embedding | |
| Fathi et al. | Research Article An Efficient Cancer Classification Model Using Microarray and High-Dimensional Data | |
| Ali et al. | Evolutionary Hybrid Machine Learning Techniques for DNA Cancer Data Classification. | |
| Mostafa | Gene expression analysis using machine learning | |
| CN116153397A (zh) | 基于蛋白质/基因序列数据的生物物种同源性分析系统 | |
| Sridharan et al. | Clustering Method Analysis for Gene Expression Data using Fire Fly Optimization and Simple K-means Algorithm with Machine Learning | |
| Ding et al. | Dynamic model inference of gene regulatory network based on hybrid parallel genetic algorithm and threshold qualification method | |
| Zeng et al. | scICML: Information-theoretic Co-clustering-based Multi-view Learning for the Integrative Analysis of Single-cell Multi-omics 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 |