(12)发明专利申请
(10)申请公布号 CN 112802543 A(43)申请公布日 2021.05.14
(21)申请号 202110048797.2(22)申请日 2021.01.14
(71)申请人 东北大学
地址 110819 辽宁省沈阳市和平区文化路3
号巷11号(72)发明人 王之琼 隋玲 曲璐渲 信俊昌
王炜祎祺 李婵 殷文强 (74)专利代理机构 沈阳东大知识产权代理有限
公司 21109
代理人 梁焱(51)Int.Cl.
G16B 5/00(2019.01)G16B 45/00(2019.01)
权利要求书2页 说明书6页 附图4页
(54)发明名称
一种基于概率图的基因调控网络分析方法(57)摘要
本发明提供了一种基于概率图的基因调控网络分析方法，属于基因调控网络分析技术领域。包括：输入概率图，其每个顶点代表一个基因，顶点间的有向边代表基因间存在的调控关系，边上的概率值表示基因间调控关系的强度；在表示基因调控网络的概率图上计算各顶点的中心性，包括度中心性、紧密中心性和中介中心性，对中介中心性计算方法进行了优化；根据概率图上各顶点的中心性，筛选出关键节点；输入源顶点集合和目的顶点，根据实际需要可选择精确计算方法或者近似计算方法计算属于源顶点集合的各源顶点到目的顶点且经过关键节点的可达概率，并根据可达概率对各源顶点进行排序。该方法可以更高效、可靠地分析基因调控网络，大大提高分析结果的准确度。
CN 112802543 ACN 112802543 A
权 利 要 求 书
1/2页
1.一种基于概率图的基因调控网络分析方法，其特征在于：包括以下步骤：步骤1：输入概率图并利用该概率图表示基因调控网络；概率图的每一个顶点代表一个基因，顶点间的有向边代表基因间存在的调控关系，边上的概率值表示基因间调控关系的强度；
步骤2：在表示基因调控网络的概率图上计算各顶点的中心性，包括度中心性、紧密中心性和中介中心性；
步骤3：根据概率图上各顶点的中心性，筛选出关键节点；步骤4：输入源顶点集合S和目的顶点t，计算属于源顶点集合S的各源顶点s到目的顶点t且经过关键节点的可达概率，并根据可达概率对各源顶点s进行排序。
2.根据权利要求1所述的基于概率图的基因调控网络分析方法，其特征在于：所述步骤2包括如下步骤：
步骤2.1：对概率图进行采样，得到多个基于概率图的可能世界；
紧密中心性与中介中心性；步骤2.2：在得到的各可能世界上计算各顶点的度中心性、
步骤2.3：根据各顶点在各可能世界上的紧密中心性与中介中心性，分别在各可能世界上对所有顶点的紧密中心性和中介中心性进行排序；
步骤2.4：计算各顶点在所有可能世界上的度中心性的平均值、紧密中心性的排序次序的平均值和中介中心性的排序次序的平均值；
步骤2.5：取各顶点在所有可能世界上的度中心性的平均值、紧密中心性的排序次序的平均值和中介中心性的排序次序的平均值，分别作为概率图上对应顶点的度中心性、紧密中心性和中介中心性。
3.根据权利要求2所述的基于概率图的基因调控网络分析方法，其特征在于：在步骤2.1中所述采样的方法为：通过线性同余发生器生成随机数，用确定性算法对概率图上的每条边生成[0,1]之间的随机数后，判断当前边的存在概率是否大于该随机数，若大于则保留该条边，否则删除该条边；当所有边都生成随机数并判断是否保留后，视为完成一次采样，得到一个可能世界。
4.根据权利要求2所述的基于概率图的基因调控网络分析方法，其特征在于：在步骤2.2中，所述中介中心性的计算方法包括：首先判断待计算中心性的顶点u的出度是否为零：若是，则顶点u的中介中心性为零；若否，则将概率图的边进行方向反转，在反转后的概率图上计算与顶点u可达的顶点集合RF(u)，在概率图上查找以RF(u)集合中的顶点为根的最短路径，并根据最短路径计算每一最短路径上其他顶点到作为根的顶点的依赖分数，最后将所有的依赖分数相加得到的结果作为顶点u的中介中心性。
5.根据权利要求1所述的基于概率图的基因调控网络分析方法，其特征在于：所述步骤3包括如下步骤：
步骤3.1：分别对概率图上顶点的各中心性进行排序并用曲线表示，得到3条曲线，计算各曲线的斜率，选取各曲线上斜率趋于0的顶点的中心性作为阈值；
步骤3.2：筛选其各中心性均大于阈值的顶点，作为该概率图的关键节点。6.根据权利要求1所述的基于概率图的基因调控网络分析方法，其特征在于：在步骤4中，所述计算各源顶点s到目的顶点t且经过关键节点的可达概率的方法包括：
步骤I.1：根据步骤3筛选到的关键节点，计算关键节点之间的可达概率，且根据关键节
2
CN 112802543 A
权 利 要 求 书
2/2页
点之间的可达概率构建概率子图，并采用紧密结构存储；
步骤I.2：查找出与源顶点可达的所有关键节点，并计算源顶点到这些关键节点的可达概率；再查找出能够到达目的顶点的关键节点，并计算这些关键节点到目的顶点的可达概率；再根据关键节点概率子图查询出这些关键节点之间的可达概率；最后根据源顶点、关键节点以及目的顶点之间的可达概率，构建概率查询子图；
步骤I.3：在概率查询子图上计算各源顶点到目的顶点的可达概率。7.根据权利要求6所述的基于概率图的基因调控网络分析方法，其特征在于：步骤I.1中所述的紧密结构是一个带权有向无环图，其中顶点为关键节点，有向边为关键节点在概率图上的指向关系，边上的权值是一个数组b，数组中只存储0和1两个值，数组的长度为采样的次数，b[i]＝0表示在第i次采样中的图中该边不存在，b[i]＝1表示第i次采样中的图中该边存在。
8.根据权利要求1所述的基于概率图的基因调控网络分析方法，其特征在于：在步骤4中，所述计算各源顶点s到目的顶点t且经过关键节点的可达概率的方法包括：
步骤J.1：根据基因调控网络上关键基因的特性，将关键节点分为多个独立的关键子图，筛选到达源顶点或者目的顶点的路径距离小于路径距离阈值且可达概率大于可达概率阈值的关键节点，作为该源顶点或者目的顶点在各关键子图的代表，称为代表关键节点；
步骤J.2：计算各源顶点经过代表关键节点与目的顶点的可达概率。
3
CN 112802543 A
说 明 书
一种基于概率图的基因调控网络分析方法
1/6页
技术领域
[0001]本发明涉及基因调控网络分析技术领域，尤其涉及一种基于概率图的基因调控网络分析方法。
背景技术
[0002]某个基因的表达水平受到其他基因的影响，这个基因的表达水平又会影响其他基因的表达水平，这种基因间相互制约的调控关系构成了复杂的基因调控网络。根据数学算法和已知的经验知识发掘数据关系信息，建立基因调控网络模型，研究网络特性，认识调控
研究分析基因表达数据之间的关系，构建合适关系和机制，对生物学发展产生深远的影响。
的基因调控网络模型来模拟生物系统的行为，从中发现生物学规律，进而认识生命现象的本质，成为了生物信息学研究的重要内容。
[0003]如何精确分析基因调控网络是一项极具挑战性的任务，现有的基因调控网络分析方法都是直接对基因调控网络进行分析，或将其定义为一个确定图。由于在数据采集过程中存在固有的噪声、不完全性、时延等问题，而且基因之间的调控关系是通过实验观察到的，具有不确定性，传统的图模型无法准确描述基因之间的调控关系，从而导致基因调控网络的分析结果准确率不够。
[0004]目前有很多基于概率图的可达查询、最短路径查询等概率图查询研究成果。但是由于基因调控网络具有网络结构特殊、调控关系复杂等特点，使用现有的对普通概率图查询的算法并不能适用于基因调控网络，从而不能有效地查询出基因间的调控关系。发明内容
[0005]本发明要解决的技术问题是针对上述现有技术的不足，提供一种基于概率图的基因调控网络分析方法，考虑基因调控网络的独有的特点，提出了适用于基因调控网络的概率图查询算法，通过使用概率图定义的基因调控网络，可得到更准确的分析结果。[0006]为解决上述技术问题，本发明所采取的技术方案是：[0007]一种基于概率图的基因调控网络分析方法，包括以下步骤：[0008]步骤1：输入概率图并利用该概率图表示基因调控网络；[0009]概率图的每一个顶点代表一个基因，顶点间的有向边代表基因间存在的调控关系，边上的概率值表示基因间调控关系的强度；[0010]步骤2：在表示基因调控网络的概率图上计算各顶点的中心性，包括度中心性、紧密中心性和中介中心性；[0011]步骤3：根据概率图上各顶点的中心性，筛选出关键节点；[0012]步骤4：输入源顶点集合S和目的顶点t，计算属于源顶点集合S的各源顶点s到目的顶点t且经过关键节点的可达概率，并根据可达概率对各源顶点s进行排序。[0013]进一步地，根据所述的基于概率图的基因调控网络分析方法，所述步骤2包括如下步骤：
4
CN 112802543 A[0014][0015]
说 明 书
2/6页
步骤2.1：对概率图进行采样，得到多个基于概率图的可能世界；步骤2.2：在得到的各可能世界上计算各顶点的度中心性、紧密中心性与中介中心
性；
步骤2.3：根据各顶点在各可能世界上的紧密中心性与中介中心性，分别在各可能
世界上对所有顶点的紧密中心性和中介中心性进行排序；[0017]步骤2.4：计算各顶点在所有可能世界上的度中心性的平均值、紧密中心性的排序次序的平均值和中介中心性的排序次序的平均值；[0018]步骤2.5：取各顶点在所有可能世界上的度中心性的平均值、紧密中心性的排序次序的平均值和中介中心性的排序次序的平均值，分别作为概率图上对应顶点的度中心性、紧密中心性和中介中心性。[0019]进一步地，根据所述的基于概率图的基因调控网络分析方法，在步骤2.1中所述采样的方法为：通过线性同余发生器生成随机数，用确定性算法对概率图上的每条边生成[0,
判断当前边的存在概率是否大于该随机数，若大于则保留该条边，否则1]之间的随机数后，
删除该条边；当所有边都生成随机数并判断是否保留后，视为完成一次采样，得到一个可能世界。
[0020]进一步地，根据所述的基于概率图的基因调控网络分析方法，在步骤2.2中，所述中介中心性的计算方法包括：首先判断待计算中心性的顶点的出度是否为零：若是，则顶点的中介中心性为零；若否，则将概率图的边进行方向反转，在反转后的概率图上计算与顶点可达的顶点集合，在概率图上查找以集合中的顶点为根的最短路径，并根据最短路径计算每一最短路径上其他顶点到作为根的顶点的依赖分数，最后将所有的依赖分数相加得到的结果作为顶点的中介中心性。[0021]进一步地，根据所述的基于概率图的基因调控网络分析方法，所述步骤3包括如下步骤：
[0022]步骤3.1：分别对概率图上顶点的各中心性进行排序并用曲线表示，得到3条曲线，计算各曲线的斜率，选取各曲线上斜率趋于0的顶点的中心性作为阈值；[0023]步骤3.2：作为该概率图的关键节点。筛选其各中心性均大于阈值的顶点，[0024]进一步地，根据所述的基于概率图的基因调控网络分析方法，在步骤4中，所述计算各源顶点s到目的顶点t且经过关键节点的可达概率的方法包括：[0025]步骤I.1：根据步骤3筛选到的关键节点，计算关键节点之间的可达概率，且根据关键节点之间的可达概率构建概率子图，并采用紧密结构存储；[0026]步骤I.2：查找出与源顶点可达的所有关键节点，并计算源顶点到这些关键节点的可达概率；再查找出能够到达目的顶点的关键节点，并计算这些关键节点到目的顶点的可达概率；再根据关键节点概率子图查询出这些关键节点之间的可达概率；最后根据源顶点、关键节点以及目的顶点之间的可达概率，构建概率查询子图；[0027]步骤I.3：在概率查询子图上计算各源顶点到目的顶点的可达概率。[0028]进一步地，根据所述的基于概率图的基因调控网络分析方法，步骤I.1中所述的紧密结构是一个带权有向无环图，其中顶点为关键节点，有向边为关键节点在概率图上的指向关系，边上的权值是一个数组，数组中只存储0和1两个值，数组的长度为采样的次数，表示在第次采样中的图中该边不存在，表示第次采样中的图中该边存在。
5
[0016]
CN 112802543 A[0029]
说 明 书
3/6页
进一步地，根据所述的基于概率图的基因调控网络分析方法，在步骤4中，所述计
算各源顶点s到目的顶点t且经过关键节点的可达概率的方法包括：[0030]步骤J.1：根据基因调控网络上关键基因的特性，将关键节点分为多个独立的关键子图，筛选到达源顶点或者目的顶点的路径距离小于路径距离阈值且可达概率大于可达概率阈值的关键节点，作为该源顶点或者目的顶点在各关键子图的代表，称为代表关键节点；[0031]步骤J.2：计算各源顶点经过代表关键节点与目的顶点的可达概率。[0032]采用上述技术方案所产生的有益效果在于：本发明提供的基于概率图的基因调控网络分析方法，结合度中心性、中介中心性和紧密中心性的计算方法，在概率图上计算各节点的中心值，查找出了概率图上的关键节点。根据查找到的关键节点，提出了经过关键节点的源节点到目的节点的可达概率的计算方法，并提出了近似方法优化了可达概率的计算效率。提高了基因调控网络分析结果的准确性，并采用更高效的方法得到查询结果，在实际的药靶定位中，可以更准确地查找到对特定基因用药效果更好的靶向基因。附图说明
[0033]图1为本发明基于概率图的基因调控网络分析方法流程图；[0034]图2为本发明方法中概率图上关键节点查找方法的流程图；[0035]图3为本发明方法中基因对调控关系精确查询方法的流程图；[0036]图4为本发明实方法中基因对调控关系近似查询方法的流程图。
具体实施方式
[0037]下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。[0038]本实施例中，基于概率图的基因调控网络分析方法，如图1所示，包括以下步骤：[0039]步骤1：输入概率图并利用该概率图表示基因调控网络；[0040]概率图的每一个顶点代表一个基因，顶点间的有向边代表基因间存在的调控关系，边上的概率值表示基因间调控关系的强度。[0041]概率图是其边上的权值表示存在概用概率图对基因调控网络建模具有重要意义。率的图。概率图可以用一个确定图并在边上增加存在概率来表示。与普通的图结构相比，概率图由于其具有概率特征的独特结构，使其可以模拟许多现实世界中的情景，通过分析这些概率图数据，可以找到许多更有价值的信息。比如推断基因之间的相似性，哪些基因对其他基因的影响特别大，某基因组对另一基因组的调控关系等等。[0042]步骤2：在表示基因调控网络的概率图上计算各顶点的中心性，如图2所示，计算过程如下：
[0043]步骤2.1：对概率图进行采样，得到多个基于概率图的可能世界。[0044]所述可能世界为概率图的一个实例，一个M条边的概率图有2M个可能世界。本实施方式采用蒙特·卡罗(Monte Carlo)方法进行采样：通过线性同余发生器生成随机数，用确定性算法对概率图上的每条边生成[0,1]之间的随机数后，判断当前边的存在概率是否大于该随机数，若大于则保留该条边，否则删除该条边。当所有边都生成随机数并判断是否保留后，则为一次采样，得到一个可能世界。实施例一中，选取乳腺癌相关基因调控网络，表示
6
CN 112802543 A
说 明 书
4/6页
该基因调控网络的概率图大小为顶点数V＝547，边数E＝3503，设定采样的次数N＝10000，采样次数由技术人员根据概率图的大小和实际经验进行确定，实施例一中经过10000次采样后，得到10000个可能世界。[0045]步骤2.2：在得到的各可能世界上计算各顶点的度中心性、紧密中心性与中介中心性，得到各可能世界上各顶点的中心值。
[0046]度中心性以顶点的直接邻居数量作为量度，认为直接邻居数量越多，顶点越重要。顶点u的度中心性归一化表示为：
[0047]
式中k(u)为顶点u的入度和出度和，n为顶点总数。
[0049]紧密中心性计量一个顶点到所有其他顶点的紧密性，即顶点距离的倒数，一个拥有高紧密性中心性的顶点拥有到其他顶点的距离最小值。常用的紧密中心性采用归一化的中心性公式进行计算，即计算顶点u到其他顶点的平均距离的倒数：
[0050][0051]
[0048]
式中v为概率图G中除u的其他任意顶点，d(u,v)表示顶点u到顶点v的最短路径距
离。
中介中心性用于寻找连接图的两个部分的桥梁顶点。本实施方式中，中介中心性的计算方法包括：首先判断待计算中心性的顶点u的出度是否为零：若是，则顶点u的中介中心性为零；若否，则将概率图的边进行方向反转，在反转后的概率图上采用BFS(Breadth First Search，广度优先搜索)或DFS(Depth First Search，深度优先搜索)方法计算与顶点u可达的顶点集合RF(u)，查找以RF(u)集合中的顶点为根的最短路径，并根据最短路径计算每一最短路径上其他顶点到作为根的顶点的依赖分数，最后将所有的依赖分数相加得到的结果作为顶点u的中介中心性。[0053]在本实施方式中，采用随机近似(Randomized‑Approximate Brandes)方法计算依赖分数。该方法定义了任意起始顶点v对另一个顶点u的依赖分数(dependency score)为：
[0054][0055][0056][0052]
任何顶点的中介中心性都可以用依赖分数来表示：
式中σσvt表示顶点v到顶点t之间最短路径数量，vt(u)是其中经过顶点u的路径数量，V(G)是概率图G的顶点集合。[0058]步骤2.3：根据各顶点在各可能世界上的紧密中心性与中介中心性，分别在各可能世界上对所有顶点的紧密中心性和中介中心性进行排序，可以根据实际需求进行降序或者升序等排序方式。[0059]步骤2.4：计算各顶点在所有可能世界上的度中心性的平均值、紧密中心性的排序次序的平均值和中介中心性的排序次序的平均值。
7
[0057]
CN 112802543 A[0060]
说 明 书
5/6页
步骤2.5：取各顶点在所有可能世界上的度中心性的平均值、紧密中心性的排序次
序的平均值和中介中心性的排序次序的平均值，分别作为概率图上对应顶点的度中心性、紧密中心性和中介中心性。[0061]步骤3：根据概率图上各顶点的中心性，筛选出关键节点；[0062]步骤3.1：分别对概率图上顶点的各中心性进行排序并用曲线表示，得到3条曲线，计算各曲线的斜率，选取各曲线上斜率趋于0的顶点的中心性作为阈值；[0063]在实施例一中，度中心性阈值(threshold‑d)为0.0135，紧密中心性阈值(threshold‑c)为0.0014，中介中心性阈值(threshold‑b)为0.2679。[0064]步骤3.2：筛选其各中心性均大于阈值的顶点，作为该概率图的关键节点；[0065]在实施例一中，筛选同时满足度中心性大于度中心性阈值0.0135、紧密中心性大于紧密中心性阈值0.0014、且中介中心性大于中介中心性阈值0.2679的顶点作为关键节点，共筛选出34个乳腺癌相关基因调控网络的关键基因，如表1所示：[0066]表1乳腺癌相关基因调控网络的关键基因表
COL1A1MMP2VCANSRPX2NEK2CEP55BUB1BFAPKIF2CCOL8A2PRRX1RACGAPTOP2ACOL10A1UBE2CKIF4ACCDC80DACT1LOXBNC2MIR100HCOL5A1COL6A3FSTL1PTTG1DLGAP5TPX2KIF11COL5A2MSRB3CCNB2CDCA3DNM3OSESPL1 [0068]步骤4：输入源顶点集合S，目的顶点t，计算属于源顶点集合S的各源顶点s到目的顶点t且经过关键节点的可达概率，并根据可达概率对各源顶点s进行排序。可以根据实际需求对各源顶点s进行降序、升序或者其他种类的排序。
[0069]本实施方式提供两种方法计算各源顶点s到目的顶点t且经过关键节点的可达概率，本领域技术人员可以选择方法一或者方法二计算各源顶点s到目的顶点t且经过关键节点的可达概率。[0070]利用方法一可以精确计算各源顶点s到目的顶点t且经过关键节点的可达概率，如图3所示，方法一包括如下步骤：[0071]步骤I.1：根据步骤3筛选到的关键节点，计算关键节点之间的可达概率，且根据关键节点之间的可达概率构建概率子图，并采用紧密结构存储；[0072]所述的紧密结构是一个带权有向无环图，其中顶点为关键节点，有向边为关键节点在概率图上的指向关系，边上的权值是一个数组b，数组中只存储0和1两个值，数组的长度为采样的次数，b[i]＝0表示在第i次采样中的图中该边不存在，b[i]＝1表示第i次采样中的图中该边存在。[0073]步骤I.2：查找出与源顶点可达的所有关键节点，并计算源顶点到这些关键节点的可达概率；再查找出能够到达目的顶点的关键节点，并计算这些关键节点到目的顶点的可达概率；最后根据关键节点概率子图查询出这些关键节点之间的可达概率。根据源顶点、关键节点以及目的顶点之间的可达概率，构建概率查询子图。[0074]步骤I.3：在概率查询子图上计算各源顶点到目的顶点的可达概率；[0075]采用递归分层抽样法(Recursive Stratified Sampling)计算概率图的可达概率：选择r条边(e1,...er)，并决定它的状态(0or 1)，对于剩下的概率图的总边数‑r条边，设
8
[0067]
CN 112802543 A
说 明 书
6/6页
置这些边的状态为*，表示“它们的状态未知”，则将整个可能世界空间Ω分为2r个子空间
使Xi＝(Xi,1,Xi,2,...,Xi,j,...,Xi,r)表示在第i层选择第r条边的状态向量，在第
i层的可能世界概率为：
[0076]
式中GP是Ωi中的一个可能世界，pj是边ej在概率图G上的存在概率。[0078]在实施例一中，根据输入的概率图的尺寸设定选取的边数r为10。[0079]利用方法二可以近似计算各源顶点s到目的顶点t且经过关键节点的可达概率，在保证计算结果准确率的情况下，提高计算效率，如图4所示，方法二包括如下步骤：[0080]步骤J.1：根据基因调控网络上关键基因的特性，将关键节点分为多个独立的关键子图，在概率图上查找与源顶点或者目的顶点可达的关键节点的过程中，根据概率图的大小选择合适的路径距离阈值和可达概率阈值，筛选到达源顶点或者目的顶点的路径距离小于路径距离阈值且可达概率大于可达概率阈值的关键节点，作为该源顶点或目的顶点在各关键子图的代表，称为代表关键节点。[0081]例如，在实施例一中，根据经验设定路径距离阈值λ为15，可达概率阈值β为0.5，则筛选到达源顶点或目的顶点的距离小于15且可达概率大于0.5的关键节点作为源顶点或目的顶点在各关键子图的代表，称为代表关键节点。[0082]步骤J.2：由于各代表关键节点之间相互独立，则直接用概率公式计算各源顶点经过代表关键节点与目的顶点的可达概率。
[0083]计算源顶点s到目的顶点t之间拥有m条独立路径的可达概率公式为：
[0084]
[0077]
式中pi为第i条独立路径的存在概率，即，组成该路径的边概率的乘积。
[0086]最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。
[0085]
9
CN 112802543 A
说 明 书 附 图
1/4页
图1
10
CN 112802543 A
说 明 书 附 图
2/4页
图2
11
CN 112802543 A
说 明 书 附 图
3/4页
图3
12
CN 112802543 A
说 明 书 附 图
4/4页
图4
13
因篇幅问题不能全部显示，请点此查看更多更全内容