生物信息学分析食管癌关键基因的表达及其临床意义

赵立然，卜 梁*

(1.厦门大学附属翔安医院胸外科，福建 厦门 361100；2.厦门大学医学院，福建厦门 361100)

食管癌是世界上常见的癌症之一[1]，是世界上第8常见的癌症，也是导致癌症死亡的第6 大原因。因其具有较强的侵袭、迁移能力和较差的预后，5 年生存率仅为15%～25%[2]。食管癌的危险因素包含环境、遗传和表观遗传因素，如吸烟和饮酒，以及改变细胞生长的表观遗传和遗传调节的变化[3]。食管癌早期通常无症状，但单独吞咽困难或伴有体质量意外减轻是最常见的症状。由于早期诊断的筛查方法不够有效，大多数人在确诊时失去了根治性切除的机会[4]。此外，长期服用抗癌药物(包括顺铂)的患者通常会出现耐药性，也会导致癌症复发。化疗、放疗和手术治疗后仍有许多患者出现疾病进展和复发[5]。新的治疗方法和肿瘤分子标志物的发现对确定治疗靶点和改善患者预后具有重要意义。因此，研究食管癌的病理机制，寻找早期检测标志物，提高患者的生存时间和预后迫在眉睫。

1 材料与方法

1.1 数据来源

选取GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)中食管癌的基因表达谱芯片GSE38129[6]和GSE20347[7]。GSE38129 包含30 例食管癌组织样本，30 例癌旁正常组织样本。GSE20347 包含17 例食管癌组织样本，17例癌旁正常组织样本。

1.2 数据处理

使用GEOquery 包在基因表达综合数据库(Gene Expression Omnibus，GEO)中下载两个食管癌芯片数据集GSE38129，GSE20347。首先对两个芯片数据集进行整合，再用R 语言中的sva 包去除批次效应，然后对标准化数据集进行主成分分析(principal component analysis，PCA)。进一步，对标准化后的基因芯片数据集进行差异分析，使用Limma 包，设置差异表达基因(differently expressed genes，DEGs)阈值为|log2(FC)|>1.5 且P<0.01 进行筛选，以降低假阳性结果的产生，FC指基因表达的差异倍数(fold change)。

1.3 DEGs的富集分析

利用clusterProfiler包对上调和下调DEGs分别进行基因本体(gene ontology，GO)功能和京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes，KEGG)通路富集分析，以P<0.05作为纳入标准。

1.4 蛋白互作网络分析

将筛选出的DEGs导入STRING(https://cn.string-db.org/)网站，进行权重算法分析并构建蛋白互作网络(protein-protein interaction，PPI)，选取combined score≥0.9(highest confidence)作为纳入标准，将结果导入Cytoscape 软件，应用MCODE 及CytoHubba 插件提取核心模块及Hub gene。

1.5 关键基因在UALCAN 数据库里的信息检索与表达差异分析

分析DEGs 在多种肿瘤中的表达情况，及其在食管癌和正常组织中的表达水平以及与分期、性别、年龄、TP53 突变的关系，并探索DGEs 表达与自身启动子甲基化水平的关系。筛选步骤为：①UALCAN 网站输入关键基因名称；②探索泛癌视图，查看食管癌和其他癌症中关键基因的表达情况；③探索关键基因在该数据库纳入的食管癌患者与正常人群中的表达，选择分期、性别、年龄、启动子甲基化水平、TP53突变选项分析。

1.6 关键基因表达与疾病预后相关性

使用GEPIA(http://gepia.cancer-pku.cn/)数据库对经UALCAN 数据库筛选出的相关基因在食管癌患者和正常人群中的表达再次验证，绘制Kaplan-Meier 生存曲线，探究关键基因表达是否影响食管癌患者的预后。

2 结果

2.1 评价批次效应去除结果、评估标准化数据集合理性与DEGs筛选结果

A：部分数据集单个样本表达量；B：部分标准化数据集单个样本表达量.

R语言中的sva包对数据集去除批次效应后得到标准化数据集，使用箱线图来评估批次效应去除效果，标准化数据集样本表达量分布情况比原始数据集更加集中，部分标准化数据集与数据集中单个样本表达量的数据如图1 所示。对标准化数据集进行主成分分析(PCA)，如图2显示该数据集中食管癌组与正常对照组组间差异大，绝大部分组内差异较小在正常误差范围内，从而可进一步分析。在标准化数据集GSE38129_GSE20347 中，设置筛选标准为|log2(FC)|>1.5且P<0.01，共筛选出390个显著差异表达基因，其中上调基因166 个，下调基因224 个。根据筛选出的390 个显著差异基因，利用Volcano Plot 包绘制火山图，对这些基因进行可视化展示如图3 所示。在上调和下调基因中挑选log2(FC)降序排名前10 基因，绘制差异基因聚类热图进行可视化展示如图4。

图2 标准化数据集的主成分分析

图3 390个显著差异基因火山图

图4 差异基因聚类热图

2.2 DEGs的GO和KEGG富集分析

对上调和下调基因进行GO 和KEGG 富集分析，GO 分析表明上调基因参与的生物进程有：有丝分裂细胞周期相变、细胞外基质组织、细胞外结构组织生成等，下调基因与表皮发育、角质细胞分化、皮肤发展等功能密切相关(均为P<0.05)，部分上、下调基因分别GO富集结果的信息见表1，部分DEGs的GO富集结果如分子功能(molecular function，MF)、生物过程(biological process，BP)和细胞组成(cellular component，CC)可视化见图5。

图5 部分差异基因的GO分析结果

KEGG 通路富集分析结果显示，上调基因主要涉及细胞周期、细胞外基质受体相互作用、阿米巴病等信号通路(均为P<0.05)，部分上调基因的KEGG 富集信号通路可视化如图6 所示。下调基因主要与视黄醇代谢、不饱和脂肪酸的生物合成、5-羟色胺能神经突触等信号通路密切相关(均为P<0.05)，部分下调基因的KEGG 富集信号通路可视化如图7 所示。部分上、下调基因分别KEGG富集结果的信息见表2。

表2 食管癌组织差异基因KEGG富集分析结果

图6 上调差异基因KEGG分析结果

图7 下调差异基因KEGG分析结果

2.3 DEGs的PPI、中心模块及核心基因筛选

将DEGs 导入STRING 网站构建PPI，得到的结果导入Cytoscape 软件，使用MCODE 插件，对所得网络进行提取分析，得到了4 个中心模块，选最主要的模块分析，使用CytoHubba 中的MCC 算法，筛选出评分最高的前20 个核心基因(CDK1、ASPM、TOP2A、AURKA、CCNB2、CDC20、AURKB、NUSAP1、KIF20A、BUB1、DLGAP5、TPX2、BUB1B、KIF2C、NDC80、PBK、TTK、NEK2、PRC1、KIF4A)，它们之间的关系如图8 所示，基因间联系密切，可能在食管癌的发生、发展中起着协同作用。

图8 20个关键基因蛋白互作网络

2.4 CDK1、TOP2A、AURKA在GEPIA与UALCAN数据库里的表达差异分析

在UALCAN数据库中对这20个关键基因进行泛癌验证，筛选出CDK1、TOP2A、AURKA在食管癌中的表达显著高于正常组织(均为P<0.05，见图9A、10A、11A)。在GEPIA数据库中，CDK1、TOP2A、AURKA在食管癌中的表达也显著升高(均为P<0.05，见图9B、10B、11B)。在UALCAN 数据库中，CDK1、TOP2A、AURKA与正常组织比较，在I～IV期表达均升高，且各期与正常组织表达水平间差异均有统计学意义(P<0.05，见图12A～C)。表达在性别、年龄、TP53 突变状态方面与正常组织比较，差异均有统计学意义(P<0.05，性别分类见图12D～F，年龄分类见图12G～I，TP53突变状态分类见图12J～L)，甲基化分析结果显示，食管癌组织中CDK1、TOP2A、AURKA启动子甲基化表达水平升高，差异具有统计学意义(P<0.05，见图13A～C)。以上指标均证明CDK1、TOP2A、AURKA与食管癌的发展有关。

图9 CDK1在食管癌肿瘤组织和正常组织的表达差异

图10 TOP2A在食管癌肿瘤组织和正常组织的表达差异

图12 CDK1、TOP2A、AURKA表达水平与临床资料相关性

图13 CDK1、TOP2A、AURKA在肿瘤组织和正常组织中的启动子甲基化水平

2.5 生存分析

把UALCAN 得到的相关基因输入GPEPIA 数据库进行验证，结果显示CDK1、AURKA与食管癌预后密切相关，CDK1高表达组总生存期(overall survival，OS)显著低于低表达组(HR=1.8，P=0.036，见图14A)，AURKA高表达组OS 亦显著低于低表达组(HR=2，P=0.033，见图14A)，利用芯片数据集GSE70409 对3 个基因进行再次验证，CDK1、TOP2A、AURKA在食管肿瘤组织中表达显著上调，与标准化数据集趋势表现一致(见表3)。

表3 食管癌中CDK1、TOP2A及AURKA在3个数据集的差异表达

图14 CDK1、TOP2A、AURKA表达水平与食管癌预后的生存曲线

3 讨论

本研究通过生物信息学方法，对食管癌基因芯片GSE38129、GSE20347 进行整合并标准化处理，共得到390个DEGs，其中166个上调，224个下调。GO和KEGG 富集分析结果显示上调基因涉及有丝分裂细胞周期相变、细胞外基质组织构成等功能以及细胞周期、细胞外基质受体相互作用等信号通路。下调基因与表皮发育、角质细胞分化等功能和5-羟色胺能神经突触等信号通路密切相关。通过PPI和cytohubba 插件选出20 个关键基因，在UALCAN 和GEPIA 数据库中进行表达差异分析及生存分析验证，显示基因CDK1、TOP2A及AURKA可能与食管癌的发生、发展密切相关。

细胞周期蛋白依赖性激酶(CDK)属于蛋白激酶家族成员一种；其中，CDK1通过G2/M期转变和激活同源重组(HR)DNA修复途径在细胞周期进程中发挥关键作用。癌细胞以永久性增殖、分裂，不受控制为特点，因此CDK1 在其生存中充当着不可替代的角色。细胞周期中有序的G2/M 转换由CDK1/CCNB 复合物控制。CDK1 参与多种癌症的发展，例如乳腺癌、胶质瘤、肝癌等。Zou 等[8]研究表明CDK1 中mRNA 表达在包括HCC 在内的多种肿瘤组织中上调。CDK1 的高表达与HCC 患者的较差预后相关，CDK1 较低的启动子甲基化可能导致HCC肿瘤组织中较高的表达水平。目前，已经开发了多种用于癌症治疗的CDK1 抑制剂，这些抑制剂可诱导延长的G2 停滞和/或使细胞对肿瘤细胞中的DNA 损伤剂敏感，从而导致细胞死亡[9]。CDK1 药物抑制剂包括：Rigosertib(III 期)、BEY-1107(II 期)等，药物临床适应症主要集中在胰腺癌和胶质瘤[10]。文献显示，食管癌中的CDK1明显上调，BIRC5可能通过影响细胞周期通路在ESCC 中发挥重要作用，而CDK1 可能是该通路的枢纽基因[11]。而本研究也发现CDK1 在食管癌中的表达明显上调，与上述文献结果一致，但有关CDK1 在食管癌中的实验及药物研究并未见报道，需要开展进一步研究予以验证。

DNA拓扑异构酶II Alpha(TOP2A)是一种蛋白质编码基因，是复制过程的核心，已发现在包括乳腺癌在内的多种恶性肿瘤中失调。TOP2A是染色质拓扑结构的保守调节剂，可催化可逆的DNA 双链断裂(DSB)，对于在转录、复制和细胞分裂等多种动态过程中维持基因组完整性至关重要。有研究表明TOP2A在乳腺癌组织中高度过表达，TOP2A的过表达与较差的总生存期和无复发生存期相关[12]。此外，TOP2A 显示与肿瘤间质高度相关，尤其是与骨髓来源的抑制细胞。Liu等[13]研究发现MDM4 和TOP2A 相互结合并在翻译后表达水平上调，导致TOP2A 蛋白稳定、p53 抑制和肿瘤细胞增殖增加，揭示了MDM4和TOP2A的新功能以及它们在肿瘤发生中的相互作用，表明抑制MDM4-TOP2A相互作用可能代表了一种新的策略，可以特异性并同时靶向TOP2A 和MDM4 进行癌症治疗。Hailati 等[14]通过生物信息学分析得出TOP2A 可能为食管癌患者诊断和预后判断的标志物。

细胞周期调节激酶(aurora kinase A，AURKA)属于丝氨酸/苏氨酸激酶家族，其激活对于通过调节有丝分裂的细胞分裂过程是必需的。Aurora 激酶家族由三种丝氨酸/苏氨酸激酶Aurora-A/B/C 组成。其中，Aurora-A和Aurora-B在有丝分裂中发挥核心作用，而Aurora-C 在减数分裂中发挥独特作用。激酶的过度表达或基因扩增已在广泛的人类恶性肿瘤中得到报道，如AURKA 的异常扩增表达与人类结直肠癌、肺癌和白血病的化疗耐药密切相关[15]，表明它们在肿瘤发生中作为有效致癌基因的作用。目前已经产生了许多激酶抑制剂(AKI)；其中一些正在接受临床评估。AURKA 的过表达已被证明会导致染色体畸变和基因组不稳定[16]。细胞实验证实，USP3和AURKA在ESCC细胞中的高表达促进了肿瘤细胞的增殖和侵袭[17]。有研究报道了AURKA在人ESCC组织和细胞系中的表达上调，这种上调导致预后不良[18]。Aurora-A 可能在多西他赛化疗敏感性中发挥重要作用，抑制其表达可能是ESCC的潜在治疗靶点[19]。

差异基因的GO 功能和KEGG 富集分析显示、有丝分裂细胞周期相变、表皮发育、细胞外基质组织、细胞周期、阿米巴病等信号通路与食管癌的发生、发展机制有着密切的联系。本文局限之处在于，未对这些基因对于食管癌的具体临床意义进行深究和验证，但本研究为食管癌发病机制的研究提供了新的思路，后续将开展进一步的实验验证。