山东科技大学，数学与系统科学学院，山东 青岛

收稿日期：2023年4月12日；录用日期：2023年5月15日；发布日期：2023年5月22日

摘要

乳腺癌是女性最常见的恶性肿瘤之一，其发展是一个长期的、多阶段、多基因改变积累的过程。临床及实验研究发现，在乳腺癌发展与转移的过程中，许多原癌基因、抑癌基因及相关蛋白质发生改变。本文提出了一种结合随机游走重启算法的综合模型，并将其应用于乳腺癌骨转移关键候选基因的鉴定中，在蛋白质相互作用网络上识别乳腺癌转移候选基因。首先利用RWR算法和置换检验规则对基因进行预选，得到候选基因集。然后构建一个候选基因间的关键子网并计算关键子网中节点的介数中心性。最后，利用交互得分规则对基因筛选，选择介数得分排名前三十的基因作为乳腺癌骨转移关键基因。相互作用分析、富集分析及文献挖掘的结果表明26个潜在的关键基因都参与了乳腺癌骨转移的起始或进程，验证了该综合模型的有效性。

关键词

肿瘤转移，复杂网络，蛋白质相互作用，随机游走重启算法，介数中心性

Network-Based Method for the Identification of Genes Related to Cancer Metastasis

Zhimeng Tian, Zihao Zhang

College of Mathematics and System Science, Shandong University of Science and Technology, Qingdao Shandong

Received: Apr. 12^th, 2023; accepted: May 15^th, 2023; published: May 22^nd, 2023

ABSTRACT

Breast cancer is one of the most common malignant tumors in women, and its development is a long-term, multi-stage, multi-genetic alteration accumulation process. Clinical and experimental studies have revealed that many proto-oncogenes, oncogenes and related proteins are altered during the process of breast cancer development and metastasis. In this paper, we propose an integrated model combining random wander restart algorithm and apply it to the identification of key candidate genes for breast cancer bone metastasis on protein interaction network to identify breast cancer metastasis candidate genes. The genes are first pre-selected using the RWR algorithm and the substitution test rule to obtain the candidate gene set. Then a key subnet among candidate genes is constructed and the betweenness centrality of nodes in the key subnet is calculated. Finally, the genes were screened using the interaction score rule, and the top thirty genes in the mediator score were selected as the key genes for breast cancer bone metastasis. The results of interaction analysis, enrichment analysis and literature mining indicated that 26 potential key genes were involved in the initiation or progression of breast cancer bone metastasis, validating the validity of this comprehensive model.

Keywords:Tumor Metastasis, Complex Networks, Protein Interactions, RWR, Betweenness Centrality

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

乳腺癌是女性最常见的恶性肿瘤之一，其发展是一个长期的、多阶段、多基因改变积累的过程 [1] [2] ，病因与发病机制一直是乳腺癌研究领域的重点。乳腺肿瘤的发展过程一般可分为以下三个阶段：起始、侵袭和转移。在乳腺癌中，常见的转移部位包括骨头、肝脏和大脑等，并且它们有完全不同的症状。其中，超过70%的乳腺癌患者死亡可归因于乳腺癌骨转移 [3] [4] 。除此之外，骨转移还同时伴随不明原因的体重下降、骨头关节疼痛、黄胆及神经症状等。因此，本研究以乳腺癌骨转移作为研究案例具有特殊的意义和价值。临床及实验研究发现，在乳腺癌发展与转移的过程中，许多原癌基因、抑癌基因及相关蛋白质发生改变 [5] 。

2011年，Ruan等人 [6] 整合蛋白质相互作用网络和基因表达数据，提出一种鉴定乳腺癌转移中的生物标志物的方法。该方法应用于两个乳腺癌微阵列数据集的生物标志物的发现，这些生物标记物中含有大量已知的癌症易感基因，并在致癌过程中涉及的生物学过程和通路中显著富集。此外，该方法选择的标记在两个数据集中具有比以往研究更高的稳定性。2016年，Tian等人 [7] 提出了一种结合蛋白质相互作用信息与表达数据的方法来检索乳腺癌转移相关基因。该论文基于随机森林MDG准则和蛋白质相互作用信息的GeneRank算法两种变量选择算法的重要得分，提出了一种新的基于随机森林的算法，来识别与乳腺癌转移高度相关的重要变量。较高的预测精度充分说明了通过该算法选择基因表的可靠性和可解释性。2017年，Cai等人 [8] 提出了一种基于蛋白质相互作用网络使用最短路径算法识别乳腺癌骨转移关键基因的方法。该方法在已确认的乳腺癌和骨癌相关基因数据的基础上，将最短路径算法应用于蛋白质相互作用网络，最终得到18个介导乳腺癌骨转移的候选基因，并进行了广泛的分析，根据研究证实了这些基因参与了乳腺癌骨转移过程。

本文在第一部分介绍了当前国内外的乳腺癌转移研究，在第二部分介绍了数据来源并提出了一种结合随机游走重启算法和介数中心性的综合模型，在第三部分将模型应用于乳腺癌骨转移关键候选基因的识别，在第四部分通过相互作用分析、富集分析及文献挖掘的分析结果从基因层面揭示了肿瘤转移机制。

2. 基于随机游走重启算法的乳腺癌转移关键基因识别

2.1. 数据与数据预处理

乳腺癌相关基因主要来源于Oncomine数据库(https://www.oncomine.org)和TCGA数据库(https://www.cancer.gov/ccg/rresearc/tcga)。Oncomine数据库是一个癌症微阵列数据库，提供了多种癌症类型的转录组数据。通过对两个数据库的检索，并对得到的基因进行整合，将它们映射到STRING数据库中的Ensembl ID上。最终得到了358个乳腺癌相关基因，其集合用 $S_1$ 表示；得到的588个骨癌相关基因，其集合用 $S_2$ 表示。

2.2. 方法介绍

在构建好的蛋白质相互作用网络上，本研究提出了一种基于随机游走重启算法的癌症转移基因识别方法，并将其运用于乳腺癌骨转移特异性关键基因的鉴定中。整个模型的流程见图1所示。

图1. 整个模型的流程图

2.2.1. 随机游走重启算法

RWR算法是一种经典的排序算法，它从一些种子节点开始，模拟其在网络中随机游走和重启，同时更新所有节点的概率得分并对节点进行排名。在本方法中，将已验证的乳腺癌与骨癌的基因作为种子节点，在PPIN上执行RWR算法，利用算法的结果对网络中其它节点进行排名，以实现对基因的预选，其主要过程如下：

其中，重启概率 $r$ 表示在随机游走过程中回到种子节点的概率，取值范围为 $\left[0,1\right]$ 。当 $r=0$ 时，种子节点的信息对于结果无影响，节点的概率得分排名由节点在网络结构中的位置决定。根据先前结合RWR算法应用于基因识别的研究，当 $r=0.8$ 时，算法的结果可以更好地融合种子节点的信息和网络的结构的影响。因此，本研究将重启概率 $r$ 设为0.8来表示种子节点的重要性。当 ${\Vert P_{i+1}-P_i\Vert }_{L_1}<{10}^{-6}$ ，迭代结果趋于稳定，算法停止并输出 $P_{i+1}$ 。在 $P_{i+1}$ 中，每个节点(基因)的概率得分表示其与种子节点在PPIN中的相似性得分。因此，具有较高概率得分的基因与已验证的骨癌和乳腺癌基因更相关，从而更有可能是转移相关基因。为了防止丢失潜在的转移相关基因，将概率得分的阈值设为 ${10}^{-5}$ ，即节点概率得分大于 ${10}^{-5}$ 对应的基因被筛选出来，这些基因被记为RWR基因。与其他算法筛选相关基因相比，RWR具有算法快速、简单及准确性高等特点。

2.2.2. 置换检验规则

上述RWR基因可能会受到PPIN结构的影响，从而存在很多与癌症转移无关的假阳性基因。例如，在网络中度中心性远高于平均度的节点即使与种子节点无关，也会获得较大的概率得分，因此使用随机选择网络中的某些节点作为种子节点这些基因仍会得到较高的概率得分排名。基于这一思想，本研究提出利用置换检验规则来筛选并删除这些假阳性基因。

首先，将总置换数(total permutations)设为1000，即随机构建了1000个Ensembl IDs集合，记为 $E_1,E_2,\cdots ,E_{1000}$ ，每个集合包含882个随机的基因Ensembl IDs。然后，通过将 $E_i\left(1\le i\le 1000\right)$ 中的882个基因设置为种子节点，在PPIN上执行RWR算法以获取每个RWR基因的概率得分。对于每个RWR基因，存在一个真实概率得分 $Ps\left(g\right)$ 和1000个随机概率得分 $P{s}_i\left(g\right)$ 。最后，对每一个RWR基因 $g$ ，计算p-value值如下：

$\text{p-value}\left(g\right)=\frac{\left\{P{s}_i\left(g\right)>Ps\left(g\right)\right\}}{\left\{\text{totalpermutations}\right\}}$ (2-2)

其中， $\left\{P{s}_i\left(g\right)>Ps\left(g\right)\right\}$ 表示1000个随机概率得分 $P{s}_i\left(g\right)$ 中大于真实概率得分 $Ps\left(g\right)$ 的个数，如果随机概率得分 $P{s}_i\left(g\right)$ 普遍大于真实概率得分 $Ps\left(g\right)$ ，说明 $g$ 更可能是因为网络结构而被选出来的假阳性基因。显然，p-value值很高的RWR基因应当被删除。由于0.05是作为被广泛接受的统计学检验传统显著性水平的阈值，本研究选择p-value值小于0.05的RWR基因作为乳腺癌骨转移的潜在候选基因做进一步分析。

2.2.3. 构建关键子网

经过RWR算法和置换检验得到的可能仅与乳腺癌基因或骨癌基因中一项高度相关，为了从候选基因中筛选出在转移过程中起到更关键作用的基因，本研究从PPIN中除候选基因对应节点以外的所有节点及连边删除，得到了一个仅保留候选基因间相互作用的关键子网。

在PPIN中，蛋白质之间连边的权重代表蛋白质间的相互作用得分，其取值范围为 $\left(0,1000\right)$ 。在基于最短路径的模型上，节点之间通常使用较低的权重来表示更强的相关性。因此，对于关键子网中蛋白质 $p_a$ 和 $p_b$ ，当相互作用得分 $S\left(p_a,p_b\right)\ne 0$ 时，利用下面的公式更新它们连边的权重。

$w\left(p_a,p_b\right)=\{\begin{array}{l}\frac{1}{S\left(p_a,p_b\right)};\text{if}S\left(p_a,p_b\right)\ne 0\hfill \\ 0\text{otherwise}\hfill \end{array}$ (2-3)

显然，当权重值 $w\left(p_a,p_b\right)$ 越小，表示关键子网中蛋白质之间的相互作用更强。根据权重更新公式函数图像，见图2可知，当相互作用得分 $S\left(p_a,p_b\right)$ 增大时，函数曲线的斜率减小，更新后的权重相差更小。因此，通过公式(2-3)对关键子网的权重进行归一化处理后，权重的分布相对之前更加集中。

图2. 权重更新公式函数图像

2.2.4. 介数中心性

通过删除网络中的无关基因并更新节点间的权重，得到了一个候选基因间的关键子网。本研究中，与乳腺癌和骨癌相关基因都有密切关联的基因更有可能是与转移有关的特异性基因。在该关键子网中，我们应用dijkstra算法计算乳腺癌关键基因集 $S_1$ 中所有基因到骨癌关键基因集 $S_2$ 中所有基因的最短路径。在每条最短路径上，连续的节点通过较低权重的边相连，表示了相对应蛋白质之间的强相关性。由于每条最短路径的起始和末端节点为已验证的乳腺癌和骨癌相关基因，这些最短路径内部节点所对应的基因可能是导致乳腺癌向骨转移的特异性基因。

对于除去882个种子节点对应基因以外的每个候选基因 $g$ ，计算它基于所有基因对 $\left(G_a,G_b\right)\left(\forall G_a\in S_1,\forall G_b\in S_2\right)$ 的介数中心性得分：

$BC\left(g\right)=\frac{B\left(g\right)}{\left|S_1\right|\cdot \left|S_2\right|}$ (2-4)

其中 $B\left(g\right)$ 代表在所有最短路径中包含基因 $g$ 对应节点的数目； $\left|S_1\right|\cdot \left|S_2\right|$ 指 $S_1$ 和 $S_2$ 中基因数目的乘积，表示可能存在的所有最短路径数目。介数中心性得分 $BC\left(g\right)$ 表示了最短路径经过基因 $g$ 的概率，因此 $BC\left(g\right)$ 值更大的候选基因更有可能是与乳腺癌骨转移有关的关键基因。

此外，基于PPI中相互作用得分较高的两个蛋白质更有可能具有相似功能的研究，利用公式(2-5)计算了每个候选基因 $g$ 的最大–最小交互得分MMIS。

$MMIS\left(g\right)=\min \left\{\max \left\{S\left(g,g^{\prime }\right):g^{\prime }\in S_1\right\}，\max \left\{S\left(g,g^{\prime }\right):g^{\prime }\in S_2\right\}\right\}$ (2-5)

本研究以STRING数据库中最高置信度值900作为MMIS的阈值，最终选择了MMIS不小于900且介数中心性得分BC排名前三十的基因作为乳腺癌骨转移的潜在关键基因，这些基因被认为在乳腺癌骨转移中发挥了重要作用。

3. 结果与分析

通过将已验证的882个乳腺癌与骨癌相关基因作为种子节点，在PPIN上执行RWR算法，选择结果中概率得分大于 ${10}^{-5}$ 的基因，得到了7105个RWR基因。此外，为了消除网络结构的假阳性影响，本研究进行了置换检验规则并选择p-value值小于0.1的2955个基因作为乳腺癌骨转移的候选基因。

为了筛选在乳腺癌骨转移过程中起关键作用的基因，本研究构建了一个候选基因间的关键子网，并在该网络上应用最短路径算法。然后，计算了每个候选基因的最大–最小交互得分MMIS和介数中心性得分BC，最后选择了MMIS不小于900且介数中心性得分BC排名前三十的基因作为乳腺癌骨转移的潜在关键基因，它们的Ensembl ID、基因名、P-value值、MMIS值及验证的参考文献在表1中展示。

表1. 通过该方法得到的30个潜在的关键基因根据BC的排名

4. 分析与讨论

根据文献中的实验证据，通过该方法得到的30个潜在关键基因中，有26个基因已被证实与乳腺癌的骨转移中的基因组整合、转移的发展过程或特定并发症有关。这些发现证明了本研究所提出的综合方法的有效性和合理性。

4.1. 相互作用分析

为了验证这些关键基因与已验证的乳腺癌和骨癌相关基因之间的关系，将这30个潜在的关键基因和集合 $S_1$ 、 $S_2$ 中重复的50个基因输入STRING数据库中，得到了它们之间的相互作用网络如图3所示。本研究将网络中边权重的阈值设为700，仅当节点对应蛋白质之间的相互作用得分大于或等于700时，它们之间的连边在网络中保留。图中左边的26个节点为已被证实与乳腺癌的骨转移有关的潜在关键基因，中间的节点表示已验证的乳腺癌和骨癌相关基因，右边的四个节点表示还未被实验证明的潜在关键基因。该蛋白质网络的节点平均度为13.5，相对于一般的网络节点间联系更紧密，表示蛋白质之间的强相关性。其中，26个已被证实与乳腺癌的骨转移有关的潜在关键基因对应节点的平均度为20.6，未验证的4个基因的平均节点度为22.8，暗示了这些潜在关键基因在乳腺癌骨转移中的重要作用。

图3. 已知的乳腺癌和骨癌相关基因与30个关键基因之间的相互作用图

4.2. 基因富集分析

为了揭示这30个潜在关键基因的分子功能及生物学背景，本研究使用了DAVID功能注释工具。DAVID是一个功能强大的注释工具，可用于对大型基因列表进行整合和系统化，并根据KEGG通路和GOterms对基因进行分析。其中GOterms主要包含三个部分：生物过程(BP)，细胞成分(CC)和分子功能(MF)。由于本研究提出的方法主要是通过蛋白质间相互作用来鉴定基因，因此基于BP进行分析更为合理。

在基于BP的富集分析中，30个潜在关键基因显著富集在15个GOterms上，图4展示了十个主要的GOterms名称及每个GOterms富集的基因数量。文献综述显示，大多数富集的GOterms与乳腺癌骨转移的生物过程有关。例如，GO:0045944-RNA聚合酶II启动子转录的正调控富集了包括AP1、HARS和CREBBP等14个基因，它主要激活了RNA聚合酶II启动子转录过程的频率、速率及程度，并与骨癌的形成和转移相关 [6] [58] 。GO:0038128-ERBB2信号通路的P-value值最低，该信号通路已被研究证明在乳腺癌的骨转移过程中被激活 [7] 。此外，GO:0008284-细胞增殖的正向调控和GO:0007219-Notch信号通路被报道促进了骨肿瘤细胞的形成 [8] [59] 。

在基于KEGG通路的分析中，30个潜在关键基因显著富集在14条通路上，图5展示了十个主要的KEGG通路名称及每条通路富集的基因数量。其中，hsa05200-癌症通路是富集最显著的通路，包含UBC、AP1、VEGF和FN1等15个基因。该通路富集的基因通常参与癌症的形成和发展，从而体现了该方法的有效性。hsa04151-PI3K-Akt信号通路富集了包括CDC37、MDM2和PIK3R1等12个基因，它是一种在调节代谢、细胞生长和血管生成中发挥重要作用的细胞内信号转导通路，已证明与乳腺癌的骨转移高度相关 [60] [61] 。此外，hsa04012-ErbB信号通路、hsa04350-TGF-β信号通路等KEGG通路已被证明有助于乳腺癌的骨转移过程 [62] [63] [64] 。

富集分析的结果表明，这些基因的功能活性和定位与乳腺癌骨转移过程直接或间接相关，进一步支持了提出方法的结果。此外，这些发现可能会拓展探索乳腺癌骨转移新机制的途径。

图4. 潜在关键基因富集的十个主要GO terms

图5. 潜在关键基因富集的十个主要KEGG通路

4.3. 基因与转移机制分析

在30个潜在的关键基因中，有26个基因已通过文献中的实验证据证实与乳腺癌骨转移的形成或发展有关。在所有关键基因中，介数中心性BC值排名最高的基因UBC (泛素C)已被证实参与了乳腺癌的转移过程。体外实验证明，UBC通过调节控制细胞周期和增殖的中枢基因在骨髓间充质干细胞的复制性衰老中起着关键作用，这一过程促进了骨骼组织中乳腺癌细胞的适应性和侵袭性 [9] 。此外，UBC末端水解酶-L3促进了EMT过程，这是癌细胞侵袭和转移的重要因素 [10] 。关键基因中另一个与泛素有关的基因UBE2I参与编码酵母泛素结合酶的人类同源物，对调节细胞周期至关重要。研究表明，UBE2I的异位表达可以促进乳腺癌细胞的侵袭和转移，包括骨转移 [39] 。UBE2C也是一个泛素相关基因，它编码破坏有丝分裂周期蛋白和细胞周期进展所需的蛋白质，据报道在乳腺癌的增殖和迁移中高度表达 [57] 。

Wnt信号通路是调控细胞增殖和分化的重要途径，大量研究表明Wnt信号通路与包括乳腺癌骨转移在内的骨恶性肿瘤密切相关 [65] [66] 。作为Wnt信号通路的一部分，CTNNB1已被证实在骨破坏过程中异常表达，这对乳腺癌骨转移的发生有重要意义 [11] 。CTNNB1编码的蛋白质是组成粘附连接(AJs)的复合蛋白的一部分，在骨肿瘤转移的骨病变阶段是必不可少的 [12] 。另一个参与Wnt信号通路的基因EGF，编码一种促进多种细胞类型生长、增殖和分化的蛋白质，EGF在患者血清中的表达对乳腺癌骨转移瘤的早期诊断具有重要价值 [30] [31] 。

在原癌基因家族中，HRAS是一种常见的致癌基因，已被报道为乳腺腺肌上皮瘤的驱动因素。体外实验证明，HRAS在乳腺癌骨转移过程中过度表达 [15] [16] 。作为乳腺癌的潜在治疗靶点，MYC编码一种核磷蛋白，在细胞周期进程、凋亡和细胞转化中发挥作用。MYC作为侵袭性和转移性乳腺癌的靶基因，其基因易位与骨微环境的修饰过程有关 [19] [20] 。另一个原癌基因MDM2参与了许多癌症的发病机制，包括乳腺癌的侵袭和转移过程。基质金属蛋白酶(MMP9)是一种增强肿瘤侵袭性的重要酶，MDM2促进了MMP9的高表达 [42] [43] 。MMP9是一种增加肿瘤侵袭性的重要酶。根据最近的研究，MDM2促进乳腺癌的侵袭和转移。对于另一个原癌基因JUN，研究表明该基因通过介导骨吸收过程中生长因子的释放增强了转移性乳腺癌细胞在骨骼中的生长和存活 [21] [22] 。

CDC37、CDC5L、CDK4、PLK1和CDC20是与细胞周期调控相关的基因。根据临床数据分析，PLK1和CDC20的表达与乳腺癌患者的无病生存率(DFS)密切相关。因此，PLK1和CDC20在促进乳腺肿瘤形成中发挥作用，并可能作为潜在的肿瘤诊断标志物和个体治疗靶标 [51] [52] 。CDK4是功能性丝氨酸/苏氨酸蛋白激酶家族的成员，参与细胞周期调节。根据其表达谱，CDK4在乳腺癌骨转移过程中被过度表达，作为乳腺癌骨转移的治疗靶点，其下调被认为是判断预后的合适指标 [17] [18] 。此外，CDC37通过稳定突变和过度表达的致癌激酶的受损结构来介导癌变和维持恶性表型，已被证实参与乳腺癌骨转移的功能途径 [13] [14] 。

对于其他一些重要的基因，CREBBP已被广泛报道与骨髓白血病的发展有关。在白血病的发病过程中，骨髓和骨组织的微环境发生了剧烈变化。作为骨组织的功能调节因子，CREBBP已被证实在乳腺癌骨转移过程中参与骨微环境的修饰过程，证实其在肿瘤转移过程中的重要作用 [40] [41] 。CD44上调了基质金属蛋白酶(MMPs)的表达，MMPs可以促进细胞从骨髓窦外渗，有利于乳腺癌细胞通过新生血管进入血液循环b并在乳腺癌骨转移中发挥特定作用 [27] [28] [29] 。BMP7编码TGF-β蛋白超家族的分泌配体，其诱导异位骨形成并促进骨转移。此外，BMP7参与了乳腺癌细胞的EMT过程，被报道直接参与乳腺癌脑转移的进程 [37] [38] 。鉴于NOTCH1在乳腺癌转移中的重要作用，它被认为是乳腺癌的核心治疗靶点。根据最近的文献研究，通过抑制NOTCH1信号通路，可以抑制乳腺癌细胞的骨转移 [48] [49] [50] 。

在尚未被验证的4个基因中，有3个存在于上述的基因家族中。其中，UBB和RPS27A是与泛素相关的基因，已发现它们在骨髓中广泛表达。通过降低UBB降低泛素水平，是一种潜在的癌症治疗方法 [67] 。RPS27A调节细胞周期进程并抑制白血病细胞的凋亡，暗示其在骨微环境的潜在作用 [68] [69] 。CDC5L是一个参与细胞周期调控的基因，最近的研究发现它是骨肉瘤和宫颈肿瘤的候选致病基因 [70] [71] 。另一个基因POLR2A是有关细胞存活的重要基因，已被确定为人类三阴性乳腺癌的精确靶标 [72] ，表明其在乳腺癌进展中的作用。

根据文献挖掘的结果及分析表明，通过该方法鉴定出的大多数潜在关键基因已被证实有助于乳腺癌骨转移的进展，这证明了所提出的综合方法可能是预测肿瘤转移相关基因和研究肿瘤转移机制的有用工具。

文章引用

田智萌,张子豪. 基于网络方法确定乳腺癌转移相关基因
Network-Based Method for the Identification of Genes Related to Cancer Metastasis[J]. 理论数学, 2023, 13(05): 1267-1280. https://doi.org/10.12677/PM.2023.135130

参考文献