胰腺癌预后相关铁死亡基因筛选及预后预测模型构建

何冲，周姝睿，李雅晴，陈少杰，练国达，陈尚祥，黄开红

胰腺癌是消化道常见恶性肿瘤，侵袭性强、恶性程度高，患者5 年生存期低于8%［1］。胰腺癌早期临床症状隐匿，病情进展迅速，多数患者确诊时已发生局部进展或远端转移，失去手术机会，即使患者早期接受手术切除，90%以上患者会出现复发［2］。探索新的预后预测指标，建立便捷、有效的预测系统，积极开展科学的临床干预，是提高胰腺癌临床预后的关键。

铁死亡是一种铁依赖性的细胞程序性死亡，主要是由线粒体内的脂质过氧化物损伤积累驱动细胞死亡［3］。据研究报道，铁死亡在胰腺癌发生发展中发挥至关重要的作用，例如天冬氨酸氨基转移酶（GOT1）驱动的代谢途径在胰腺癌用于维持氧化还原平衡可负性调控铁死亡，促进肿瘤细胞增殖、肿瘤生长［4］。近年来的研究还发现，铁死亡与MAPK 通路密切相关。在胰腺癌细胞中，铁死亡诱导剂诱导stat3 的激活需要MAPK/ERK 通路的激活［5］。研究还发现ASK1⁃p38 MAPK 通路是脂质过氧化物下游的调节因子之一［6］。但铁死亡相关基因与胰腺癌患者的预后是否相关仍未明确。为探讨铁死亡相关基因对胰腺癌患者预后预测的价值，并进一步探索其分子机制，本研究通过分析公共数据库中胰腺癌患者的转录组表达谱和临床资料，筛选差异表达的铁死亡相关基因，验证其在胰腺癌患者预后预测中的准确性，以便为胰腺癌患者预后评估及治疗提供依据。

1 材料与方法

1.1 数据收集

从 GTEx（Genotype⁃Tissue Expression Project）（https：//www.gtexportal.org/）数据库中下载正常胰腺组织转录组表达谱数据，从TCGA（The Cancer Genome Atlas）（https：//portal.gdc.cancer.gov/）数 据库及ICGC（International Cancer Genome Consortium）（https：//dcc.icgc.org）数据库中下载胰腺癌患者的转录组表达谱及相应的临床病理数据，排除TCGA和ICGC 中无生存状态和生存时间的样本，保留同时具有转录组和临床数据的样本，最终在GTEx 数据库中纳入328 例正常胰腺组织转录组表达谱数据，在TCGA 数据库中纳入了176 名胰腺癌患者作为训练队列，在ICGC 数据库中纳入了81 名胰腺癌患者作为独立外部验证队列。

从铁死亡数据库（FerrDb）（http：//www.zhounan.org/ferrdb/legacy/index.html）中下载铁死亡相关基因，包括123 个铁死亡标记基因、150 个铁死亡驱动基因和109 个铁死亡抑制基因，并与TCGA 和ICGC 数据库中的基因表达谱取交集，最终纳入了382 个铁死亡相关基因进行研究。

1.2 方法

1.2.1 铁死亡相关差异表达基因的鉴定和功能分析 使用 R 软件（Version 4.0.2）“limma”包研究mRNA 的差异表达。在TCGA 或GTEx 中分析了调整后的P值以纠正假阳性结果。使用Fold change和校正后P值绘制火山图。筛选条件为“|log2（Fold change）|≥2 且AdjustedP<0.05”，检测肿瘤与癌旁正常组织之间的差异表达基因。使用R 软件“pheatmap”包可视化肿瘤与癌旁正常组织之间与癌旁正常组织之间的差异程度范围。然后，将铁死亡相关基因与差异表达基因进行交集，得到铁死亡相关差异表达基因。

使用 Metascape Online（https：//metascape.org/gp/index.html#/main/step1）进行功能分析［7］。将铁死亡相关基因上传到Metascape 进行功能分析，构建PPI 网络。我们使用MCODE 插件（筛选条件为P<0.05）进行模块化分析。

1.2.2 铁死亡相关差异表达基因预后模型的构建与验证 TCGA 数据集作为训练队列，用于建立铁死亡相关差异表达基因预后模型，ICGC 数据集作为验证队列验证预后模型的预测效果。首先对总生存期（overall survival，OS）进行单变量COX 回归分析以确定有与预后相关的铁死亡相关差异表达基因（以P<0.05 表示差异具有统计学意义）。我们进一步使用 R 软件（Version 4.0.2）“glmnet”包，整合生存时间、生存状态和与预后相关的铁死亡相关差异基因表达数据，利用LASSO COX 方法进行回归分析，构建了一个基于铁死亡相关基因的预后模型。该模型根据各基因的FPKM 标准化表达水平及其相应的回归系数计算患者的风险评分。根据中位值将队列中所有患者分为高风险组与低风险组，分析患者风险评分与患者的总生存期、生存状态及各个基因的表达变化的关系。输入患者的总生存期、生存状态及铁死亡相关差异表达基因预后模型风险评分，使用R 软件（Version 4.0.2）“pROC”包进行了ROC（receiver operating character⁃istic curve，受试者工作特征曲线）分析以评估预测模型准确性。并根据风险评分分组，使用R 软件“survival ROC”包绘制生存曲线，进行生存分析。

1.2.3 训练队列和验证队列高低风险组差异基因功能富集分析 根据预后预测模型计算TCGA 训练队列和ICGC 验证队列各患者风险评分，以每个队列风险评分中位数划分高低风险组，分别对TCGA 训练队列和ICGC 验证队列数据进行差异基因重分析。采用 R 软件（Version 4.0.2）“limma”包进行分析，筛选条件为“0.25

1.2.4 训练队列肿瘤微环境免疫浸润分析 根据TCGA 训练队列所有肿瘤组织的转录组表达谱数据，我们采用R 软件（Version 4.0.2）“CiberSort”包计算TCGA 训练队列所有肿瘤组织中浸润的免疫细胞，并对同一类免疫细胞在高低风险组肿瘤组织中的浸润情况进行统计学分析。

1.2.5 统计学分析 所有统计分析均采用R（Ver⁃sion 4.0.2）或 SPSS（Version 25.0）进行。如果未特别注明，以P<0.05 表示差异具有统计学意义。

2 结 果

2.1 胰腺癌中差异表达的铁死亡基因

从TCGA 数据库和GTEx 数据库中下载胰腺癌转录组数据，共获得510 个样本（肿瘤样本178 个，正常组织样本332 个）。使用“limma”程辑包研究mRNA 的差异表达，共筛选出4149 个差异基因（上调基因有4011 个，下调基因有138 个）（图1A⁃B），与铁死亡数据库（FerrDb）中铁死亡相关基因取交集，维恩图显示最终获得103 个差异基因（图1C）。为了探索胰腺癌铁死亡的潜在机制，我们使用Metascape Online 进行基因功能富集分析。如图1D⁃F 所示，GO（Gene Ontology，基因本体论）分析结果表明，这些差异表达的铁死亡相关基因主要在活性氧应激、核受体通路及细胞应激等通路富集。如图1F 所示，基于Metascape Online 的PPI（protein⁃protein interaction，蛋白质⁃蛋白质相互作用）网络和MCODE 插件确定了这些铁死亡基因中的5 个重要模块。因此，铁死亡相关差异基因可能通过这些机制影响胰腺癌患者预后。

图1 TCGA 及GTEx 中胰腺癌内差异表达的铁死亡相关基因分析

2.2 在TCGA 队列中构建预后模型

首先下载 TCGA⁃PAAD 归一化的 mRNA 表达数据和相应的患者信息，提取103 个铁死亡相关差异表达基因的表达数据，对总生存期（overall survival，OS）进行单因素COX 回归分析。如图2A显示了单变量Cox 回归分析结果的森林图，37 个铁死亡相关基因被鉴定为预后相关基因。将得到的37 个预后相关的铁死亡相关基因进行LASSO COX 回归分析以构建预后模型，根据最优的λ 值（λ=0.07）确定了一个15基因的预后模型（图2B⁃C）。根据风险评分的中位数将TCGA⁃PAAD队列分为高风险组（n=86）和低风险组（n=86）。Kaplan⁃Meier生存曲线显示，高风险组胰腺癌患者的总生存期明显低于低风险组患者（P<0.05）（图2E）。随后，我们进行了ROC 分析来评估模型的预后预测效率，发现AUC（area under curve，曲线下面积）分别为 0.74（1 年）、0.82（3 年）和 0.88（5 年），表明预后预测模型的预测准确性良好（图2F）。

图2 构建铁死亡相关差异表达基因预后模型

2.3 在ICGC 队列中验证模型

此外，为了评估预测模型的可靠性和准确性，我们在ICGC⁃PAAD 队列中验证了预测模型的功能。我们使用相同计算公式计算出了ICGC⁃PAAD队列中所有胰腺癌患者的风险评分，根据风险评分的中位数将TCGA⁃PAAD 队列分为高风险组（n=41）和低风险组（n=40）。Kaplan⁃Meier 生存曲线显示，高风险组胰腺癌患者的总生存期明显低于低风险组患者（P<0.05）（图3B）。我们进一步建立了受试者操作特征曲线（ROC）来评估模型的预后预测效率，发现AUC（area under curve，曲线下面积）分别为0.61（1 年）、0.69（3 年）和 0.77（5 年），表明在ICGC 队列中预后预测模型的预测准确性良好（图2F）。

图3 验证铁死亡相关差异表达基因预后模型

2.4 在TCGA 队列和ICGC 队列中进行差异表达基因功能再分析

为了研究高低风险组之间的分子异质性和潜在的生物学过程与途径，我们使用R 软件的“limma”软件包研究 TCGA⁃PAAD 队列及 ICGC⁃PAAD 队列高低风险组的差异表达基因，TCGA⁃PAAD 列中发现了6189 个差异表达基因（|log2（Fold change）|≥2且AdjustedP<0.05），ICGC⁃PAAD队列发现了1804个差异表达基因（|log2（Fold change）|≥2 且AdjustedP<0.05），并分别对两组队列的差异基因进行了GO和KEGG分析。TCGA⁃PAAD队列的GO和KEGG分析显示，这些差异表达基因主要富集在以下几个方面：适应性免疫反应、细胞间信号传递、细胞间黏附及配体-受体反应等（图5A⁃B）。ICGC⁃PAAD队列的GO 和KEGG 分析显示，这些差异表达基因主要富集在以下几个方面：组织生长、细胞外基质形成、免疫系统反应及细胞黏附等（图5C⁃D）。免疫相关通路同时出现在TCGA 和ICGC 队列的通路富集结果中，提示我们应继续研究铁死亡基因集与肿瘤免疫微环境之间的关系。

2.5 免疫微环境分析

图4 TCGA 和ICGC 队列高低风险组差异表达基因的通路富集分析

为了进一步探讨风险评分与免疫之间的关系，我们基于TCGA⁃PAAD 队列基因表达数据，应用“CiberSort”算法探索了TCGA⁃PAAD 队列中免疫细胞的比例。如图5A 所示，CD4+静息T 细胞、巨噬细胞（包括M0、M1 和M2 亚群），在TCGA⁃PAAD中占浸润性免疫细胞的很大比例。此外，我们采用箱式图来可视化低风险和高风险得分组之间的免疫细胞亚群分布（图5B）。我们发现高风险组中幼稚B 细胞、浆细胞、CD8+T 细胞浸润程度较低风险组低，但M0 和M1 巨噬细胞的浸润程度高（P<0.05）。此外，我们还研究了PD⁃1 和CTLA4在TCGA⁃PAAD 队列低风险组和高风险组中的表达，发现在TCGA⁃PAAD 队列低风险组患者中PD⁃1和CTLA4 的表达较高风险组患者更高（图5C⁃D）。

图5 TCGA 训练队列中高低风险组免疫细胞分布

3 讨 论

铁死亡是一种铁离子依赖的细胞程序性死亡［3］。近年来的研究发现，选择性诱导肿瘤细胞铁死亡，可增强肿瘤细胞化疗效果［8，9］，但在胰腺癌中相关机制尚不明确。在胰腺癌中，胰腺癌具有半胱氨酸依赖的特性，通过半胱氨酸合成谷胱甘肽降低自身铁死亡水平［10］。研究胰腺癌铁死亡相关基因表达和预后的关系，寻找胰腺癌铁死亡治疗靶点，对于增强胰腺癌患者肿瘤治疗疗效、改善患者预后，具有重要临床意义。本研究构建了胰腺癌铁死亡相关基因预后预测模型，并发现TCGA 队列中高低风险组的免疫微环境存在差异。本研究提示，胰腺癌患者的肿瘤微环境、临床预后与铁死亡相关基因表达存在一定联系，可为评估患者预后、制定化疗方案提供线索。研究表明，铁死亡核心机制如下：Fe2+依赖（为氧化过程提供电子）、GPX4 失活（抗氧化剂失衡）和ROS 产物（主要来源于芬顿反应）［11］。铁死亡是多种酶参与的级联反应，铁死亡相关基因在肿瘤进展中发挥重要作用［12，13］。本研究筛选出胰腺癌中肿瘤组织和正常癌旁组织的差异表达基因，有103个是铁死亡相关基因。进一步分析发现，其中31 个差异表达的铁死亡相关基因与胰腺癌患者的预后相关。这一发现提示铁死亡相关基因会影响胰腺癌患者的预后。

通过LASSO 回归分析，我们构建了一个基于15 个铁死亡相关基因（ANO6、AURKA、CAV1、CD⁃KN1A、CISD1、ENPP2、HSPB1、MAP1LC3A、NRAS、PML、PROM2、PTGS2、RRM2、SLC2A6 和SP1）的预后预测模型。ANO6 作为Ca2+活化的氯离子通道，可以破坏质膜的磷脂双分子层，介导肿瘤细胞铁死亡［14］。AURKA 是一种丝氨酸/苏氨酸激酶，AURKA 受抑制会减少GPX4 的表达并促进细胞铁死亡［15］。CAV1 可正向调节胰腺癌细胞的增殖、迁移侵袭和耐药性［16］。CDKN1A 可通过抑制细胞凋亡来应对氧化应激，研究表明，p53 介导的CDKN1A 表达延缓了肿瘤细胞中半胱氨酸剥夺之后的铁死亡［17］。CISD1是一种含铁的线粒体膜外膜蛋白，负向调节肿瘤细胞铁死亡，促进肿瘤增殖和转移［18，19］。ENPP2 作为脂质激酶参与脂质代谢，对Erastin 诱导的心肌细胞铁死亡有保护作用［20］。热休克蛋白（HSPs）在真核生物中高度保守，作为分子伴侣参与蛋白质的组装、折叠、输出和翻转的调节［21］，其中HSPB1 是铁死亡的负性调节因子［22］。MAP1LC3A 通过 ATG5⁃ATG7⁃NCOA4 途径激活细胞自噬，导致铁蛋白降解，从而导致细胞内游离态铁积聚，最终促进肿瘤细胞铁死亡［23］。NRAS 已被证实与肝细胞癌中的erastin耐药有关［24］。PML可通过调节线粒体代谢增强人卵巢癌化疗敏感性［25］。PROM2 在乳腺癌中能通过促进胞内铁向胞外转运，增强肿瘤细胞对铁死亡耐受性［26］。PTGS2 参与铁死亡过程当中的脂质代谢，可正向调节细胞铁死亡［27］。RRM2 在肝细胞癌的谷胱甘肽合成中发挥作用，下调肝癌细胞铁死亡水平［28］。SLC2A6在NRF2 通路中发挥作用，抑制胰腺癌细胞铁死亡［29］。SP1是调节ACSL4表达的一个关键的转录因子，通过与ACSL4启动子区域结合来增加ACSL4 的转录，参与调节铁死亡过程中的脂质代谢［30］。这些基因都可以通过调节肿瘤铁死亡影响肿瘤发展。这一预测模型的预后预测价值在ICGC 队列中得到了验证，在胰腺癌的临床诊疗中具有应用价值。

为了探索这一预测模型内在分子机制，本研究分析了TCGA 队列中高低风险组的差异表达基因并进行了功能富集分析，结果提示差异基因的功能通路富集在适应性免疫反应上，并且高低风险组的免疫微环境存在差异。通过分析TCGA 队列中高低风险组的免疫细胞浸润情况，我们发现高风险组患者的CD8+T 细胞浸润程度明显低于低风险组。研究表明，CD8+T 细胞可以调节癌细胞铁死亡，进而可以影响癌症免疫治疗的疗效［31］。免疫检查点抑制剂已被用于治疗膀胱癌和非小细胞肺癌，具有显著的疗效，但只有一小部分胰腺癌患者可以从免疫检查点抑制剂治疗中受益［32］。分析 TCGA 队列中高低风险组 PD⁃1 和 CTLA4 表达情况发现，高风险组PD⁃1 和CTLA4 表达明显低于低风险组。这一结果提示高风险组患者或许不能从免疫检查点抑制剂治疗中受益，我们需要为高风险组胰腺癌患者寻找新的治疗靶点。

然而，本研究依旧存在一些局限性。首先，本研究通过TCGA 数据集构建模型，通过ICGC 数据集验证该模型，但这一模型并没有在临床队列中得到验证；其次，本模型中纳入的铁死亡相关基因，没有通过基础实验验证它们在胰腺癌中的生物学功能；第三，本研究只是初步证明胰腺癌中铁死亡相关基因表达和患者预后及患者肿瘤免疫微环境存在一定联系，但内在分子机制仍需要更多的基础研究来证明。