杨天濠 王新赠
(1.中国石油大学(华东)计算机科学与技术学院 青岛 266580)
(2.山东科技大学数学与系统科学学院 青岛 266580)
肿瘤转移是指恶性肿瘤细胞从原发部位,经淋巴道、血管或体腔等途径,到达其他部位继续生长的过程[1]。恶性肿瘤的转移通常发生在癌症晚期,是导致癌症患者死亡的主要原因之一。其中,肺癌的转移是一个较为复杂、由多基因参与的过程,它严重影响肺癌患者治疗的疗效和预后。骨是肺癌远处转移常见的靶部位之一,临床发现,约40%的晚期肺癌患者会发生骨转移,同时会引发多种并发症[2~3]。生物学研究证明,特定的遗传背景对癌症转移有重要的影响,有些基因可能参与骨转移发展的进程。因此,对肺癌骨转移相关候选基因的鉴定和筛选对于肺癌患者的诊断和治疗具有迫切而重要的意义。
目前,肿瘤转移相关基因的鉴定和验证主要依赖于临床医学及生物学实验,需要花费大量的时间和成本,限制了发现的能力。随着生物信息学的发展,一些计算方法被应用于识别疾病相关基因及其驱动因子[4~6]。相对于临床医学和生物学实验,计算方法具有高效、低成本的优点。
针对恶性肿瘤转移关键基因的发现问题,我们提出了一种基于蛋白质-蛋白质相互作用网络(PPIN)的癌症转移基因识别方法,将其应用于肺癌骨转移关键候选基因的鉴定。首先利用随机游动重启(RWR)算法对基因进行分析和预选,然后通过置换检验规则消除网络结构的影响,并利用交互得分规则和富集分析对基因进一步筛选,最终获得了12 个可能与肺癌骨转移有关的关键基因。根据文献挖掘的结果,这些基因中有9 个基因已被证实与肺癌骨转移的形成或发展有关,并揭示了这些基因可能参与的潜在分子过程,为利用计算方法研究肿瘤转移机制提供了新的思路。
肺癌和骨癌相关基因主要来源于Oncomine 数据库和TCGA 数据库。Oncomine 数据库是一个整合了264个独立的数据集,涉及35种癌症类型的综合型癌症数据库。TCGA 数据库是目前为止可以获得的公开数据库里面数据相对全面的一个,在各个领域得到了广泛的应用。通过对两个数据库的检索,我们最终得到了412 个肺癌相关基因,其集合用S1表示;以及348 个骨癌相关基因,其集合用S2表示。
通过对STRING 数据库(版本11.0)的检索,我们得到了5,879,727 个涵盖19,354 种蛋白质的人类PPI(蛋白质相互作用)。研究证明这些PPI反映了蛋白质之间的直接(物理)和间接(功能)关联。其中每个PPI 包含两个Ensembl ID,分别代表蛋白质pa和pb,以及一个范围在150 和999 的得分S( )pa,pb,代表它们的相互作用强度。基于这些数据,我们构造了一个无向加权的PPIN,包含19,354个节点和5,879,727条边。
在本研究中,我们提出了一种基于PPIN 的癌症转移基因识别方法,以鉴定肺癌骨转移特异性关键基因。首先,结合收集的肺癌和骨癌基因数据,在PPIN 上执行RWR 算法,对基因进行预选。然后,通过置换检验消除网络结构的影响,得到候选基因集。最后,利用交互得分规则和富集分析对基因筛选,增强结果的准确性,得到肺癌骨转移关键基因集。整个方法的过程如图1所示。
图1 方法流程图
RWR 算法是一种经典的排序算法,它从一些种子节点开始,模拟其在网络中随机游走和重启,同时更新所有节点的概率得分并对节点进行排名[7]。它已被用于解决疾病基因的发现和药物重定位等问题[5~6]。RWR算法的主要过程如下:
输入:PPIN 的列归一化的邻接矩阵A,初始概率得分向量P0={Ps1,Ps2,…,Psn} (n=19354)
初始化:将S1与S2中基因整合并删去重复的基因,得到682 个节点作为种子节点,它们在P0中的概率得分设为1/682,其他节点的初始得分设为0;令重启概率r=0.8
过程:Fori=0 do
十月怀胎,真的不容易。尽管小心翼翼,在怀孕期间还是出现了高血压和其它并发症,经过保胎治疗,两个孩子在子宫内生长发育着,这让我饱含憧憬。
执行迭代Pi+1=(1 -r)APi+rP0(1)
直到‖Pi+1-Pi‖L1<10-6
End
输出:Pi+1中概率得分大于阈值10-5的节点对应的基因集合
算法的最终结果表示种子节点在网络中随机游走到其他节点的概率,体现了种子节点与其他节点在PPIN 中的相似性。因此,具有较高概率得分的基因与已验证的骨癌和肺癌基因更相关,从而更有可能是转移相关基因。概率得分大于阈值10-5的基因最终被筛选出来,这些基因统称为RWR基因。
通过RWR 算法得到的基因可能会受到PPIN结构的影响,从而存在很多与癌症转移无关的基因。为了尽可能排除这些基因,我们提出了置换检验规则。
首先,我们将总置换数设为1000,即随机构建了1000 个Ensembl IDs 集合,记为E1,E2,…,E1000,每个集合包含682 个随机的基因Ensembl IDs。然后,通过将Ei(1 ≤i≤1000 )中的682 个基因设置为种子节点,在PPIN 上执行RWR 算法以获取每个RWR 基因的概率得分。对于每个RWR 基因,存在一个真实概率得分Ps(g)和1000 个随机概率得分Psi(g)。最后,对每一个RWR 基因g,计算p-value值如下:
如果随机概率得分Psi(g)普遍大于真实概率得分Ps(g),说明g更可能是因为网络结构而被选出来的假阳性基因。显然,p-value 值很高的RWR基因并不是与肺癌骨转移特异性相关的基因,应当被删除。由于0.05 是作为被广泛接受的统计学检验传统显著性水平的阈值,我们选择p-value 值小于0.05的RWR基因作为肺癌骨转移的潜在候选基因做进一步分析。
根据研究证明,PPI 中交互得分高的两个蛋白质更有可能具有相似功能[8]。我们可以利用这一信息筛选出同时与肺癌和骨癌基因在功能上相似的候选基因。对于每个候选基因g,计算它的最大-最小交互得分MMIS:
其中,S1与S2分别表示2.1节中的肺癌相关基因集合与骨癌相关基因集合,因此MMIS 较高的候选基因至少同时与一个已验证的肺癌相关基因和骨癌相关基因密切相关。在STRING 数据库中,900 是蛋白质之间的最高置信度值,因此选择MMIS 得分不小于900的候选基因做进一步研究。
基因本体论(GO)可以从分子功能、生物学过程和细胞成分三个方面描述给定的基因及其产物;京都基因与基因组百科全书(KEGG)数据库提供了多个基因之间的生物学代谢途径。与已知肺癌和骨癌基因共享相同或相似的GO terms 和KEGG通路的候选基因更有可能是与转移相关的基因[9]。首先,根据富集分析的结果计算每个候选基因g与所有GO terms 和KEGG 通路的关系值,得到向量ES(g)。对 于 两 个 基 因g与g′ 在GO terms 和KEGG 通路上的富集分析相似性得分可以通过余弦定理计算:
具有更高Δ(g,g' )值的两个基因通常在分子功能和生物学过程等方面有很强的相关性。对于每个候选基因g,再计算最大-最小富集得分MMES:
在本研究中,我们尝试将0.9 作为MMES 的阈值,即筛选出MMES 大于0.9 的候选基因作为最终的转移关键基因。在整个方法中,对于由RWR 算法和置换检验规则产生肺癌骨转移候选基因,通过交互得分规则和富集分析进行评估,选择MMIS 不小于900 并且MMES 大于0.9 的基因作为肺癌骨转移关键基因,这些基因被认为在肺癌骨转移中发挥了重要作用。
如3.1 节所述,我们将与肺癌和骨癌相关的682 个基因作为种子节点,在PPIN 上执行RWR 算法,筛选概率得分大于10-5的基因后,得到了6850个RWR 基因。其次,我们采取了置换检验规则来消除网络结构对结果的影响,得到了964个p-value值小于0.05的候选基因做进一步研究。
为了更准确地识别肺癌骨转移相关基因,我们通过交互得分规则和富集分析测试对候选基因进行了评估与筛选。通过计算,对于每个侯选基因得到了一个MMIS 和MMES,我们选择MMIS 不小于900 并且MMES 大于0.9 的12 个基因作为肺癌骨转移关键基因,如表1 所示。文献挖掘的结果证明这些关键基因基因大部分参与了肺癌骨转移的发展过程,与肺癌骨转移特异性显著相关。
表1 12个肺癌骨转移关键基因及其概率得分、P-value、MMIS及MMES值
在获得的12个肺癌骨转移关键基因中,有9个基因已被证实与肺癌骨转移的形成或发展有关,其中包括骨髓毛细血管的侵袭和外渗,对趋化因子的反应以及对骨细胞外基质的粘附等。根据以往的研究,肺癌细胞的上皮细胞-间质细胞转化(EMT)过程和骨微环境的改变被认为是肺癌形成骨转移的关键因素[10~11]。大多数潜在的关键基因都直接或间接地参与了这两个过程,体现了它们在肺癌骨转移中的特殊作用。
根据12 个肺癌骨转移关键基因的基因家族,我们将它们分为5个簇,如图2所示,并进行了相应的分析。其中,MDM2 与许多癌症的发病机制有关,它刺激基质金属蛋白酶(MMPs)的表达,促进骨髓窦细胞外渗,有利于肺癌细胞通过新生血管进入血液循环,对肺癌骨转移有特异性作用[12~13]。此外,CD44 同样上调了MMPs 的表达,对肺癌细胞在骨骼组织中的适应性和侵袭性起着重要作用[14]。
图2 12个关键基因的基因家族分布
EMT过程发生在肺癌骨转移的初始阶段,有利于降低细胞间黏附力,加速相邻细胞脱落。其中,BMP7、CTBP1 基因调控并参与肺癌细胞的EMT 过程[15~16],表明了它们在肺癌骨转移过程中的影响。骨微环境是调节骨组织并维持其动态平衡的重要环境,肺癌骨微环境的改变在骨转移的进展中起着重要的作用。其中,APC、PROCR 及IL6 基因参与了骨微环境的改变过程[17~18],为骨转移瘤提供生长所需营养物质,有利于肺癌细胞在骨组织中的生长和扩散。此外,原癌基因MET 可以激活骨微环境中RANK 信号通路,诱导破骨细胞的活化,最终导致溶骨性转移的发生[19]。RAF1 是一种参与RAS信号通路的功能性原癌基因,被广泛报道参与癌症转移过程[20]。
在尚未确定的3 个基因中,UBE2C 是与泛素相关的基因,编码细胞周期进展所需的蛋白质,在骨髓中广泛表达[21]。NOTCH3 在调控肿瘤细胞的凋亡、增殖的分化中起着重要作用,是多种肿瘤治疗的潜在靶标[22]。PAK1 参与细胞粘附、迁移、增殖、凋亡、有丝分裂等多种细胞生物学过程,促进肺癌细胞的增殖及侵袭能力[23]。这些基因可能是潜在的肺癌骨转移相关基因,值得进一步研究。
肿瘤转移是一个复杂的过程,通常是促进肿瘤加重的主要原因。肿瘤转移相关基因的鉴定可为肿瘤转移的治疗提供分子靶点,有助于癌症患者的治疗和预后。在本研究中,基于两种相互作用更强的蛋白质更可能具有相似功能的假设,我们在PPIN 上设计了一种综合方法来识别癌症转移相关的基因。我们将该方法运用于肺癌骨转移相关基因的鉴定,最终获得了12 个潜在的肺癌骨转移关键基因并进行了广泛分析。结果表明,大多数鉴定的基因已被证实有助于肺癌骨转移的进程,体现了该方法的有效性和合理性。我们希望这一贡献将有助于识别肿瘤转移特异性基因,并为肿瘤转移的机理研究提供启示。