Identifying Candidate Disease Genes in Multilayer Heterogeneous Biological Networks
-
摘要: 现有大多数用于识别候选疾病基因的随机游走方法通常优先访问高度连接的基因, 而可能与已知疾病有关的不知名或连接性差的基因易被忽略或难以识别. 此外, 这些方法仅访问单个基因网络或各种基因数据的聚合网络, 导致偏差和不完整性. 因此, 设计一种能控制随机游走运动方向和整合多种数据源的候选疾病基因识别方法将是一个迫切需要解决的问题. 为此, 首先构建多层网络和多层异构基因网络. 然后, 提出一种游走于多层网络和多层异构网络的拓扑偏置重启随机游走(Biased random walk with restart, BRWR)算法来识别疾病基因. 实验结果表明, 游走于不同类型网络上的识别候选疾病基因的BRWR算法优于现有的算法. 最后, 应用于多层异构网络上的BRWR算法能预测未诊断的新生儿类早衰综合征中涉及的疾病基因.Abstract: Most of existing random walk methods to identify candidate disease genes preferentially visit highly-connected genes, while unwell-known or poorly-connected genes probably relevant to known diseases are more easily ignored or complicated to identify. Moreover, these methods access only a single gene network or an aggregated network of various gene data, leading to bias and incompleteness. Therefore, it is a pressing challenge for controlling the motion direction of random walk and for integrating multiple data sources involving different information for disease-gene identification. To this end, we first construct a multilayer network and multilayer heterogeneous genetic network. Then, we propose a topologically biased random walk with restart (BRWR) algorithm applicable to multilayer and multilayer heterogeneous networks for the identification of candidate disease genes. Experimental results show that the BRWR algorithm to identify candidate disease genes outperforms the state-of-the-art ones on different types of networks. Finally, the BRWR algorithm on multilayer heterogeneous networks is used to predict disease genes implicated in the undiagnosed neonatal progeroid syndrome.
-
近年来, 生物大分子之间物理和功能相互作用的识别是生物信息学中最活跃的研究领域之一[1]. 例如, 研究者所筛选出的多种生物体中蛋白质−蛋白质相互作用(Protein-protein interaction, PPI)[2-3]可揭示数千种蛋白质之间功能或物理相互作用. 这些生物分子间的相互作用通常可以表示为其节点是蛋白质或基因, 边是物理或功能相互作用的生物网络. 大量生物网络研究表明, 在物理或功能上彼此接近的基因往往参与相似的生物学路径, 并对候选疾病基因的识别起关键作用. 基于此, 研究者提出了许多识别生物网络中候选疾病基因的方法[4-9], 旨在能够有效准确地识别出与疾病基因相关的候选基因, 进而预测基因和疾病间的关系, 找出疾病的发病机制. 例如, Chen等[4]提出一种最大化网络信息流(Maximize the information flow, MAXIF)的计算方法以发现与疾病关联的基因. Zhang等[5]提出一种改进的双标签传播算法(Improved dual label propagation, IDLP)来序列化候选疾病基因. Chen等[6]提出整合人类组学数据的基于回归模型基因识别方法(Based on regression to identify disease genes, BRIDGE). Lee等[7]开发一种新的生物学系统来揭示转移性乳腺癌的突变驱动因素. Yang等[8]提出一种药物−疾病关联的异构网络嵌入方法来预测药物与疾病间可能存在的关联. Yang等[9]提出网络增强方法对变体进行基因标注. 此外, 生物信息学中用于生物计算的关联推断最常用的重启随机游走(Random walk with restart, RWR)方法[10-15]已广泛应用于新的疾病关联基因的识别.
实际上, 在候选疾病基因的识别方法中, 探索生物网络的重启随机游走方法致力于获得与疾病基因关联的节点中心性排名以识别出新的致病基因, 而中心性排名是复杂网络中关键或重要性成份的排名(本文为基因或疾病节点的排名), 所谓重要性节点是指相比网络中的其他成份, 对网络结构和功能有更大程度影响的节点. 例如, Köhler等[10]在一个PPI网络上使用RWR来获得候选疾病基因中心性排名进而识别出关键的候选疾病基因. Li等[11]扩展RWR算法到由一个PPI网络和一个表型网络(Phenotype network)构成的异构网络上以推断基因与疾病的关系. Valdeolivas等[16]将RWR扩展到多重和异构网络上旨在探索基因和蛋白质之间的物理和功能相互作用.
目前, RWR方法被认为是识别疾病关联基因中心性排名中常用方法之一, 它能够度量出种子节点和网络中所有其他节点之间的距离或接近度(即与种子节点关联的节点中心性排名), 通常认为该方法是生物信息学中用于生物计算的关联推断最常用的方法[16-17]. 在生物网络中, 根据与已知疾病相关联的种子节点的接近度, RWR算法可对网络中所有的节点进行中心性排名. 各种扩展的RWR方法得以提出, 其目的在于提高候选基因序列化预测的准确性. 例如, Li等[13]提出一种融合多个异构基因数据和表型数据的多图网络的RWR方法, 该方法具有较强的噪声容错功能, 从而提高疾病基因排序的性能.
然而, 现有RWR方法存在两个明显的缺陷: 1) RWR方法通常偏向高度连接的基因或蛋白质[18](基因即为蛋白质). 在这种情况下, 可能与已知疾病有关的不知名或连接度不高的蛋白质易被忽略或难以被检测. 2) RWR方法常用于单个基因网络或一个由不同基因数据构成的聚合网络. 由于单个或聚合数据源容易产生噪声、偏差和不完整性, 因此需要设计一种能集成各种基因组数据源的RWR方法, 使之能够可靠地识别出与给定疾病关联的候选基因. 即便RWR已经用于探索聚合了基因和表型数据的多图网络, 但它却忽视了蛋白质间边的类型、网络的拓扑结构和网络间相互连接等重要信息.
为解决前一问题, 本文将拓扑偏置随机游走[19-20]引入RWR方法中, 以便于灵活地识别出与种子疾病弱连接或高度连接的关联基因. 为此, 本文提出一个通用的拓扑偏置随机游走方法, 该方法允许随机游走偏置目标节点的各种拓扑属性(如节点的度, 强度或聚类系数). 特别是, 拓扑偏置的游走可表示为马尔科夫过程, 其转移概率是偏置游走参数的函数. 在此情况下, 通过调整偏置函数的参数, 可以使游走者优先访问或避开具有高拓扑属性值的节点. 为解决后一问题, 受文献[21]启发, 一个较好的策略是将基于拓扑偏置随机游走的RWR方法应用于能描述多种基因相互作用关系的多网络[22-26]和多层异构生物网络. 借助于基因多网络, 基于拓扑偏置随机游走的RWR方法能更准确可靠地识别出与种子疾病关联的候选基因序列, 从而缓解单个或聚合的数据源上候选疾病基因的识别所产生噪声、偏差和不完整性等问题.
为此, 本文构造了一个由PPI、共表达(Co-expression, COEX)和路径(Pathway, PATH)网络组成的多层基因网络(如图1(a)), 由一个表型(即疾病)网络和基因网络通过基因−表型关联关系组成的二分关系网络(即异构网络, 如图1(b)). 实际上, 一个异构网络是具有特定类型节点和边的两个单网络组成的多层网络, 其中两个单网络通过二分相互关系连接(二分相互关系表示为两种不同类型网络中节点的连接关系[27], 例如基因网、表型网络和基因−表型二分相互关系网络可以构成一个异构网络, 如图1(b)). 借助这些网络, 本文还构造了一个多层异构网络(如图1(c)), 它是由一个多层基因网络、表型网络和基因−表型二分关系网络组成. 多层基因网络中每一层通过基因−表型二分关联关系与一个表型网络连接起来构成一个异构网络, 这样多层异构网络就是一个由异构网络组成的多层网络. 然后, 本文提出一个通用的拓扑偏置随机游走方法, 并将其引入RWR, 称为拓扑偏置重启随机游走(Biased random walk with restart, BRWR). 通过调整随机游走中的偏置参数, BRWR方法可以灵活地获得与种子基因关联的候选基因中心性排名, 从而可靠地识别出一组候选疾病基因, 这些基因要么是高度连接, 要么是连接不好或两种情况兼之. 此外, 本文分别将BRWR方法应用于基因多网络(BRWR applied to gene multiplex network, BRWR-M)和多层异构网络(BRWR applied to gene multilayer heterogeneous network, BRWR-MH)来识别出更可信的候选疾病基因. 实验中, 本文设计留一交叉验证策略来评估BRWR方法在不同类型网络(包括单层网络、聚合网络、异构网络、多层网络和多层异构网络)上的性能. 与目前候选疾病基因中心性排名方法相比, 选择合适的偏置参数且同时考虑各种相互作用源的BRWR方法能更准确地识别出候选疾病基因. 最后, 通过使用不同的偏置参数, 本文应用BRWR-MH方法来探索隐藏在新生儿早衰症(Neonatal progeroid syndrome, NPS)中的候选致病基因, 并展示其网络近邻表示.
图 1 多层网络、异构网络、多层异构网络以及探索它们的随机游走路径(箭头的实线)的示意图Fig. 1 Schematic of multilayer, heterogeneous and multilayer heterogeneous networks, together with paths of random walks (arrow solid lines)本文的贡献如下:
$ 1) $ 构建了一个由多层基因网络、表型网络和基因−表型网络组成的多层异构网络.
$ 2) $ 提出了一个融入拓扑随机游走的BRWR方法, 并将BRWR方法分别应用到多层网络和多层异构网络中来确定候选疾病基因中心性排名.
$ 3) $ 设计留一交叉验证方法来评估BRWR方法在不同类型网络上的性能, 揭示了BRWR-MH方法胜过现有方法, 且是一个随参数适度变化的健壮方法.
$ 4) $ 应用BRWR-MH方法来预测未确诊的新生儿早衰症所涉及的未知疾病基因, 并展示其网络表示.
1. 研究方法
本节首先介绍拓扑偏置随机游走方法; 然后, 详述单层网络、多层网络和多层异构网络上拓扑偏置重启随机游走BRWR方法.
1.1 拓扑偏置随机游走
受文献[19, 28]启发, 本文提出游走在基因网络上的拓扑偏置随机游走的一般表达: 位于基因节点$ i $上的随机游走者(即粒子)在每个时间步移向它的一个直接邻居节点$ j $, 其移动概率比例于节点$ j $的拓扑属性$ f_j $的$ b\in{{\bf{R}}} $次幂. 在此情景下, 随机游走的转移概率定义为
$$ \begin{equation} \begin{aligned} \pi_{ij} = \frac{g_{ij}(f_j)^{b}}{\sum\limits _{k = 1}^{N}g_{ik}(f_k)^{b}} \end{aligned} \end{equation} $$ (1) 其中, $ g_{ij} $表示具有$ N $个节点的基因网络所对应邻接矩阵的节点$ i $和$ j $的一个元素, $ \pi_{ij} $表示位于基因节点$ i $上的粒子跳跃到其邻居$ j $的概率. 注意, 概率$ \pi_{ij} $既依赖于网络拓扑也依赖于随机游走中偏置的类型. 这是个通用的定义, 因为它适用于节点的各种拓扑属性, 如节点的度、强度、聚类系数、多参与系数或其他拓扑属性. 本文$ f_j $表示到达基因节点$ j $的度. 偏置指数$ b $是一个控制参数, 可以调整随机游走者对目的节点拓扑属性的依赖性. 当$ b>0 $时, 随机游走者将会优先访问具有较大值$ f $的节点. 当$ b<0 $时, 它倾向于更频繁地移向具有较小值$ f $的节点. 当$ b = 0 $时, 它将恢复到无偏(即经典)的情况. 随后, 本文将拓扑偏置随机游走引入到重启随机游走RWR中, 旨在定义新颖的RWR以度量种子节点和网络中所有其他节点之间的中心性排名.
1.2 基于BRWR的单层网络基因识别
RWR是一种常见的排名方法, 源于经典的PageRank[29], 它引入重启概率. 随机游走者在每个时间步骤上或者从当前节点迭代地移动到其直接邻居, 或者重启在指定的种子节点上. 对于它的每次迭代, 随机游走者以预先设定的概率$ r\in(0,1) $[30]在种子节点上重启. 这确保随机游走者不陷入孤立节点而能到达稳态概率分布. 在此情景下, 游走者将探索整个网络, 从而获得种子近邻的中心性排名. 稳态概率分布给出种子与网络中所有其他节点之间的中心性排名.
事实上, 经常用来获得候选基因中心性排名的RWR方法通常优先访问高度连接节点[18]. 生物网络中高度连接节点(即中心)可以是细胞中高度连接的蛋白质或基因. 在这种情况下, 那些与已知疾病相关的不知名或连接度差的蛋白质, 比高度连接的蛋白质更难识别. 为解决此问题, 将拓扑偏置随机游走引入到RWR中, 这样的随机游走即拓扑偏置重启随机游走(BRWR). 在此情况下, 通过调整BRWR中的偏置参数, 游走者也可以跳跃到较少连接或不知名的节点. 这里设定$ p_i(t) $是给定时间步骤$ t \in{\bf{R}} $下的节点$ i $上的时间相关概率变量, 并且$ {\boldsymbol{\pi}}_g = [\pi_{ij}]_{i,j = 1}^N $ ($ \pi_{ij} $已经在式(1)中定义)表示与基因网络相关的转移矩阵. 因此, 时间$ t $~$ t+1 $的基因网络下的BRWR主方程定义为
$$ \begin{equation} {\boldsymbol{p}}(t+1) = (1-r){\boldsymbol{\pi}}_g {\boldsymbol{p}}(t)+r {\boldsymbol{p}}_0 \end{equation} $$ (2) 其中, ${\boldsymbol{p}}{{(t)}} = [ p_1(t), p_2(t), \cdots, p_N(t)]$表示在时间步骤$ t $下的$ N $个元素的向量, $ {\boldsymbol{p}}_0 $是种子节点不等于零的初始概率向量, $ r\in(0,1) $是随机游走重启概率.
经过多次迭代后, $ {\boldsymbol{p}}(t+1) $和$ {\boldsymbol{p}}(t) $之间的差异变得可以忽略不计, 概率达到稳定状态, 本文假定当$ {\boldsymbol{p}}(t+1) $和$ {\boldsymbol{p}}(t) $间直接差异低于$ 10^{-6} $时, 即达到稳定概率. 稳定概率$ {\boldsymbol{p}}(\infty) $表示节点$ i $到给定种子节点的中心性排名. 注意, 如果$ {\boldsymbol{p}}_i(\infty)>{\boldsymbol{p}}_j(\infty) $, 那么节点$ i $比节点$ j $更接近种子节点.
对于全局重启概率, 根据文献[15-16, 31], 本文将同样采用$ r = 0.7 $的经验值. 在下面的BRWR方法的所有版本中, 重启概率将保持不变.
1.3 基于BRWR的多层网络基因识别
在本节中, 将BRWR引入到多层基因网络中以确定候选疾病基因中心性排名. 通常, 在下一非重启步骤中, 相互连接多层网络上的随机游走者可有两种行走情况: 1)随机游走者从当前节点跳转到同一层的一个直接邻接节点; 2)从当前节点跳跃至其他层中对应的节点, 如图1(a)所示. 事实上, 在多层网络上的随机游走者并不局限于单一层上的运动, 它可以利用所有层中指向当前节点的所有连接进行随意移动.
这里考虑一种特定类型的相互连接的多层网络, 称为节点对齐和连接的多网络[32]. 基因多网络是由$ L $个单层基因网络组成, 其中每层分享着一样的基因节点集$ V = \{1,2,\cdots,N\} $, 而边集则取决于具体的层. 每一单层基因网络$ \alpha \in \{1,2,\cdots,L\} $对应的邻接矩阵为$ {\boldsymbol{A}}^\alpha = [a_{ij}^\alpha]_{N\times N} $. 此外, 层间连接所对应的矩阵定义为${\boldsymbol{D}}^{\alpha\beta} = [d_{ii}^{\alpha\beta}]_{N\times N}$, 它表示层$ \alpha $中节点和层$ \beta $中同样节点之间的连接, 其中$ d_{ii}^{\alpha \beta}\in [0,+\infty) $量化了从层$ \alpha $中节点$ i $游走到层$ \beta $中节点$ i $的“代价”. 即一个层间连接基因多网络可表示为一个超级邻接矩阵$ {\boldsymbol{A}}_M = [A_M(ij)] \in {\bf{R}}^{NL\times NL} $, 即表示为
$$ \left[\begin{array}{cccc} (1-\delta){\boldsymbol{A}}^1 &\delta {\boldsymbol{D}}^{12} &\cdots &\delta {\boldsymbol{D}}^{1L}\\ \delta {\boldsymbol{D}}^{21} &(1-\delta){\boldsymbol{A}}^{2} &\cdots &\delta {\boldsymbol{D}}^{2L}\\ \vdots &\vdots &\ddots &\vdots\\ \delta {\boldsymbol{D}}^{L1} &\delta {\boldsymbol{D}}^{L2} &\cdots &(1-\delta){\boldsymbol{A}}^{L} \end{array}\right] $$ (3) 其中, $ \delta\in[0,1] $控制随机游走在不同层之间跳跃或停留在同一层中的概率. 如果$ \delta = 0 $, 那么随机游走者将会仅在同一层内移动. 通过用$ A_M(ij) $取代式(1)中的成份$ g_{ij} $, 这样就能得到多层网络上随机游走者的转移概率矩阵$ {{\boldsymbol{\pi}}}_M = [\pi_{ij}]_{i,j = 1}^{NL} $. 因此, 可以自然地将偏置重启随机游走BRWR方法引入多层网络(BRWR-M), 那么式(2)将变为
$$ \begin{equation} {\boldsymbol{p}}_M(t+1) = (1-r){{\boldsymbol{\pi}}}_M {\boldsymbol{p}}_M(t)+r {\boldsymbol{p}}_{M0} \end{equation} $$ (4) 其中, $ {\boldsymbol{p}}_M(t) = [{\boldsymbol{p}}_1(t), {\boldsymbol{p}}_2(t),\cdots,{\boldsymbol{p}}_L(t)] $是一个具有$ NL $成份的超向量, 它描述多层网络随机游走的概率分布, 这些向量由每一层上节点的概率分布组成. 重启向量$ {\boldsymbol{p}}_{M0} $描述了每个层的初始概率分布, 这里设置为$ {\boldsymbol{p}}_{M0} = \tau {\boldsymbol{p}}_0 $, 其中参数$ {{\boldsymbol{\tau}}} = [\tau_1,\tau_2,\cdots,\tau_L] $量化了多层网络中每一层种子节点的重启概率. 重要的是, 可以通过修改参数$ {\boldsymbol{\tau}} $来调整每一层的重要性.
本文考虑为每一层构建相同的重新启动概率$ {{\boldsymbol{\tau}}} = (1/L,1/L,\cdots,1/L) $, 也为随机游走在不同的层间或在同一层内跳跃设置相等的概率$ \delta = 0.5 $. 对于处于层$ \alpha $中节点$ i $和层$ \beta $中节点$ i $之间的层间连接参数$d_{ii}^{\alpha\beta}$, 将同样的值$ d_x = 1 $赋值给所有的$d_{ii}^{\alpha\beta}$. 参数$ d_x $对BRWR性能的影响将在实验部分进行分析.
1.4 基于BRWR的多层异构网络基因识别
受文献[11]方法的启发, 本文将拓扑偏置重启随机游走BRWR方法引入到多层异构网络, 初始假定随机游走者位于多层网络中某一层的一个节点上. 那么, 在下一非重启步骤, 随机游走者可以: 1)在同一层内移动; 2)或者跳到其他层中其对应节点; 3)或者借助基因−表型关联关系跳到另一个网络, 其具体动态过程如图1(c)所示.
这里考虑一个由$ L $层和$ N $个基因节点组成的基因多网络, 其对应的超级邻接矩阵$ {\boldsymbol{A}}_{M(NL\times NL)} $在式(3)中已经定义. 本文也考虑具有$ M $个表型的表型网络, 其对应邻接矩阵为$ {\boldsymbol{A}}_{P(M\times M)} $, 并且构建了超级二分关联矩阵. 根据基因−表型关联关系, 对应邻接矩阵为$ {\boldsymbol{B}}_{(N\times M)}^{1,\cdots,L} $的多层二分网络将基因多层网络每一层中的基因节点与表型节点关联起来. 事实上, 多层二分网络是由$ L $个不相同的二分网络组成. 因此, 可以用$ {\boldsymbol{B}}_{(N\times M)} $来描述所有的网络, 并且通过堆叠$ L $个二分网络$ {\boldsymbol{B}}_{(N\times M)} $来创建一个与多层二分网络相关联的超级二分关联矩阵$ {\boldsymbol{B}}_{MP} = [{\boldsymbol{B}}_{(N\times M)}, {\boldsymbol{B}}_{(N\times M)}, \cdots, {\boldsymbol{B}}_{(N\times M)}]^{\rm{T}} $.
接下来, 本文构建多层异构网络所对应的全局邻接矩阵, ${{\boldsymbol{A}}_{MH}} = \left[ \begin{aligned} {{\boldsymbol{A}}_M},{{\boldsymbol{B}}_{MP}}\\ {{\boldsymbol{B}}_{PM}},{{\boldsymbol{A}}_P} \end{aligned} \right]$, ${\boldsymbol{B}}_{PM}$是矩阵$ {\boldsymbol{B}}_{MP} $的转置. 这里同样需要计算BRWR在不同网络上的转移概率, 以便将其引入到多层异构网络中. 让${\boldsymbol{R}} = \left[ \begin{aligned} {{\boldsymbol{R}}_{MM}},{{\boldsymbol{R}}_{MP}}\\ {{\boldsymbol{R}}_{PM}},{{\boldsymbol{R}}_{PP}} \end{aligned} \right]$ 表示多层异构网络所对应的转移矩阵, 其中$ {\boldsymbol{R}}_{MM} $和$ {\boldsymbol{R}}_{PP} $分别表示位于基因多网络和表型网络内的随机游走的转移概率; $ {\boldsymbol{R}}_{MP} $和$ {\boldsymbol{R}}_{PM} $描述基因多网络和表型网络之间的随机游走的跳跃概率. 对于一个给定的节点, 如果存在二分基因−表型关联关系, 那么随机游走者将在基因多网络和表型网络之间跳跃, 或者以跳跃概率$ \lambda \in [0,1] $停留在当前网络. $ \lambda $控制着多层异构网络中两个子网络的强化作用. 因此, $ \lambda $值越大, 网络间跳跃概率就越高; 如果$ \lambda = 0 $, 基因和表型独立排名.
假定随机游走者位于多层异构基因网络中某一层的一个基因节点$ g_i $, 下一时间步骤随机游走者或者以如下的转移概率跳跃到基因节点$ g_j $:
$$ \begin{split} &{\boldsymbol{R}}_{MM}(i,j) = \\ &\;\;\;\left\{ \begin{aligned} &{\frac{{\boldsymbol{A}}_{M}(i,j){(f_j)^{b_{M}}}}{\sum \limits_{k = 1}^{NL}\boldsymbol{A}_{M}(i,k)(f_k)^{b_{M}}}},&\;\;\sum\limits_{k = 1}^{M}{\boldsymbol{B}}_{MP} (i,k) = 0 \\ &{\frac{(1-\lambda){\boldsymbol{A}}_{M}(i,j){(f_j)^{b_{M}}}}{\sum\limits_{k = 1}^{NL}{\boldsymbol{A}}_{M}(i,k)(f_k)^{b_{M}}}}, &\quad{{\rm{否则}}}\qquad\qquad \qquad \end{aligned} \right. \end{split} $$ (5) 或者通过一个基因−表型关联关系以概率$ {\boldsymbol{R}}_{MP}(i,j) $跳跃到一个表型节点$ p_j $上:
$$ \begin{split} &{\boldsymbol{R}}_{MP}(i,j) =\\ &\;\;\;\left\{ \begin{aligned} &{\frac{\lambda {\boldsymbol{B}}_{MP}(i,j){(f_j)^{b_{MP}}}}{\sum\limits_{k = 1}^{M}{\boldsymbol{B}}_{MP}(i,k)(f_k)^{b_{MP}}}}, &{{\rm{}} \sum\limits_{k = 1}^{M} {\boldsymbol{B}}_{MP} (i,k)\neq0} \\ &{0}, &\quad{{\rm{否则}}} \;\qquad\qquad\quad\quad \end{aligned} \right. \end{split} $$ (6) 类似地, 当随机游走者位于表型网络中一个表型节点$ p_i $时, 它能以$ {\boldsymbol{R}}_{PP}(i,j) $概率移动到表型网络中另一节点$ p_j $:
$$ \begin{split} & {\boldsymbol{R}}_{PP}(i,j) = \\ &\;\;\;\left\{ \begin{aligned} & {\frac{{\boldsymbol{A}}_{P}(i,j){(f_j)^{b_{P}}}}{\sum\limits_{k = 1}^{M}{\boldsymbol{A}}_{P}(i,k)(f_k)^{b_{P}}}},&{}\;\;\sum\limits_{k = 1}^{NL} {\boldsymbol{B}}_{PM} (i,k) = 0 \\ &{\frac{(1-\lambda){\boldsymbol{A}}_{P}(i,j){(f_j)^{b_{P}}}}{\sum\limits_{k = 1}^{M}{\boldsymbol{A}}_{P}(i,k)(f_k)^{b_{P}}}}, &\;\;{{\rm{否则}}} \qquad\qquad\quad\quad \end{aligned} \right. \end{split} $$ (7) 或者以$ {\boldsymbol{R}}_{PM}(i,j) $概率跳跃到一个基因节点$ g_j $:
$$ \begin{split} &{\boldsymbol{R}}_{PM}(i,j) =\\ &\;\;\;\left\{ \begin{aligned} &{\frac{\lambda {\boldsymbol{B}}_{PM}(i,j){(f_j)^{b_{PM}}}}{\sum\limits_{k = 1}^{NL}{\boldsymbol{B}}_{PM}(i,k)(f_k)^{b_{PM}}}},&\;\;\sum\limits_{k = 1}^{NL} {\boldsymbol{B}}_{PM} (i,k)\neq0 \\ &{0}, &\quad{{\rm{否则}}} \;\;\;\qquad\qquad\quad \end{aligned} \right. \end{split} $$ (8) 其中, 与每个子网络(即多层、基因和基因−表型网络)相关的$ b_M $, $ b_P $, $ b_{MP} $和$ b_{PM} $是控制参数, 允许调整随机游走者对目标节点拓扑属性的依赖性. $ f_j $可表示目标节点的度、强度或聚类系数等拓扑属性.
因此, 在多层异构网络上BRWR-MH方法定义为
$$ \begin{equation} {\boldsymbol{p}}_{MH}(t+1) = (1-r){\boldsymbol{R}} {\boldsymbol{p}}_{MH}(t)+r {\boldsymbol{p}}_{MH0} \end{equation} $$ (9) 由于BRWR-MH方法同时对$ L $个不同层的$ N $个基因节点和表型网络的$ M $个表型进行排序, 因此向量$ {\boldsymbol{p}}_{MH}(t+1) $, $ {\boldsymbol{p}}_{MH}(t) $和 $ {\boldsymbol{p}}_{MH0} $的维数是$N\;\times L +M$. 特别地, 在重启之后, 随机游走者可以返回到表型种子或基因种子. 值得指出的是, 它可以通过预先定义$ {{\boldsymbol{p}}_{MH0}} = {\left[ {(1 - \eta ){{\boldsymbol{\mu }}_0}\;\;\eta {{\boldsymbol{\upsilon }}_0}} \right]^{\rm{T}}}$来调整每一子网络的重要性, 其中, $ {\boldsymbol{\mu}}_0 $和$ {\boldsymbol{\upsilon}} _0 $分别表示多层网络和表型网络的初始概率分布. 参数$ \eta $控制着每个子网络的重要性(即多层或表型). 本文考虑两个参数的取值都等于0.5, 即$ \lambda = 0.5 $和$ \eta = 0.5 $, 也事先设置所有的偏置参数等于1, 即$b_M = b_P = b_{MP} = b_{PM} = 1$.
需要指出的是, 除另有说明, 上述BRWR方法的所有版本中的所有参数的默认值将始终用于以下实验.
2. 实验及结果分析
本节首先介绍各种基因网络、异构网络和多层异构网络的构造. 接着, 描述如何设计一个留一交叉验证方法来测试各种方法的性能. 然后, 将BRWR方法与同类的方法进行比较, 以验证其在识别关键候选疾病基因方面上的性能. 最后, 测试各种参数对BRWR-MH方法性能的影响.
2.1 数据集
根据文献[16, 23]和一些基因数据库网站提供的数据集, 本文构造了实验所需的全部数据集.
1)基因网络. 借助于人类蛋白质或基因(基因即为蛋白质)之间的不同相互作用数据库, 本文构建了3个网络: a)由人类蛋白质图谱数据库
1 获得的RNA-Seq数据构建的共表达(COEX)网络; b)由CCSB相互作用组数据库[33]和 PSICQUIC数据库[34]创建的蛋白质−蛋白质相互作用(PPI)网络; c)由Biocarta2 、Kegg[35]、Reactome[36]、Panther[37]和PID[38]数据库的数据构造的路径(PATH)网络. 表1显示了这些网络的节点、边和平均度的统计信息.表 1 表型、基因和聚合网络的统计属性Table 1 Statistical properties of phenotype, gene and aggregated networks网络 节点数 边数 平均度 COEX 10 415 998 712 47.44 PPI 12 893 70 141 7.73 PATH 10 966 274 051 13.47 聚合网络 17 611 1 342 703 25.79 表型网络 7 324 29 853 4.38 这些网络具有以下某一特性: a)作为一个多网络, 其层分别对应于PPI、COEX和PATH网络; b)作为一个聚合(Aggregated)网络, 其节点和边对应于3个单层网络的并, 即17 611个节点和1 342 703条边(见表1); c)独立地作为单层网络. 需要指出的是, 多网络的每一层都共享相同的节点集, 对应于3个单层网络节点的并, 未在层中出现的基因将被手动地添加到层中. 此外, 假设层中节点与其他层中同样节点的连接是存在的. 图1(a)中展示了一个生物多层网络.
2)异构网络. 为构造一个异构生物网络, 需要3种类型的网络: 基因网络、表型网络和基因−表型二分关系网络. 其中, 基因网络是上面描述的3个单网络中的任意一个网络. 对于表型网络数据, 首先从人类表型本体(Human phenotype ontology, HPO)
3 库下载最新的HPO本体图结构[39]和标注文件phenotype annotation.tab, 其中包括疾病与它们相关联的表型数据. 这里只保留在线人类孟德尔遗传(Online mendelian inheritance in man, OMIM)4 中的疾病记录, 并且提取HPO项中每个表型的最小集合. 表型的最小集合就是描述没有冗余的疾病[40]. 接下来, 依据共享表型的数量, 计算所有疾病对之间的表型相似性. 最后, 正如Li等[11]的描述, 依据相似性得分, 将每种疾病与其最近的5个邻居连接起来, 构建表型−表型相似性网络(即表型网络). 生成的表型网络包括7 324种表型和29 853条边(见表1). 基因−表型网络表征表型实体与相关疾病直接的二分关联关系. 通过使用BiomaRt[40], 可以从OMIM的基因−表型关联中抽取出它们的二分关联关系. 这里下载BiomaRt的最新数据. 由此产生的基因−表型网络由基因和表型间5 483条边组成.因此, 通过使用基因−表型间的二分关联关系, 本文构建了由基因网络和表型网络组成的异构网络, 还构建了由基因多层网络、表型网络和基因−表型网络组成的多层异构网络, 其中每层的基因根据已知的二分关联关系与其相关表型连接, 可以构建3个完全不同的异构网络, 即基因−表型异构网络、共表达−表型异构网络和路径−表型异构网络. 图1(c)展示了由层为异构网络构成的多层异构网络.
2.2 留一交叉验证
本文设计留一交叉验证(Leave-one-out cross validation, LOOCV)方法[41]来测试各种方法的性能. 简言之, 从文献[42]中获得疾病/表型和相关基因, 依赖于被测试的不同方法, 可提取这些基因−表型数据集的不同子集, 并且仅考虑那些至少关联到两个基因的表型. 接着, 每个留出(left-out)基因依次地被删除, 剩余的基因成为方法中的种子基因. 然后, 所有的网络节点都根据它们与种子的接近程度进行评分和排序, 当前运行中疾病关联的留出基因的排名被记录下来. 此外, 当运行LOOCV方法时, 随着阈值的变化, 计算出已知疾病的真阳性率(True positive rate, TPR)和假阳性率(False positive rate, FPR)的值. TPR表示实际带病基因且按该LOOCV实验标准高于某阈值被正确地判为与疾病相关的疾病基因的百分比, FPR表示实际无病但根据实验高于阈值的疾病基因的百分比. 最后绘制出TPR对FPR的受试者工作特征曲线(Receiver operating characteristic, ROC)[43], 并利用该曲线下面积(Area under curve, AUC)作为衡量性能的指标.
2.3 非异构网络上的性能比较分析
本节将提出的BRWR方法与经典的RWR[10]、$ N $维次序统计(N-dimensional order statistics, NDOS)[41]、贴现评定系统(Discounted rating system, DRS)[9]、序列化疾病基因(Prioritization of disease genes, ProDige)[44]和基于多元回归模型的BRIDGE[6]方法进行对比. 这些方法将会在没有集成基因−表型网络的非异构网络上进行. 这些非异构网络包括3个基因网络(PPI、COEX和PATH), 对应的聚合网络(Aggregated)和构成的多层网络(Multilayer). 相应地, 这些方法中使用的种子仅为与已知疾病相关的基因节点. 需要指出的是, 每个疾病至少涉及到两个相关基因的测试集. 测试集相关的基因被包含在具有4 529个共享节点的3个基因网络中. 为最大化这个测试集, LOOCV方法用于测试从DisGeNET[42]中提取的基因−表型关联关系. 共享于3个基因网络中DisGeNET数据包括1 525种基因−表型关联关系, 其中每种关系对应319种不同的表型. 对于每种疾病, 将其关联的基因逐一删除, 而剩下的基因则用作方法的种子. 然后, 通过检索剩余基因以比较不同方法的性能, 所有方法的AUC值如表2所示. 图2(a)和图2(b)分别绘制BRWR方法应用于PPI、COEX、PATH、Aggregated和Multilayer上获得的ROC曲线以及所有方法仅仅应用于多层网络上的ROC曲线.
表 2 不同的非异构网络上的不同方法的AUC值(%)Table 2 AUC values of different algorithms on different non-heterogeneous networks (%)PPI COEX PATH Aggregated Multilayer RWR 73.35 72.84 74.43 76.53 77.98 ProDige 79.12 73.63 80.29 83.27 84.12 NDOS 78.27 74.78 79.86 84.49 87.95 DRS 78.93 74.94 80.87 84.78 88.45 BRIDGE 79.91 74.26 81.51 85.13 89.33 BRWR 81.15 75.20 84.18 86.73 90.17 
图 2 非异构基因网络上不同方法的ROC曲线及其对应的AUC值Fig. 2 ROC curves and AUC values of different algorithms on the non-heterogeneous gene networks在表2中, 对于共表达网络, 所有的方法都表现出最差的性能. 而运行在3个单层网络上的BRWR方法也显示了一般的性能, 如图2(a)所示. 这些结果表明, 仅从mRNA表达数据所生成的网络中检索与已知疾病关联的候选疾病基因似乎是困难的. 然而, 在这3个单层网络中, 运行于从管理信号路径数据库中所提取的PATH网络上的方法显示出最佳的性能, 这可能是路径数据库通常是根据管理生物学知识创建所致. 在所有非异构网络上的ProDige、NDOS、DRS和BRIDGE方法相对于传统的RWR方法在排序候选基因上具有更好的性能(见表2), 这些结果完全与文献[12, 44]中讨论的结论一致. 此外, 所有的方法在聚合网络上比在3个单层网络上能取得更好的性能, 这可能是因为不同数据源的集成比单个数据源更能容忍数据的偏差和不完整. 然而, 在多层网络上, 所有方法产生的疾病关联的留出基因的排名优于其对应聚合网络上的排名(见表2). 尽管多层网络和聚合网络都具有一样的相互作用生物数据源, 但是多层网络却能完全地追踪每层网络的拓扑结构, 并且也考虑了层间相互作用的影响, 这些因素都导致运行在多层网络的方法能取得较好的性能.
通过调整偏置参数, BRWR方法在探索多层网络中高度连接或连接较差的节点, 相比于所有其他方法性能更佳(表2和图2(b)所示). 结果表明, 当要检索候选疾病基因时, 通过适当调整偏置参数的BRWR方法是更为有效的方法. 此外, 所有的方法在利用多个相互作用源的多层框架时表现最佳. 尤其是, 即使单独的共表达网络在序列化候选疾病基因表现出最差性能, 然而它作为多层网络的一层明显地能提高所有方法的性能. 总之, 候选疾病相关基因中心性排名的最佳结果是通过采用多层网络框架实现的(表2和图2).
2.4 异构网络上的性能比较分析
如前所述, 可以将提出的BRWR方法应用于各种异构网络, 称异构网络上偏置重启随机游走(BRWR on the heterogeneous network, BRWRH). 这样, 可以进一步将提出的BRWRH方法与推断基因−表型关系的异构网络的RWR方法(RWR on the heterogeneous network, RWRH)[11]、相互作用表型网络预测疾病基因(Correlating interaction phenotype network to predict disease genes, CIPHER)[45]方法、异构网络上基于RWR的标准拉普拉斯方法(Laplacian based RWR on the heterogeneous networks, LapRWRH)[15]、异构网络上基于网络RWR方法(Network-based RWR on the heterogeneous network, NRWRH)[46]、基于最大信息流MAXIF方法[4]和改进的双标签传播方法IDLP[5]进行比较. 通过使用二分基因−表型关联关系构建的异构网络包括: PPI异构网络(PPI heterogeneous network, PPIH)、共表达异构网络(COEX heterogeneous network, COEXH)、路径异构网络(PATH heterogeneous network, PATHH)、三者的聚合异构网络(Aggregated heterogeneous network, AggregatedH)和多层异构网络(Multilayer heterogeneous network, MultilayerH). 对于异构网络, 方法中使用的种子节点或是基因节点, 或是疾病节点, 或二者兼之. 为执行LOOCV方法, 这里构建了一个由已知疾病组成的测试集, 其中每种疾病至少有两个相关的基因, 由4 529种基因被3个基因网络所共享. 从OMIM中提取的测试集包含243种疾病和884种基因−表型关联关系. 类似于上述的LOOCV方法应用, 与已知疾病相关联的基因迭代地被删除, 疾病本身和剩余疾病相关的基因作为种子节点. 相应地, 表3展示所有方法的AUC值. 同样地, 图3(a)和图3(b)分别绘制了BRWRH方法应用于所有异构网络上获得的ROC曲线以及所有方法仅仅应用于多层异构网络上的ROC曲线.
表 3 不同异构网络上的不同方法的AUC值(%)Table 3 AUC values of different algorithms on different heterogeneous networks (%)PPIH COEXH PATHH AggregatedH MultilayerH CIPHER 74.52 73.51 78.30 77.89 78.31 RWRH 80.37 75.34 79.47 83.67 86.53 MAXIF 80.91 76.56 80.15 84.02 88.43 LapRWRH 81.91 77.80 80.90 84.93 88.78 NRWRH 81.36 78.38 82.70 86.56 89.36 IDLP 82.08 79.25 83.37 87.79 90.16 BRWRH 82.36 80.91 85.17 89.65 91.09 
图 3 异构基因网络上不同方法的ROC曲线及其对应的AUC值Fig. 3 ROC curves and AUC values of different algorithms on the heterogeneous gene networks在这种情况下, 异构网络上执行拓扑偏置重启随机游走方法BRWRH被期望具有更好的性能. 通过对比表2和表3中的数据, 应用于异构网络的这些方法的性能与应用于它们对应的非异构网络的情况类似. 但是, 异构网络上每种方法的AUC值稍高于非异构网络上其对应方法的值. 这些结果表明, 异构网络上的方法比那些非异构网络上的方法在序列化候选基因时效果更好或更准确, 这完全与文献[11, 15]中结论一致. 此外, 这些结果突出了表型网络在加强疾病基因序列化的重要性. 重要的是, 多层异构网络上的BRWRH方法通过取合适的偏置参数, 在序列化任务中(表3和图3)展示出显著的性能改善.
综上所述, 游走于多层异构网络上的拓扑偏置重启随机游走BRWRH方法预测候选疾病关联基因相比于在标准的异构网络或多层网络上取得更准确的性能.
2.5 参数对BRWR-MH方法的影响
除偏置参数外, BRWR-MH方法还包括很多参数, 如$ d_x $, $ r $, $ \delta $, $ \tau $, $ \lambda $, $ \eta $. 这里将使用LOOCV策略来评估这些参数对该方法的影响. 如前所述, 同样构建已知疾病的测试集, 每种疾病至少关联两个基因, 共17 611种基因包含在由3个基因网络构成的基因并集中. 这个测试集由389种疾病和1 120个基因−表型二分关联关系组成. 正如早期的研究[38], 本文绘制了留出基因排名的累积分布函数(Cumulative distribution function, CDF). CDF能有效地展示前$ k $位排名的留出基因的百分比, 并可用于评估和比较不同方法的性能. 这主要集中绘制前$ k = 80 $节点排名的累积分布. 与ROC相比, CDF可以对网络中的所有节点进行排名, 即使目前的研究只关注排名前80的节点. 相反地, ROC方法仅仅排名留出节点相关联基因的一组子集, 例如, 前100种最接近的疾病基因[10-11, 13, 15], 相比于ROC, CDF能获得一个更为普适的验证.
1)偏置参数的影响. 为捕获不同的偏置参数对候选疾病基因排名的影响, 本部分论述了偏置参数$ b_M $, $ b_P $, $ b_{MP} $, $ b_{PM} $对BRWR-MH方法性能的影响, 这些参数允许调整随机游走者对到达节点拓扑属性的依赖性. 通过调整偏置参数, 随机游走者可优先从当前节点跳到具有高或低的拓扑属性值的邻居节点. 为简便起见, 本文设$ b_M $, $ b_P $, $ b_{MP} $, $ b_{PM} $取相同的$ b $值.
图4展示了在BRWR-MH方法中随着偏置参数$ b $的变化而获得的候选疾病留出基因排名的CDF, 当参数$ b $变化时, 所有其他参数保持其默认值. 从图4可知, 当偏置参数从$ b = -4 $增加到$ b = 4 $时, LOOCV性能的显著差异说明了参数的变化对其性能有相对高的敏感性. 此外, 当偏置参数取值越来越接近$ b = -4 $或$ b = 4 $时, 节点排名有一个轻微的变化. 这可从两个方面解释所出现的现象: a)当偏置参数$ b\rightarrow+\infty $时, 偏置的随机游走者将会陷入具有最大度的节点(可能为中心节点)和它的直接邻居节点. 结果使得中心节点和连接中心节点的链将成为网络的中心部分. b)当$ b\rightarrow -\infty $时, 随机游走更频繁地访问连接度较差的节点或悬挂节点. 因此, 在BRWR-MH方法中选择合适的偏置参数以获得更好的疾病基因排名是很重要的. 为简单起见, 依据实验结果, 本文将取偏置参数$ b = 4 $来讨论其他参数对BRWR-MH方法性能的影响.
图 4 排名随偏置参数$ b $变化的累积分布Fig. 4 The cumulative distributions of the ranking with change of the biased parameter $ b $上面已经讨论了偏置参数$ b $的变化对候选疾病关联基因排名的影响, 在固定其值的情况下, 受文献[16, 44]启发, 接下来将要深入地讨论BRWR-MH方法中参数$ d_x $, $ r $, $ \delta $, $ \tau $, $ \lambda $, $ \eta $如何影响候选疾病关联的基因排名. 图5显示了LOOCV中由BRWR-MH方法所产生的留出基因的累积分布函数. 当一个参数变化时, 所有其他参数都保持其默认值, 图5(a) ~ 5(f)分别展示了不同参数变化时的验证结果.
图 5 排名随参数变化的累积分布Fig. 5 The cumulative distributions of the ranking with change of the parameters2)重启参数$ r $的影响. 已经证明[10-11, 47]重启参数$ r $的变化对候选疾病关联的基因仅有较小的影响(见图5(b)). 这里固定$ r $的值为0.7.
3)层间连接的影响. 层间连接参数$d_{ii}^{\alpha\beta} \in [0,+\infty)$量化了从层$ \alpha $中的节点$ i $到层$ \beta $中的节点$ i $游走的“成本”, 并且在不同层间的相同节点处起到桥接作用. 为研究层间连接强度如何影响候选疾病基因的排名, 在固定所有其他参数值的情况下改变层间连接的值. 根据文献[48]的研究, 将所有的层间连接$ d_{ii} $分配一样的值$ d_x \in [0,+\infty) $, 即处于任何层间同样节点对间的交换有一样的“成本”. 由图5(a)所示, 对应较小的层间连接$ (d_x = 0.1 )$, 运行在多层异构网络上的偏置重启随机游走BRWR-MH方法产生了最差的性能, 然而对于小的层间连接$( d_x = 1) $, 该方法在排序任务中表现出显著的性能改善. 由图5(a)中所示, $ d_x $针对BRWR-MH方法性能的影响可大致划分为3种情况.
a)$ d_x\ll 1 $. 随机游走者移动到其他层是困难的, 即随机游走者可能一直在同样的层内移动, 节点的排名可能恢复到单层网络的情况.
b)$ d_x = 1 $. 随机游走者更容易移动到其他层. 没有偏置随机游走的“小世界”效应, 多层异构网络上的BRWR-MH方法能高效地执行. 特别地, 由于$ d_x = 1 $, 其他层中对应节点相对于当前层中该节点的一个直接邻居节点. 由此说明这样的事实: 处于不同层对间的转移概率相比于同一层的概率而言更大, 从而能显著提高排名方法性能. 这阐明了一层中的节点对其他层中其对应节点的重要性, 因为它是随机游走者访问任何其他层的唯一途径.
c)$ d_x\gg 1 $. 类似$ d_x = 1 $的情况, 对于层中一个节点, 其他任何层中其对应的节点应该视为一个普通的邻居, 以致获得较差的节点排名. 对于$ d_x\gg 1 $, 随机游走者的移动主要局限于节点与其对应的副本之间的链上, 忽略了整个多层网络中层内节点间的连接, 从而导致不能高效地检索整个网络中与已知疾病关联的基因, 自然地取得较差的候选疾病基因序列化排名.
因此, 为在疾病关联基因序列化排名中获得最佳性能, 在下面的所有实验中, 所有的$ d_{ii}^{\alpha\beta} $分配一样的值$ d_x = 1 $, 即在不同层之间的任何节点上, 交换层具有相同的“1”的成本.
4)参数$ \delta $和$ \tau $的影响. 这里研究了多层异构网络上偏置随机游走相关的参数$ \delta $和$ \tau $对候选基因排名的影响. 参数$ \delta $控制着在一个非重启步骤下随机游走者从当前层的节点跳跃到其他层的相同节点上的概率. 如果$ \delta = 0 $, 那么随机游走者仅仅在同一层内移动; 如果$ \delta = 1 $, 随机游走者在每一个时间步骤下将会一直在层间跳跃. 由图5(c)所示, 该参数的微小变化对BRWR-MH的性能没有显著的影响, 然而较大$ \delta $变化能明显影响BRWR-MH的性能. 例如$\delta = 0.999\;9$时BRWR-MH的性能明显高于$\delta = 0.000\;1$的性能, 主要因为当$\delta = 0.000\,1$时, 随机游走者几乎在同一层内移动, 而层间移动很少, 这样BRWR-MH的性能将可能恢复到单层网络的情况; 当$\delta = 0.999\;9$时, 随机游走者可以在节点与其他层间对应节点间移动, 也能在层内节点间游走, 从而可以探索整个多层网络, 自然地取得较好排名性能. 参数$ \tau $控制多层网络不同层的重启概率. 理论上, 这将允许利用单层网络上BRWR性能的知识. 例如, 在PATH网络中支持重启, 同时在COEX网络中阻碍重启似乎是合理的. 然而, 随着这个参数的改变, 图5(d)并没有显示出LOOCV性能的显著差异. 特别需要指出的是, 即使此参数发生根本变化, 性能也不会急剧下降. 因为虽然某一层中作为种子的粒子没有重启, 但由于跳跃, 它将会继续探索不同的网络层, 并且仍然可以利用多层异构网络信息.
5)参数$ \lambda $和$ \eta $的影响. BRWR-MH方法在多层异构网络上涉及的参数是$ \lambda $和$ \eta $. 依据基因−表型二分关联关系, 参数$ \lambda $量化了在多层网络和表型网络之间跳跃的概率. 参数$ \lambda $值越大, 网络间跳跃的概率就越高. 如果$ \lambda = 0 $, 那么随机游走者将不会利用表型网络和多层网络之间的二分关联关系; 相反, 如果$ \lambda = 1 $, 那么二分基因−疾病关联则主导着游走者, 并且不允许随机游走者深入探索多层或表型网络的拓扑结构. 但是这个参数的微变化只显示了在性能上的细微变化(见图5(e)). 然而, 当$ \lambda $值趋于0或1时, BRWR-MH性能下降. 因为$ \lambda = 0 $时, 粒子仅在表型网络或多层网络上移动; $ \lambda = 1 $时, 粒子可能仅在二分基因−疾病网络上移动, 这样粒子不能游走于整个多层异构网络, 导致其性能下降. 参数$ \eta $量化了随机游走者在多层网络或表型网络中重新启动的概率. 如果$ \eta = 0 $, 则随机游走者将一直在多层网络中重启; 相反, 如果$ \eta = 1 $, 则随机游走者将一直在表型网络中重启. 在这种情况下, 参数的变化也仅仅略微影响BRWR-MH方法的性能(见图5(f)). 事实上, 不管参数$ \eta $如何变化(处于0 ~ 1间), 即使粒子不在二者中重启, 由于跳跃, 它仍然可以探索多层异构网络, 从而导致BRWR-MH性能不发生大的变化.
总之, 当不考虑层间连接和偏置参数时, BRWR-MH将是一个稳健的方法, 因为适度的参数变化不会导致排名的较大变化.
3. 应用
为探索BRWR-MH方法在疾病相关基因预测中的应用, 本节用新生儿早衰症(NPS)[49]作为该方法的种子节点来预测以前从未诊断的疾病基因, 并展示其网络表示. 新生儿早衰症, 其特征是宫内发育迟缓, 导致随后发育不良和身材矮小. 患者还表现出早衰表现, 皮下脂肪减少, 毛发稀少和大头畸形[50]. 目前只有少数已发表的病例被记录, 但没有基因被确定为这种综合征的致病基因.
为在实验中识别那些高或弱连接的候选疾病基因, 本部分考虑偏置参数的取值分别为−5, −1, 0, 1, 5, 并且所有其他参数都保持默认值. 图6 ~ 10展示了在这些偏置参数(即$ b = -5,-1,0,1,5 $)时NPS (粗体大圆形)作为BRWR-MH方法的种子所产生的排名前30位的基因和疾病的网络表示. 特别地, 本部分详述了偏置参数$ b = 5 $时与NPS高度关联的疾病及基因情况, 图6显示了NPS作为BRWR-MH方法的疾病种子(粗体大圆形)时所获得排名前30的基因和疾病的网络表示[51]. 在图6 ~ 10中, 圆形表示疾病, 正方形表示基因. 为了检索连接度高的节点, 所有偏置的参数取等于5的值, 其他参数保持其默认值. 由图6所示, 许多排名靠前的预测的候选基因(如LMNA、LETM1和POLD1), 都与表型相似的NPS疾病密切关联. LMNA突变导致了综合症早年衰老症候群和其他一些过早衰老的症状, 如A型脂肪代谢障碍引起的下颌骨发育不良. 然而, 少数NPS患者的靶向LMNA测序仍不能识别这种突变[52-53]. BRWR-MH还预测ZMPSTE24基因, 它会导致严重的前期综合症如限制性皮肤病[54]. 但是, 在5个NPS患者中却没有发现这种基因的突变[55].
图 7 所有偏置参数为 −5时的网络表示Fig. 7 Network representation when all the biased parameters are −5
图 8 所有偏置参数为 −1时的网络表示Fig. 8 Network representation when all the biased parameters are −1
图 10 所有偏置参数为1时的网络表示Fig. 10 Network representation when all the biased parameters are 1另一组有趣的候选是由4个基因组成的子网络, 包括IGF2, INS, INSR和RPS6KA3 (见图6). 所有这些基因都参与胰岛素路径, 并且与NPS具有相同表型的疾病相关(即多诺休综合症、高胰岛素原血症和严重生长受限). 胰岛素路径被怀疑在NPS中起了关键性的作用[56]. 类似地, 一组与细胞周期和DNA修复有关的蛋白质通过Wolf-Hirschhorn综合症连接到NPS, 而DNA修复缺陷也被怀疑参与了NPS[55].
在图6 ~ 10中, 与偏置参数为5对比发现, 本文提出的方法不仅可以识别与NPS相关的高连接基因, 还可以识别与NPS相关的一些低连接基因(如图8中ORC4、ORC6和ERCC6). 此外发现一些疾病的症状与NPS直接相关. 重要的是, 当偏置参数$ b = -5 $时, 也观察到一些基因(如图7中AREDYLD ERCC6和ERCC8)更与NPS密切相关. 一些疾病(如图7中CHROMOSOME Xq26.3 DUPLICATION SYNDROME, COCKAYNE SYNDROME A, DYSKERATOSIS CONGENITA, AUTOSOMAL RECESSIVE 6, HYPOSPADIAS-MENTAL RETARDATION SYNDROME, MEIER-GORLIN SYNDROME 2, MEIER-GORLIN SYNDROME 4, 和MEIER-GORLIN SYNDROME 5)也与NPS相关联. 特别地, 当偏置参数$ b = 0 $时, 只有3个基因或疾病直接关联NPS; 当偏置参数$ b = -5 $时, 多达10个基因与NPS直接相关. 这些发现可能为今后NPS的治疗靶点提供一些新的证据. 这证明了提出方法的灵活性、适用性及优越性.
4. 结束语
最近, 针对候选疾病基因中心性排名的各种方法已经被研究, 其中包括基于序列、功能、信息流和RWR的方法. 特别地, RWR方法被认为是排序候选疾病基因的最成功的方法之一. 然而, RWR方法有两个限制: 1)通常是基于单一数据源; 2)通常偏向于高度连接的基因或蛋白质. 为了解决这些问题, 本文提出应用于多个相互作用数据源的BRWR方法, 以便能够提高重要的候选疾病基因识别的性能. BRWR方法是拓扑偏置随机游走方法的扩展版本, 可以应用于多层或多层异构网络. 借助于LOOCV策略, 通过调整BRWR-M方法中的偏置参数, 基因序列化排名结果比非异构网络中的现有方法要好. 本文还构建了多层异构网络, 并提出了BRWR-MH方法. BRWR-MH方法相对于最新的方法能明显地提高序列化基因排名的性能. 此外, 除层间连接和偏置参数之外, 针对不同参数的变化, 本文已经证明BRWR-MH方法是一种较为稳定的方法. 这与BRWR方法参数改变的结论是一致的[11, 15]. 然而, 值得指出的是, 尽管当参数改变时, LOOCV的CDF全局曲线不会发生显著变化, 但在实际应用中, 对前30种疾病和基因的集中分析及网络邻近表示已经揭示出变化规律. 最后, 本文将BRWR-MH方法应用于多层异构网络中来预测与NPS综合症相关的候选基因, 从而展示了此方法在研究疾病病因、揭示靶向治疗和帮助诊断患者方面的有用性.
本文重点关注了由PPI、PATH和COEX网络组成的多层网络. 考虑到疾病具有一定的组织特异性[57]和常见病间的多病关系[58-59], 本文将会从更多生物网络数据源收集数据并将其整合到多层异构网络框架, 旨在进一步探索与疾病相关的未知候选疾病基因. 例如, 为了识别基于基因表达谱与癌症相关的miRNAs[57], 本文将从NCBI gene expression Omnibus
5 收集不同类型癌症的基因表达数据集构建基因表达网络, 以探索癌症与miRNAs间的关联关系. 此外, 如何集成一些包括转录因子、靶向基因、非编码RNAs及药物和治疗靶标的网络也是未来的一个研究方向.
-
图 1 多层网络、异构网络、多层异构网络以及探索它们的随机游走路径(箭头的实线)的示意图
Fig. 1 Schematic of multilayer, heterogeneous and multilayer heterogeneous networks, together with paths of random walks (arrow solid lines)
图 2 非异构基因网络上不同方法的ROC曲线及其对应的AUC值
Fig. 2 ROC curves and AUC values of different algorithms on the non-heterogeneous gene networks
图 3 异构基因网络上不同方法的ROC曲线及其对应的AUC值
Fig. 3 ROC curves and AUC values of different algorithms on the heterogeneous gene networks
图 4 排名随偏置参数$ b $变化的累积分布
Fig. 4 The cumulative distributions of the ranking with change of the biased parameter $ b $
图 5 排名随参数变化的累积分布
Fig. 5 The cumulative distributions of the ranking with change of the parameters
图 7 所有偏置参数为 −5时的网络表示
Fig. 7 Network representation when all the biased parameters are −5
图 8 所有偏置参数为 −1时的网络表示
Fig. 8 Network representation when all the biased parameters are −1
图 10 所有偏置参数为1时的网络表示
Fig. 10 Network representation when all the biased parameters are 1
表 1 表型、基因和聚合网络的统计属性
Table 1 Statistical properties of phenotype, gene and aggregated networks
网络 节点数 边数 平均度 COEX 10 415 998 712 47.44 PPI 12 893 70 141 7.73 PATH 10 966 274 051 13.47 聚合网络 17 611 1 342 703 25.79 表型网络 7 324 29 853 4.38 
下载: 导出CSV
表 2 不同的非异构网络上的不同方法的AUC值(%)
Table 2 AUC values of different algorithms on different non-heterogeneous networks (%)
PPI COEX PATH Aggregated Multilayer RWR 73.35 72.84 74.43 76.53 77.98 ProDige 79.12 73.63 80.29 83.27 84.12 NDOS 78.27 74.78 79.86 84.49 87.95 DRS 78.93 74.94 80.87 84.78 88.45 BRIDGE 79.91 74.26 81.51 85.13 89.33 BRWR 81.15 75.20 84.18 86.73 90.17
下载: 导出CSV
表 3 不同异构网络上的不同方法的AUC值(%)
Table 3 AUC values of different algorithms on different heterogeneous networks (%)
PPIH COEXH PATHH AggregatedH MultilayerH CIPHER 74.52 73.51 78.30 77.89 78.31 RWRH 80.37 75.34 79.47 83.67 86.53 MAXIF 80.91 76.56 80.15 84.02 88.43 LapRWRH 81.91 77.80 80.90 84.93 88.78 NRWRH 81.36 78.38 82.70 86.56 89.36 IDLP 82.08 79.25 83.37 87.79 90.16 BRWRH 82.36 80.91 85.17 89.65 91.09
下载: 导出CSV
-
[1] Guala D, Sonnhammer E L L. A large-scale benchmark of gene prioritization methods. Scientific Reports, 2017, 7: Article No. 46598 doi: 10.1038/srep46598 [2] Schwikowski B, Uetz P, Fields S. A network of protein-protein interactions in yeast. Nature Biotechnology, 2000, 18(12): 1257−1261 doi: 10.1038/82360 [3] Sharma V, Ranjan T, Kumar P, Pal A K, Jha V K, Sahni S, et al. Protein-protein interaction detection: Methods and analysis. Plant Biotechnology. New York: Apple Academic Press, 2018. 391−411 [4] Chen Y, Jiang T, Jiang R. Uncover disease genes by maximizing information flow in the phenome-interactome network. Bioinformatics, 2011, 27(13): i167−i176 doi: 10.1093/bioinformatics/btr213 [5] Zhang Y G, Wang Y, Liu J H, Liu X H, Hong Y X, Fan X, et al. IDLP: A novel label propagation framework for disease gene prioritization. In: Proceedings of the 22nd Pacific-Asia Conference on Knowledge Discovery and Data Mining (PAKDD). Melbourne, Australia: Springer, 2018. 261−272 [6] Chen Y, Wu X B, Jiang R. Integrating human omics data to prioritize candidate genes. BMC Medical Genomics, 2013, 6: Article No. 57 doi: 10.1186/1755-8794-6-57 [7] Lee J H, Zhao X M, Yoon I, Lee J Y, Kwon N H, Wang Y Y, et al. Integrative analysis of mutational and transcriptional profiles reveals driver mutations of metastatic breast cancers. Cell Discovery, 2016, 2: Article No. 16025 [8] Yang K, Zhao X Z, Waxman D, Zhao X M. Predicting drug-disease associations with heterogeneous network embedding. Chaos, 2019, 29(12): Article No. 123109 doi: 10.1063/1.5121900 [9] Yang A Y, Chen J Q, Zhao X M. nMAGMA: A network-enhanced method for inferring risk genes from GWAS summary statistics and its application to schizophrenia. Briefings in Bioinformatics, 2021, 22(4): Article No. bbaa298 doi: 10.1093/bib/bbaa298 [10] Köhler S, Bauer S, Horn D, Robinson P N. Walking the interactome for prioritization of candidate disease genes. The American Journal of Human Genetics, 2008, 82(4): 949−958 doi: 10.1016/j.ajhg.2008.02.013 [11] Li Y J, Patra J C. Genome-wide inferring gene-phenotype relationship by walking on the heterogeneous network. Bioinformatics, 2010, 26(9): 1219−1224 doi: 10.1093/bioinformatics/btq108 [12] Li Y J, Patra J C. Integration of multiple data sources to prioritize candidate genes using discounted rating system. BMC Bioinformatics, 2010, 11(S1): Article No. S20 doi: 10.1186/1471-2105-11-S1-S20 [13] Li Y J, Li J Y. Disease gene identification by random walk on multigraphs merging heterogeneous genomic and phenotype data. BMC Genomics, 2012, 13(S7): Article No. S27 [14] Xie M Q, Xu Y J, Zhang Y G, Hwang T, Kuang R. Network-based phenome-genome association prediction by bi-random walk. PLoS One, 2015, 10(5): Article No. e0125138 doi: 10.1371/journal.pone.0125138 [15] Zhao Z Q, Han G S, Yu Z G, Li J Y. Laplacian normalization and random walk on heterogeneous networks for disease-gene prioritization. Computational Biology and Chemistry, 2015, 57: 21−28 doi: 10.1016/j.compbiolchem.2015.02.008 [16] Valdeolivas A, Tichit L, Navarro C, Perrin S, Odelin G, Levy N, et al. Random walk with restart on multiplex and heterogeneous biological networks. Bioinformatics, 2019, 35(3): 497−505 doi: 10.1093/bioinformatics/bty637 [17] Doncheva N T, Kacprowski T, Albrecht M. Recent approaches to the prioritization of candidate disease genes. WIREs Systems Biology and Medicine, 2012, 4(5): 429−442 doi: 10.1002/wsbm.1177 [18] Yang K, Lu K Z, Wu Y, Yu J, Liu B Y, Zhao Y, et al. A network-based machine-learning framework to identify both functional modules and disease genes. Human Genetics, 2021, 140(6): 897−913 doi: 10.1007/s00439-020-02253-0 [19] Bonaventura M, Nicosia V, Latora V. Characteristic times of biased random walks on complex networks. Physical Review E, 2014, 89(1): Article No. 012803 [20] Ding C F, Li K. Centrality ranking in multiplex networks using topologically biased random walks. Neurocomputing, 2018, 312: 263−275 doi: 10.1016/j.neucom.2018.05.109 [21] Pio-Lopez L, Valdeolivas A, Tichit L, Remy É, Baudot A. MultiVERSE: A multiplex and multiplex-heterogeneous network embedding approach. Scientific Reports, 2021, 11(1): Article No. 8794 doi: 10.1038/s41598-021-87987-1 [22] Peng J, Zhou Y Y, Wang K. Multiplex gene and phenotype network to characterize shared genetic pathways of epilepsy and autism. Scientific Reports, 2021, 11(1): Article No. 952 doi: 10.1038/s41598-020-78654-y [23] Novoa-del-Toro E M, Mezura-Montes E, Vignes M, Térézol M, Magdinier F, Tichit L, et al. A multi-objective genetic algorithm to find active modules in multiplex biological networks. PLoS Computational Biology, 2021, 17(8): Article No. e1009263 doi: 10.1371/journal.pcbi.1009263 [24] Zhao B H, Hu S, Liu X E, Xiong H J, Han X, Zhang Z H, et al. A novel computational approach for identifying essential proteins from multiplex biological networks. Frontiers in Genetics, 2020, 11: Article No. 343 doi: 10.3389/fgene.2020.00343 [25] Dursun C, Smith J R, Hayman G T, Kwitek A E, Bozdag S. NECo: A node embedding algorithm for multiplex heterogeneous networks. In: Proceedings of the IEEE International Conference on Bioinformatics and Biomedicine (BIBM). Seoul, South Korea: IEEE, 2020. 146−149 [26] Bentley B, Branicky R, Barnes C L, Chew Y L, Yemini E, Bullmore E T, et al. The multilayer connectome of Caenorhabditis elegans. PLoS Computational Biology, 2016, 12(12): Article No. e1005283 doi: 10.1371/journal.pcbi.1005283 [27] Shi C, Li Y T, Zhang J W, Sun Y Z, Yu P S. A survey of heterogeneous information network analysis. IEEE Transactions on Knowledge and Data Engineering, 2017, 29(1): 17−37 doi: 10.1109/TKDE.2016.2598561 [28] Gómez-Gardeñes J, Latora V. Entropy rate of diffusion processes on complex networks. Physical Review E, 2008, 78(6): Article No. 065102 [29] Brin S, Page L. Reprint of: The anatomy of a large-scale hypertextual web search engine. Computer Networks, 2012, 56(18): 3825−3833 doi: 10.1016/j.comnet.2012.10.007 [30] Pan J Y, Yang H J, Faloutsos C, Duygulu P. Automatic multimedia cross-modal correlation discovery. In: Proceedings of the 10th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. Seattle, USA: ACM, 2004. 653−658 [31] Smedley D, Köhler S, Czeschik J C, Amberger J, Bocchini C, Hamosh A, et al. Walking the interactome for candidate prioritization in exome sequencing studies of Mendelian diseases. Bioinformatics, 2014, 30(22): 3215−3222 doi: 10.1093/bioinformatics/btu508 [32] Kivelä M, Arenas A, Barthelemy M, Gleeson J P, Moreno Y, Porter M A. Multilayer networks. Journal of Complex Networks, 2014, 2(3): 203−271 doi: 10.1093/comnet/cnu016 [33] Rolland T, Taşan M, Charloteaux B, Pevzner S, Zhong Q, Sahni N, et al. A proteome-scale map of the human interactome network. Cell, 2014, 159(5): 1212−1226 doi: 10.1016/j.cell.2014.10.050 [34] Del-Toro N, Dumousseau M, Orchard S, Jimenez R C, Galeota E, Launay G, et al. A new reference implementation of the PSICQUIC web service. Nucleic Acids Research, 2013, 41(W1): W601−W606 doi: 10.1093/nar/gkt392 [35] Kanehisa M, Sato Y, Kawashima M. KEGG mapping tools for uncovering hidden features in biological data. Protein Science, 2022, 31(1): 47−53 doi: 10.1002/pro.4172 [36] Gillespie M, Jassal B, Stephan R, Milacic M, Rothfels K, Senff-Ribeiro A, et al. The reactome pathway knowledgebase 2022. Nucleic Acids Research, 2022, 50(D1): D687−D692 doi: 10.1093/nar/gkab1028 [37] Mi H Y, Muruganujan A, Casagrande J T, Thomas P D. Large-scale gene function analysis with the PANTHER classification system. Nature Protocols, 2013, 8(8): 1551−1566 doi: 10.1038/nprot.2013.092 [38] Schaefer C F, Anthony K, Krupa S, Buchoff J, Day M, Hannay T, et al. PID: The pathway interaction database. Nucleic Acids Research, 2009, 37(S1): D674−D679 [39] Köhler S, Gargano M, Matentzoglu N, Carmody L C, Lewis-Smith D, Vasilevsky N A, et al. The human phenotype ontology in 2021. Nucleic Acids Research, 2021, 49(D1): D1207−D1217 doi: 10.1093/nar/gkaa1043 [40] Greene D, BioResource N I H R, Richardson S, Turro E. Phenotype similarity regression for identifying the genetic determinants of rare diseases. The American Journal of Human Genetics, 2016, 98(3): 490−499 doi: 10.1016/j.ajhg.2016.01.008 [41] Aerts S, Lambrechts D, Maity S, Van Loo P, Coessens B, De Smet F, et al. Gene prioritization through genomic data fusion. Nature Biotechnology, 2006, 24(5): 537−544 doi: 10.1038/nbt1203 [42] Piñero J, Bravo À, Queralt-Rosinach N, Gutiérrez-Sacristán A, Deu-Pons J, Centeno E, et al. DisGeNET: A comprehensive platform integrating information on human disease-associated genes and variants. Nucleic Acids Research, 2017, 45(D1): D833−D839 doi: 10.1093/nar/gkw943 [43] Hanley J A, McNeil B J. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology, 1982, 143(1): 29−36 doi: 10.1148/radiology.143.1.7063747 [44] Mordelet F, Vert J P. ProDiGe: Prioritization of disease genes with multitask machine learning from positive and unlabeled examples. BMC Bioinformatics, 2011, 12: Article No. 389 doi: 10.1186/1471-2105-12-389 [45] Wu X B, Jiang R, Zhang M Q, Li S. Network-based global inference of human disease genes. Molecular Systems Biology, 2008, 4: Article No. 189 doi: 10.1038/msb.2008.27 [46] Chen X, Liu M X, Yan G Y. Drug-target interaction prediction by random walk on the heterogeneous network. Molecular BioSystems, 2012, 8(7): 1970−1978 doi: 10.1039/c2mb00002d [47] Blatti C, Sinha S. Characterizing gene sets using discriminative random walks with restart on heterogeneous biological networks. Bioinformatics, 2016, 32(14): 2167−2175 doi: 10.1093/bioinformatics/btw151 [48] De Domenico M, Solé-Ribalta A, Gómez S, Arenas A. Navigability of interconnected networks under random failures. Proceedings of the National Academy of Sciences of the United States of America, 2014, 111(23): 8351−8356 [49] Pivnick E K, Angle B, Kaufman R A, Hall B D, Pitukcheewanont P, Hersh J H, et al. Neonatal progeroid (Wiedemann-Rautenstrauch) syndrome: Report of five new cases and review. American Journal of Medical Genetics, 2000, 90(2): 131−140 doi: 10.1002/(SICI)1096-8628(20000117)90:2<131::AID-AJMG9>3.0.CO;2-E [50] Kiraz A, Ozen S, Tubas F, Usta Y, Aldemir O, Alanay Y. Wiedemann-Rautenstrauch syndrome: Report of a variant case. American Journal of Medical Genetics Part A, 2012, 158A(6): 1434−1436 doi: 10.1002/ajmg.a.35336 [51] Kohl M, Wiese S, Warscheid B. Cytoscape: Software for visualization and analysis of biological networks. Data Mining in Proteomics: From Standards to Applications. Totowa: Humana, 2011. 291−303 [52] Becerra C H, Contreras-García G A, Perez Vera L A, Díaz-Martínez L A, Beltran Avendaño M A, Salazar Martínez H A. Wiedemann-Rautenstrauch syndrome prenatal diagnosis. Journal of Perinatology, 2014, 34(12): 954−956 doi: 10.1038/jp.2014.156 [53] Paolacci S, Bertola D, Franco J, Mohammed S, Tartaglia M, Wollnik B, et al. Wiedemann-Rautenstrauch syndrome: A phenotype analysis. American Journal of Medical Genetics Part A, 2017, 173(7): 1763−1772 doi: 10.1002/ajmg.a.38246 [54] Navarro C L, Esteves-Vieira V, Courrier S, Boyer A, Duong Nguyen T, Huong L T T, et al. New ZMPSTE24 (FACE1) mutations in patients affected with restrictive dermopathy or related progeroid syndromes and mutation update. European Journal of Human Genetics, 2014, 22(8): 1002−1011 doi: 10.1038/ejhg.2013.258 [55] Beauregard-Lacroix E, Salian S, Kim H, Ehresmann S, D'Amours G, Gauthier J, et al. A variant of neonatal progeroid syndrome, or Wiedemann-Rautenstrauch syndrome, is associated with a nonsense variant in POLR3GL. European Journal of Human Genetics, 2020, 28(4): 461−468 doi: 10.1038/s41431-019-0539-6 [56] Arboleda G, Ramírez N, Arboleda H. The neonatal progeroid syndrome (Wiedemann-Rautenstrauch): A model for the study of human aging? Experimental Gerontology, 2007, 42(10): 939−943 doi: 10.1016/j.exger.2007.07.004 [57] Zhao X M, Liu K Q, Zhu G H, He F, Duval B, Richer J M, et al. Identifying cancer-related microRNAs based on gene expression data. Bioinformatics, 2015, 31(8): 1226−1234 doi: 10.1093/bioinformatics/btu811 [58] He F, Zhu G H, Wang Y Y, Zhao X M, Huang D S. PCID: A novel approach for predicting disease comorbidity by integrating multi-scale data. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 2017, 14(3): 678−686 doi: 10.1109/TCBB.2016.2550443 [59] Dong G Y, Feng J F, Sun F Z, Chen J Q, Zhao X M. A global overview of genetically interpretable multimorbidities among common diseases in the UK Biobank. Genome Medicine, 2021, 13(1): Article No. 110 doi: 10.1186/s13073-021-00927-6 期刊类型引用(1)
1. 邱馨锐，赵文，吴长汶，张斐，杨朝阳，李灿东，周常恩. 刍议高阶关系在中医隐性知识显性化研究中的关键作用. 福建中医药. 2024(08): 24-26+30 . 
百度学术其他类型引用(1)
-
下载:
计量
- 文章访问数: 814
- HTML全文浏览量: 422
- PDF下载量: 127
- 被引次数: 2
