胰腺癌基质细胞评分相关疾病预后特征的研究

Study on the Prognostic Characteristics of Pancreatic Cancer Related to Stromal Cell Score

摘要: 胰腺癌是预后最差的癌症之一。改善胰腺癌患者的预后是当下研究的目标。在本研究中，我们使用“estimate”软件包分析了胰腺癌的表达数据，计算了肿瘤微环境中基质细胞浸润的含量，并根据基质细胞得分的差异总结获得了差异表达基因(differentially expressed genes, DEGs)。这些DEGs被用于构建胰腺癌预后模型。本研究发现，模型中的风险分数能很好地区分不同特征胰腺癌患者的预后。风险分数与一些免疫细胞和免疫功能密切相关。我们的模型在胰腺癌研究中具有广阔的应用前景。

Abstract: Pancreatic cancer is one of the worst prognosis cancers. It is the goal of researchers to improve the prognosis of patients with pancreatic cancer. In our study, we used “estimate” package to analyze the expression data of pancreatic adenocarcinoma, calculated the content of stromal cell infiltration in the tumor microenvironment, and obtained the differentially expressed genes according to the difference of stromal cell score. These differentially expressed genes are used to construct a prognostic model for pancreatic cancer. We found that the riskScore in our model can well distinguish the prognosis of patients with different characteristics. The riskScore was closely related to some immune cells and immune function. Our model has broad application prospects in the study of pancreatic cancer.

1. 引言

胰腺癌的发病率在癌症中并不算高，远低于肺癌和肝癌 [1] 。然而，由于胰腺癌早期症状不明显，没有准确的诊断标志物，往往导致胰腺癌被发现时已处于中晚期，癌细胞恶性程度高，使胰腺癌患者的预后较其他肿瘤差 [2] [3] 。因此，寻找与胰腺癌相关的预后标志物尤为重要。

胰腺癌的生物学行为在很大程度上取决于肿瘤细胞与周围微环境之间的相互作用 [4] [5] 。胰腺癌中存在大量细胞外基质致密性沉积，使胰腺癌成为发病率最高的实体瘤之一。以往的研究侧重于胰腺肿瘤细胞的治疗，但结果并不乐观 [6] 。近年来，人们发现肿瘤微环境在调控肿瘤组织的发展中发挥着重要作用。因此，肿瘤微环境在胰腺癌中被广泛研究。Hegde等人发现，树突状细胞缺乏会导致胰腺癌患者体内免疫监视功能失调 [7] 。Yang等人发现，胰腺癌的神经周围侵袭对免疫微环境进行了重编程 [8] 。Choueiry等人发现，CD200促进了胰腺癌微环境中的免疫抑制过程 [9] 。

在我们的研究中，我们使用“estimate”软件包对胰腺癌基质细胞的含量进行评分，并发现了高分组和低分组之间的DEGs。我们利用这些基因构建了一个预后模型，并在外部数据集上进行了验证。我们发现，我们的模型与许多指标高度相关。希望我们的研究能为胰腺癌的治疗提供新思路。

2. 材料与方法

2.1. 数据下载

我们从癌症基因组图谱(The Cancer Genome Atlas, TCGA, https://tcga-data.nci.nih.gov/tcga/)中下载了胰腺导管腺癌(Pancreatic ductal adenocarcinoma, PDAC)项目的表达数据，并将其校正为片段每千碱基百万片段(Fragments perkilobase million, FPKM)，作为训练队列和PAAD的临床数据。GSE57495的表达数据和临床数据从基因表达总库(GEO, https://www.ncbi.nlm.nih.gov/geo/)下载。我们通过从Ensembl (http://asia.ensembl.org)获得的基因转移格式(gene transfer format, GTF)文件对数据进行了注释。

2.2. 模型的构建与验证

用R软件(4.0.0)中的“limma”软件包(https://bioconductor.org/packages/limma/)筛选DEGs。使用“survival”(https://cran.r-project.org/package=survival)、“glmnet”(https://cran.r-project.org/package=glmnet)和“survminer”(https://cran.r-project.org/package=survminer)软件包对数据进行Cox危险分析和Lasso回归分析。“survival ROC”软件包用于绘制接收者操作特征曲线，“survival”软件包用于绘制生存曲线。

2.3. 基因组富集分析(Gene Set Enrichment Analysis, GSEA)

本研究利用GSEA来发现TCGA队列中不同风险组之间的差异。研究选择了一个注释基因组文件(c2.cp.kegg.v7.0.symbols.gmt)作为参考。确认阈值为FDR q-val < 0.05。

2.4. 免疫分析

通过CIBERSORT分析，TCGA PAAD患者的免疫浸润解旋结果显著。TCGA队列中每个样本的“StromalScore”、“ImmuneScore”和“ESTIMATEScore”均由“estimate”软件包完成。“GSVA”和“GSEABase”软件包用于对每位患者进行ssGSEA分析。各指标的相关性分析通过spearman检验完成。

3. 结果

3.1. 在TCGA队列中构建预后模型

对不同PAAD患者的免疫微环境进行评分后，我们将StromalScore从低到高排序。对第一季度(43例)和最后一季度(44例)患者的DEGs进行分析，共发现2322个DEGs。结合临床预后，我们在TCGA队列中通过单变量Cox危险分析筛选出256个基因。我们使用Lasso回归和多变量Cox危险分析来缩小基因数量，最终得到14个基因来优化模型(图1)，并计算出每个样本的riskScore (riskScore = PIP5K1B* − 0.1440 + B3GNT3*0.0154 + PDE4C*0.1717 + IDO1*0.0339 + SPTBN2*0.0767 + NR2F1* − 0.1421 + TLE6* − 0.6760 + BICC1*0.1070 + VASH1* − 0.1907 + LY6G6C*0.3764 + ASPHD1* − 0.0326 + LLGL2* − 0.0458 + FOXA1* − 0.1337 + MYEOV*0.0267)，并用风险分数的中位数来区分高、低风险组。2、3、5年的ROC曲线下面积AUC值分别为0.840、0.893、0.961 (图2A)。高危组和低危组的预后存在明显差异(图2B)。热图显示，PIP5K1B、B3GNT3、PDE4C、IDO1、SPTBN2、NR2F1 TLE6、BICC1、VASH1、LY6G6C、ASPHD1、LLGL2、FOXA1、MYEOV在不同风险组的表达水平不同(图3)，肝细胞癌(Hepatocellular carcinoma, HCC)患者的死亡风险随riskScore的增加而增加(图4)。

Figure 1. Multivariate Cox hazard analysis results

图1. 多变量Cox危险分析结果

Figure 2. (A) Comparison of survival status between high-risk and low-risk groups. (B) ROC curves for different years in the TCGA cohort

图2. (A) 高危组和低危组的生存状况比较。(B) TCGA队列中不同年份的ROC曲线

Figure 3. Expression levels of 14 genes in different groups

图3. 不同组别中14个基因的表达水平

Figure 4. Survival status of patients with different risk scores

图4. 不同风险分数患者的生存状况

Figure 5. Confirmatory prognostic model (A) ROC curves for different years in the GSE 57495 group. (B) Comparison of the survival status of different groups in the GSE 57495 cohort

图5. 验证预后模型(A) GSE 57495组不同年份的ROC曲线。(B) GSE 57495队列中不同组别的生存状况比较

Figure 6. Relationship between risk scores and clinicopathological features. Univariate Cox risk factor analysis (A) and multivariate Cox risk analysis (B) for patient characteristics. (C) Relationship between quantity and riskScore

图6. 风险分数与临床病理特征之间的关系。患者特征的单变量Cox危险因素分析(A)和多变量Cox危险分析(B)。(C)数量和riskScore的关系

3.2. 验证预后模型

我们在GSE57495队列中验证了该模型。在2、3、5年中，ROC曲线下的面积AUC值分别为0.661、0.709、0.768 (图5A)。该模型能明显区分GSE57495队列中高风险组和低风险组患者的预后(图5B)。

3.3. 风险分数作为独立的预后指标

我们分析了riskScore与临床病理特征(年龄、性别、组织学分级、临床分期和TNM)之间的关系。临床病理特征的单变量Cox危险分析表明，年龄、N和riskScore的p值均小于0.05，危险比大于1 (图6A)。临床病理特征的多变量Cox危险分析表明，N和riskScore的p值均小于0.05，危险比大于1 (图6B)。我们绘制了N和riskScore的提名图(图6C)。不同等级、M、N和T组的riskScore有显著差异(图7)。不同年龄、性别、组织学分级、M0、N0、分期、T组中不同riskScore组的预后有明显差异(图8)。

3.4. 不同风险分数组的GSEA

高风险组发现0个基因组(NOM p-val)。在低风险组中，我们发现了3个基因组(PHOSPHATIDYLINOSITOL_ SIGNALING_SYSTEM、CALCIUM_SIGNALING_PATHWAY、NEUROACTIVE_LIGAND_RECEPTOR_ INTERACTION) (NOM p-val < 0.05) (图9)。

Figure 7. (A) The distribution of risk scores for different grades. (B) The distribution of risk scores for different M. (C) The distribution of risk scores for different N. (D) The distribution of risk scores for different T

图7. (A) 不同等级的风险分数分布。(B) 不同M的风险分数分布。(C) 不同N的风险分数分布。(D) 不同T的风险分数分布

Figure 8. Risk scores can predict survival of patients with different age, sex, histological grade, M0, N0, stage, and T

图8. 风险分数预测不同年龄、性别、组织学分级、M0、N0、分期、T患者的生存率

Figure 9. The genome of the low-risk group

图9. 低风险组的基因组

3.5. 风险分数与免疫

我们发现，StromalScore可以很好地区分不同的风险分数组。在不同的风险分数组中，一些免疫细胞的含量也存在明显差异(图10)。风险分数与T细胞CD4记忆激活、巨噬细胞M1存在明显的正相关(图11)。ssGSEA分析显示，高风险组和低风险组的B细胞、pDCs、Trh和TIL的细胞含量存在显著差异，细胞溶解活性、T细胞成本刺激、II型IFN应答在两组间存在显著差异(图12)。我们还发现，风险分数与免疫检查点之间的相关性并不明显。

Figure 10. Differences in immune cells of different risk scores

图10. 不同风险分数组免疫细胞的差异

Figure 11. Association of immune cells with risk scores

图11. 免疫细胞与风险分数的相关性

Figure 12. ssGSEA analysis of immune cells and immune function in different risk scores

图12. 不同风险分数组免疫细胞和免疫功能的ssGSEA分析

4. 讨论

如何改善胰腺癌的预后一直困扰着研究人员。但近年来基因组学的飞速发展加深了人们对胰腺癌的认识。基因组测序信息被广泛应用于胰腺癌的诊断和治疗。我们的研究基于转录组数据描述了胰腺癌免疫微环境中浸润的基质细胞数量，分析了高分组和低分组胰腺癌基因表达的差异，并构建了预后模型。

构成本研究模型的14个基因在肿瘤中发挥着重要作用。Xiong等人发现，PIP5K1B是头颈部鳞状细胞癌的预后相关基因 [10] 。Ferrucci等发现PIP5K1B在三阴性乳腺癌中驱动肿瘤相关巨噬细胞的极化 [11] 。Li等人发现，B3GNT3在靶向糖基化PD-L1根治三阴性乳腺癌中具有重要作用 [12] 。Sun等人发现B3GNT3会促进肺癌的发展，B3GNT3的敲低可以抑制肺癌细胞的生长和体外侵袭 [13] 。Samuel等人比较了研究各PDE4亚型在不同恶性肿瘤中的功能作用的证据 [14] 。Bao等发现PDE4C是人类胶质瘤的预后生物标志物 [15] 。Zhai等人发现，IDO1靶向策略可提高癌症患者免疫疗法的疗效 [16] 。Cheong等发现IDO1/TDO2-KYN-AhR通路在癌症免疫治疗中非常重要 [17] 。Yang等发现SPTBN2与结直肠癌患者的预后有关 [18] 。Zhang等发现NR2F1促进乳腺癌血管生成 [19] 。Zhang等发现TLE6是结直肠癌中与m6A相关的生物标志物 [20] 。Zhao等发现BICC1是胃癌的新型预后生物标志物 [21] 。Wang等发现VASH1可调控肺癌的脑转移 [22] 。Mette发现LY6G6C可能是溃疡性粘膜炎的潜在保护性生物标志物 [23] 。Saito等人发现，LLGL2能缓解ER(+)乳腺癌的营养压力 [24] 。Adams等发现FOXA1基因突变会改变前列腺癌的表型 [25] 。Liang等人发现，MYEOV会促进胰腺癌的进展 [26] 。

经过外部数据集的验证，我们研究构建的预后模型具有很高的可信度。该模型生成的风险评分能显著区分不同患者组的预后，并与一些免疫细胞和免疫功能有很好的相关性。然而，我们的模型还需要更多的生物学实验和更大规模的多中心数据验证。