基于生物信息学的食管鳞癌相关基因筛选及其预后生存分析

杨 莉，杨东红，冯丽萍，谢 忠，李海文（广东医科大学附属医院肿瘤中心，广东湛江 524000）

中国是食管癌的高发地区，超过90%为鳞癌[1]，近年来治疗手段虽取得进展，但未显著提高疗效[2]，5 a生存率仅为29.7%[3]。由于肿瘤的异质性，疗效差异显著[4]，因此，筛选预测食管癌生存预后的标志物，对选择治疗策略和改善生存预后意义重大[5]。本研究使用生物信息分析筛选与食管鳞癌（ESCC）生存预后相关基因，旨为其提供预测生存预后的标志物。

1 资料和方法

1.1 数据信息的获取

基因芯片表达谱数据集GSE17351、GSE26886、GSE 77861、GSE161533、GSE20347、GSE23400 及GSE75241数据集均来自于Gene Expression Omnibus（GEO）数据库（https://www.ncbi.nlm.nih.gov/geo/），其中GSE17351数据集中包括5 例食管癌旁组织和5 例食管鳞癌组织的样本，GSE26886 数据集中包括19 例食管癌旁组织和9 例食管鳞癌组织的样本，GSE77861 数据集中包括7 例食管癌旁组织和7 例食管鳞癌组织的样本，GSE161533 数据集中包括28 例食管癌旁组织和28 例食管鳞癌组织的样本，GSE20347 数据集中包括17 例食管癌旁组织和17 例食管鳞癌组织的样本，GSE23400 数据集中包括53 例食管癌旁组织和53 例食管鳞癌组织的样本，GSE75241 数据集中包括15 例食管癌旁组织和15 例食管鳞癌组织的样本。基因芯片表达谱及患者的临床信息特征从TCGA（http://tcga‐data.nci.nih.gov）数据库下载，其中包括1 例食管正常组织患者和80 例食管鳞癌患者的相关数据。

1.2 癌旁组织中差异表达基因（DEGs）的筛选

使用R 软件中limma 包对食管癌旁组织与食管鳞癌组织之间的DEGs 进行筛选，选条件为变化倍数|log2（FC）|>1 且校正后的P<0.05。单基因芯片筛选出的DEGs 使用R 软件中ggplot2 软件包进行绘制火山图，多基因芯片联合筛选DEGs 使用R 软件中Robust RankAggreg 及pheatmap 软件包进行绘制热图。

1.3 DEGs 及核心基因基因本体（GO）功能及KEGG信号通路分析

采用R 软件中clusterProfiler、org.Hs.eg.db、DOSE、enrichplot、colorspace 等程序包对DEGs 和核心基因进行GO 功能富集分析和KEGG 代谢通路富集分析[6]。GO类别富集分析包括基因参与的生物学过程（biological process，BP）、基因所处的细胞组分（cellular component，CC）和基因执行的分子功能（molecular function，MF）。在富集分析过程中，使用费舍尔精确检验来检验差异表达基因在某个网络中是否富集，设置校正后的P<0.05作为筛选条件，使用R 软件中ggplot2 软件包进行绘制柱状图。

1.4 DEGs 蛋白互作网络分析和筛选核心基因

DEGs 提交至STRING11.0 数据库（https://string‐db.org/）进行蛋白互作网络（PPI）网络分析，并用Cytoscape软件对蛋白质互作网络进行可视化，设置最低互作分数为 0.7。用Cytoscape3.6.0 生物信息学软件平台的图论聚类算法MCODE（molecular complex detection）工具对构建的生物学网络进行关联度分析，筛选出核心蛋白质簇和关键节点蛋白，MCODE 选择标准：cut‐off值为评分>10 及节点数>20。使用该软件中CytoHubba 插件（含10 种算法）计算位于共表达网络核心的核心关键基因，使用R 软件中UpSetR 软件包进行统计10 中算法中均出现且频率最高的核心基因进行绘制网络结构中的核心差异基因（度值最高的基因）柱状图。

1.5 核心基因的生存分析

从TCGA 数据库下载的ESCC 相关数据，并进行转录组数据整理，提取患者生存数据，并根据核心基因表达的中位数分高表达组和低表达组，分析比较两组间的生存时间进行判断目的基因与食管鳞癌患者预后的关系，使用R 软件中survival 和suvminer 软件包进行绘制核心基因的生存曲线。

1.6 统计学处理

所有统计数据均采用R3.6.3 对数据进行统计学分析，R 软件包从Bioconductor 网站下载（https:// www.bioconductor.org）。采用Wilcox 检验进行统计筛选DEGs，两组间比较采用t检验进行分析；费舍尔精确检验关联性；采用卡方检验比较分类变量；生存分析采用Kaplan‐Meier 法log‐rank 检验。检验水准α=0.05（双尾）。P<0.05 为差异有统计学意义。

2 结果

2.1 DEGs 筛选结果

7 个基因芯片表达谱共包括144 个食管癌旁组织及134 个食管鳞癌组织样本，每个基因芯片表达谱数据集及7 个基因芯片表达谱数据集联合分析DEGs，结果每个基因芯片表达谱数据集DEGs 如图1A 所示，对7 个基因芯片DEGs 进行映射，发现775 个共差异表达基因，其中上调表达基因325 个、下调表达基因450个，上调与下调差异表达基因中变化最显著的前20 个基因如图1B 所示。

图1 DEGs 筛选结果

2.2 GO 和KEGG 通路富集分析

把775 个DEGs 进行GO 类别分析，筛选P<0.05 的结果有377 条，其中BP、CC、MF 每个类别最显著的前5 条GO 分析结果如图2 所示。生物过程最具代表性的是：细胞外结构组织、细胞外基质组织、皮肤发育、表皮发育、角质化；最具有显著意义的细胞组分：含胶原蛋白的细胞外基质、细胞外基质成分、角质包膜、胶原三聚体复合物、细胞顶端；最具有代表性的分子功能：细胞外基质结构成分、丝氨酸型内肽酶活性、丝氨酸型肽酶活性、丝氨酸水解酶活性、具有拉伸强度的细胞外基质结构成分。KEGG 通路富集分析结果筛选出13 条信号通路有显著性差异（P<0.05）。其中最显著前5 条信号通路分别是：IL‐17 信号通路、细胞外基质受体相互作用、病毒蛋白与细胞因子及细胞因子受体的相互作用、阿米巴病、细胞周期信号通路。结果如图2。

图2 GO 功能富集分析及KEEG 分析信号通路最显著的前5 个结果

2.3 PPI 的构建和筛选核心基因

采用STRING 数据库分析构建775 个DEGs 的PPI，共有452 个DEGs 被过滤到PPI 网络中，其中包含452个节点和2 226 个交互作用，有 94 个属于核心基因，分别属于5 个PPI 子网络（得分≥10）。PPI 子网络A（得分35.4）包含41 个节点、728 个交互作用；PPI 子网络B（得分14.7）包含15 个节点、103 个交互作用；PPI子网络C（得分14.0）包含14 个节点、91 个交互作用；PPI 子网络D（得分13.0）包含13 个节点、78 个交互作用；PPI 子网络E（得分10.0）包含11 个节点、50 个交互作用。见图3A‐E。使用R 软件中UpSetR 软件包统计10 种算法中出现频率最高的核心基因并进行绘制核心基因集分布柱状图，结果10 种算法均出现的核心基因数目有12 个，结果如图3F 所示，12 个核心基因分别是：泛素结合酶E2C（ubiquitin‐conjugating enzyme E2C，UBE2C）、微小染色体维持蛋白7（mini‐chromosome maintenance protein 7，MCM7）、III 型胶原α1 链（Collagen Type III Alpha 1 Chain，COL3A1）、VII 型胶原α1 链（Collagen Type VII Alpha 1 Chain，COL7A1）、V 型胶原α2 链（Collagen Type V Alpha 2 Chain，COL5A2）、IV 型胶原α1 链（Collagen Type IV Alpha 1 Chain，COL4A1）、C‐X‐C 基序趋化因子配体10（C‐X‐C Motif Chemokine Ligand 10，CXCL10）、膜联蛋白A1（Annexin A1，ANXA1）、神经调节肽受体（Neuromedin U，NMU）、ISG15 泛素样修饰物（ISG15 Ubiquitin Like Modifier，ISG15）、重组人光蛋白聚糖（Lumican，LUM）、细胞周期蛋白A1（Cyclin A1，CCNA1）。

图3 PPI 的构建和核心基因的筛选

2.4 12 个核心差异基因进行GO 功能分析和KEGG通路富集分析

GO 功能富集BP 主要包括细胞外基质组织、细胞外结构组织、胶原原纤维组织等；CC 主要包括：胶原三聚体复合物、细胞外基质成分、胶原三聚体；MF 主要包括：具有拉伸强度的细胞外基质结构成分、细胞外基质结构成分、血小板衍生生长因子结合；其中，ANXA1/MCM7/CCNA1 主要BP 功能富集于有丝分裂细胞周期的G1/S 转换和细胞周期G1/S 相变；ANXA1/MCM7 主要BP 功能富集于DNA 双链展开、DNA 几何变化和DNA 构象变化；KEGG 信号通路富集主要包括：蛋白质消化吸收、EBV 病毒感染、RIG‐I 样受体信号通路和细胞周期信号通路（MCM7/CCNA1）等。其中BP、CC、MF 每类别最显著的前3 条GO 注释术语及最显著前5 条信号通路结果如图4。

图4 12 个核心差异基因的 GO 功能富集分析最显著的前3个结果和KEEG 信号通路富集分析最显著的前5 个信号通路

2.5 核心基因的生存分析

用TCGA 中的80 例ESCC 病例临床数据来验证核心基因的表达情况与患者生存的关系。将筛选出与食管癌相关的12 个核心基因进行验证分析，根据每个核心基因的表达中位数划分成低表达组和高表达组，结果发现MCM7、CCNA1、ANXA1 共3 个基因与食管癌病人生存预后显著相关，其中MCM7 高表达组患者生存率比低表达组患者高（高、低表达组的临界点是38.368），差异有统计学意义（P=0.047）；CCNA1 和ANXA1 高表达组患者生存率比低表达组患者低（高、低表达组的临界点分别是3.678 和325.510），差异有统计学意义（分别P=0.025、P=0.035）。见图5。

图5 Kaplan‐Meier 预测ESCC 患者生存率

3 讨论

食管癌是我国常见的恶性肿瘤之一，目前全球食管癌5 a 生存率仅为10%～30%[3]，因此挖掘应用于食管癌预测生存预后的分子标志物显得尤为重要。食管癌的发生是一个涉及多基因变化的复杂过程，采用基因芯片从分子水平研究肿瘤的发生、发展过程，寻找与食管癌预后密切相关的基因备受关注。本研究整合了GEO 数据库中ESCC 的7 个数据集，利用生物信息学方法对这些数据集进行深入分析，挖掘出12 个核心基因，其中MCM7、CCNA1、ANXA1 与食管鳞癌患者生存预后密切相关。

本研究对12 个核心基因进行GO 功能富集分析发现ANXA1/MCM7 主要富集于DNA 双链展开、DNA几何变化和DNA 构象变化，调控DNA 复制的进程等生物过程。DNA 复制是细胞增殖的核心过程，而异常DNA 复制被认为是肿瘤发生的驱动力，MCM7 在启动DNA 复制中起重要作用，MCM 复合体可通过保护复制叉结构参与检查点控制，而且可以通过募集捕获后的检查点蛋白以协助重启复制，促进肿瘤的发生、发展。Qiu 等[7]发现大约12%的ESCC 以及超过4%的头颈部鳞癌和胃癌中MCM7 有异常扩增，siRNA 敲除MCM7 可显著抑制了体外KYSE510 和EC9706 食管癌细胞系的细胞增殖、集落形成和迁移。Qu 等[8]发现，将肝细胞癌细胞敲除MCM7 后可显著抑制体外细胞增殖和降低在动物体内移植瘤的致瘤能力，同时发现MCM7 高表达的肝细胞癌患者的总体生存率较差，MCM7 高表达可作为预测肝细胞癌患者生存预后的标志物，且具有较高的敏感性和特异性。本研究对12 个核心基因进行KEGG 富集分析亦发现MCM7/CCNA1富集于细胞周期信号通路，且ESCC 患者MCM7 表达水平与患者的生存预后密切相关，提示MCM7 可潜在作为预测ESCC 患者的生存预后标志物。类似的，张琼等[9]发现MCM7、肿瘤蛋白P53 蛋白联合检测可作为食管癌早期诊断和预测生存预后的标志物。以上相关研究结果提示，MCM7 表达水平与食管癌预后密切相关。然而，王利等[10]却发现MCMs 家族成员中仅MCM4、MCM6 与食管癌患者生存预后具有相关性。综上所述，MCM7 是否能作为预测ESCC 患者生存预后因素尚有争议。本研究为筛查预测ESCC 患者生存预后标志物，从TCGA 数据库下载ESCC 患者的相关临床数据分析发现MCM7 高表达组患者生存预后好于较低表达组患者（P=0.047），本研究结果与文献报道结果MCM7 高表达是肿瘤患者预后不良因素结论不一致，一方面考虑可能是P=0.047 接近统计有显著性差异临界值0.05，未能有效检验MCM7 表达水平与ESCC 患者预后的显著相关性；另一方面，由于患者例数仅有80 例，病例数有限，需纳入更多的病例数来提高结果的可靠性和进一步进行实验验证。下一步我们将收集相关临床样本，研究MCM7 表达水平与ESCC 患者预后的相关性。

CCNA1 的异常表达与多种实体瘤的癌变、进展和转移密切相关。本研究发现CCNA1 高表达组患者生存率比低表达组患者低（P=0.025），提示CCNA1 表达水平与ESCC 患者预后呈负相关性。类似的，He 等[11]发现ESCC 中CCNA1 mRNA 水平比邻近非癌组织中增加了10 倍以上，约 70%的ESCC 患者CCNA1 表达水平上调，CCNA1 mRNA 过表达的患者往往临床分期更晚期、淋巴结转移率增加、侵袭性增加和生存预后（无病生存率和总生存率）比低表达组差，研究结果提示CCNA1 在ESCC 的发展中具有重要作用，并且可能作为ESCC 的预后生物标志物和治疗靶点。Chang 等[12]研究发现，在530 例结直肠肿瘤患者中CCNA1 基因拷贝数显著增加，且CCNA1 过表达的患者无进展生存期显著缩短；当使用奥沙利铂化疗时，CCNA1 过表达的近端结直肠肿瘤患者的总体生存率和无进展生存期比未进行化疗的结直肠癌更差，研究结果表明CCNA1 过表达可能有助于优化近端结直肠癌的治疗和监测，并有可能开发成分子检测标志物。本研究对核心基因进行KEGG 信号通路富集发现MCM7/CCNA1 主要涉及细胞周期信号通路，类似的，Zhang 等[13]发现lncRNA LINC00304 在前列腺癌细胞和组织中的表达高于正常前列腺细胞和组织，其过表达可以显著促进前列腺癌细胞中CCNA1 的过表达，并且显著促进前列腺癌细胞的细胞增殖和细胞周期进程。Chen 等[14]发现上调miR‐1271 和下调CCNA1 可抑制肝细胞癌增殖、迁移和侵袭，同时通过激活AMPK 信号通路加速细胞凋亡，从而抑制肝细胞癌的进展。有研究报道CCNA1 基因不仅可作为甲状腺乳头状癌的潜在诊断标志物[15]，而且也可作为预测卵巢癌对紫杉醇、多柔比星和5‐氟尿嘧啶的耐药性的生物标志物[16]，化疗药物耐药是肿瘤治疗失败的主要原因，导致患者预后差。本研究结果发现CCNA1 高表达ESCC 患者生存率降低，与相关文献报道结果一致，提示CCNA1 的表达水平可能可作为预测ESCC 患者生存预后的标志物。由于本研究从TCGA下载的ESCC 患者临床数据主要来自非中国患者的临床数据且患者样本例数仅80 例，尚需在国内收集我国患者的临床数据，并进一步扩大病例数进行验证。

本研究发现核心基因ANXA1 高表达组患者生存率比低表达组患者低（P=0.035），提示ANXA1 表达水平与ESCC 患者生存预后呈负相关。ANXA1 是一种内源性免疫调节蛋白，是一种调节转录因子和miRNA的稳态蛋白，与肿瘤的多种特性有关，包括细胞增殖、凋亡、化学敏感性、转移和侵袭[17]。而吴晓鹏等[18]发现miR‐196a 能够作为食管癌预后的预测指标之一，其作用机制可能是miR‐196a 通过抑制ANXA1 等基因的表达参与调控食管癌细胞的增殖、侵袭和迁移等。类似的，Han 等[19]将构建ANXA1 过表达的质粒转染到ESCC 细胞系EC109 中，发现转染ANXA1 过表达质粒的细胞中ANXA1 表达水平上调，细胞增殖、迁移和侵袭性显著增加。Hu 等[20]研究发现miR‐196a 可能通过靶向ANXA1 调控ESCC 细胞的增殖、侵袭和迁移，miR‐196a/ANXA1 轴可能代表 ESCC 的潜在治疗靶点。Takaoka 等[21]发现 ANXA1 和COX‐2 在食管癌和胃癌的癌前病变中不受调控的表达，从而促进胃癌、食管癌发生。也有研究报道ANXA1 可作为ESCC 肿瘤三级预防策略中的早期诊断ESCC 的潜在工具[22]。Han 等[23]研究发现ESCC 患者血清ANXA1 水平显著高于健康人群，若放化疗后ANXA1 水平升高，尤其是较基线水平升高2 倍以上的无进展生存期显著降低，结果表明放化疗后ESCC 患者血清ANXA1 的表达水平增加是患者预后不良的因素。综上所述，ANXA1 表达水平升高不仅可促进肿瘤细胞的增殖、侵袭和迁移，而且还可作为肿瘤潜在治疗的靶点和潜在早期诊断标志物，甚至可预测肿瘤患者的生存预后。本研究亦发现ANXA1高表达ESCC 患者生存率降低，提示ANXA1 表达水平高低有可能可作为ESCC 的预后生物标志物及早期诊断的潜在工具，与既往相关研究结果一致，但仍需更进一步收集临床ESCC 患者的标本及相关临床数据进行验证，为其在临床应用提供更可靠的证据。

此外，本文也存在如下不足之处。首先，该研究从GEO 选用的基因芯片表达谱数据集并非是同一个平台数据，不同平台间数据可能有差异；其次，由于是基于生物信息学方法进行的研究，未对筛选出的基因进行进一步实验验证，为之后的临床应用提供更可靠的依据；最后，由于ESCC 疾病本身发生是受环境等多方面因素影响、发生机制的异质性和TCGA 临床数据样本数量的有限性及缺乏亚洲来源的样本进行验证，未来仍需纳入更多的外部数据如我国食管癌患者数据提高结果的可靠性和进一步行实验验证。

综上所述，本研究采用基因芯片数据与生物信息学技术相结合，挖掘出与ESCC 生存预后关系紧密相关的核心基因MCM7、CCNA1、ANXA1，发现MCM7低表达及CCNA1、ANXA1 高表达是ESCC 生存预后不良因素，提示其可能可作为预测ESCC 患者生存预后的标志物，但需要进一步利用我国ESCC 患者数据及扩大样本进行实验验证，为后续其作为ESCC 的治疗靶点和预测生存预后的标志物提供依据。