胰腺癌是一种恶性程度极高的肿瘤，在全球范围内死亡率极高。本研究的目的是分析胰腺癌的肿瘤免疫微环境，确定预后相关免疫生物标志物和潜在靶点。
利用癌症基因组图谱(The Cancer Genome Atlas，TCGA)胰腺癌数据库和基因表达综合数据库(Gene Expression Omnibus，GEO数据库)，进行单样本基因集富集分析以评估肿瘤组织浸润24种免疫细胞类型的比例。通过最小绝对收缩和选择算子分析确定关键免疫细胞类型，并用于建立风险评分模型(risk scoring model，RSM)。
RSM在TCGA队列中被验证为一种强有力的总体生存预后工具(P<0.05)，1年、3年和5年总生存期的曲线下面积分别为0.789、0.764和0.887。此外，基于患者临床信息构建的列线图生成了拟合良好的校准曲线。通过分析免疫相关基因和肿瘤突变负荷(tumor mutational burden，TMB)，发现胰腺癌TMB与RSM评分高度相关，RSM评分可以在一定程度上对免疫治疗的预后进行预测。RSM的可靠性通过GEO数据库中的GSE62452和GSE28735数据集成功验证。
通过单样本基因集富集方法可预测胰腺癌浸润免疫细胞图谱，基于浸润免疫细胞建立了一个有效评估胰腺癌预后和免疫治疗反应的风险评分模型。
本刊版权归人民卫生出版社有限公司所有。任何机构或个人欲转发本刊图文，请与本刊联系。凡未经授权而转载、摘编本刊文章者，本刊将保留追究其法律责任的权利。
由于早期诊断的困难和目前晚期癌症治疗的局限性，胰腺癌是全球癌症相关死亡的主要原因之一[1]。鉴于这种不良预后，鉴定用于评估预后的有效生物标志物至关重要。肿瘤免疫微环境已被证明在胰腺癌的发展中很重要[2,3]，因此，肿瘤免疫微环境可能为预后分析和提高风险评分模型(risk score models，RSM)[4]的精度和可靠性提供潜在的生物标志物。然而，由于其动态和复杂的性质，我们对肿瘤免疫微环境的理解仍然不足。在肿瘤进展中起重要作用的肿瘤浸润性免疫细胞(tumor-infiltrating immune cells，TIC)的特征是器官特异性的，并且根据肿瘤阶段而有所不同。最近有学者分析了多种癌症类型中TIC与预后的相关性[4,5]，但无法获得胰腺癌的全面免疫概况。
TIC的子集通常通过免疫组织化学或流式细胞术测量。这些技术取决于表型标记的可用性，并且实验程序的质量控制和标准化具有挑战性[6]。随着生物信息学的发展，现在可以通过RNA表达数据挖掘来表征TIC的组成。单样本基因集富集分析(gene set enrichment analysis，GSEA)是一种根据RNA测序(RNA-seq)表达矩阵分析免疫细胞组分的新方法，比其他基于RNA-seq的方法更准确[7,8]。使用GSEA，可以通过标记基因的特征来确定异质样品中TIC的多样性，该方法已在泛癌免疫基因组分析中得到验证[9]。
TIC的表征对癌症的治疗具有重要意义。免疫疗法已被公认为是多种晚期和侵袭性癌症的有效治疗方法[10,11]，并且免疫疗法的功效与多种癌症中的特定基因和肿瘤突变负荷(tumor mutational burden，TMB)密切相关[12,13]。因此，了解肿瘤微环境中免疫治疗反应性的分子机制非常重要。
本研究根据胰腺腺癌(pancreatic adenocarcinoma，PAAD)的RNA-seq数据进行了单样本GSEA(single-sample gene set enrichment analysis，ssGSEA)分析，以评估24种免疫细胞类型的比例。进行LASSO-Cox分析以识别关键的免疫预测因子并构建RSM，该RSM可以作为评估胰腺癌患者生存率的有力工具。此外，RSM与免疫相关基因和TMB的表达高度相关。综上，本研究开发的RSM可用于评估胰腺癌治疗中的预后和对免疫治疗的反应。
RNA-seq数据的表达谱、体细胞突变数据和相应的临床信息从癌症基因组图谱网站(The Cancer Genome Atlas，TCGA)(https：//tcga-data.nci.nih.gov/tcga/)下载。GEO数据集GSE62452[14]和GSE28735[15]来自NCBI基因表达综合数据库(Gene Expression Omnibus，GEO数据库)(https：//www.ncbi.nlm.nih.gov/gds/)。
为了量化免疫细胞的比例，TIC的基因列表从Bindea等[7]的补充表S1下载。使用R软件[16]的GSVA包分析每百万转录本(transcripts per million，TPM)标准化表达数据。研究了24种浸润性免疫细胞：活化树突细胞(activate dendritic cell，aDC)、B细胞、CD8+T细胞、细胞毒性细胞、树突细胞(dendritic cell，DC)、嗜酸性粒细胞、未成熟树突细胞(immature dendritic cell，iDC)、巨噬细胞、肥大细胞、中性粒细胞、自然杀伤(natural killer，NK)CD56高表达细胞、NK CD56低表达细胞、NK细胞、浆细胞样树突状细胞(plasmacytoid dendritic cell，pDC)、T细胞、T辅助细胞、中央记忆T细胞(central memory T cell，Tcm)、效应记忆T细胞(effector memory T cell，Tem)、T滤泡辅助细胞(T follicular helper cell，TFH)、T gamma delta(Tgd)细胞、Th1细胞、Th17细胞、Th2细胞和调节性T细胞(regulatory T cell，TReg)。使用Wilcoxon检验评估每个样品的不同免疫细胞类型。
将PAAD患者的总生存(overall survival，OS)期作为生存分析的终点。根据纳入标准收集样本：①原发性肿瘤组织；②完整的OS信息可用性。然后进行LASSO分析以识别TIC的预测性生物标志物[17]。RSM基于相应的系数构建。进行了Kaplan-Meier生存曲线和受试者操作特征曲线(receiver operating characteristic curve，ROC曲线)分析以评估RSM。最后，我们基于风险模型和患者的临床信息构建了预后列线图，为预后评估提供了有用的工具。建立校准曲线以评估预期和观察到的存活率之间的差异。
TMB定义为每百万个碱基中检测到的体细胞基因编码错误总数。我们使用R软件中的"maftools"包和默认参数计算了TMB。
我们根据以下标准选择验证队列[18]：①患者诊断为PAAD；②样本量大于100例；③生存信息的可用性。基于这些标准，选择GSE62452和GSE28735数据集来验证我们的RSM。
从TCGA队列中提取了177例PAAD患者的数据(表1)。根据图1所示的步骤分析数据。我们最初通过ssGSEA分析了177个癌组织和4个癌旁组织中的免疫细胞组分。与CD8+T细胞、细胞毒性细胞、中性粒细胞、Tem细胞和Tgd细胞相对应的基因特征在邻近组织中的表达比在癌组织中更高。相比之下，与癌组织相比，癌旁组织中与NK CD56高表达细胞和Th2细胞相对应的特有基因的表达显著降低(图2A，表2)。免疫细胞类型的差异值如图2B所示。为了验证这些差异，分析了包括61个癌旁组织和69个癌组织在内的GSE62452队列。细胞毒性细胞、嗜酸性粒细胞、TFH细胞和Th17细胞的基因特征在癌旁组织中的表达显著高于癌组织，而对aDC、DC、巨噬细胞、肥大细胞、NK CD56低表达细胞、T辅助细胞、Th1细胞和Th2细胞癌组织表达显著高于癌旁组织(图2C)。免疫细胞类型的差异值如图2D所示。与之前的研究一致[20]，TCGA和GEO队列之间的免疫细胞组成不同，这可能归因于样本量的差异，因为PAAD队列中仅包括4个对照样本。
注：A，TCGA队列中癌旁组织和癌组织中表达水平的小提琴图；B，TCGA队列中癌旁组织和癌组织之间组成的差异值；C，GSE62452中癌旁组织和癌组织中表达水平的小提琴图；D，GSE62452中癌旁组织和癌组织之间组成的差异值。Wilcoxon检验用于计算P值；NK，自然杀伤。
胰腺腺癌患者的临床资料
胰腺腺癌患者的临床资料
| 参数 | 计数(例) | 比例(%) | |
|---|---|---|---|
| 年龄 | |||
| >65岁 | 84 | 47.46 | |
| ≤65岁 | 93 | 52.54 | |
| 性别 | |||
| 男性 | 97 | 54.80 | |
| 女性 | 80 | 45.20 | |
| 种族 | |||
| 亚洲人 | 11 | 6.21 | |
| 黑人或非裔美国人 | 6 | 3.39 | |
| 白种人 | 156 | 88.14 | |
| 无记录 | 4 | 2.26 | |
| 病理分期 | |||
| Ⅰ期 | 21 | 11.86 | |
| Ⅱ期 | 146 | 82.49 | |
| Ⅲ期 | 3 | 1.69 | |
| Ⅳ期 | 4 | 2.26 | |
| 无记录 | 3 | 1.69 | |
| 病理分期T | |||
| T1 | 7 | 3.95 | |
| T2 | 24 | 13.56 | |
| T3 | 141 | 79.66 | |
| T4 | 3 | 1.69 | |
| TX* | 2 | 1.13 | |
| 病理分期N | |||
| N0 | 49 | 27.68 | |
| N1 | 123 | 69.49 | |
| NX* | 5 | 2.82 | |
| 病理分期M | |||
| M0 | 79 | 44.63 | |
| M1 | 4 | 2.26 | |
| MX* | 94 | 53.11 | |
注：*，表示无相关记录数据。
Wilcoxon检测胰腺腺癌中癌旁组织和癌组织的免疫细胞类型结果
Wilcoxon检测胰腺腺癌中癌旁组织和癌组织的免疫细胞类型结果
| 细胞类型 | TCGA Wilcoxon检测P值 | GSE62452 Wilcoxon检测P值 |
|---|---|---|
| 活化树突细胞 | 0.337 | <0.001 |
| B细胞 | 0.141 | 0.068 |
| CD8+T细胞 | 0.022 | 0.083 |
| 细胞毒性细胞 | 0.028 | 0.005 |
| 树突细胞 | 0.858 | 0.001 |
| 嗜酸性粒细胞 | 0.717 | 0.001 |
| 未成熟树突细胞 | 0.204 | 0.794 |
| 巨噬细胞 | 0.265 | <0.001 |
| 肥大细胞 | 0.172 | 0.006 |
| 中性粒细胞 | 0.025 | 0.061 |
| NK CD56高表达细胞 | 0.011 | 0.056 |
| NK CD56低表达细胞 | 0.065 | 0.018 |
| NK细胞 | 0.141 | 0.063 |
| 浆细胞样树突状细胞 | 0.139 | 0.547 |
| T细胞 | 0.063 | 0.249 |
| T辅助细胞 | 0.703 | 0.010 |
| 中央记忆T细胞 | 0.050 | 0.441 |
| 效应记忆T细胞 | 0.005 | 0.341 |
| T滤泡辅助细胞 | 0.163 | <0.001 |
| Tgd细胞 | 0.044 | 0.063 |
| Th1细胞 | 0.965 | <0.001 |
| Th17细胞 | 0.572 | <0.001 |
| Th2细胞 | 0.012 | <0.001 |
| 调节性T细胞 | 0.100 | 0.532 |
我们对入选肿瘤样本进行生存分析。免疫细胞组成的热图如图3A所示。构建LASSO回归模型，并将7种免疫细胞类型确定为潜在预测因子(图3A~C)。这7种免疫细胞类型用于根据以下公式构建Cox模型：RSM风险评分＝(aDC细胞水平×1.844)+(巨噬细胞水平×2.972)+(NK CD56低表达细胞水平×0.770)+[Tem水平×(－8.336)]+[TFH水平×(－2.086)]+(Th2细胞水平×8.625)+[TReg水平×(－1.582)]。根据个体预后风险评分，PAAD患者以中位数为临界值分为低风险组和高风险组。Kaplan-Meier曲线证实，高风险组的OS显著较差(P<0.001，图4A)。此外，风险评分可以估计患者长期的OS，因为1年、3年和5年OS的曲线下面积(area under the curve，AUC)分别为0.789、0.764和0.887(图4B)。图4C给出了患者风险评分的分布。
注：A，PAAD癌组织中浸润免疫细胞的热图；B，LASSO回归的参数选择(虚线之间的区域表示适当的值)；C，LASSO回归中每种浸润性免疫细胞类型的系数图；PAAD，胰腺腺癌；L1范数，面板中每个元素的绝对值之和。
注：A，总生存期的Kaplan-Meier生存曲线；B，风险评分模型的受试者操作特征曲线分析；C，胰腺腺癌患者的风险评分、生存状态和基因表达水平的分布；AUC，曲线下面积；ROC曲线，受试者操作特征曲线。
为了为PAAD患者预后的临床评估提供有用的工具，我们构建了基于风险评分和临床信息的预后列线图，包括年龄、性别、TNM分期、分级和新事件(图5A)。随后的校准曲线分析用于测试列线图的可靠性，在TCGA队列中的预测和观察之间产生了良好的拟合(图5B)。
注：A，用于预测TCGA数据集中胰腺腺癌患者1年、3年和5年生存概率的列线图；B，列线图的校准曲线；C，胰腺腺癌中的基因突变图谱；D，TMB分析与风险评分的关系；TCGA，癌症基因组图谱。
我们使用"maftools"包可视化PAAD中的突变图谱(图5C)。综合分析显示TMB分值与风险评分呈正相关(图5D)，表明RSM模型中的高风险评分患者，理论上对免疫治疗有更高敏感性。提示我们，在RSM高风险评分患者中，经传统治疗后预后较差的PAAD患者更有可能对免疫治疗产生更好的响应。
我们分析了风险评分与先前研究中7个免疫治疗相关标记基因(CD274、CTLA4、LAG3、LIF、PDCD1、PGAM1、TIGIT)的相关性。然而，这些标记基因的表达水平与我们的RSM无关，并且对PAAD队列中的患者生存没有显著影响(图6)。
注：左下角的数字代表皮尔逊相关系数，P>0.05用"×"表示；risk score，风险评分。
我们使用单变量和多变量Cox分析来确认风险评分作为独立于其他临床参数的临床指标的能力(表3)。使用来自GEO数据库的GSE62452和GSE28735数据集对我们构建的RSM模型进行了验证。与TCGA队列的结果一致，验证数据集中风险评分高的患者预后较差(P<0.001，图7A、图7C)。此外，验证数据集中基于RSM的患者3年和5年OS的AUC具有显著性(图7B、图7D)。
注：A为GES62452的Kaplan-Meier风险生存分析；B为GES62452的受试者操作特征曲线；C为GSE78299的Kaplan-Meier风险生存分析；D为GSE78299的受试者操作特征曲线。
相关临床特征和风险评分的预测值
相关临床特征和风险评分的预测值
| 特征 | 单因素Cox分析 | 多因素Cox分析 | ||
|---|---|---|---|---|
| HR（95%CI） | P值 | HR（95%CI） | P值 | |
| 年龄 | 1.027(1.005~1.049) | 0.018 | 1.017(0.995~1.039) | 0.140 |
| 性别(女＝0，男＝1) | 0.781(0.514~1.187) | 0.247 | ||
| 病理分期T | 1.490(0.924~2.405) | 0.102 | ||
| 病理分期N | 2.180(1.283~3.706) | 0.004 | 1.931(1.129~3.303) | 0.016 |
| 病理分期M | 0.992(0.243~4.050) | 0.991 | ||
| 病理分期 | 1.233(0.820~1.855) | 0.314 | ||
| 分级 | 1.338(0.998~1.795) | 0.052 | ||
| 新事件 | 2.237(1.330~3.765) | 0.002 | 1.742(1.029~2.950) | 0.039 |
| 风险评分 | 4.787(3.051~7.511) | <0.001 | 4.152(2.650~6.507) | <0.001 |
TIC与癌症的发生和发展密切相关[21]。最近的研究表明，肿瘤进展是由肿瘤免疫状态和宿主免疫反应之间的不平衡引起的[22]。因此，了解肿瘤微环境中的免疫状态并研究其在肿瘤进展中的功能和机制至关重要。
目前，已经建立了几种基于免疫组织化学的免疫评分模型来评估癌症患者的预后[23,24]。然而，这些免疫组织化学数据的信息价值受限于少数细胞类型。在本研究中，我们使用ssGSEA综合分析了PAAD中24种TIC的分布。与癌旁组织相比，癌组织中的免疫细胞组成存在显著差异。NK CD56高表达细胞和Th2细胞特异性基因在肿瘤组织中的表达较高，而CD8+T细胞、细胞毒性细胞、中性粒细胞、Tem细胞和Tgd细胞的基因在邻近组织中的表达较高。使用LASSO-Cox方法分析这些TIC数据，以选择关键免疫生物标志物来构建RSM。由此产生的RSM由7个成分组成：Tem、TFH和Treg细胞，它们是与预后良好相关的保护因素，以及aDC细胞、巨噬细胞、NK CD56低表达细胞和Th2细胞，它们与预后不良有关。
在RSM中的保护因子中，TFH细胞是T辅助细胞的一个子集，可调节抗原依赖性B细胞免疫的发展，据报道在乳腺癌中发挥防御作用[25]。Tem细胞是肿瘤免疫中重要的效应细胞类型，与卵巢癌患者的长期无复发生存期有关[26]。关于Treg细胞在癌症预后中的作用存在互相矛盾的报道：Miracoo等[27]报道，Treg细胞的高浸润与黑色素瘤的不良预后有关，而Treg浸润与头颈癌的良好预后有关[6]。在与预后不良相关的因素中，Th2细胞与肿瘤侵袭性有关[28]。aDC和NK CD56低表达细胞均在抗肿瘤途径中发挥重要作用，许多研究表明这些细胞类型与长期存活有关[29]。然而，由于肿瘤环境中的免疫细胞动力学，效应免疫细胞可能在恶性转化过程中反应更强烈，导致在某些情况下较高的表达与较差的存活率相关[30]。巨噬细胞在肿瘤进展中也起着关键作用，不同的巨噬细胞亚型(M0、M1和M2)因肿瘤分期而异，并在癌症进展中发挥不同的免疫功能[31,32]。由于标记基因的可用性有限，本研究无法分析巨噬细胞亚群，在未来的工作中需要对巨噬细胞进行更详细的研究。
在本研究中，RSM的高风险评分与TCGA队列中PAAD患者的不良预后显著相关。为了更好地理解RSM，我们根据风险评分、年龄、性别、TNM阶段、等级和新事件构建了列线图。校准曲线显示列线图预测和观察到的OS之间的拟合良好。RSM还使用从GEO数据库下载的2个独立数据集GSE62452和GSE28735进行了验证。Kaplan-Meier风险分析证实了RSM在这2个数据集上的良好性能。3年和5年OS曲线的AUC值显著；1年OS缺乏显著结果可能归因于数据集大小和数据集之间的患者平均寿命的差异。
