基于生物信息学筛选的胰腺癌差异基因及免疫浸润机制预测靶向中药

龙思丹，季双双，陈美池，贠张君，马银杰，薛 鹏，朱世杰*

1中国中医科学院望京医院肿瘤科，北京 100102；2北京中医药大学研究生院，北京 100029

胰腺导管腺癌(pancreatic ductal adenocarcinoma，PDAC)的全球负担在过去25年中增加了一倍以上，预计2030年，PDAC将成为癌症相关死亡的第2大原因[1]。PDAC患者五年的平均生存率仅8%，并且高达71%患者即使手术完全切除后也面临复发的风险[2]。这种高致死率和高复发率的肿瘤特性是由PDAC紧密的基质成分(物理屏障)决定的[3]，其防止免疫效应细胞浸润到肿瘤中，使胰腺癌细胞能够逃避免疫监视。目前，许多临床试验尝试通过免疫治疗改善PDAC患者的整体疗效，包括免疫检查点抑制剂、癌症疫苗、过继细胞转移等[4]，单抗治疗的临床效果不佳，仍未取得突破性进展。

随着高通量遗传分析的出现，基因表达谱分析成为了鉴定各种疾病差异表达基因的有效方法。CIBERSORT是一种分析工具，它使用微阵列数据或RNA测序数据来评估样本中免疫细胞的表达，并获得各种免疫细胞比率[5]。虽然目前已有多篇生物信息学针对PDAC差异基因进行了分析[6,7]，但缺乏进一步对于免疫机制和治疗靶向预测的深入研究，尤其是在植物药的治疗潜力方面，中医已对于PDAC常伴见的黄疸、腹水等病症已形成系统认识，同时中药以多靶点的天然优势能有效干预复杂的免疫微环境，值得深入挖掘。因此，本研究利用基因芯片数据库(gene expression omnibus，GEO)，选择PDAC患者配对组织样本进行相关分析，从基因水平挖掘疾病机制，进行生存差异和免疫相关分析，通过CIBERSORT反卷积算法描绘PDAC组织中的免疫细胞浸润模式，最后综合预测和筛选潜在的有效中药，为临床药物的选择、研发提供依据。

1 资料与方法

1.1 差异基因筛选

以“Pancreatic Cancer”为关键词检索GEO基因表达综合数据库，筛选条件：(1)“Expression profiling by array”(2)“Homo sapien”(3)PDAC和癌旁组织。(4)组织来自同一患者的成对匹配样本。(5)全基因表达谱。(6)数据详细，包括Gene Symbol等信息。最终筛选得到GSE15471基因芯片为本次研究分析的对象。GSE15471数据集的处理过程：在Affymetrix U133 plus 2.0全基因组微阵列上分析，并且进行了重复的微阵列杂交，以测量技术测量误差，因此总共进行了78次基因芯片杂交。其中一个样本不符合质量控制被排除，随后使用RMA算法对微阵列数据进行标准化。应用GEO2R(http://www.ncbi.nlm.nih.gov/geo/geo2r)在线工具及Excel软件筛选芯片胰腺癌组织与正常组织的差异表达基因，筛选条件|log2FC| ≥ 1和校正P< 0.05进行差异基因分析。

1.2 蛋白网络互作及功能注释

核心差异基因通过STRING数据库获取蛋白网络互作(protein-protein interaction，PPI)网络模型，并将结果导入Cytoscape软件中，采用Cytoscape插件工具cyto Hubba中的MCC方法识别重要基因，将重要基因进行可视化运用R语言程序中的“cluster Profiler”插件包进行差异基因的基因本体论(gene ontology，GO)与京都基因和基因组百科全书(kyoto encyclopedia of genesand genomes，KEGG)分析，包括生物过程(biological process，BP)、细胞组分(cellular component，CC)、分子功能(molecular function，MF)。

1.3 核心基因表达水平差异

利用GEPIA数据库(http://gepia.cancer-pku.cn/)明确核心差异基因在胰腺癌中的表达水平，匹配的正常数据来自TCGA和GTEx数据。免疫组化(immunohistochemistry，IHC)通过利用抗体 - 抗原结合的高特异性，可以揭示蛋白质的相对分布和丰度。免疫组化数据来自人类蛋白质图谱(https://www.proteinatlas.org/)，比较了正常和癌症组织中目标基因在蛋白水平的表达差异。

1.4 核心差异基因与胰腺癌患者生存时间、组织免疫浸润的相关性

Kaplan-Meier Plotter数据库(http://kmplot.com/analysis/)中Kaplan-Meier曲线对排名前五位的核心差异基因进行总生存率(overall survival，OS)和无复发生存率(relapse-free survival，RFS)分析；利用TIMER数据库(https://cistrome.shinyapps.io/timer/)，调整PC的肿瘤纯度后，探索前五位核心差异基因与B细胞、CD4+T细胞、CD8+T细胞、树突状细胞、巨噬细胞等免疫细胞浸润程度的相关性。

1.5 获取免疫浸润矩阵、免疫细胞间相关性及组间差异分析

采用R语言程序并链接CIBERSORT反卷积算法进行22种免疫细胞转录特征的模拟计算。设定模拟次数为1000次，采用Kruskal-Wallis秩和检验对P< 0.05的数据进行后续分析。

对CIBERSORT反卷积算法以P< 0.05基础上筛选的可信样本数据中进行不同免疫细胞间Pearson相关系数计算，并采用秩和检验比较胰腺癌组与对照组差异。

1.6 免疫细胞特征分析

对纳入样本中PC患者不同胰腺组织区域来源的肿瘤组织及正常组织样本进行免疫细胞差异分析，获得癌组织的免疫细胞浸润情况及特征。

1.7 潜在有效中药的预测

Coremine Medical(https://coremine.com/medical/)数据库是一个开放的生物医学数据分析平台，记录了大量生物医学术语间的关系。将两组差异的核心靶基因和GO富集与免疫相关BP导入Coremine Medical数据库，分别筛选具有潜在效应机制的中药。设置筛选条件P< 0.01，筛选可用于胰腺癌治疗的中药。

2 结果

2.1 胰腺癌差异基因

选取GEO数据库中的GSE15471基因芯片进行分析，该芯片来自36位PC患者，病理类型为最常见的PDAC，以其肿瘤组织和自身匹配的正常胰腺组织为样本，并且为减少技术误差共进行了78次基因芯片杂交。筛选获得差异基因2 653个，其中上调基因2 305个，下调基因348个，前20个上下调基因的表达分布情况(见图1A)，差异基因整体表达情况(见图1B)。

图1 PDAC基因差异表达热图及火山图分析Fig.1 Heat map and volcano map of differential gene expression in PDAC注：A：Group A为正常组织；Group B为肿瘤组织；B：差异基因分布，蓝色、灰色、红色分别表示下调、无差异、上调。Note:A,Group A is normal tissue;Group B is tumor tissue;B:Blue,grey and red indicate down-regulation,no difference,and up-regulation,respectively.

2.2 核心基因、PPI及功能注释

将2 653个差异表达基因导入STRING，选择最高置信值0.9,去除无连接的基因后得到PPI。将分析结果文件导入Cytoscape软件进行可视化显示共有94个节点，160条边(见图2)。经cyto Hubba计算列取前10位作为核心基因，表达均呈现上调(见图3)，按重要程度依次排序为I型胶原蛋白a-2(collagen type I alpha 2 chain，COL1A2)、I型胶原蛋白a-1(collagen type I alpha 1 chain，COL1A1)、III型胶原蛋白a-1(collagen tpe III alpha 1 chain，COL3A1)、V型胶原蛋白a-1(collagen type V alpha 2 chain，COL5A1)、V型胶原蛋白a-2(collagen type V alpha 2 chain，COL5A2)、具有血小板反应蛋白基序的去整合素和金属蛋白酶2(ADAM metallopeptidase with thrombospondin type 1 motif 2，ADAMTS2)、基膜聚糖(lumican，LUM)、整合素Α2(integrin alpha 2，ITGA2)、纤维粘黏蛋白1(matrix metallopeptidase 2，MMP2)。

图2 差异基因PPI Fig.2 Differential gene PPI注：颜色越深代表该节点的度值越大，重要程度越高。Note:The darker the color means the greater the degree value of the node,the higher the importance.

图3 核心蛋白基因模块Fig.3 Core protein gene module注：颜色越深代表该节点的度值越大，重要程度越高。Note:The darker the color means the greater the degree value of the node,the higher the importance.

2.3 核心基因表达水平差异

为了验证以上筛选出来的核心基因是否具有代表性和是否具有作为治疗目标的潜力，我们分析了核心基因在癌症和癌症邻近组织中的差异性表达。GEPIA数据库的分析结果显示，前五位核心基因COL1A2、COL1A1、COL3A1、COL5A1、COL5A2在PDAC组织中的表达水平均明显高于正常胰腺组织(见图4A～4E)，差异具有统计学意义(P<0.05)。HPA数据库的IHC结果显示，在蛋白表达水平，COL1A2、COL1A1、COL3A1、COL5A1在胰腺癌中的表达高于正常胰腺组织(见图4F～4I)，COL5A2缺乏相关数据。该部分结果一定程度上验证了以上生信分析的结果。

图4 关键基因在胰腺癌组织和正常胰腺组织中的差异性表达Fig.4 Differential expression of key genes in PDAC tissues and normal pancreatic tissues

2.4 核心差异基因介导的生存分析、免疫浸润

将前5个核心差异基因进行胰腺癌生存分析，结果显示(见图5)COL1A2、COL1A1、COL3A1、COL5A1、COL5A2基因上调表达的患者均表现出较差的OS(排列在前)和RFS(排列在后)，并且以上基因对于患者RFS均具有显著的影响(P< 0.05)，但是对于患者OS并没有统计学意义，说明针对PDAC患者OS的影响，核心差异基因表达水平的改变并非独立的影响因素。

图5 核心差异基因与生存分析Fig.5 Core differential genes and survival analysis

进一步探讨核心差异基因是否通过介导免疫浸润模型改变肿瘤微环境，从而影响最终生存率，结果显示(见图6)，该5个核心差异基因与肿瘤纯度均呈现负相关(P< 0.05)，然而除了COL5A1和COL5A2与CD4+T细胞无显著相关，其余基因与常见六种免疫浸润细胞均呈显著正相关。

图6 核心差异基因与免疫细胞浸润的相关性分析Fig.6 Correlation analysis of core difference genes and immune cell infiltration

2.5 GO及KEGG富集分析结果

为了更加全面了解疾病发病机制中差异基因的功能及其相关免疫过程，利用R语言程序对差异基因进行生物学功能及通路富集分析。差异基因的GO分析(见图7A)。其中，获取BP 2 950条，与免疫相关的生物学过程主要涉及中性粒细胞激活与脱颗粒、血小板脱颗粒、趋化因子介导的信号通路、T辅助细胞的分化、免疫系统过程的负向调节、体液免疫反应等，详见表1。KEGG富集的相关通路主要涉及蛋白质消化与吸收、细胞外基质受体相互作用、PI3K/Akt信号传导、衰老等相关信号通路(见图7B)。

表1 胰腺癌差异基因免疫浸润相关GO富集分析列表Table 1 List of GO enrichment analysis associated with immune infiltration of differential genes in pancreatic cancer

图7 GO及KEGG富集分析Fig.7 GO and KEGG enrichment analysis注：圆形越大占比越高，颜色越深则差异越大。Note:The larger the circle the higher the percentage,and the darker the color the greater the difference.

2.6 免疫浸润细胞分布特征及关联性分析

采用CIBERSORT反卷积法以P< 0.05为筛选条件对芯片进行筛选，共得到77个可信样本，热图左侧35个为正常胰腺组织，右侧42个为PDAC组织。结果显示(见图8)，与正常胰腺组织相比，胰腺癌组织中CD8+T细胞浸润程度低、而CD4+记忆T细胞、M2巨噬细胞、树突状细胞等浸润程度较高。进一步分析了22种免疫细胞在样本中的分布变化特征(见图9)。

图8 正常胰腺组织与胰腺癌组织的免疫细胞分布热图Fig.8 Heat map of immune cell distribution between normal pancreatic tissue and pancreatic cancer tissue

图9 22种免疫细胞浸润比例柱状图Fig.9 Histogram of the proportion of 22 types of immune cell infiltration

NK细胞与CD8+T细胞(r= 0.73)、活化树突状细胞与CD8+T细胞(r= 0.75)、调节T细胞与活化树突状细胞(r= 0.71)、肥大细胞与CD8+T细胞(r=0.64)等呈现较强正相关；M1巨噬细胞和记忆T细胞分别与活化树突状细胞(r= -0.59；r= -0.64)、活化肥大细胞与未活化肥大细胞(r=-0.62)之间呈现较强负相关(见图10)。

图10 22种免疫细胞相关性热图Fig.10 Heat map of 22 immune cell correlations

2.7 免疫细胞浸润差异性分析

通过小提琴图(见图11)对不同组织样本免疫浸润细胞差异分析进行可视化，发现肿瘤组织中活化的记忆CD4+T细胞、未活化NK细胞、M0/M1/M2巨噬细胞、未活化树突状细胞、中性粒细胞的浸润程度显著增多(P< 0.05)；同时，初始B淋巴细胞、血小板、CD8+T细胞、未活化记忆CD4+T细胞、辅助T细胞、调节T细胞、活化NK细胞、单核细胞、活化树突状细胞的浸润程度显著减少(P< 0.05)。

图11 正常组织与肿瘤组织之间的免疫浸润差异Fig.11 Differences in immune infiltration between normal and tumor tissues注：蓝色代表正常胰腺组织，红色代表胰腺癌组织。Note:Blue represents normal pancreatic tissue,red represents pancreatic cancer tissue.

进一步筛选获得38组符合要求的配对样本，样本中各种免疫细胞改变趋势进行分析的结果表明(见图12)，在PDAC发病过程中个体免疫细胞变化趋势不尽完全一致，但可以观察到M0巨噬细胞多呈上升趋势(P= 0.001)、活化的NK和DC细胞多呈下降趋势(P= 0)。

图12 胰腺癌配对样本免疫浸润特征组合Fig.12 Immune infiltration characteristics of paired samples of pancreatic cancer

2.8 中药预测分析

通过Coremine Medical预测具有潜在治疗PDAC的中药，输入前十名核心差异基因及筛选得到的免疫相关生物过程名称。若以P< 0.05为筛选标准获得不同靶点的中药227味，删除重复的53味中药后获得174味中药。为更加精准选择用药，最终确定以P< 0.01为筛选标准，规范中药名称后，进行整理共获得77种中药，其中某些基因及生物学过程缺乏符合要求的中药，最终结果详见表2。将最终得到的中药输入Cytoscape，出现2次以上为高频中药，葶苈子、姜黄、皂荚、白附子、桂枝、杜仲、人参、雷丸、香椿、鹿角、三七的度值较高(见图13)。整体预测中药的功效分布以清热类(22味)、补益类(10味)、祛风湿(9味)为主，多归肝经(见图14A、14B)。

表2 中药靶向核心差异基因及免疫相关生物学过程的预测结果Table 2 Predicted results of Chinese medicine targeting core differential genes and immune-related biological processes

图13 靶向关键基因和免疫相关通路的中药预测Fig.13 Traditional Chinese medicines predicted by key genes and immune-related pathways注：图13中黄色菱形、三角形、圆形节点分别表示免疫相关生物学过程、核心基因、预测中药；颜色越深代表关联性越强；Note:The yellow prism,triangle and circle nodes in Figure 13 indicate immune-related biological processes,core genes,and predicted herbal medicines,respectively;the darker the color represents the stronger the association.

图14 中药预测特征可视化Fig.14 Visualization of Chinese medicine prediction features注：A为中药归经雷达图；B为中药功效分布图。Note：A is the radar map of herbal medicine attribution;B is the distribution of herbal medicine efficacy.

3 讨论与结论

目前，胰腺癌的免疫治疗取得了一定进展，给临床治疗胰腺癌带来了希望，然而，由于胰腺癌的微环境通常包括免疫抑制细胞增多、免疫细胞失活和肿瘤突变负荷低等状态，多数免疫治疗结果并不令人满意。随着高通量技术的发展，Lv等[8]通过多个数据集来挖掘胰腺癌的差异核心基因，Xiao等[9]通过“炎-癌”差异基因挖掘中药等等，都加深了对胰腺癌发病分子机制的理解，但是对于该病免疫相关机制的工作开展不足，不利于改善患者的长期生存预后。因此，本研究首次在差异基因的基础上，进一步挖掘相关免疫浸润机制和潜在有效中药，旨在完善机制和促进药物开发。首先，本研究通过GEO数据库筛选合适的基因芯片，比较同一个体患者的胰腺肿瘤与正常组织样本的基因表达谱，得到2 653个差异基因，并通过PPI得到核心差异基因，分别是细胞外基质组织和胶原纤维组织相关基因(COL1A2、COL1A1、COL3A1、COL5A1、COL5A2)、ADAM金属肽酶(ADAMTS2)、基膜聚糖(LUM)、整合素亚基α2(ITGA2)、纤维连接蛋白(FN1)、基质金属蛋白酶(MMP2)。以上差异基因经证明在肿瘤组织中均表达上调，与胶原蛋白代谢相关，多属于COL1、COL3和COL5家族[10]，其中COL1A2的表达通过qRT-PCR测定验证在组织中表达最高[11]，反应此次结果可信度高。与其他实体瘤不同，PDAC肿瘤微环境是一个包含胰腺星状细胞、癌相关成纤维细胞、免疫细胞和细胞外基质过度纤维化的致密基质环境，基质中密集纤维化反应和肿瘤在免疫环境中的改变被认为是目前PDAC治疗失败的主要原因。进一步通过生存分析研究揭示该类核心差异基因与PDAC患者RFS显著相关，并且与多种肿瘤细胞浸润程度呈正相关，提示该类基因可能通过改变肿瘤免疫微环境介导PDAC患者的临床预后，如何借助该点改变临床结局成为可能的突破口之一。

进一步对差异基因进行GO及KEGG富集分析，提示差异基因的生物学过程涉及中性粒细胞脱颗粒、迁移、趋化因子、T辅助细胞分化、体液免疫反应等过程。PDAC微环境本质上是炎症的，趋化因子等促进肿瘤炎性微环境的事实也已经被证明，除细胞因子和趋化因子外，中性粒细胞的激活往往是第一个分子事件，随后发生迁移向肿瘤组织趋化。已经明确的是PDAC癌细胞可以将多形核中性粒细胞募集到肿瘤附近，但大多无法产生抗肿瘤反应，通常高丰度的中性粒细胞数量会促进PDAC的进展，与预后不良有关[12]；最近体外实验表明，胰腺癌细胞诱导中性粒细胞胞外陷阱形成，从而促进PDAC肝脏微转移的发生[13]；同时高度纤维化导致的局部缺氧变化决定了免疫环境，活性氧的持续释放导致局部缺氧，从而增强中性粒细胞浸润，因此，中性粒细胞在胰腺癌微环境中至关重要[14]。并且，本次研究结果提示差异基因直接参与了T细胞的分化和体液免疫的调节。

为深入探索免疫浸润在PDAC发病中的作用，本研究利用CIBERSORT反卷积法对总体肿瘤样本进行分析，与正常胰腺组织相比，肿瘤组织中M2巨噬细胞、静息记忆CD4+T细胞显著增高，而CD8+T细胞显著下降，在配对样本中观察到了同样趋势。M2巨噬细胞产生促进肿瘤进展的抗炎信号，经确认是胰腺癌预后不良的标志之一。2013年，一项通过IHC观察PDAC患者免疫细胞浸润程度的研究表明，较高水平的肿瘤浸润性泛巨噬细胞M2、Tregs与CD4+T与较短的生存期显著相关[15]。在PDAC小鼠模型中已经证明肿瘤相关巨噬细胞通过释放细胞因子、蛋白酶和生长因子介导免疫抑制和血管生成并促进肿瘤进展[16]，并且，抑制CCR2已被证明可通过阻断单核细胞募集到肿瘤微环境中来提高化疗效果、抑制转移、增强放疗效果并增加T细胞免疫浸润[17,18]。因此，巨噬细胞极可能成为人类PDAC新治疗策略的潜在目标。既往认为由成纤维细胞和促纤维母细胞组成的PDAC免疫抑制微环境会限制T细胞浸润，但最近的研究揭示了T细胞的特定空间分布，正如相关性热图显示胰腺癌组织中静息的记忆CD4+T细胞与CD8+T细胞之间存在显著负相关(r = -0.64)，这已被证明是预后不佳的标志之一，尤其是PDAC，CD8+T细胞在肿瘤组织中的浸润与肿瘤细胞紧密接近和患者存活率紧密相关[19]。最近有学者强调在长期幸存者中建立成功的免疫更需要高丰度的肿瘤内CD8+T细胞，了解CD8+T细胞的驱动因素以及CD8+T细胞的浸润机制，对于进一步开发PDAC治疗方案至关重要[14]。目前，PDAC免疫浸润细胞之间的关联尚缺乏大规模的验证，本预测具有一定参考意义。

PDAC在中医属于“脾积”“癥瘕”“积聚”“痞气”“伏梁”等范畴。病例系列文献研究和Meta分析结果提示联合中医药治疗PDAC能够有效降低患者死亡率，改善远期预后[20]。最近基础研究揭示中草药治疗PDAC的机制多通过改善肿瘤微环境而间接改善医疗结局[21]，故探索肿瘤免疫微环境对寻找有效中药十分重要。目前，PDAC中医病机多认为是“本虚标实”，以湿热蕴毒为著，而本研究预测中药的功效多属于清热类(黄芩、山慈姑、牡丹皮、紫草等)，其次为补益药(人参、杜仲、仙茅、龟甲、阿胶等)，符合中医病机。从病位而言，对比中医古籍中关于解剖位置和生理功能的描述，多数学者支持将胰归属于五脏中的脾，但本次研究显示绝大多数有效预测中药归属肝经，其次是肺、肾经，提示我们或许从肝论治湿热，结合扶助脾虚能够进一步提高治疗有效率，经检索发现该理论已有散在经验[22]和综述报道[23]，值得深入探讨形成系统认识。此外，本研究利用反向预测的方法获得了数味中药，在关联度较高的中药队列中，研究发现人参皂苷能够通过下调MMP2等表达从而减少血管生成，降低肿瘤细胞营养供给[24]；姜黄素通过Wnt/β-Catenin信号通路活性抑制PDAC干细胞的自我更新[25]，并可改善PDAC皮下移植瘤模型组织的局部微循环,增加肿瘤微环境中的供血供氧,抑制肿瘤细胞和间质细胞的低氧应激反应[26]；以上文献支持了本研究的部分预测结果，但有关中药治疗PDAC相关免疫机制的研究开展仍十分不足。

本研究在筛选基因芯片时考虑到免疫机制分析的个体差异，严谨筛选了采用配对样本分析的PDAC实验数据，同时考虑到本研究本质上是对既往单个GSE数据的二次挖掘，采用了不同水平来验证核心差异基因的表达，降低了风险。但是本研究仍存在以下不足：(1)由于样本量过少导致数据分析具有一定偏倚，后期纳入更多符合要求的数据集来进行分析；(2)CIBERSORT反卷积算法分析是基于有限的遗传数据，这些数据可能会偏离细胞异型相互作用、疾病的诱发因素以及疾病表型的可塑性，临床和基础试验可以更加全面、客观地反映个体免疫情况。总之，本研究首次利用生物信息学技术结合CIBERSORT反卷积算法探索了PDAC的免疫浸润机制，并预测潜在可能有效的中药，对未来该领域的研究具有一定价值。