基于内质网应激基因识别胃癌分子亚型并建立 预后模型

朱 畅, 王红兵, 吕爱红, 徐 方

（1. 徐州医科大学, 徐州 221000；2. 徐州医科大学附属医院肿瘤科, 徐州 221000）

胃癌（stomach cancer, STAD）是发病率和死亡率均较高的恶性肿瘤, 其许多治疗策略包括靶向治疗和各种阻断性免疫检查点疗法, 如程序性死亡蛋白1（programmed cell death protein 1, PD-L1）已显示出显著的治疗疗效[1]。KEYNOTE-061和KEYNOTE-062试验显示：与单纯接受化疗的患者相比, 接受帕博利珠单抗（一种人源化的抗PD-1单克隆抗体）单药治疗的患者的总生存期（overall survival, OS）和无进展生存期（progression fres survival, PFS）均延长, 并且帕博利珠单抗在一线治疗中比化疗更有效[2]。针对HER2阳性患者, 以KEYNOTE-811（帕博利珠单抗联合曲妥珠单抗和化疗治疗HER2阳性胃癌III期临床试验）为代表的研究发现, 与曲妥珠单抗+化疗组相比, 联合应用帕博利珠单抗可以显著提高患者的中位OS。然而, 由于STAD的异质性, 一部分STAD患者对化疗、靶向治疗或免疫治疗产生耐药性, 导致癌症进展、复发或死亡。内质网是蛋白质合成、折叠、加工、运输的重要场所。蛋白质折叠和分泌受损导致未折叠或错误折叠蛋白质在内质网腔内积累, 这种现象称为内质网应激（endoplasmic reticulum stress, ERS）, 并引发一系列保护性级联信号反应, 被称为未折叠蛋白反应（unfolded protein response, UPR）[3]。肿瘤细胞经常侵入或转移到外部环境中, 诸如缺氧、葡萄糖缺乏、氧化应激等不利条件会损害内质网中的蛋白质折叠[4]。许多研究发现, 包括多发性骨髓瘤和乳腺癌在内的原发性人类肿瘤类型中, UPR的三个分支（PERK、ATF6、IRE1α）持续高水平激活[5]。此外, ERS通路的激活与Hp阳性胃癌显著相关[6]。GRP78是ERS通路中的启动因子, 在STAD中过表达, 并且与浸润深度、疾病分期、淋巴结转移、分化、淋巴浸润和血管侵犯正相关, 而与STAD患者总生存率负相关[7]。上皮-间充质转化（epithelial-mesenchy mal transition, EMT）是恶性肿瘤细胞迁移和侵袭的重要步骤。已有研究发现, ERS可以促进不同细胞系的EMT进程, 表明ERS也是引发EMT的重要因素之一[8]。因此, 细胞内质网应激是一种具有极大潜力的抗癌方法。本研究中, 我们基于内质网应激基因对STAD进行分型并构建新的预后模型, 同时探讨了分型间差异基因的功能通路、基于亚型构建的预后模型、免疫浸润差异等, 为STAD患者的个性化治疗提供了启示。

1 材料与方法

1.1 数据库

从癌症基因组图谱数据库（the cancer genome atlas program, TCGA）和基因表达综合数据库（gene expression omnibus database, GEO）分别获得了375例和433例STAD样本的RNA测序数据和相应的临床信息。

1.2 分子亚型的构建

从MSigDB数据库（http://gsea-msigdb.org）下载了61个内质网应激基因。通过“ConsensusCluster-Plus”R包基于内质网应激基因对STAD样本进行共识聚类且区分成不同的分子亚型[9], 使用经验累积分布函数（cumulative distribution function, CDF）图来帮助选择多个集群并获得集群稳定性。

1.3 STAD分子亚型的表征

通过使用“pheatmap”R包, 设置筛选标准为|log2FC| ＞0.585且错误发现率（flase discovery rate, FDR）=0.05, 分析不同分子亚型的基因表达水平和临床病理特征之间的相关性。此外, 在TCGA队列中, 使用“survival”R包通过Kaplan-Meier方法分析总生存率, 并通过对数秩检验评估生存分布之间的差异, 在GEO队列中用同样的方法进行验证。

1.4 功能富集分析

使用“limma”R包识别分型间的差异表达基因（differential expressed genes, DEGs）, 阈值为|log2FC|＞1且FDR=0.05。然后基于DEGs, 使用“clusterProfiler”R包进行基因本体论（gene ontology, GO）和京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）通路富集分析。调整后的P值＜0.05的GO通路和KEGG通路被认为具有统计学意义。

1.5 构建预后模型

利用单因素Cox分析来评估两种分子亚型之间的DEGs与STAD患者OS之间的关系。筛选出P＜0.05的DEGs, 进一步进行LASSO回归分析。风险评分等于每个mRNA的标准化表达水平和对应LASSO回归系数乘积的总和。根据风险评分中位值将STAD样本分为高低风险组, 并通过ROC曲线和Kaplan-Meier分析验证风险评分系统。利用参与风险评分的基因与临床病理参数通过“rms”R包构建列线图。

1.6 免疫细胞浸润概况分析

通过单样本基因集富集分析（single sample gene set enrichment analysis, ssGSEA）以及免疫差异分析来评估高、低风险组之间的免疫状态差异[10]。

2 结果与分析

2.1 STAD分子亚型的鉴定与验证

对内质网应激基因进行单回归Cox分析, 并应用“ConsensusClusterPlus”R包将所有肿瘤样本分为k（k=2～9）个不同的亚型。聚类分析结果（图1a）显示,k=2时, 组内相关性最高, 组间相关性最低, 表明STAD患者可以被准确分为两个亚型（C1型和C2型）。CDF曲线在k=2时分割效率较为优异, 表明STAD患者分为两个亚型时分布较为稳定（图1b）。对两个亚型间基因进行差异分析, 筛选出影响STAD患者分型的主要内质网应激基因LRRK2、DDIT3、MAGEA3、RNF186、TRIB3、BCL2（表1）。

图1 内质网应激相关STAD亚型的构建Fig. 1 Construction of endoplasmic reticulum stress-related STAD molecular subtypes

表1 两种分型间的内质网应激差异基因Tab. 1 Endoplasmic reticulum stress differential genes between the two types

2.2 分型的差异基因表达与临床相关性

基于TCGA数据库中STAD的mRNA表达数据和临床信息, 对 C1、C2分型进行进一步的差异分析, 可以看到基因表达谱与临床病理参数之间的相关性, 包括年龄（≤65或＞65岁）、性别（男性或女性）、临床分期（I、II 或III、IV )以及病理分级（G1、G2、G3）在两个亚型中的分布（图2）, 且病理分级（P＜0.01）与STAD分型显著相关。

Kaplan-Meier分析（图3a）显示, C1分型的生存结局较C2分型更好（P=0.029）。同样, 我们使用来自GEO数据库的433名STAD患者作为验证集以进一步评估分型的稳健性和可靠性。类似的结果在GEO队列得到验证（图3b）, C1分型的生存结局较C2分型更好（P=0.008）。

2.3 两种分型间的功能富集分析

GO富集和KEGG通路分析结果表明, 分型间差异基因参与了多种细胞和生物学功能（图4）。GO富集显示, DEGs主要参与细胞因子-细胞因子受体相互作用, 细胞黏附分子生成, Th1、Th2、Th17细胞分化, 中性粒细胞胞外陷阱形成, 补体与凝血级联等。KEGG通路分析显示, DEGs主要参与T细胞活化、细胞间黏附调节、T细胞受体复合物和质膜信号受体复合物等, 以此, 可以推测STAD的分型差异与细胞黏附、细胞分化、细胞因子及受体相互作用的相关性较大。

图4 两种分型间差异基因的功能富集分析Fig. 4 Functional enrichment analysis of differential genes between two molecular subtypes

2.4 预后模型的构建

基于两种分型间DEGs单因素Cox分析, 得到预后相关基因的森林图（图5a）。通过LASSOCox回归分析, 构建了一个预测STAD患者预后的8基因风险评分系统, 即风险评分=0.058×GPC3+0.181×DCLK1+0.229×ADRA1B+0.230×ANKRD1+0.051×SLC7A2+0.067×GAMT+0.057×RAB3B+0.068×FKBP10。单因素生存分析结果（图5b～5g）显示,GPC3、DCLK1、ADRA1B、SLC7A2、GAMT和RAB3B的较低表达与更好的OS 相关（P＜0.05）。

2.5 风险评分系统的验证与独立预后价值

根据风险评分中位值将STAD患者分为低风险组和高风险组, Kaplan-Meier分析（图 6a）显示了高风险组与低风险组之间的OS差异（P＜0.001）。主成分分析（principal component analysis, PCA）（图 6b）表明, 基于8基因特征的风险评分可将高低风险进行良好的区分。

将风险评分与临床病理参数纳入单因素和多因素Cox回归模型中。单因素Cox回归分析（图7a）表明, STAD患者的年龄、T分期、M分期、N分期、风险评分与 OS 相关（P＜0.05）。多因素Cox回归分析（图7b）揭示, 风险评分可作为患者OS评估的独立预测因子（P＜0.001）。此外, 年龄、M分期、N分期也可作为患者OS评估的独立预测因子（P＜0.05）。基于多因素Cox回归分析结果, 构建了一个可以预测STAD患者1、3和5年生存率的列线图模型（图7c）。ROC曲线下面积（area under curve, AUC）（图7d）分别为0.702（1年）、0.741（3年）、0.762（5年）, 这表明该模型在预测STAD患者的生存方面具有良好的性能。

2.6 免疫细胞浸润差异分析

ssGSEA以及免疫差异分析发现, 共有8个免疫细胞亚群及免疫相关功能在高低风险组分布差异具有统计学意义（图8）。其中, 低风险组的NK cells（P＜0.01）、Th2 cells（P＜0.001）、APC co-inhibition（P＜0.001）、Cytolytic activity（P＜0.05）、MHC class I（P＜0.001）评分显著高于高风险组, Mast cells（P＜0.001）、T helper cells（P＜0.05）、Type II IFN response（P＜0.001）在高风险组中得分更高。

图8 免疫细胞浸润差异分析Fig. 8 Differential analysis of immune cell infiltration

3 讨论

由于早期诊断困难、肿瘤易复发性和预后不良等原因, STAD在临床治疗上具有挑战性。研究表明, ERS是多种疾病的潜在治疗靶点, 包括糖尿病、神经退行性疾病和癌症[3]。此外, ERS影响肿瘤微环境和肿瘤免疫治疗[11]。研究ERS在STAD中的作用有助于确定不同分型中的潜在治疗靶点。目前, 已鉴定出一种新的内质网应激基因来预测胶质瘤患者预后[12]。然而, 尚未有研究阐明内质网应激基因在STAD中的作用, 我们的研究旨在阐明这一作用, 为STAD患者的分析提供另一种研究方向。

两种分型间的内质网应激差异基因包括LRRK2、DDIT3、MAGEA3、RNF186、TRIB3、BCL2。LRRK2是一种多结构域细胞质蛋白[13], 其与强大的抗肿瘤活性有关, 如抑制肿瘤细胞的增殖、迁移和侵袭, 以及诱导细胞凋亡和细胞周期停滞[14]。DDIT3是一种内质网应激基因, 可介导DNA损伤、ER应激和细胞生长的反应, 其蛋白属于CCAAT增强子结合蛋白b（CEBPb）家族, 可调节细胞分化[15]。此外, 有研究发现,DDIT3可以促进STAD中肿瘤进展和CSC干性[16]。MAGEA3是睾丸癌抗原（cancer testis antigen, CTA）的主要成员, 位于染色体Xq28上, 其表达受DNA甲基化或组蛋白乙酰化的调节。已有研究证明,MAGEA3的高表达与STAD的淋巴结转移有关[17]。RNF蛋白家族是一组包含一个RNF结构域的复杂蛋白质, 该家族的一些成员与结肠癌的发展和进展有关。已有研究表明, RNF186可通过调节结肠上皮细胞中的内质网应激来维持肠道稳态[18]。Tribbles假激酶家族由TRIB1、TRIB2、TRIB3、STK40组成, 它们是在真核生物的大多数组织中发现的进化保守的假激酶[19]。TRIB3由位于人类第20号染色体上的6个外显子编码, 是一种约65 kD的蛋白质, 含有358个氨基酸[20]。TRIB3参与多种信号通路的激活, 如MAPK通路。内质网应激、缺氧和营养缺乏会上调TRIB3表达[21]。BCL2具有明显抑制细胞凋亡的作用, 其过表达或异常激活直接参与维持癌细胞的存活和生长以及促进治疗耐药性。此外, 异常的BCL2激活与多种癌症的发生、进展、转移和复发有关, 因此,BCL2已成为评估临床治疗效果和预后的重要指标[22]。

在参与构建风险评分系统的基因中,GPC3是一种膜结合硫酸乙酰肝素蛋白多糖, 在大多数肝细胞癌病例中过表达。一项研究发现,GPC3的高表达可用于提示胃腺癌的不良预后[23]。DCLK1是胃肠道簇细胞的标志物, 从机制上来讲, 它可以从NOTCH、WNT分子信号通路等促进癌细胞的生长和进展。先前的研究表明,DCLK1在STAD中过度表达, 并在肿瘤的EMT和癌症干细胞中发挥功能性作用[24]。ADRA1B是肾上腺素能受体α1（ADRA1）亚家族的成员, ADRA1可以促进癌症的转移, 并且持续激活的ADRA可以通过p53诱导细胞凋亡[25]。有研究发现,ADRA1B的高水平表达与STAD较差的OS相关[25]。SLC7A2对L-精氨酸、赖氨酸和鸟氨酸的转运至关重要, 其遗传多态性与结直肠癌进展相关,SLC7A2还被发现在非小细胞肺癌的放射抗性中起作用[26]。GAMT编码的蛋白质是一种甲基转移酶, 它使用S-腺苷甲硫氨酸作为甲基供体将胍基乙酸转化为肌酸,GAMT的高表达被发现与STAD的预后不良相关[27]。RAB3B是RAS癌基因家族的成员, 被证明是miR-200b的靶点, 被认为是STAD中的肿瘤抑制因子[26]。

在本研究中, 我们基于内质网应激基因对STAD患者进行分型并建立预后模型, 探索了分型间差异基因所涉及的相关信号通路。GO富集和KEGG通路分析显示, DEGs主要参与细胞因子-细胞因子受体相互作用, 细胞黏附分子生成, Th1、Th2、Th17细胞分化, 中性粒细胞胞外陷阱形成、补体与凝血级联等, 以及T细胞活化、细胞间黏附调节、T细胞受体复合物、质膜信号受体复合物等方面。

此外, 低风险组的NK cells、Th2 cells、APC-co inhibition、Cytolytic activity、MHC class I评分显著高于高风险组, Mast cells、T helper cells、Type II IFN response在高风险组中得分更高。通过对高、低风险组患者免疫细胞亚群及免疫相关功能差异的对比, 我们不难推测出低风险组患者的预后优于高风险组患者。此外, 我们的研究表明, 这些免疫细胞及相关功能或许在STAD的发生、发展中发挥着作用, 这为寻找新的免疫检查点起到一定的指导作用。除了这些分析外, 我们还通过建立预测患者1、3、5年的列线图, 有效地预测STAD患者1、3、5年的生存概率, 有助于临床医生进一步了解患者的预后。

然而, 本研究也存在一些局限性, 仅使用了在线数据库TCGA和GEO的数据集进行分析, 后续需要更多来自不同地区的患者数据进行验证。此外, 由于数据集来源于数据库, 缺乏治疗反应的数据, 需要进一步分析分型和风险组间的治疗反应, 以帮助临床治疗。最后, 确定了的差异基因所涉及的功能和途径, 以及它们如何调节免疫浸润, 需要进一步的生物学试验来验证。