炎症对胰腺癌患者预后的评估及免疫状态的影响

刘 基，孙振纲，邓 岩，唐子扬

（1.长江大学医学部，湖北 荆州 434000；2.长江大学附属荆州医院，湖北 荆州 434000）

胰腺癌（pancreatic adenocarcinoma，PAAD）是发生于胰腺组织的恶性肿瘤，早期症状隐匿，恶性度高，进展迅速，预后差，5 年生存率仅为2%～9%[1，2]。PAAD 风险因素包括吸烟、慢性胰腺炎家族史、糖尿病、肥胖等[3，4]，在众多风险因素中炎症过程已成为PAAD 发生和发展的关键介质[5]。炎症和癌症之间的功能关系是当前临床研究的热点之一，多项研究表明[6-8]，肿瘤微环境在很大程度是由炎性细胞协调，炎症是肿瘤过程不可缺少的参与者，可促进其增殖、存活与转移。本研究基于炎症相关的PAAD 预后模型，以期为PAAD 分子机制以及预后情况提供新的思路。

1 资料与方法

1.1 数据收集 从GTEx 数据库（https://xena.ucsc.edu/）中下载167 例正常胰腺样本的RNA 测序数据，从TCGA 数据库（https://portal.gdc.cancer.gov/）中下载4例正常胰腺样本和178 例PAAD 样本的RNA 测序数据和随访信息，以及从GEO（https://www.ncbi.nlm.nih.gov/geo/）数据库GSE57495 和GSE62452 共下载132 例PAAD 样本的RNA 测序数据和随访信息，使用制造商提供的注释文件将探针与基因匹配，如果有多个探针与单个基因匹配，则取RNA 基因测序数据的中值，根据GTEx、TCGA 和GEO 的数据访问政策和发行的指南（数据均为公开数据），在Molecular Signatures 数据库（https://www.gsea-msigdb.org/gsea/msigdb/）[9]中下载炎症反应相关基因。

1.2 炎症反应相关基因的预后模型构建 肿瘤样本和正常样本之间差异表达基因（DGEs）通过R 语言的“limma”包以logFC 大于2 和P值小于0.05 进行差异分析和筛选，单变量Cox 分析用于筛选具有预后价值的炎症反应相关基因，并通过Benjamini&Hochberg（BH）校正方法调整P值。Lasso Cox 分析用于构建预后模型，同时缩小过度拟合的风险，使用R语言“glmnet”包用于去收缩变量，使得回归系数相当于零，建立预后模型。用于预后分析的表达差异基因的标准化表达矩阵是自变量，在TCGA 中患者的生存期和生存状态是因变量，根据每个炎症反应相关基因的表达水平及其相应的回归系数计算患者的风险评分。公式如下：分数=sum（每个基因的表达量×对应系数）。根据风险分数的中间Cut-off 值，将患者分为高风险组和低风险组。采用R 语言“scatterplot”包进行PCA 分析探讨不同风险组的分布情况，采R 语言“survminer”包进行生存分析，采用“survival”R 包和“timeROC”R 包进行时间依赖性ROC 曲线分析，以评估预后特征的预测价值。此外，行单变量和多变量Cox 分析，以探讨基因特征的独立预后价值。

1.3 功能富集分析和基因富集分析 采R 语言对差异表达基因行京都基因与基因组分析（KEGG）和基因本体分析（GO）分析，使用R 语言“bioconductor”包行基因名和富集通路名称的匹配，以差异基因表达谱为变量，以KEGG 和GO 分子特征和通路的表达谱为参考，同时计算logFC 值，P值根据BH 方法进行调整。使用“DOSE”“clusterProfiler”“pathview”“enrichplot”包进行KEGG 和GO 富集分析，使用R语言“digest”“GOplot”包进行基因和富集通路联系分析，以logFC 值为表示两者之间的关系度。

1.4 特征基因验证 为了测试从TCGA 队列构建的模型稳定性和胰腺癌生存状态和时间的影响因素，GEO 数据库中GDE57495 和GSE62452 的患者也根据TCGA 队列的风险系数取中间Cut-off 值分为高风险组或低风险组。

1.5 肿瘤微环境和免疫反应分析 将200 个炎症基因用R 语言“limma”包进行数据矫正，使用R 语言“estimate”包给每个样本行肿瘤环境、免疫环境、肿瘤纯度进行评分，使用“pheatmap”绘制高、低风险组在肿瘤免疫微环境的表达情况，使用R 语言“ggpubr”包分析高、低风险组在不同肿瘤免疫微环境的表达差异。

1.6 生存预测列线图的建立与评估 为了准确预测1、2、3 年的总生存率（0S），通过整合年龄、性别、风险分数构建预后列线图，Harrell’s 一致性指数（C 指数）评估预测准确性，C 指数范围从0.5（无预测能力）到1（完美预测）。校准图用于评估列线图的性能特征。

1.7 数据分析 采用Wilcoxon 检验比较癌症组织和正常组织之间的差异表达基因，采用Kaplan-Meier分析比较各组OS 差异。采用单变量和多变量Cox分析筛选OS 的独立预测因子。使用estimate 包对各组肿瘤进行基质细胞评分、免疫细胞评分、综合评分（基质细胞评分+免疫细胞评分）、肿瘤纯度评分。R软件（版本4.1.1）和软件包ggplot2、igrph、pheatmap、ggpubr 和survminer 被用来创建平面图。双尾P值小于0.05 表示差异有统计学意义。

2 结果

2.1 炎症相关差异表达基因的识别 共得到29 个差异表达基因，见图1A、1B。除HPN 基因在正常组织中表达高于肿瘤组织，其余基因在肿瘤组织表达均高于正常组织（P＜0.05）；单变量COX 分析显示，17个基因与OS 相关，MET 基因的风险率为1.960（95%CI1.555～2.470，P＜0.001），见图1C、1D。

图1 TCGA 和GTEx 数据序列中识别炎症相关差异表达基因

2.2 炎症相关基因预后模型的构建 将差异基因用Lasso Cox 方法进行分析，构建预后模型，根据最佳λ 值确定了6 个标记基因（图2A、2B），风险分数根据score=0.047×CXCL9 的基因表达量+0.028×LY6E的基因表达量+0.038×HBEGF 的基因表达量+0.524×MET 的基因表达量+0.096×CXCL10 的基因表达量+0.051×RTP4 的基因表达量进行计算，根据中间Cut-off 值划分为高、低风险组（图2C）；散点图显示，高风险组患者比低风险的患者预后差（图2D）；PCA 分析显示，患者在高、低风险组中可以被分为两个不同方向（图2E）；Kaplan-Meier 曲线显示，高风险组总体生存率OS 低于低风险组（P＜0.05，图2F）；ROC 曲线分析显示，曲线下面积（AUC）等于0.774，表明预后模型准确性较高（图2G）。使用GEO数据库测试在TCGA 数据库中构建模型的稳定性，从GEO数据库抽取GSE57495、GSE62452 根据TCGA 数据风险分数的计算方法和中间值将患者划分为高风险组和低风险组，得到的结果和TCGA 相类似（图2I），PCA 将患者划分为两个方向的分布（图2H），患者在高风险组较低风险组提前死亡（图2J），并且相对低风险组生存期更短（图2K）；此外，6个基因AUC 曲线是0.632（图2L）。

图2 6 个预后基因在GEO 和TCGA 预后模型的分析

2.3 胰腺癌的预后因素 单因素Cox 分析显示，风险评分和OS 相关（HR=3.19，95%CI=2.077～4.911，P＜0.001）；在纠正其它混合因素后，多因素Cox 分析显示，风险评分也和OS 相关（HR=3.302，95%CI=2.074～5.260，P＜0.001）；ROC 曲线分析显示，风险分数具有良好预后预测准确度，且PAAD 具有良好的预后价值，见图3。

图3 胰腺癌的预后因素

2.4 生物功能和通路富集分析 通过GO 富集分析显示，29 个差异表达基因在PAAD 患者的生物进程（BP）主要在细胞趋化作用、细胞对脂多糖的反应、细胞对细菌起源分子的反应、细胞对生物刺激的反应、骨髓白细胞游走等5 个进程占主导作用，癌症细胞主要有膜锚定构件、分泌颗粒膜、膜筏、膜微区等主要细胞组分（CC），分子功能（MF）主要是受体配体活动、信号受体激活物活性、细胞因子活性、G 蛋白偶联受体结合等分子功能参与活动（图4A）。通过KEGG 富集分析显示，29 个差异表达基因主要参与细胞因子受体相互作用、病毒蛋白与细胞因子及细胞因子受体的相互作用生物功能及肿瘤坏死因子信号通路、趋化因子信号通路、Toll 样受体信号通路（图4B）；同时，建立这些通路和基因之间的联系（图4C～图4F）。

图4 生物功能和富集通路分析

图4 生物功能和富集通路分析（续）

2.5 免疫活性和肿瘤微环境分析 低风险组免疫通路较高风险组更加活跃（图5A），另低风险组免疫评分、肿瘤评分和免疫及肿瘤微环境综合评分高于高风险组，相反肿瘤纯度呈现相反趋势，肿瘤纯度从低风险到高风险组呈上升趋势（Kruskal-Wallis test，P＜0.01，图5B～图5E）。

图5 免疫活性和肿瘤微环境

图5 免疫活性和肿瘤微环境（续）

2.6 建立列线图对总生存率进行预测 在TCGA 数据中，基于年龄、性别、风险分数预后因素建立关于PAAD 患者总生存率的列线图（图6A），校准图进一步证实列线图能够较好的预测患者2、3 年的总体生存率，是一个较好的预测模型（图6B、图6C）。

图6 列线图的构建

3 讨论

随着精准医学和公共卫生的发展，PAAD 得到了较好的治疗以及控制，然而因为实用标记物较少，现在仍然无法有效的进行早期诊断以及PAAD 患者治疗评估。CA199 作为PAAD 中的一种诊断性肿瘤标志物的作用已被证实，由于受到宿主炎症反应和梗阻性黄疸的干扰，从良性疾病中准确诊断PAAD 的个体能力较低[10，11]。然而，作为PAAD 预后因子的炎症反应相关基因标记物极少报导。既往研究表明[12-14]，铁死亡相关基因、免疫相关基因、m6A相关基因显著影响PAAD 预后。本研究中构建的炎症反应相关基因与上述基因相比具有更多优势，如本研究结合了GTEx 数据库中更多的正常样本，使得差异基因更加具有差异性，同时探索了免疫活性和肿瘤微环境在高、低风险组的表达状态，并且构建列线图对2、3 年PAAD 总体生存率进行预测。

本研究建立的预后模型由6 个炎症反应相关差异表达基因组成（CXCL9、LY6E、HBEGF、MET、CXCL10、RTP4），这些基因在PAAD 中均上调。CXCL9、CXCL10 属于趋化因子，与肿瘤微环境中免疫环境变化有关，其中CXCL9 可以通过STAT3 信号通路影响CD8+T 细胞，CXCL9 表达量的提升会抑制CD8+T 细胞的抗肿瘤细胞因子的增殖、活化、分泌[15-17]，而CXCL10 可以通过其受体CCR7 和CXCR3 促进PAAD 细胞的迁移[18]，HBEGF 可以诱导EGFR 核转位和组蛋白H4 甲基化，进而介导EGFR 的代谢激活，而EGFR 配体的转录上调引起EGFR 信号的传导，导致谷氨酰胺饥饿，从而促进体外胰腺腺泡细胞DNA 损伤的消退。在刺激YES 相关蛋白（YAP）的活性时，则影响PAAD 的生存预后[19-21]。RTP4 由Ⅰ型IFN（IFN-Ⅰ）诱导，并与TANK 结合激酶（TBK1）复合物结合，通过干扰TBK1 和IFN 调节因子3 的表达和磷酸化来负调节TBK1 信号传导[22]。本研究中通过对PAAD 预后模型的构建，进一步验证此6 个PAAD 相关炎症基因可以作为PAAD 预测分子。

已有研究表明[23]，Toll 样受体（TLR）在炎症和免疫反应中起关键作用，其主要参与病原体相关分子模式（PAMP）和危险相关分子模式（DAMP）。本研究在免疫微环境分析中，高风险组患者的免疫细胞浸润、抗肿瘤免疫活动及两者综合活动均较多，表明高风险组患者的免疫功能整体受损，免疫活动的增加可以解释低风险组患者具有良好的临床结果，然而高风险组的肿瘤纯度高于低风险组，这表明不良预后与肿瘤微环境无关。同时，有研究表明先天免疫系统由已经存在于体内的免疫细胞组成，这些细胞可以立即被招募到感染部位。天然免疫细胞包括粒细胞、巨噬细胞、肥大细胞、自然杀伤细胞（NK）和树突状细胞（DC）。中性粒细胞通常会触发对炎症的快速反应并分泌细胞因子，而成熟后，巨噬细胞可以分化为M1 和M2 极化细胞，进而抑制炎症反应和促进肿瘤的生长[24]，这同本研究中低风险人群在免疫浸润活动是一致的。

综上所述，炎症反应相关基因的预后特征是预测PAAD 患者预后的良好工具，但仍需要进一步的研究证实。