基于大数据的前列腺癌生物信息学分析

李志标，唐福才，卢泽潮，黄伟娜，何朝辉

（1.广州医科大学第三临床学院，广东广州 511436；2.中山大学附属第八医院泌尿外科，广东深圳 518033；3.广州医科大学第一临床学院，广东广州 511436）

前列腺癌（prostate cancer，PC）是泌尿系统中的常见肿瘤之一。2019 年美国预计有174 650 例新诊断前列腺癌病例和31 620 例死亡病例，其发病率位居男性恶性肿瘤发病率的首位，患者死亡率高居男性恶性肿瘤的第二位［1］。而在中国，前列腺癌的发病率也在逐年增高，已成为男性最常见的十大肿瘤之一，其死亡率也在随着年龄的增长而上升［2］。前列腺癌起病隐匿，临床发现多为晚期［3］且目前治疗方式主要以手术治疗和雄激素阻断治疗为主，但激素依赖性前列腺癌大部分会恶化发展为去势抵抗性前列腺癌，目前尚不能有效地治疗这种病情［4］。因此，前列腺癌的早期诊断，有助于提升患者的生存率，筛选出前列腺癌的生物标志物对其预防和前期的诊治具有重要的意义。人类肿瘤相关的基因表达汇编（gene expression omnibus，GEO）和癌症基因组图谱（the cancer genome alas，TCGA）是目前最大的两个公共组学数据库，近年来，国内外学者通过生物信息学的方法对这两个数据库的数据进行了多项研究［5-7］。生物信息学是目前自然科学的研究重点之一，它以计算机作为工具，从核酸和蛋白质序列进行研究，分析序列当中蕴含的结构功能信息。本研究结合两大数据库的基因数据集，筛选共同的差异表达基因，并进一步进行生物信息学分析，为前列腺癌的发生机制与预后预测提供依据。

1 材料与方法

1.1 材 料

以“prostate cancer；prostate carcinoma；prostate neoplasms；prostatic neoplasm；cancer of the prostate；prostatic cancer”作为关键词，从GEO 数据库（www.ncbi.nlm.nih.gov/geo）中搜索基因芯片数据集，筛选标准：①数据集为mRNA 基因组芯片；②具有癌组织与癌旁组织；③以Affymetrix Human Genome U133 Plus 2.0 Array 平台作为研究平台；④样本数＞20。基于以上筛选标准，基因芯片数据集GSE55945 和GSE46602 纳入研究。GSE46602 有36 个癌组织样本，14 个癌旁组织样本；GSE55945 有13 例癌组织样本，8 例癌旁组织样本。为了进一步利用不同基因组库的大数据集进行验证，从TCGA 数据库（http：//carcinomagenome.nih.gov/）下载前列腺癌的RNA-seq 数据及相应的随访数据，剔除TCGA 数据库中临床随访数据缺失或不全的样本数据，筛选出495 例前列腺癌基因数据与98 例癌旁数据及其对应的临床随访数据。

1.2 数据处理与基因筛选

利用GEO2R（https：//www.ncbi.nlm.nih.gov/geo/geo2r/）在线分析网站对GEO数据库中的GSE46602和GSE55945 基因芯片进行数据处理，筛选差异基因；同时采用R 语言的edgeR 包对TCGA 数据库下载的RNA-seq 数据进行差异表达分析。差异基因的阈值设为log2（fold change）＞1，P.adjust＜0.05。三个数据集差异分析结果绘制火山图，并采用韦恩图取并集。

1.3 通路与功能富集分析

利用R 软件的clusterProfiler 包［8-10］对上述共同差异基因进行GO（gene ontology，GO）功能分析及KEGG（kyoto encyclopedia of genes and genomes，KEGG）通路分析，以P.adjust＜0.05 为阈值筛查差异基因的主要富集功能和通路。

1.4 蛋白互作网络分析

将交集的差异基因导入string（https：//stringdb.org/）在线分析网站进行蛋白质互作网络分析（protein-protein interaction，PPI），将所得结果导入Cytoscape 软件，筛选出关键的蛋白表达模块与关键节点基因。

1.5 生存分析

利用R 软件survival 包，将1.4 中获得的前列腺癌关键节点基因表达值按中位数分为高表达和低表达组，结合TCGA 数据库中下载的前列腺癌患者预后随访信息，包括总生存时间和生存状态，对两组绘制K-M 生存曲线，并进行log-rank 统计检验，分析关键节点基因的预后关系，P＜0.05 为差异有统计学意义。

2 结果

2.1 前列腺癌差异表达基因的筛选

利用GEO2R 在线网站对GEO 筛选的数据集进行数据处理分析，GSE55945 基因芯片筛选结果显示，表达上调的基因622 个，表达下调的基因237 个；GSE46602 基因芯片筛选结果显示，表达上调的基因934 个，表达下调的基因542 个。利用R软件的edgeR 包对TCGA 下载的RNA-seq 数据进行基因差异表达分析，结果显示，表达上调的基因1 958个，表达下调的基因1 759个。利用R语言的Volcano Plot 包做火山图，并利用Venn Diagram 包进行基因合集处理（图1）。

2.2 差异基因的通路与功能富集分析

利用R 语言的clusterProfiler 包对筛选出来的278 个共同差异表达基因进行GO 功能分析与KEGG 通路富集分析，以P.adjust 从小到大进行排列，展示GO与KEGG 富集分析中的前10富集结果（图2）。结果显示差异基因主要与上皮细胞的增殖和含苯化合物的代谢过程等功能相关（表1），主要富集表达于谷胱甘肽代谢和粘着斑等信号通路（表2）。

表1 前10 的GO 富集分析结果Table 1 Top 10 GO enrichment analysis result

表2 KEGG 通路富集分析结果Table 2 KEGG pathways enrichment analysis result

图1 差异表达基因火山图和韦恩图Fig.1 Volcano Map and Venn Diagram of differentially expressed genes

图2 功能与通路富集分析结果Fig.2 Function and pathway enrichment analysis results

2.3 蛋白互作网络分析

将交集的差异表达基因导入string 在线分析网站进行蛋白互作网络分析，结果进一步导入Cytoscape 软件进行PPI 网络可视化处理（图3），利用cytoHubba 插件对PPI 网络中的基因进行节点得分计算，以Degree≥6 作为阈值，筛选得到前12个关键基因CAV1、ANXA1、MME、EDNRB、PTGS2、BIRC5、AMACR、PGR、EDN3、GSTP1、CYP3A5及CCK基因，并利用MCODE 插件进行重点蛋白表达模块的筛选，以默认参数：Degree≥6，Node Score Cutoff≥0.2，K-core≥2 及Max.Depth=100 作为筛选标准（图4）。

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

图4 重点蛋白表达模块Fig.4 Keyprotein expression modules

2.4 生存分析

为了探究关键节点基因与前列腺癌的预后关系，利用TCGA随访数据结合关键节点基因表达值进行生存分析，分析结果发现EDNRB、EDN3以及AMACR与前列腺癌预后有紧密联系（P＜0.05）。前列腺癌患者中低表达EDN3、EDNRB以及AMACR的患者生存率明显降低（图5）。

图5 生存曲线Fig.5 Survival curves

3 讨论

前列腺癌是泌尿系统好发的癌症之一，全球发病率每年皆呈上升趋势。据文献报道，2015 年我国前列腺癌新增人数已超过6 万人，死亡人数已达26 600 人［2］。前列腺癌有易转移、预后差的特点，因此，筛选前列腺癌的生物标志物对其预防和前期的诊治具有重要的意义，目前，前列腺癌的早期诊治的生物标志物前列腺特异抗原PSA（prostate specific antigen，PSA）已广泛应用于早期筛查中［11］。而在近年来，随着肿瘤分子生物学的发展，基因的靶向治疗也已逐渐得到人们的重视，它具有特异性抗肿瘤的作用，且副作用小，是肿瘤早期筛查与诊治的新技术新方向，具有很好的发展前景［12］。

本研究立足于GEO 和TCGA 两大组学数据库的基因数据集，采用生物信息学分析方法，对前列腺癌的基因数据进行分析和处理，筛选出共同差异基因278 个，差异基因主要汇集于上皮细胞的增殖和含苯化合物的代谢过程等相关功能，并且主要富集于谷胱甘肽代谢和粘着斑等信号通路。同时对交集的差异基因进行PPI 网络分析，结果显示CAV1、ANXA1、MME、EDNRB、PTGS2、BIRC5、AMACR、PGR、EDN3、GSTP1、CYP3A5及CCK这12 个基因可能在前列腺癌发生的蛋白表达调控机制发挥重要的作用。PPI 网络筛选以及生存分析结果显示，EDN3、EDNRB与AMACR很可能在前列腺癌的发生发展及预后中发挥着重要作用。文献报道，目前美国的前列腺癌患者的5 年生存率接近100%［13］，因此，对前列腺癌远期预后的评估尤为重要。AMACR是CAIB-BALF 辅酶A 转换酶家族成员之一，编码基因为P504s，是目前前列腺癌诊断的阳性标志物之一［14-15］。Eryilmaz 等［16］的结果显示，前列腺癌中AMACR表达明显上调。研究表明，AMACR还可作为前列腺癌的治疗靶点，Festuccia 等［17］证明抑制AMACR能有抗前列腺癌的作用，AMACR可作为一个有效的药物作用靶点；Honma 等［18］证明AMACR可作为前列腺癌细胞毒性T 淋巴细胞免疫疗法的新靶点。本研究发现AMACR在前列腺癌中表达明显上调，与上述文献结果一致，但有关AMACR的表达情况在前列腺癌患者的预后评估未见报道，需要进一步的研究验证。EDN，全称为人内皮素（human endothelin，EDN），据调查，EDN 主要由大脑神经细胞、肾小管上皮细胞、前列腺细胞和结肠细胞等细胞分泌得出，主要作用于内皮素受体EDNRA 和EDNRB 受体［19-20］，EDN在多种癌细胞中有着一定的调节作用，在细胞增殖、侵袭以及扩散等过程中起着关键的作用。EDN3与EDNRB的亲和力较高，易结合发挥细胞信号传导的作用［21］。在癌症中，EDNRB可通过促进细胞的凋亡作用及清除ET-1（endothelin-1，ET-1）的效用来拮抗肿瘤的进程［22］，文献结果显示，ET-1在癌症中的明显高表达，参与癌症的自分泌和旁分泌途径的调控［23］，这表明在前列腺癌中EDN3/EDNRB可能抑制癌症进展。并且文献报道显示，EDNRB在胃癌、结肠癌及前列腺癌具有高甲基化表达的现象，与癌症的进展密切相关，具有重要的研究价值［24-26］。上述报道表明EDN3/EDNRB信号通路能对多种癌症的发生进行调控，EDN3与EDNRB均在前列腺癌中低表达，而生存分析结果也显示，EDN3与EDNRB在前列腺癌中发挥着抑癌基因的作用。因此，综合文献挖掘的结果显示，EDN3与EDNRB很可能是前列腺癌的潜在靶基因，对前列腺癌的发生发展及患者的预后评估具有重要的意义，有望成为前列腺癌的早期筛查标记物或基因治疗靶点。

通过对交集的差异表达基因进行通路与功能富集分析结果可知，谷胱甘肽代谢和粘着斑等信号通路与前列腺癌的发生发展机制有密切的联系。研究表明，谷胱甘肽代谢过程中的关键酶谷胱甘肽S-转移酶（glutathione S-transferases，GST）可通过表达沉默参与前列腺癌的发生。Zhao 等［27］和Zhou 等［28］发现前列腺患者中GSTP1的启动子甲基化发生率明显高于非前列腺癌患者，可能是诊断前列腺癌的潜在生物标记。有学者认为GSTP1的高甲基化可导致其表达沉默，可作为前列腺癌发生早期诊断的重要靶点之一，可调节前列腺癌的发生发展进程［29-31］，且在多种癌症中都发挥着抑制癌症发展的作用，可用作早期癌症的检测或诊断的生物标志物［32-33］。虽然目前AMACR、EDN3与EDNRB与GST的关系尚未有报道，但我们的研究结果显示，EDN3/EDNRB在肿瘤的作用与其有着一定的相关性，可以有效抑制肿瘤的发展过程且在肿瘤中一般呈高甲基化现象，而AMACR是前列腺癌中的重要标志物，其高表达可能与GSTP1的高甲基化抑制有关，具体机制仍待进一步的实验探索。粘着斑激酶（focal adhension kinase，FAK）与多种癌症的发生发展功能相关，已有研究证明FAK 的表达强度与乳腺癌患者的远端转移程度有密切相关［34］，在胃癌中FAK信号通路传导的增强也会增强其侵袭能力［35］。据文献报道，FAK 通路的激活会导致多种激酶的激活，增强肿瘤的存活和增殖功能［36］，而FAK 抑制剂有着潜在的抗癌剂的功效［37］，文献结果显示，EDN3/EDNRB的抑制也可激活FAK 信号通路，增强肿瘤的发生发展进程［38-39］。而AMACR与FAK 间的关系尚未见报道，推测AMACR有可能在FAK 信号通路的激活上也存在着调控作用，导致癌症的进展加快，但仍需进一步的研究探讨。

综上所述，本研究通过对前列腺癌基因芯片数据集的收集，利用生物信息学分析，筛选出前列腺癌与癌旁组织表达差异的基因，并发现上皮细胞的增殖和含苯化合物的代谢过程等相关功能以及谷胱甘肽代谢和粘着斑等信号通路可能对前列腺癌的发生机制有着密切的关系，同时通过PPI网络分析与生存分析结果表明，EDN3、EDNRB与AMACR作为关键节点基因，很可能在前列腺癌发生发展与预后评估中起重要作用，为前列腺癌发生机制的研究以及新靶点的发现提供理论依据。