基于共表达网络分析鉴定膀胱癌预后相关核心基因

金敏 纪珊珊 徐珊珊 王淞 严林 张帆 刘佳伟

1华北理工大学研究生学院 河北唐山 063210；2唐山工人医院；3天津医科大学总医院滨海医院

膀胱癌(Bladder cancer，BLCA)是十大恶性肿瘤之一。据相关统计报道，全世界每年有42.9万例新发病例和16.5万例死亡病例[1]。早期膀胱癌患者的5年生存率高达到95.7%，而转移性患者的5年生存率仅为5%[2]。膀胱癌[3-4]按临床分期可分为两大类:非肌肉浸润性膀胱癌(non-muscle invasive bladder cancer, NMIBC)和肌肉浸润性膀胱癌(muscleinvasive bladder cancer, MIBC)。NMIBC具有成纤维细胞生长因子受体3(fibroblast growth factor receptor, FGFR3)突变的共激活、高复发率的特点，且约20%的患者可进展为MBIC[5]。MIBC的特征为肿瘤蛋白53(tumor protein53, TP53)的高突变率及转移率，患者一经发展为MIBC，则预示着不良预后[5]。目前膀胱癌缺乏能够准确预测其疾病进展和预后的有效生物标志物，因此迫切需要研究新的生物标志物以提高膀胱癌检测的诊断或预后性能。

1 材料与方法

1.1数据下载与处理 通过美国国立生物技术信息中心(NCBI)平台下的基因表达综合数据库(GEO)下载亚洲人群的芯片数据集GSE13507，包括58例正常组织和165例原发膀胱癌组织。通过TCGA(https://portal.gdc.cancer.gov/)在线下载BLCA的基因表达谱数据，包括患者的年龄、分期和生存时间等相应的临床信息。利用R软件的Limma包对数据集的差异基因(differentially expressed genes，DEGs)进行筛选，筛选标准：P<0.05，差异倍数>2.0。

利用R软件的WGCNA软件包对TCGA的表达数据进行共表达网络的构建，确定合适的软域值。将邻接矩阵转化为拓扑重叠矩阵(topological overlap measure，TOM )使用层次聚类法，对基因进行聚类，确定关键模块以筛选关键基因。将与临床特征具有高相关系数的模块筛选出来用于后续分析。

使用R包VennDiagram将筛选出的DEGs和从共表达网络中提取的共表达基因进行可视化展示，三者的重叠基因被鉴定为BLCA中具有潜在预后影响的关键基因。

1.2基因功能(GO)和京都基因与基因组百科全书(KEGG)富集分析 使用R包ClusterProfiler包来探索具有潜在预后影响的基因之间的功能，并以调整的P<0.05为标准，进行基因功能(gene ontology, GO)、京都基因和基因组百科全书(kyoto encyclopedia of genes and genomics, KEGG)分析，预测其生物学功能。

1.3PPI的构建和HubGenes的筛选 使用了STRING在线检索工具预测蛋白质-蛋白质相互作用(protein-protein interaction, PPI)，构建选定基因的PPI网络。使用STRING数据库，选择得分≥0.9的基因，将数据导入Cytoscape(v3.7.2)进行可视化展示。使用Cytoscape的插件CytoHubba计算MCC值最高的10个基因，作为中心基因进行后续分析。

1.4基因的预后价值验证 根据TCGA数据库中的数据，通过使用R软件中的生存软件包对患者的总体生存(OS)进行了Kaplan-Meier单变量生存分析。在此次研究中，仅选择随访时间完整的患者进行生存分析，然后根据中枢基因的中位表达值将其分为两组。Log-rankP<0.05被认为差异具有统计学意义。

2 结果

2.1加权基因共表达模块的构建 利用WGCNA软件包分析TCGA-BLCA患者基因的功能簇，模块的排序是基于TOM矩阵的基因分层聚类树状图，在本研究中总共确定了14个模块，绘制了模块与癌组织、癌旁组织之间关系的热图，见图1、2。表明TCGA-BLCA中的棕色模块与正常组织的相关性最高(模块：r=0.59，P=5e-2)。

图1 共表达网络的聚类树状图

图2 各个模块与临床的关系

差异基因和共表达模块之间的基因鉴定通过利用R软件Limma包对差异基因的筛选得到，在TCGA数据集图中共有3457个DEGs，见图3A；在GSE13507数据集图中有411个DEGs肿瘤组织中表达失调，见图3B。在TCGA数据集的蓝绿色模块和TCGA差异基因与GSE13507的差异基因中发现了共表达基因。总共提取了178个重叠基因用于后续分析，见图4。

图3 TCGA和GSE13507差异基因的火山图(P<0.05)

图4 TCGA-DEGs和GEO-DEGs与TCGA中共表达模块基因的维恩图

2.2178个基因的功能富集分析 为进一步了解与DEGs和TCGA共表达模块重叠的178个基因的潜在功能，使用ClusterProfiler包进行了基因富集分析。在筛选GO富集分析之后，观察到富集基因集，见图5。178个基因的生物过程(biological process，BP)主要富集在肌动蛋白纤维组织分化；对细胞成分(cell component，CC)的结果分析发现这些基因主要与肌原纤维、粘着斑有关；糖胺聚糖绑定主要富集在分子功能(molecular munction，MF)中。178个基因明显富集在丝裂原活化蛋白激酶(MAPK)信号通路与磷脂酰肌醇3激酶-蛋白激酶B(PI3K-Akt)信号通路上，见图6。

图5 重叠关键基因GO富集分析

图6 重叠关键基因KEGG分析

2.3PPI网络的构建和枢纽基因鉴定 利用STRING数据库构建了重叠基因的PPI网络，使用CytoHubba插件的最大相关原则算法(maximum correlation criterion，MCC)从PPI网络中选择出枢纽基因，见图7。根据MCC计算后，得分最高的前十个基因分别为：TPM1、TPM2、ACTG2、ACTA2、MYL9、MYLK、LMOD1、SORBS1、MYH11和GNG11。

图7 构建可视化PPI并筛选枢纽基因

3 讨论

免疫治疗在膀胱癌不断开展，但由于缺乏精确的分子靶标，患者的预后通常较差。因此，需要开发特定的生物标志物用于膀胱癌的预后和进展。WGCNA最初于2005年由Zhang和Horvath学者提出[6]，在高通量基因表达数据中构建具有高度相关性基因的表达模块，进而发现具有潜在生物学意义的基因，是鉴定恶性肿瘤核心模块和关键基因的理想方法；其强大的分析效能在生物医学研究中得到了广泛应用，尤其在肿瘤相关基因的高通量数据挖掘中。目前已有多项研究采用此方法对恶性肿瘤进行分析，如非小细胞肺癌[7]、低级别胶质瘤[8]、胰腺导管腺癌[9]等。

既往的研究多单独应用WGCNA寻找关键基因，本研究采用R-Llimma软件包筛选膀胱癌组织和癌旁组织的差异基因，并与WGCNA筛选共表达模块基因相结合的方式探索膀胱癌发生的核心基因，使结果更具稳健性。在这项研究中，在TCGA和GSE13507数据库中鉴定出178个关键基因。通过分析发现，这些基因与MAPK通路、PI3K-Akt通路高度相关。通过CytoHubba筛选出了与膀胱癌发生发展密切相关前10个基因：TPM1、TPM2、ACTG2、ACTA2、MYL9、MYLK、LMOD1、SORBS1、MYH11和GNG11。查阅文献后发现，TPM1、MYL9、MYH11、LMOD1可能是膀胱癌发生发展的关键基因。

原肌凝蛋白1(Tropomyosin1，TPM1)是高度保守的肌动蛋白结合蛋白家族成员，是编码高分子量TMs的亚型，可调控肿瘤细胞的增殖、侵袭、转移和运动[10-11]。研究表明，多种恶性肿瘤中TPM1的表达水平与生存相关，TPM1可能是肿瘤发生的潜在生物标志物[12]。肌球蛋白轻链9(Myosin Light Chain9，MYL9)是肌球蛋白超家族中的重要成员，其磷酸化会导致平滑肌细胞收缩。肌球蛋白重链11(Myosin Heavy Chain 11，MYH11)是由MYH11基因编码的平滑肌肌凝蛋白，属于肌凝蛋白重链家族[13]。MYH11是一种收缩蛋白，在主动脉组织中，MYH11的破坏可导致血管平滑肌细胞的丢失、紊乱和增生，是导致胸主动脉瘤和夹层的机制之一[14-15]。近年研究表明，MYH11基因的突变与癌症的发生密切相关，这可能是与MYH11参与肿瘤细胞迁移、细胞粘附蛋白的相互作用过程有关[16-17]。平滑肌细胞限制性基因1(leiomodin1，LMOD1)是LMOD家族的一员。此表型是肌动蛋白丝组装中断、细胞骨架排列和平滑肌收缩性受损的结果[18]。目前为止对LMOD1在实体肿瘤中的表达及生存分析较少，仅在平滑肌肉瘤中有过报道。本研究发现TPM1、TPM2、ACTG2、ACTA2、MYL9、MYLK、LMOD1、SORBS1、MYH11和GNG11可能与膀胱癌的发生发展相关，但作为膀胱癌诊断和治疗的潜在靶点，仍需要进一步实验验证。