R语言下食管鳞状细胞癌关键驱动基因富集分析的生物信息学研究

王子明，李新阳，姚 歌，王钰鲲，王新帅，原 翔，张广平

食管癌是常见的消化系统恶性肿瘤之一，其发病率和病死率在全球恶性肿瘤中分别居第 8 位和第6位，每年死亡人数可达400 000人[1-2]。食管癌在世界范围内呈现出显著的地域性和一定的民族差异，在欧美国家以腺癌居多，在亚洲国家以鳞状细胞癌为主，食管鳞状细胞癌(esophageal squamous cell carcinoma,ESCC)是食管癌的主要组织学亚型，约占90%[3]。尽管对食管癌的多种治疗方法已取得一定的进展，但食管癌仍然是一种恶性疾病，在西方国家总体5 a生存率为12%～20%[4]。近年来，不断有研究尝试食管癌生物靶向治疗手段，探索相关分子作用的靶点及相关信号，如针对表皮生长因子受体(epidermalgrowth factor receptor,EGFR)、针对血管内皮细胞生长因子(vascularendothelial growth factor,VEG)及VEGF受体(vascularendothelial growth factor receptor,VEGFR)、针对原癌基因人类表皮生长因子受体2(human epidermal growth factor receptor 2,HER2)等。但是，目前根据已发现的靶点而研究出的靶向药物较少。因此，本研究通过运用R语言及其相关软件包，在基因组水平上探索食管癌患者发生、发展的潜在分子靶点及相关信号途径，为进一步发现新的食管癌治疗靶标提供重要的理论依据。

1 材料与方法

1.1 GEO数据库及数据来源和处理GEO数据(http://www.ncbi.nlm.nih.gov/geo/)库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据。本研究运用R语言平台(3.6.0版本)下的Rstudio(3.6.0版本)和相关的软件包，从GEO数据库中获得数据GSE38129。GSE38129平台注释信息：GPL571[HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array。对微表达矩阵芯片数据进行去除缺失值、log2处理、标准化、取平均表达量最高的探针用于探针ID和基因名匹配，将筛选得到的基因表达矩阵数据用作本研究的数据来进行分析。

1.2 差异表达基因矩阵选定运用R语言下limma包(3.40.2版本)对肿瘤和正常样本组织筛选得到的微表达矩阵芯片数据进行进一步筛选，以“FDR<0.05及logFC≥2或log2FC≤-2”为筛选阈值，得到显著差异表达基因矩阵，并用R语言下ggplot2包(3.1.1版本)绘制火山图和热图来验证已得到的差异表达基因矩阵。

1.3 富集分析分析差异表达基因本研究对差异表达基因分别进行包括分子功能(molecular function)、细胞组成(cellular component)及生物过程(biological process)的GO功能富集分析和KEGG通路富集分析，并将所得数据分别通过ggplot2包和pathview包(1.24.0版本)进行GO功能富集分析可视化展示和KEGG通路富集分析可视化图示。

1.4 显著表达差异基因蛋白—蛋白相互作用(PPI)网络构建和分析STRING数据库(https://string-db.org/)是一个用于构建蛋白之间的相互作用网络并分析和探索预测的蛋白相互作用网络在线数据库，本研究将筛选得到的差异基因上传到STRING数据库(10.5版本)和Cytoscape(3.7.2版本)及MCODE插件用于差异基因PPI网络可视化分析，并以“combined score>0.4、MCODE>2为筛选阈值，得到相应的功能模块网络图，为了进一步分析关键基因，选取其中MCODE Score>5，得到2个功能模块网络图，从而得到潜在的关键基因。

1.5 生存预后分析在本研究中，所得关键基因均通过运用以The Cancer Genome Atlas(TCGA)数据库中81例ESCC患者数据为参考数据的Kaplan-Meier plotter(http://kmplot.com/analysis/)ESCC数据库做生存预后分析。

1.6 统计方法P<0.05被认为具有显著的统计学差异，多检验假设采用false discovery rate (FDR)方法对P值进行调整，且FDR<0.05作为筛选差异表达基因的阈值，在进行富集分析时运用超几何分布方法，同时本研究使用ESCC Kaplan-Meier曲线并用log-rank检查进行了比较。

2 结果

2.1 ESCC组织样本和正常细胞组织样本数据清洗本研究运用R语言相关函数从GEO数据库获得GSE38129数据(30个正常样本组织，30个ESCC样本组织)，共得到22 277个微表达矩阵芯片数据，通过对数据进行初步清洗，得到12 403个基因表达数据矩阵。

2.2 表达基因矩阵选定对所得12 403个基因表达矩阵通过运用R语言下limma包进行差异表达分析，得到12 403个差异基因表达矩阵(上调基因或者下调基因)。以“FDR<0.05及logFC≥2或log2FC≤-2”为阈值，共得到132个显著差异表达基因，其中51个上调基因、81个下调基因(图1)，同时将差异表达基因进行可视化做图(图2)。

红色部分代表51个上调基因，蓝色部分代表81个下调基因。图1 12 403个差异表达的火山图

图2 差异表达基因的热图

2.3 富集分析差异表达基因通过clusterProfiler包对筛选得到的差异表达基因进行富集分析。首先分别对差异表达基因进行分子功能、细胞组成、生物过程富集分析，发现这些基因表达产物分别表现肽链内切酶丝氨酸型肽酶活性、细胞外基质结构成分、丝氨酸型肽酶活性、丝氨酸水解酶活性等分子功能上(图3、表1)，细胞外基质、纤维胶原三聚体、编码角质化包膜、纺锤体等细胞组成中(图4、表2)，显著富集在表皮生长、细胞外基质组织、胶原蛋白分解代谢过程、女性减数分裂核分裂等生物过程内(图5、表3)。其次，对差异基因进行KEGG通路富集分析，发现这些基因显著富集在蛋白质消化吸收通路、Amoebiasis通路、类风湿性关节炎通路上(表4)。

图3 正常组织和ESCC组织的差异表达基因分子功能富集分析可视化热图

表1 正常组织和ESCC组织的差异表达基因(FDR<0.05及logFC≥2或log2FC≤-2)分子功能富集分析结果

注：FDR：false discovery rate。

图4 正常组织和ESCC组织的差异表达基因细胞组成富集分析可视化网络图

GO Path编号细胞组成PFDR基因数GO:0031012细胞外基质1.50E-133.00E-1123GO:0005583纤维胶原三聚体6.75E-094.50E-075GO:0098643带状胶原原纤维6.75E-094.50E-075GO:0062023含胶原细胞外基质1.84E-089.20E-0716GO:0044420细胞外基质成分3.35E-081.34E-067GO:0098644胶原三聚体复合物1.20E-074.02E-065GO:0016324顶端质膜1.26E-063.60E-0512GO:0045177细胞顶端部分1.63E-064.08E-0513GO:0001533编码角质化包膜1.32E-052.93E-045GO:0030496中间体2.64E-055.27E-048GO:0005581胶原蛋白三聚体3.02E-055.49E-046GO:0005819纺锤体4.67E-047.78E-039

注：FDR：false discovery rate。

表3 正常组织和ESCC组织的差异表达基因(FDR<0.05及logFC≥2或log2FC≤-2)生物过程富集分析结果

注：FDR：false discovery rate。

图5 正常组织和ESCC组织的差异表达基因生物过程富集分析可视化柱状图

编 号通路名称PFDR基因数hsa04974蛋白质消化吸收7.40E-068.88E-047hsa05323类风湿性关节炎7.12E-044.27E-025hsa05146Amoebiasis1.08E-034.33E-025

注：FDR：false discovery rate。

2.4 差异表达基因蛋白—蛋白相互作用网络构建和分析本研究将得到的132个差异表达基因上传到STRING数据库进行差异基因蛋白—蛋白相互作用网络构建(图7A)。为了进一步识别关键基因，本研究将STRING数据库得到的网络图数据上传至Cytoscape软件。将combined score>0.4、MCODE>2设定为筛选阈值，得到4个功能模块网络图，选取其中MCODE Score>5模块，得到2个功能模块网络图(图7B、C)，根据连接度和MCODE评分大小，得到DLGAP5、MMP1、TOP2A、MMP13、CXCL12、SPP1、ATAD2、PBK、VCAN共9个关键基因。

A：131个差异基因蛋白相互作用网络图；B：由15个节点、105个边组成的蛋白相互作用网络图；C：由10个节点、45个边组成的蛋白相互作用网络图。图7 差异表达基因蛋白—蛋白相互作用图

2.5 生存预后分析所筛选出来的差异表达基因均以TCGA数据库中ESCC生存分析数据为参考数据，通过Kaplan-Meier plotter进行生存预后分析，根据连接度和MCODE评分大小，得到DLGAP5、MMP1、TOP2A、MMP13、CXCL12、SPP1、ATAD2、PBK、VCAN共9个关键基因与生存预后关系(图8)。

3 讨论

ESCC恶性程度较高，预后较差，目前ESCC发生、发展的确切分子作用机制及信号通路尚不清楚[4]。本研究利用ESCC和正常组织微表达矩阵芯片的表达数据，进行ESCC中的特异性表达基因的生物学分析。

在本研究中，通过以“FDR<0.05及logFC≥2或log2FC≤-2”为选定阈值，筛选出来51个上调基因和81个下调基因，并把这些基因进行生物信息学分析。首先，对这些基因进行富集分析、分子功能富集分析表明在细胞外基质结构成分、丝氨酸型肽酶活性等处显著富集(表1)，VCAN、MMP1、CXCL12、ATAD2等差异基因参与了这些分子功能。ATAD2基因可能是核受体ESR 1的转录辅助激活因子，可诱导雌激素靶基因的表达，如CCND1、MYC和E2F1，可参与雌激素诱导的乳腺癌细胞增殖和细胞周期进程，可能在某些ESR1靶基因启动子中参与CREBBP的合成或占据[5]，可能导致ESCC的发生、发展。细胞组成富集分析显示在细胞外基质、纺锤体、染色体分离等处显著富集(表2)，VCAN、TPO2A、MMP1、SSP1等差异基因参与了这些细胞组成，SPP1基因编码的蛋白是一种细胞因子，可上调干扰素-γ和白细胞介素-12的表达。该基因对细胞-基质相互作用很重要，且作为细胞因子参与促进干扰素-γ和白细胞介素-12的产生和减少、白细胞介素-10的产生，在导致Ⅰ型免疫的途径中扮演重要角色[6-7]。生物过程富集分析表明在细胞外基质组织、胶原蛋白分解代谢过程等处显著富集(表3)，DLGAP5、VCAN、TPO2A、MMP1、MMP13、SSP1等差异基因参与了这些生物过程，DLGAP 5是一个蛋白质编码基因。该基因在泛素-蛋白酶体通路调控有丝分裂磷蛋白，是细胞黏附、连接、完整性和分化的关键调节因子，可能是在细胞癌变过程中起作用的潜在细胞周期调节因子[8-9]，从而可能导致ESCC细胞过度增殖分化。MMP1基因是11号染色体上MMP基因簇的一部分，其编码基质金属蛋白酶(MMPs)M10家族的一个成员。在正常的生理过程中，如胚胎发育、繁殖和组织重塑，以及关节炎和转移等疾病过程中，该家族的蛋白质参与细胞外基质的破坏。先前研究表明MMP1可能与膀胱癌、口腔鳞状细胞癌等恶性肿瘤相关[10]，也可能参与ESCC发生过程。

同时，把这些基因进行KEGG通路分析，得到蛋白质消化吸收通路、Amoebiasis通路、类风湿性关节炎通路3条通路，其中COL11A1、COL1A1、COL1A2、COL5A2等COL基因家族参与了蛋白质消化吸收通路(表4)，COL1A1是一种蛋白质编码基因，该基因编码Ⅰ型胶原蛋白的亲α1链，其三螺旋结构由两条α1链和一条α2链组成。该基因和血小板衍生生长因子β基因所在的第17和22号染色体之间的相互易位与一种称为突起皮肤纤维肉瘤的特殊类型的皮肤肿瘤有关，这种肿瘤是由于生长因子的不受管制的表达而引起的，其可能是ESCC发生的潜在基因[11-12]。

A：DLGAP5，P=0.000 19；B：MMP1，P=0.12；C：TOP2A，P=0.018；D：MMP13，P=0.12；E：CXCL12，P=0.22；F：SPP1，P=0.089；G：ATAD2，P=0.001 9；H：PBK，P=0.000 62；I：VCAN，P=0.003 7。红色代表基因高表达，黑色代表基因低表达。图8 所有差异表达基因在以TCGA数据库中81例ESCC患者做Kaplan-Meier总生存分析

其次，把这些差异表达基因上传至STRING数据库和 Cytoscape软件(图7A、B、C)，根据连接度和MOCDE评分大小，得到DLGAP5、MMP1、TOP2A、MMP13、CXCL12、SPP1、ATAD2、PBK、VCAN共9个关键基因。最后，9个关键基因在Kaplan-Meier plotter数据库做生存分析(图8)。PBK、VCAN基因与食管癌患者生存期明显相关(P<0.01)(图8H、I)，到目前为止，没有文章或者报道出PBK、VCAN和ESCC发生、发展相关，PBK基因编码丝氨酸/苏氨酸蛋白激酶，与双特异性丝裂原活化蛋白激酶(MAPKK)家族相关。这一基因的过度表达与肿瘤的发生有关，有关研究已证实PBK参与了血液学肿瘤发生[13-14]，可能参与了ESCC发生。VCAN基因编码的蛋白是一种大型硫酸软骨素蛋白多糖，是细胞外基质的主要组成部分。该基因编码的蛋白参与细胞黏附、增殖、迁移和血管生成，在组织形态发生和维持中起核心作用。此基因可能在细胞间信号传递和连接细胞与细胞外基质中发挥作用，可能参与细胞运动、生长和分化的调控[15-16]，可能对ESCC的发生起到了调控作用。

除此之外，DLGAP5、TPO2A、ATAD2也与ESCC患者生存期显著相关(P<0.01)(图8A、C、G)，既往研究表明，DLGAP5与前列腺癌、结直肠癌生存预后相关[17-18]，可能参与了恶性胶质母细胞瘤、肝细胞癌等恶性肿瘤发生、发展过程[19-20]。ATAD2可以调控众多经典的肿瘤相关信号通路，其中包括诱导MYC、E2F和甲状腺激素受体和视黄醛激活因子，ATAD2表达与肿瘤的组织学分级以及肿瘤转移、疾病复发和患者生存期等相关，可作为肿瘤预后评价分子，具有重要的临床应用价值，且在前列腺癌、乳腺癌、肺腺癌等多种恶性肿瘤的发生和发展过程中发挥重要作用[21]。TOP2A基因编码一种DNA拓扑异构酶，这种核酶参与了染色体的浓缩、染色单体的分离以及DNA转录和复制过程中扭转应激的缓解。它催化了两条双链DNA的瞬间断裂和重新连接，使两股DNA相互通过，从而改变了DNA的拓扑结构[22-23]。TOP2A参与了结肠癌、恶性胶质母细胞瘤、胆管癌、乳腺癌等多种恶性肿瘤发生过程[24-25]。MMP1、MMP13、CXCL12、SPP1基因高表达与ESCC患者生存期呈负相关(P>0.05)(图8B、D、E、F)，无统计学意义，可能与本研究所选ESCC患者样本数量较少有关。

总之，本研究通过ESCC与正常组织的基因表达谱综合生物信息学分析，表明PBK、VCAN、DLGAP5、ADAT2、TOP2A可能是ESCC发生、发展和预后的关键基因。虽然数据挖掘和整合是预测恶性肿瘤生物标志物的一种很有前途的工具，但是本研究只是一个初步的结论，且本研究中的样本数量是有限的。由于肿瘤生物标志物只有与临床数据相结合才有意义，因此需要进行进一步的实验来证实本研究的结论。