基于生物信息学对食管鳞状细胞癌新辅助放化疗敏感性生物标志物的发掘*

杜芝霖，伍健，徐鹏，郎锦义

646000 四川 泸州，西南医科大学附属医院 肿瘤科（杜芝霖、伍健、郎锦义）；610041 成都，四川省肿瘤临床医学研究中心，四川省肿瘤医院·研究所，四川省癌症防治中心，电子科技大学附属肿瘤医院放疗中心/放射肿瘤学四川省重点实验室（徐鹏、郎锦义）

根据全球癌症统计，2020 年有超过604 000 例食管癌新发病例和544 000 例死亡病例，全球发病率第七，死亡率第六［1］。中国是食管癌的高发地区，2015 年中国癌症统计，食管癌的发病率排名第三，死亡率排名第四［2］。其中，食管鳞状细胞癌（esophageal squamous cell carcinoma，ESCC）是食管癌中最常见的组织学类型。

ESCC 的早期症状并不明显，许多患者直到癌症晚期才被明确诊断。尽管近几十年出现了新的治疗策略和方案，但食管癌的预后仍然很差，5 年生存率约为20%［3-4］。而如今，新辅助放化疗（neoadjuvant chemoradiotherapy，NCRT）联合根治性食管切除术已逐渐成为食管癌的主要治疗方法［5］。然而，不同患者接受NCRT 的结果不一样。获得病理完全缓解（pathologic complete response，PCR）的患者比未获得PCR 的患者有更好的生存率［6-7］。但一些临床参数例如TNM 分期、肿瘤位置和肿瘤大小等并不能预测食管鳞癌患者对NCRT 是否敏感［8］。如何有效地预测食管鳞癌患者对NCRT 治疗的敏感性是本研究的目的。

随着高通量测序技术和生物信息学的发展，对肿瘤发生机制的认识有了很大的提升［9］。有相关研究提出，一些基因水平的变化可以预测NCRT 敏感性，这有助于ESCC 患者的个体化治疗［10-11］。如SOX17基因在食管鳞癌患者中的过表达对NCRT更为敏感［12］。研究人员也通过建立35 个基因突变谱的预测模型，来预测食管癌患者接受NCRT 治疗的敏感性［13］。还有研究发现了能预测食管癌患者NCRT 治疗敏感性的5 个基因［14］。因此基因的异常表达或突变会影响到患者接受NCRT 的敏感性。在本研究中，通过使用生物信息学方法，筛选出与食管鳞癌NCRT 治疗敏感性相关的具有预测价值的生物标志物，并采用免疫组化方法进行初步验证，以期获得新的治疗靶点，为食管癌患者的个体化治疗带来希望。

1 材料和方法

1.1 数据源和数据处理

从GEO 数据库中筛选与食管鳞癌NCRT 敏感性相关的数据集，纳入标准：①提供食管癌患者RNA-Seq 的数据集；②包含食管癌NCRT 后不敏感患者的组织与敏感患者的组织。最终筛选出数据集GSE45670［15］。

回顾性收集了四川省肿瘤医院2017～2021 年食管鳞癌患者的相关信息。纳入标准：① 患者在疾病确诊前未进行过任何抗肿瘤治疗；② 患者年龄18～75 岁，KPS 评分 ＞ 80；③ 我院初治的中晚期胸段食管鳞癌并且无远处转移的患者 （cT1b-cT2N+或cT3-cT4任何N），并于我院行NCRT 及手术治疗；④ 所有患者治疗前均需完善相关基线检查，明确病灶情况。所有患者新辅助同步放化疗治疗后第28天复查增强CT 进行初次筛选NCRT 敏感性。影像学提示肿瘤病灶没有消退，称之为疾病稳定（stable disease，SD）或是肿瘤病灶较前增大，称之为疾病进展（progression disease，PD）的患者归类于放化疗不敏感；影像学提示肿瘤病灶完全消失，称之为完全缓解（complete response，CR）或者肿瘤病灶较治疗前缓解 ＞ 50%的患者归类于放化疗敏感。放化疗后约6 周，患者接受食管鳞癌根治性手术切除原发肿瘤以及区域淋巴结。术后病理是否获得PCR 作为最终划分食管鳞癌患者接受NCRT 敏感性的标准。

1.2 差异表达基因（differentially expressed genes，DEGs）的获得及蛋白互作（protein-protein interaction，PPI）分析

使用R 软件的“limma”包对GSE45670 数据集的NCRT 敏感组和NCRT 不敏感组的 DEGs 进行筛选，阈值设置为P＜ 0.05 和LogFC ＞ 1，并使用“ggplot2”包进行了热图和火山图的绘制［16-17］。

将筛选出的DEGs 导入String 数据库（https://string-db.org）进行PPI 网络分析。通过设置PPI 最低互作得分（0.4），去除该网络中不相关的节点，得到初次筛选的PPI 网络［18］。将其导入Cytoscape 软件中将PPI 网络可视化。通过Cytoscape 软件的MCODE 插件进一步过滤，设置MCODE 分数 ＞ 4，节点数 ＞ 5［19］。使用该软件中Cytohubba 和CytoNCA 两个插件，运用Degree 算法，筛选出核心差异基因［20-21］。

1.3 构建加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA）

利用WGCNA 分析方法对数据集GSE45670构建共表达网络［22］。首先，筛选出方差前25%的5 837 个基因。其次，计算每对基因之间的皮尔逊相关性，得到一个相似性矩阵。第三，使用“pick Soft Threshold”函数来选择理想的软阈值（β = 5）。第四，利用幂函数将相似度矩阵转换为邻接矩阵。然后，将邻接矩阵转换为拓扑重叠矩阵（topoligical overlap matrix，TOM）。最后，建立了基于TOM 的不相似度（1-TOM）的层次聚类。每个模块树状图至少包含50 个基因，将相似表达模块聚类为相同模块的阈值设为0.25，最终得到模块。

1.4 筛选候选基因进行基因富集分析

使用R 软件的“VennDiagram”包，将DEGs 和WGCNA 中相关度最高的模块内的基因进行交叉重叠后，结合PPI 网络中获得核心基因和候选基因。

对候选基因进行基因通路富集分析，使用R 包“clusterProfiler”和“enrichplot”包对其进行注释，对可能参与的生物过程和信号通路进行相应的分析［23-24］。

1.5 对候选基因进行生存分析

从TCGA 数据库中筛选食管鳞癌患者行Kaplan-Meier 生存分析，获得食管鳞癌患者基因表达水平与生存时间之间的关系。在Log-rank 检验中，P＜ 0.05 为差异有统计学意义。

1.6 基质金属蛋白酶（matrix metalloproteinase，MMP）9 和醛固酮还原酶（aldo-keto reductase，AKR）1C3 的免疫组化分析

对已明确诊断为ESCC 患者的活检样本进行免疫组化分析。采用抗体MMP9（1∶100）和抗体AKR1C3（1∶200）进行苏木精染色。目标蛋白的表达量通过使用Image-Pro Plus 图像软件进行分析。MMP9 阳性蛋白主要在细胞质和细胞膜中表达为褐色。AKR1C3 阳性蛋白表现出中度到较强的细胞质反应，颜色为褐色。通过测量出每张图片的累积光密度值 （integrated option density， IOD）以及区域面积值（area），再计算出平均光密度值（mean density），即 mean density = IOD/area，此值反映了目标蛋白的单位面积浓度。最后取每个样本的 3 个随机区域mean density 的平均值即为此样本的值。

1.7 统计学方法

利用R 软件（4.2.1）、Cytoscape 软件（3.8.0）、ImageProPlus 软件、String 在线网站（https: //string-db.org/）进行数据分析。采用Kaplan-Meier 法进行生存分析，并采用Log-rank 检验。利用TCGA 数据集，来获取相关候选基因的表达水平的高低与食管鳞癌患者总生存期之间的关系。并且对其进行了相应曲线的绘制。免疫组化的结果采用GraphPad Prism 7.0 进行统计分析。采用Mann-WhitneyU检验或双尾非配对检验来确定显著性水平。本研究中使用的所有数据均满足统计检验的假设。所有定量数据均以均值±标准差呈现。P＜ 0.05 为差异有统计学意义。

2 结 果

2.1 数据收集

从GEO 数据库中筛选得到的数据集GSE45670，包含了28 个食管鳞癌患者的样本，其中有17 个食管鳞癌NCRT 治疗不敏感的患者和11 个食管鳞癌NCRT 治疗敏感的患者。而临床样本共收集了26 例活检组织标本，包括11 例NCRT 治疗敏感者，也是影像学上提示肿瘤病灶CR 的患者，和15 例NCRT治疗不敏感者，其中11 例患者影像学提示肿瘤病灶SD，4 例患者影像学提示肿瘤病灶PD。

2.2 DEGs 的筛选及PPI 网络分析

以LogFC ＞ 1，P＜ 0.05 为阈值，共筛选得到350个差异基因，其中上调基因166 个，下调基因184 个（图1）。

图1 差异基因筛选结果Figure 1. Volcano Plot of DEGs

运用String 在线数据库构建了350 个DEGs 的PPI 网络，再导入Cytoscape 软件中，运用MCODE 插件进行筛选，再使用Degree 算法，运用该软件中的CytoHubba 和CytoNCA 插件最终确定MMP9是该网络的核心基因（图2）。

图2 由差异基因构成的PPI 网络Figure 2. Protein-Protein Interaction Networks Composed of Differentially Expressed Genes

2.3 加权基因共表达网络的构建

WGCNA 是近年来广泛使用的一种技术，旨在避免仅仅分析DGEs 而造成忽视核心基因的可能性［25-26］。通过选取相关系数为0.85 时对应的软阈值为 5（图3A）构建无尺度网络进行 WGCNA，再通过动态杂化剪切合并相似的模块，获得13 个模块（图3B）。其中没有被分配到任何集群中的灰色模块，在本研究中没有进行进一步的分析。结果显示，紫色模块与食管鳞癌患者接受NCRT 敏感性的相关性最高，共159个基因，差异具有统计学意义（P= 0.03；图3C）。

图3 构建WGCNA 模块Figure 3. Construction of WGCNA Modules

2.4 获得候选基因

将紫色模块中的159 个基因与350 个DEGs 进行交集，提取到5 个重叠基因（图4）。将这5 个重叠基因与PPI 网络中的核心基因结合，得到6 个候选基因，包括AKR1C3、GPX2、RRBP1、GSTM3、AKR1C1和MMP9。

图4 DEGs 和紫色模块基因的韦恩图Figure 4．Screening of Candidate Genes: Venn Diagram of Intersections between DEGs and the Purple Module

2.5 候选基因的功能富集分析

这些候选基因共得到335 个相关的GO 条目。选出其中10 个GO 条目（图5A）。同时，分别从生物过程、细胞组分和分子功能这三个方面选取了相关的条目（图5B）。基因和基因组的KEGG 信号通路分析结果发现，这些基因主要富集于“氧化应激反应”中，也可能与“细胞解毒”“细胞间桥”等途径有关（图5C）。

图5 候选基因的功能富集分析Figure 5. Functional Enrichment Analysis of the Candidate Genes

2.6 候选基因的生存分析

从TCGA 的食管癌数据库中筛选出81 个食管鳞癌患者的临床数据来对这6 个候选基因进行生存分析，通过生成最佳截止点，形成基因高表达和低表达两组。从图6 可知，AKR1C3和AKR1C1基因高表达组患者与低表达组患者相比，总生存率（overall survival，OS）明显较差，差异具有统计学意义（P= 0.011，P= 0.035），而食管鳞癌患者的OS 并不受其它4 个基因的表达影响，差异不具有统计学意义（均P＞ 0.05）。

图6 TCGA 数据库中ESCC 中候选基因的Kaplan-Meier 生存曲线Figure 6. Kaplan-Meier Survival Curves of Candidate Genes in ESCC from TCGA Database

2.7 通过免疫组化检测食管鳞癌标本中AKR1C3和MMP9 的表达

对GSE45670 数据集中6 个基因的表达进一步分析，P＜ 0.01 被认为是基因表达的显著性差异，发现AKR1C3、AKR1C1、MMP9和RRBP1在NCRT 不敏感的组织和NCRT 敏感组织中的表达差异更为显著（图7）。通过对本研究收集的活检标本进行免疫组化分析，与NCRT 敏感的组织相比，NCRT 不敏感组织中AKR1C3 蛋白表达上调（图8）；MMP9 蛋白则在NCRT 敏感组织中表达上调（图9）。

图7 GSE45670 数据集中6 个候选基因的接受NCRT 食管鳞癌不敏感组织和接受NCRT 食管鳞癌敏感组织之间的基因表达差异Figure 7．Expressions of 6 Candidate Genes in Non-responder and Responder to NCRT in GSE45670 Dataset

图8 AKR1C3 在NCRT 敏感组织和 NCRT 不敏感组织表达水平的比较（DAB 法，×20）Figure 8．Expression Levels of AKR1C3 in Responder and Non-responder to NCRT (DAB, ×20)***P ＜ 0.001.

图9 MMP9 在 NCRT 敏感组织和 NCRT 不敏感组织表达水平的比较（DAB 法，×20）Figure 9．Expression Levels of MMP9 in Responder and Non-responder to NCRT（DAB，×20）

3 讨 论

本研究利用生物信息学对基因集GSE45670 进行分析。首先通过DEGs 分析筛选出 350 个差异基因。其次，基于WGCNA 生物信息学方法，得出与食管鳞癌放化疗治疗中异质性显著相关的159 个基因，与350 个DEGs 进行交叉，获得了5 个重叠基因（AKR1C1、AKR1C3、GPX2、GSTM3、RRBP1）。 同时，基于350 个DEGs，生成了相应的PPIs 网络。利用MCODE、CytoHubba 和CytoNCA 这三个插件，将其导入Cytoscape 软件进行可视化，得到了1 个核心基因（MMP9）。最后，得到了6 个候选基因，包括AKR1C1、AKR1C3、GPX2、GSTM3、RRBP1和MMP9。通过进一步对这6 个候选基因进行功能富集分析发现，“氧化应激反应”途径在其中起着非常重要的作用。放射治疗通过多种方式诱发肿瘤细胞的死亡，例如产生活性氧（reactive oxygen species，ROS）、DNA损伤和亚细胞细胞器中的应激反应。ROS 是由电离辐射所诱发的，不仅与蛋白质、DNA 损伤、细胞膜表面酶和基因组不稳定相关，还引发损伤修复，可以破坏氧化还原稳态，从而导致氧化应激反应，最终导致细胞死亡［27-29］。并且化疗会诱导产生大量的ROS，从而导致癌细胞死亡。该通路对于癌症的发展、治疗的耐药性和基因组的不稳定性等方面都有着不可忽视的作用。

同时，通过设置P＜ 0.01 为阈值，得出AKR1C3/1、MMP9和RRBP1在GSE45670 数据集中NCRT 敏感组和NCRT 不敏感组的表达差异更为显著。

RRBP1 是一种内质网膜蛋白，在核糖体结合、新生蛋白易位以及内质网应激和未折叠蛋白反应中起重要作用。对于新生蛋白核糖体的结合和转运，RRBP1 是必不可少［30］。许多研究证实，RRBP1 在许多恶性肿瘤中过表达，例如宫颈癌、肺癌、膀胱癌和结直肠癌，并且影响患者的预后［31-34］。同时，也有研究证实，RRBP1 可能通过调节YAP1 的表达而成为顺铂耐药的主要驱动因素［35］。这些证据都说明，RRBP1基因水平的变化与癌症的发生发展具有一定的相关性。

MMP9 是一种参与多种生物活性的基质蛋白酶，是金属蛋白酶类中分子量最大的一种。MMP9在作用于细胞外基质（extracellular matrix，ECM）的构成和膜蛋白裂解中扮演着重要的角色，它可以切割多种ECM 蛋白来使得ECM 重构。MMP9 对基底膜的降解过程也起着非常重要的作用，基底膜中含有胶原，包括IV 型胶原，而这些胶原可被MMP9 降解，而基底膜破裂通常是肿瘤形成过程中肿瘤发生侵袭和转移的关键阶段［36-37］。目前，有部分研究已经证实，MMP9 在一些恶性肿瘤中过表达，例如非小细胞肺癌、宫颈癌和胰腺癌［38-40］。同时，一些研究提出，MMP9 与放射治疗有关，部分原因是ROS 的形成可以破坏锌离子和Cys99 之间的相互作用，从而激活MMP9，增加MMP9 的表达，影响放射治疗效果［41-42］。血浆MMP9 水平在肺癌、乳腺癌和肝癌放疗期间会升高［43-44］。在晚期非小细胞肺癌患者中，血清MMP9 水平在放射治疗后显著降低，但在SD和PD 的患者中，放射治疗后的血清MMP9 水平没有变化［45］。ESCC 患者NCRT 后，健康食管组织近端甚至远端MMP9 水平升高［46］。在最近的一项研究中，研究发现MMP9 的表达与小细胞肺癌患者的顺铂耐药有关［47］。因此，MMP9 与癌症病理广泛相关，不仅在肿瘤侵袭、转移和血管生成中发挥作用，同时还介导了肿瘤的微环境。

AKR1C3 是细胞内的一种多功能酶，是AKRs的一员。AKR1C1/3基因具有较高的序列同源性。AKR1C3 主要在内分泌器官中表达，参与肾上腺和肿瘤中类固醇的新生物合成，通过调节各种激素的合成以及其与相应受体结合的量，影响了肿瘤的发生发展［48］。AKR1C3 表达也会影响到与上 皮 间 充 质转 化 （epithelial-mesenchymal transitions， EMT）相关的许多转录因子，而EMT 也是促进肿瘤转移、肿瘤侵袭的主要作用机制之一［49］。有相关研究提出，基因AKR1C3 在乳腺癌和前列腺癌中表达上调，是其预后不良的标志［50-51］。在恶性肿瘤细胞中，AKR1C3基因表达上调可消除ROS，积累PGF2，进一步刺激MAPK 信号，防止G2/M 期阻滞或减少凋亡，促进癌细胞增殖，而导致其抗辐射［52-53］。AKR1C3 在非小细胞肺癌中表达较高，而在小细胞肺癌中表达则相反。而小细胞肺癌在临床中通常对放射治疗非常敏感［54］。同时，AKR1C3 表达上调可维持肿瘤细胞内的还原环境，不仅导致ROS 浓度降低，还可能与多西他赛耐药和蒽环类药物耐药有关［55-56］。因此，AKR1C3 不仅与多种恶性肿瘤的发生和发展密切相关，并且与多种抗肿瘤治疗也有一定的联系。

同时，本研究通过在食管鳞癌组织标本进行免疫组化分析，发现MMP9 和AKR1C3 蛋白表达水平差异与ESCC 患者对NCRT 敏感性具有一定的相关性。AKR1C3 阳性蛋白在食管鳞癌接受NCRT治疗不敏感的组织中有着中度到较强的细胞质反应，如图8 所示。MMP9 阳性蛋白在食管鳞癌接受NCRT 敏感的组织中的细胞质和细胞膜中有着更高的表达，如图9 所示。生存分析显示，AKR1C1 和AKR1C3 高表达的食管癌患者的OS 低于AKR1C1和AKR1C3 低表达的食管癌患者。但是，我们的研究也有相应的缺陷。首先，与ESCC 患者NCRT 敏感性相关的样本数量是有限的，需要在更大的队列中评估这些基因的疾病预测效用。其次，本研究提出的生物标志物应进行进一步实验验证，以明确其确切的生物学行为。

综上所述，本研究结果表明，基因MMP9和AKR1C3的异常表达与ESCC 患者接受NCRT 敏感性的生物学特性有着密切的关系，这为实现食管癌患者个体化治疗提供了线索。

作者声明：本文全部作者对于研究和撰写的论文出现的不端行为承担相应责任；并承诺论文中涉及的原始图片、数据资料等已按照有关规定保存，可接受核查。

学术不端：本文在初审、返修及出版前均通过中国知网（CNKI）科技期刊学术不端文献检测系统的学术不端检测。

同行评议：经同行专家双盲外审，达到刊发要求。

利益冲突：所有作者均声明不存在利益冲突。

文章版权：本文出版前已与全体作者签署了论文授权书等协议。