基于多重生物信息学分析对膀胱癌关键基因及其调控网络的研究①

卢俅 陈兵海（江苏大学附属医院，镇江212000）

膀胱癌是泌尿系统最常见的恶性肿瘤之一，2018年全世界约有55万新发病例，占所有恶性肿瘤的3.0%，发病率居恶性肿瘤第十位，全世界每年约有20万人死于膀胱癌［1］。在我国，膀胱癌发病率位于全身恶性肿瘤的第十二位，每年约有8万新诊断膀胱癌患者，约占全部恶性肿瘤的1.87%，死亡率位于恶性肿瘤的第十二位，约占所有恶性肿瘤死亡率的1.1%［2］。膀胱癌主要起源于上皮组织，约75%的膀胱癌属于非肌层浸润性膀胱癌，同时具有较高的复发率，患者预后较差［3-4］。目前膀胱癌的诊断主要依靠膀胱镜检查，膀胱镜检查虽然存在诊断阳性率较高的优势，但仍然有着巨大局限，如属于有创检查，存在一定的手术禁忌证等［5］。因此，目前迫切需要找到膀胱癌的新型分子标志物，用于膀胱癌的诊断及预后判断。基因表达综合数据库（Gene Expression Omnibus，GEO）和癌症基因组图谱（The Cancer Genome Atlas，TCGA）是目前世界上最大的两个癌症基因信息数据库，近年来，已有众多学者通过生物信息学的方法对这两个数据库展开研究［6-7］。本研究立足于GEO和TCGA两大数据库，筛选差异交集基因集，利用多重生物信息学分析工具，发现CDC20、TPM1、ACTA2、MYH11、MYLK以及CALD1可能是膀胱癌关键的核心基因及治疗靶点。

1 材料与方法

1.1 膀胱癌差异基因筛选从GEO数据库（www.ncbi.nlm.nih.gov/geo）下载2个基因芯片数据，分别为GSE13507和GSE7476。其中GSE13507的基因平台是GPL6102 Illumina human-6 v2.0 expression beadchip，它包括188个膀胱癌组织和68个正常膀胱组织，GSE7476的基因平台是GPL570 Affymetrix Human Genome U133 Plus 2.0 Array，它包括9个癌组织和3个正常组织。为了进一步利用不同数据库进行验证，又从TCGA数据库（https：//tcga.xena‑hubs.net/）下载膀胱癌的RNA-seq数据及相应临床随访数据。其中包括407个膀胱癌样本，19个正常膀胱组织样本。剔除临床随访数据缺失的样本，共获得405个有完整临床数据的膀胱癌样本。利用SangerBox数据分析软件中DECeter内置limma包进行差异基因的筛选，设置筛选条件：P<0.05，log2（FC）>1，得到3个差异基因集组成的数据集。利用SangerBox数据分析软件中火山图绘制工具绘制每个数据集的火山图，采用维恩图（http：//bioinformatics.psb.ugent.be/webtools/Venn/）获得交集基因。

1.2 GO和KEGG通路富集分析使用Metascape在线分析网站（https：//metascape.org/）对301个差异基因进行基因的本体论（Gene Ontology，GO）功能分析和京都基因与基因组百科全书（Kyoto Encyclo‑pedia of Genes and Genomes，KEGG）通路分析，设置P<0.05为差异具有统计学意义。

1.3 蛋白质互相作用网络（protein-protein interac‑tions network，PPI）分析为了了解膀胱癌的新陈代谢和分子机制的重要途径，将得到的301个差异基因集导入String在线分析网站（https：//string-db.org/）进行PPI分析，为了获得更直观的可视化，将分析结果中TSV格式的文件导入Cytoscape软件，筛选出关键蛋白表达模块与关键核心基因。

1.4 生存分析采用从TCGA下载的临床数据及基因表达矩阵，对筛选获得的核心基因，利用Graph‑Pad Prism 8.0.1软件进行生存分析，绘制K-M生存曲线（P<0.05表示差异具有统计学意义）。

2 结果

2.1 膀胱癌差异基因筛选利用SangerBox数据分析软件中DECeter内置limma包，设置筛选条件：P<0.05，log2（FC）>1。在GSE13507基因芯片获得在膀胱癌组织中75个表达上调基因，424个表达下调基因，GSE7476中膀胱癌组织中得到584个上调基因，1 520个下调基因，TCGA数据库中获得1 634个上调基因，2 671个下调基因。利用每个数据集的表达差异基因绘制火山图及维恩图（图1）。

图1 差异基因的火山图和维恩图Fig.1 Volcano Maps and Venn Diagram of differentially expressed genes

2.2 GO和KEGG通路富集分析使用Metascape在线分析网站进行GO和KEGG通路富集分析，结果显示301个差异基因主要与含胶原的细胞外基质、超分子纤维组织、肌肉结构发展、血管的发育、细胞基质黏附、生长发育等功能相关，主要富集表达于黏着斑和Wnt信号等通路，富集结果以P值从小到大排列（表1、2）。

表1 前10的GO富集分析结果Tab.1 Top 10 GO enrichment analysis results

2.3 PPI和核心基因筛选将301个差异基因集导入String在线分析网站，为了获得更直观的可视化，再将结果导入Cytoscape软件中（图2）。结果发现在PPI关系中，上调基因（红色标记）和下调基因（蓝色标记）相互作用关系很少，为了更好获得蛋白质相互作用的核心基因，先是利用Cytoscape软件的Mcode插件根据Degree得分筛选得到排名前六的蛋白质表达模块（图3）。再分别对排名前二重点模块进行核心基因的筛选，得到2组核心基因（表3、4）。

表2 KEGG通路富集分析结果Tab.2 KEGG pathways enrichment analysis results

表3 前10蛋白表达模块1的核心基因Tab.3 Top 10 hub genes of protein expression module 1

图2 蛋白质互相作用网络Fig.2 Protein-protein interaction network

图3 重点蛋白表达模块Fig.3 Key protein expression modules

2.4 生存分析为了探究2.3获得2组核心基因与膀胱癌预后的关系，采用TCGA随访临床数据及相应基因表达量进行生存分析，使用GraphPad Prism 8.0.1绘制K-M曲线，研究结果发现CDC20、TPM1、ACTA2、MYH11、MYLK以及CALD1与膀胱癌患者预后明显相关（图4），差异有统计学意义，其他核心基因差异无统计学意义。膀胱癌组织中CDC20、TPM1、ACTA2、MYH11、MYLK及CALD1表达量越高，患者总生存期越短。

图4 生存曲线Fig.4 Survival curves

2.5 核心基因对膀胱癌重要的诊断价值为了探究与生存明显相关的6个核心基因对膀胱癌的诊断价值，利用TCGA数据库基因的表达量对407例膀胱癌组织及19例正常膀胱组织进行分析，受试者操作特征（receiver operating characteristic，ROC）曲线使用SPSS Statistics 22.0软件绘制，结果提示：CDC20、TPM1、ACTA2、MYH11、MYLK以及CALD1对膀胱癌诊断的特异度和灵敏度均较高（图5）。

图5 ROC曲线分析显示CDC20、TPM1、ACTA2、MYH11、MYLK和CALD1可能作为膀胱癌诊断的分子标志物Fig.5 ROC curve analysis demonstrated that CDC20，TPM1，ACTA2，MYH11，MYLK，and CALD1 may be diagnostic biomarkers in patients with BC

表4 前8蛋白表达模块2的核心基因Tab.4 Top 8 hub genes of protein expression module 2

3 讨论

本研究立足于GEO和TCGA两大数据库，共获得301个差异的交集基因。这些基因在膀胱癌组织中表达上调的共有41个，下调的共有260个，随后将301个差异基因进行GO功能和KEGG通路富集分析。GO功能富集分析发现差异基因与整合素、细胞基质的黏附相关。猜测在膀胱癌中整合素通过对细胞和细胞外基质的黏附发挥调节作用，进而影响肿瘤的进展。在膀胱癌细胞中，MARTINO等［8］研究发现膀胱癌细胞基质黏附的改变可以促进膀胱癌细胞的侵袭、转移。通过KEGG通路富集分析，发现黏着斑信号通路可能与膀胱癌相关，已有研究证明FAK信号通路与多种肿瘤发生发展相关，在结直肠癌中FAK信号通路的激活会加强癌细胞的增殖、侵袭能力［9］。国外众多学者研究发现FAK途径的激活可以促进膀胱癌侵袭、增殖、转移［10］。通过对膀胱癌差异基因的GO和KEGG分析，猜测细胞基质黏附的改变及FAK信号通路激活可能是膀胱癌进展的潜在机制，可能为膀胱癌研究新型的治疗方案提供参考。随后将获得的301个差异基因进行PPI及生存分析，共得到6个与生存预后密切相关的核心基因。

在这些核心基因中，ACTA2的表达上调可以激活FAK信号通路，进而促进肿瘤的进展［11］，CDC20、TPM1、MYH11、MYLK、CALD1与FAK信号通路在肿瘤中关系未见相关报道。进一步研究发现AC‑TA2在膀胱癌组织中明显下降（P=1.68E-13），且通过生存分析发现膀胱癌组织中ACTA2表达量与患者的预后密切相关（P=0.035）。本研究得出ACTA2在膀胱癌中具有一定特异性，推测ACTA2在膀胱癌中通过FAK信号通路发挥一定调节作用，导致膀胱癌进展加快。GOULET等［12］研究发现ACTA2作为癌相关纤维母细胞的生物标志物之一，其在膀胱癌中表达量增高，患者预后明显变差，这与本研究结果相符。

CDC20是本课题组获得的另外一个核心基因，全称为细胞分裂周期20，是细胞分裂周期检查点的重要调控分子，在细胞有丝分裂的后期发挥重要的调控作用，其表达水平在恶性脑胶质瘤、肺癌是明显增多的［13-14］。CHOI等［15］研究发现，CDC20在膀胱癌组织中明显高表达，CDC20表达水平与患者高龄、更高级别的病理分期、更高的临床分期、近处转移等呈正相关。同时他们发现CDC20高表达的膀胱癌患者也预示着存在更短的无复发生存期（recur‑rence-free survival，RFS）和总生存期（overall survival，OS）。本次研究发现，CDC20在膀胱癌组织中明显高表达［log2（FC）=3.84，P=2.36E-29］，ROC曲线分析得出CDC20对膀胱癌诊断特异度和灵敏度特别高（特异度：0.947；灵敏度：0.779），生存分析发现膀胱患者CDC20表达量越高，患者总生存期越短，这与以往的研究结论相符。因此，认为CDC20可能是膀胱癌的一个良好诊断和预后判定的分子。

本研究发现TPM1、ACTA2、MYH11、MYLK以及CALD1在膀胱癌组织中均明显下降，为了进一步研究这些基因在膀胱癌中发挥的作用，将以上基因导入String在线分析网站。结果发现这些基因均富集在平滑肌收缩信号通路上（P=5.50E-14），我们猜测TPM1、ACTA2、MYH11、MYLK以及CALD1可能与血管平滑肌细胞的功能有关，通过影响细胞的黏附、迁移能力进而促进膀胱癌的进展。以TPM1为例，TPM1全称是原肌球蛋白1，是原肌球蛋白家族的一员，该家族高度保守，主要参与横纹肌、平滑肌的收缩系统及非肌肉细胞的骨架组成。近年来研究发现TPM1在肿瘤发生、发展中发挥重要的调节作用［16-17］。WANG等［16］发现TPM1的表达与肾透明细胞癌患者的肿瘤大小、吸烟状态、Fuhrman分级及预 后 相 关。LIN等［18］研 究 发 现TPM1通 过 调 节miR1835p.1促进胃癌细胞的增殖、转移和侵袭。LIU等［19］发现TPM1在膀胱癌组织中是低表达的，且可以抑制膀胱癌细胞的增殖，促进癌细胞的凋亡，但TPM1影响膀胱癌机制还不清楚。本研究发现，TPM1在膀胱癌组织中表达量明显减少（P=3.07E-15），对膀胱癌诊断的灵敏度较高（灵敏度：0.914，AUC=0.857），生存分析发现膀胱癌患者高表达的TPM1与较低总生存期呈正相关，GO富集分析发现TPM1与细胞基质黏附相关，因此推测TPM1通过降低膀胱平滑肌细胞或癌细胞间黏附力，加速了膀胱癌细胞的恶化和迁移，进而促进膀胱癌的进展。此外，已有相关研究提示MYH11和MYLK可能是膀胱癌的生物标志物［20-21］，但其调控膀胱癌机制目前尚不清楚。通过本次研究，猜测MYH11和MYLK可能是通过影响细胞的黏附、迁移能力进而促进膀胱癌的进展。CALD1全称为钙调素结合蛋白，是茶碱家族的靶点，参与细胞增殖和迁移的调控，CALD1与CDC20的GO富集分析结果均与姐妹染色单体分离有关［22］。因此，研究结果显示CALD1与CDC20可能通过调控细胞周期进而发挥对膀胱癌的调控作用。

总之，立足于GEO和TCGA两大数据库，本次研究运用多重生物信息学分析工具对膀胱癌进行分析，从多个角度全面探讨了膀胱癌中的关键基因及其相关调控网络，结果发现CDC20、TPM1、ACTA2、MYH11、MYLK和CALD1这6个核心基因对膀胱癌诊断的特异度和灵敏度均较高，可能是膀胱癌潜在的生物标志物，当然这些还需要进一步的体内外实验证实；膀胱癌中FAK信号通路激活、细胞的黏附、细胞周期的改变等可能是这些核心基因调控膀胱癌的分子机制，这为膀胱癌未来治疗提供了新的思路。