基于多平台数据识别胰腺导管腺癌预后相关肿瘤微环境基因

蒲垠全 马雨凡 彭莉 汤小伟 彭燕

西南医科大学附属医院消化内科，泸州 646000

PDAC是最常见的恶性肿瘤之一，其预后不良，目前居癌症死亡原因第7位，预计将在未来上升到第3位[1]。目前研究认为，胰腺周围复杂的区域解剖结构、高侵袭性和独特的肿瘤微环境是导致其预后不良的主要原因。尽管在过去的几十年里，PDAC在外科手术、放疗和化疗方面都取得了长足的进展，但其5年总生存率仍不到5%。免疫治疗的临床试验也未显示出对大多数PDAC患者的临床益处[2]，同时也缺乏可靠的生物学标志物来选择有效的分子亚型。肿瘤微环境通过对T细胞浸润形成物理和化学屏障、阻止药物释放等方式降低治疗效果[3]，PDAC独特的肿瘤微环境被认为是免疫治疗抵抗的主要原因之一[4]。本研究旨在探索PDAC相关肿瘤微环境，以确定影响患者预后的生物学标志物和潜在的免疫治疗靶点。

材料与方法

一、数据集的选择

从肿瘤基因组图谱(The Cancer Genome Atlas, TCGA)数据库下载142例PDAC患者的基因表达谱和临床随访数据，从基因表达综合(Gene Expression Omnibus，GEO)数据库的两个微阵列数据集(GSE21501和GSE62452)下载168例具有完整临床随访信息的PDAC患者的mRNA数据。

二、细胞类型富集分析及比较

xCell是一个用来对不同平台的转录组数据进行肿瘤微环境数字剖析的工具[5]。它整合了反卷积基因富集分析的优点，能根据基因表达谱计算64种细胞类型(包括免疫细胞和基质细胞)富集分数，其富集分数范围为-1～1。本研究运用xCell网络工具对基因表达数据进行细胞类型富集分析。TCGA数据库中的数据分析结果可以从xCell门户网站上直接下载。对于GEO数据则需要使用R/Bioconductor软件中的sva包来移除数据集中的批间差，然后再将规范化后的数据传输到xCell网站进行分析。以细胞富集评分中位数将TCGA的142例患者分为高分、低分2组，使用R/Bioconductor软件中的survival包进行单因素生存分析确定有预后价值的细胞类型，再使用GEO收集的数据集进行验证，P<0.05为差异有统计学意义。

三、基因表达差异分析及比较

四、差异表达基因的功能富集分析

为了进一步了解与预后相关DEGs的功能，使用R/Bioconductor软件中的clusterprofiler包对预后相关DEGs进行功能富集分析，包括基因本体(gene ontology, GO)分析和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes, KEGG) 通路富集分析，其中GO包括分子功能(molecular function, MF)、生物学过程(biological process, BP)和细胞成分(cellular component, CC)。FDR<0.05为差异有统计学意义。

五、蛋白质-蛋白质交互作用网络的建立

运用STRING(Search Tool for the Retrieval of Interacting proteins)数据库建立预后相关DEGs的蛋白质-蛋白质相互作用(protein-protein interaction, PPI)网络，然后通过Cytoscape软件进行PPI网络的重建，并使用Cytoscape MCODE插件查找集群，以定位密集连接的区域，同时计算网络各节点的连接度。

六、TISIDB数据库分析共同预后相关DEGs

TISIDB数据库是一个用于肿瘤和免疫系统相互作用的集成存储门户网站。本研究先确定TCGA数据集及GEO数据集两个平台数据的共同预后相关DEGs，再使用TISIDB数据库来评估共同预后相关DEGs在肿瘤与免疫系统的交互作用。

结 果

一、细胞类型富集评分及其与PDAC患者生存率的相关性

通过xCell网络工具对细胞类型富集评分，并根据富集评分中位数将PDAC患者分成高分、低分2组，结果显示，在64种细胞类型中只有Th1细胞和角质形成细胞在TCGA和GEO数据集中具有相同的预后价值(图1)。Th1细胞富集评分中位数在TCGA数据集和GEO数据集中分别为0.04155(0.00000～0.22100)、0.2183(0.00000～0.64560)，角质形成细胞富集评分中位数在TCGA数据集和GEO数据集中分别为0.12850(0.01030～0.39400)、0.2700(0.00000～0.50350)，Th1细胞和角质形成细胞高分组患者的总生存率中位数均显著低于低分组，差异均有统计学意义(P值均<0.05)。

二、与PDAC患者预后相关的DEGs筛选

TCGA数据库的142例患者Th1细胞富集高分组与低分组比较，高分组有379个基因上调，294个基因下调(图2A)；角质形成细胞富集高分组与低分组比较，高分组有481个基因上调，298个基因下调(图2B)。 通过单变量生存分析并做韦恩图(图2C)，共筛选出216个与PDAC患者预后相关的DEGs。

三、216个DEGs的功能富集分析

对216个DEGs进行了GO分析和KEGG通路富集分析，结果显示，以-log(矫正P值)>3设为阈值，21个BP、3个CC和5个MF被鉴定为差异有统计学意义(P值均<0.05)。21个BP中有9个(淋巴细胞分化、抗原受体介导的信号通路调控、淋巴细胞活化的调控、免疫系统过程的负调节、趋化因子介导的信号通路、B细胞活化、T细胞活化、对趋化因子应答、细胞对趋化因子的应答)与免疫过程密切相关(图3A)。5条KEGG通路中4条通路(病毒蛋白与细胞因子和细胞因子受体的相互作用、细胞因子-细胞因子受体相互作用、趋化因子信号通路、原发性免疫缺陷)与免疫过程密切相关(图3B)。

图1 从GEO数据集筛选的Th1细胞(1A)和角质形成细胞(1B )富集评分以及从TCGA数据集筛选的Th1细胞(1C)和角质形成细胞(1D)富集评分的患者生存曲线

图2 Th1细胞高低富集评分组(2A)和角质形成细胞高低富集评分组(2B)的DEGs火山图以及预后DEGs的韦恩图(2C)

四、216个DEGs的PPI网络和模块化分析

为了更好地理解216个DEGs之间的相互作用关系，使用STRING网络工具建立了一个包含158个节点和431个边的PPI网络(图4A)，并使用Cytoscape MCODE插件鉴定了一个MCODE得分>5分的模块(图4B)，结果显示，CCR7、CD27、CD5、CXCL13、ZAP70、MS4A1和CCL19基因的连通度显著高于该模块中的其他基因。

五、216个DEGs在GEO数据集中的验证

Kaplan-Meier单变量分析结果显示，15个基因在GEO和TCGA数据集中具有相似的预后价值，其中6个基因(CNTNAP4、DLGAP1、GIMAP7、LGI1、SFRP1、TPH1)的高表达与生存期延长有关，9个基因(FAM83A、HMGA2、ITGA3、ITGB6、KRT7、KRT13、PTGES、S100A2、UBE2C)高表达与生存期缩短有关(图5)。

图3 参与免疫相关的9个生物过程(3A)以及参与5条KEGG通路(3B)的基因

图4 216个与预后相关DEGs的PPI网络(4A)及MCODE为11.333分的模块图(4B)。从蓝色到红色的梯度代表基因从下调到上调的变化过程，节点大小代表指定蛋白的连通度

图5 TCGA(5A)和GEO(5B)数据集具有相似预后价值的15个基因的Kaplan-Meier生存曲线图

七、共同差异表达基因与免疫相关因素的相关性

在TISIDB数据库检索上述15个共同基因，结果发现GIMAP7与胰腺癌的免疫相关因素淋巴细胞、趋化因子受体、趋化因子、免疫抑制因子、免疫刺激因子、主要组织相容性复合体(major histocompatibility complex, MHC)均呈显著正相关(图6)。

讨 论

本研究通过TCGA数据库中142例PDAC患者和GEO的168例PDAC患者资料筛选出与预后相关的Th1细胞及角质形成细胞，依据富集评分高低分组比较来确定与PDAC预后相关的肿瘤微环境基因。通过DEGs分析和进一步的生存分析明确与PDAC预后有关的216个 DEGs。对这些DEGs的功能分析结果显示在生物学过程GO条目和KEGG通路中有近一半的GO条目和大多数KEGG通路与免疫反应密切相关。216个DEGs的PPI网络显示CCR7、CD27、CD5、CXCL13、ZAP70、MS4A1和CCL19具有更高的连通度，被鉴定为与PDAC密切相关的中枢基因。这7个基因参与了上述9个免疫相关的生物过程和4条免疫相关的KEGG通路，表明这些基因不仅与PDAC的预后密切相关，而且广泛参与了免疫调节。

上述7个基因在肿瘤分化、转移和预后中的作用均已有相关研究报道[5-11]。CCR7、CCL19、CD27和CXCL13是参与功能富集最多的4个基因，最有可能成为治疗靶点。趋化因子CCL19及其受体CCR7可趋化树突状细胞进入肿瘤组织[12]，诱导T细胞杀伤肿瘤细胞[13]及免疫细胞浸润肿瘤组织并释放细胞因子，还能通过ERK和PI3K/AKT信号转导上调Twist基因的表达，从而调控上皮-间充质转化(epithelial-mesenchymal transition, EMT)过程，抑制肿瘤的增殖、侵袭和转移。因此，CCR7/CCL19被广泛推荐为癌症治疗的前瞻性和实用性靶点[6]。CD27基因是肿瘤坏死因子超家族的一员，其与CD28协同作用，促进T细胞的初始激活和扩增，还能维持外周效应细胞的扩张和存活，并增强其细胞溶解和细胞因子功能[9]。近年来有报道表明，抗人CD27单克隆抗体与PD-1阻滞剂联合应用可提高肿瘤疫苗和小鼠免疫治疗的效果[14-15]。CXCL13基因通过其受体CXCR5发挥生物学功能，CXCL13/CXCR5轴通过整合多个信号事件包括PI3K/AKT途径的激活、基质金属蛋白酶表达的增加等来促进癌细胞的增殖、存活和迁移。CXCL13/CXCR5轴还可通过下调T效应细胞直接帮助肿瘤细胞逃避宿主免疫监视，及诱导IL-10途径间接帮助肿瘤细胞逃避T效应细胞免疫[16]。最近，针对CXCL13或CXCR5的敲除或抑制基因疗法、抑制性抗体和载药纳米粒在其他癌症中取得了一些进展[16-18]，然而针对上述基因的PDAC的治疗还鲜有报道。

本研究从GEO数据库下载的168例PDAC患者的基因表达谱数据对上述216个DEGs 进行交叉验证，最终确定了15个与预后相关的肿瘤微环境相关基因，包括CNTNAP4、DLGAP1、FAM83A、GIMAP7、HMGA2、ITGA3、ITGB6、KRT7、KRT13、LGI1、PTGES、S100A2、SFRP1、TPH1和UBE2C。这些基因已被多项研究证明与肿瘤的发生发展有关，或被用作预测不同肿瘤预后的生物学标志物。

EMT是一个由多种肿瘤微环境因子触发的生理过程，导致癌细胞从上皮到间充质表型发生形态学和遗传学改变，它是胰腺癌细胞高转移潜能和耐药性的基础[19]。先前研究表明，HMGA2、ITGA3、KRT13、S100A2和UBE2C可在许多癌症中诱发EMT，而SFRP1可抑制EMT[20-23]。本研究结果显示，SFRP1高表达与较长的生存期相关，而HMGA2、ITGA3、KRT13、S100A2和UBE2C基因高表达与较短的生存期相关，该结果与它们在EMT中的作用机制是一致的。

本研究运用TISIDB数据库分析上述15个共同DEGs与免疫系统的关系，显示GIMAP7与多种免疫相关因子包括淋巴细胞、趋化因子、趋化因子受体、免疫抑制因子、免疫刺激因子和MHC分子存在更高的相关性。进一步观察发现，GIMAP7还与鉴定的中枢基因包括CCR7、CCL19、CXCL13和CD27也密切相关。GIMAP7蛋白是一种GTP酶，属于免疫相关蛋白(group of immunity-associated proteins，GIMAP)家族，定位于内质网和高尔基体。GIMAP家族蛋白在调节T细胞的发育、选择、存活、自身免疫和内环境稳定方面起着关键作用[24-25]。Schwefel等[25]报道，GIMAP7可刺激GIMAP2的GTP水解。Dion等[26]报道，在不同胸腺亚群中，双阴性(double negative，DN)和双阳性(double positive，DP)个体之间的大鼠GIMAP7表达显著增加，并且在从CD4+CD8+DP胸腺细胞过渡到CD4+或CD8+单阳性胸腺细胞期间增加。但GIMAP7参与免疫调节的机制尚不清楚，需要进一步研究。

总之，本研究鉴定了一组与PDAC预后相关的肿瘤微环境基因及其7个中枢基因，并提供了CCR7、CCL19、CD27、CXCL13和GIMAP7 5个新的PDAC免疫治疗潜在靶点。

利益冲突所有作者均声明不存在利益冲突