基于GEO 数据库的胰腺癌循环肿瘤细胞代谢通路及关键基因研究

林杏仪，田振烽，潘乐乐，苏铭昕，陈茵婷

胰腺癌是一种恶性程度极高的消化道肿瘤，5 年生存率仍低于10%。由于缺乏典型的症状和体征，大多数胰腺癌初诊时，已处于肿瘤晚期，甚至伴有淋巴结、肝脏等远处转移［1］。转移是胰腺癌患者死亡的主要因素。胰腺癌特殊的肿瘤微环境可通过多种方式促进肿瘤转移，如上皮间质转化、转移前生态位形成（pre⁃metastatic niche forma⁃tion）。其中，转移前生态位是指原发肿瘤可以在转移部位形成支持性微环境。循环肿瘤细胞（cir⁃culating tumor cells，CTCs）是指从原发性或转移性肿瘤脱落进入血液循环的肿瘤细胞群，常被认为是肿瘤转移的“种子”［2］。CTCs 可能优先种植于转移前生态位已形成的器官，因为这些部位能提供适合的“土壤”供“种子”生长［3，4］。CTCs 还可以与其他类型的细胞，如免疫细胞、血小板、肿瘤相关成纤维细胞形成特殊的相互作用。例如，在乳腺癌中，中性粒细胞可以通过这种相互作用，提高CTCs 的转移潜力［5］。CTCs 可以单个细胞或细胞簇的形式存在，有研究报告，细胞簇比单个细胞更易于转移［6］。

代谢重编程被认为是肿瘤的标志［7］，为了应对肿瘤微环境中低氧、营养物质有限等限制，肿瘤细胞通过代谢重编程来支持细胞的生长和增殖。胰腺癌细胞具有广泛代谢重编程，例如糖代谢、谷氨酰胺代谢等改变。该重编程由癌基因介导的自主通路、独特的肿瘤微环境以及肿瘤与周围非癌细胞的相互作用共同驱动［8］。然而胰腺癌CTCs 发生的代谢重编程、机制及其作用仍未明确。

本研究旨在通过对CTCs 的生物信息学分析，探索CTC 的代谢重编程可能的机制，并研究其中关键基因的临床意义，为后续研究CTCs 与胰腺癌转移的研究提供重要的分子基础。

1 资料与方法

1.1 数据资料搜集

于 GEO 基因表达数据库（https：//www.ncbi.nlm.nih.gov/geo/）检索胰腺癌循环肿瘤细胞（circu⁃lating tumor cells，CTCs）相关基因芯片信息。纳入标准：①种属为人类（Homo sapiens）；②同时具备原位肿瘤组织样本和CTCs样本信息；③未经化疗或其他药物处理；④具有全基因组RNA表达谱。经检索后，符合的数据集有：GSE118556［9］和 GSE18670［10］。其中前者来源于胰腺癌细胞系PANC⁃1，后者来源于临床胰腺癌患者的样本，本研究中，后者为验证数据集。

1.2 数据集PCA 主成分分析、差异基因（differen⁃tially expressed gene，DEGs）筛选

在R（version 4.1.2）中利用R 包“GEOquery”进行数据集下载、探针⁃基因名转换以及总体表达谱的可视化。采用“FactoMineR”及“factoextra”包进行PCA 主成分分析。利用“limma”包进行DEGs 分析，可获得差异倍数（fold change，FC）、P值、校正P值等数据。两数据集均以|FC|≥ 1.5 且校正P值<0.05 为阈值筛选DEGs。基因在数据集中的表达情况通过“pheatmap”包可视化。对GSE118556 基因集体内组和体外组DEGs 的上调基因、下调基因分别取交集，获得共同上调及共同下调基因，共同上、下调基因统称为共同DEGs。共同上/下调基因的火山图和韦恩图可视化分别通过R 包“EnhancedVolcano”和“VennDiagram”实现。

1.3 GO 和 KEGG 通路富集

利用 DAVID 网站（https：//david.ncifcrf.gov/）进行GO 富集分析，并以“ggplot2”可视化。为研究DEGs 的代谢重编程情况，从KEGG 数据库的615条通路中提取了92 条代谢通路，并从中提取了1703 个代谢相关基因。对代谢相关基因与共同DEGs取交集，获取代谢相关DEGs。在R中，对代谢相关 DEGs 利用“clusterProfiler”R 包［11］进行 KEGG通路富集并可视化。

1.4 蛋白互作（protein⁃protein interaction，PPI）网络构建

将代谢 相 关 DEGs 导 入 STRING（https：//cn.string⁃db.org/）进行PPI 网络构建，隐藏无连接节点后，利用Cytoscape 软件可视化并计算节点的连线（degree）。

1.5 TKT 生存分析和与临床病理分期的关系

从 UCSC Xena（https：//xenabrowser.net/datapag⁃es/）中下载TCGA 胰腺癌临床信息及基因表达矩阵，使用“survival”和“survminer”R 包，Kaplan⁃Meier方法来绘制OS 曲线。使用“ggpubr”和“ggplot2”对TKT 和临床病理分期、肿瘤分级等进行分析并绘制箱式图，TKT 表达水平以中位数（四分位数）来表示。两组间的比较采用t检验。

1.6 TIMER 数据库分析

使用TIMER（Tumor IMmune Estimation Resource）数据库（https：//cistrome.shinyapps.io/timer/）对TKT、TKTL1、TKTL2 与肿瘤免疫细胞浸润水平相关性进行分析。研究的肿瘤免疫细胞包括：B 细胞（B cell）、CD8+T 细 胞（CD8+T cell）、CD4+T 细 胞（CD4+T cell）、巨噬细胞（macrophage）、中性粒细胞（neutrophil）和树突状细胞（dendritic cell）。P<0.05时有统计学意义。

2 结 果

2.1 筛选差异基因

GSE118556 数据集中，作者把亲本胰腺癌细胞 PANC⁃1（即 PANC⁃1⁃P，in vitro）注射到 SCID 小鼠皮下。一月后，心脏穿刺获取约1mL 小鼠血液，处理并鉴定后获得循环肿瘤细胞PANC⁃1⁃CTC（in vitro）。再将上述两种细胞分别注射至小鼠腹壁使其成瘤，分别获得PANC⁃1⁃P（in vivo）与PANC⁃1⁃CTC（in vivo）样本。数据集所有样本的总体表达谱情况如图1A 箱式图所示，并对其进行PCA 主成分分析（图1B），可见P 组（亲本组）与CTCs 组两个分组有明显的差异。进一步对四组样本的基因表达谱分析，分别分析体外组（in vitro）和体内组（in vivo）：以各自的亲本组（PANC⁃1⁃P）为对照组，CTCs 组为实验组进行差异分析，从而获得两组差异基因（differentially expressed genes，DEGs）。以差异倍数≥ 1.5，校正P值<0.05 作为筛选 DEGs 的阈值。详见图1C，其中体外组DEGs 共7215 个，体内组DEGs共3337个。如图1D的韦恩图所示，体外组上调基因5245 个，下调基因1970 个；体内组上调基因1594 个，下调基因1743 个。在两组均上调的基因有536 个，298 个基因均下调。将共同上调和下调的基因合并，形成本数据集的共同DEGs，共834 个基因。

图1 GSE118556 差异基因分析

2.2 基因本体（gene ontology，GO）分析

为进一步了解共同DEGs 在细胞内的生物过程（biological process，BP）、细胞组分（cellular com⁃ponent，CC）、分子功能（molecular function，MF）情况，利用GO 数据库富集DEGs，并对P值小于0.05的组分进行分析。如图2 所示，在生物过程（BP）方面，主要与RNA 聚合酶Ⅱ启动子的正性转录调控（positive regulation of transcription from RNA poly⁃meraseⅡpromoter）、信号转导（signal transduction）、和DNA 模板的转录负调控（negative regulation of transcription，DNA⁃templated）相关。在细胞组分（CC），DEGs 主要与胞质（cytosol）、胞外外泌体（extracellular exosome）、和胞膜（membrane）相关。而DEGs 的分子功能（MF）主要集中在：蛋白结合（protein binding）、钙离子结合（calcium ion bind⁃ing）、和序列特异的 DNA 结合（sequence⁃specific DNA binding）。

图2 GSE118556的基因本体（GO）富集分析

2.3 CTCs 的代谢重编程

为了进一步探索CTCs 中代谢重编程的情况，笔者进一步分析了KEGG 通路数据库和共同DEGs。从KEGG 数据库的92 条代谢物富集途径中，提取1703 个代谢相关基因。而同时出现在代谢相关基因和共同DEGs 的基因有66 个，所有代谢相关DEGs 在数据集中的表达情况如图3A 所示，对代谢相关DEGs 进行KEGG 通路富集分析，以校正P值＜0.05 为界限，通路富集情况见表1。其中，如图3B 所示，基因富集前三位的代谢通路为半胱氨酸与蛋氨酸代谢（cysteine and methionine metabolism）、一碳代谢（carbon metabolism）和辅酶因子的生物合成（biosynthesis of cofactors）。说明这三条代谢通路可能在PANC⁃1 的CTCs 形成上起关键作用。为了解代谢相关DEGs 所编码蛋白的相互作用，利用STRING 构建蛋白互作网络，并以Cytoscape 软件可视化。如图3C 所示，其中GOT2（glutamic⁃oxaloacetic transaminase 2，谷草转氨酶2）、TKT（transketolase，转酮醇酶）与较多节点（基因）连接，说明两基因处于该调控网络的中心地位。

图3 GSE118556 代谢相关DEGs

2.4 临床数据集验证代谢重编程

由于数据集GSE118556 的样本基于胰腺癌细胞株PANC⁃1，因此笔者利用来源于临床样本的GSE18670 数据集里面的部分数据进行验证。提取数据集中的肿瘤样本和CTCs 样本数据，以肿瘤样本为对照组，CTCs 样本为实验组。提取的两组样本的PCA 主成分分析见图4A，可见两组样本的差异较大，可进行后续的差异分析。利用“limma”包在R 中进行差异分析，仍以差异倍数≥1.5，校正P值＜0.05 为阈值，筛选获得 1119 个 DEGs，其中539 个上调基因，580 个下调基因，详见图4B。与代谢相关基因比对，发现DEGs 中有106 个代谢相关基因，其在数据集的表达情况如图4C。对这些基因进行KEGG 通路富集，发现前三位的富集通路分别为一碳代谢（carbon metabolism）、嘌呤代谢（purine metabolism）、甘油磷脂代谢（glycerophos⁃pholipid metabolism），如图4D。如表1 所示，一碳代谢通路同时在GSE118556 数据集的代谢通路富集中处于第二位，说明一碳代谢通路可能在CTCs形成和维持中可能起关键作用。在两数据集的一碳代谢通路中，TKT 在两者中均出现差异表达并上调，因此TKT 可能在一碳代谢通路中起关键的调控作用。

表1 GSE118556 及GSE18670 差异基因KEGG 代谢通路富集分析

图4 GSE18670 临床数据集验证

2.5 TKT 与预后和肿瘤分期的关系

经上述生物信息学分析，TKT 可能发挥关键调控作用，因此利用TCGA 胰腺癌临床数据库，对TKT的临床意义进行分析。如图5A所示，TKT高表达时患者的总体生存时间（overall survival，OS）显著下降（P=0.00082）。在分期方面（图5B），分别与I 期对比，IIA 期与IIB 期的患者的TKT 水平显著升高（P=0.00091，P=0.0031）。在组织学分级中（图5C），G2 和G3 期患者的TKT 表达均较G1 期明显上调（P=0.02，P=0.0078）。而对于胰腺癌TNM 分期，描述原发肿瘤的T 分期中（图5D），T1、T2、T3 期患者的TKT 水平呈逐渐上升的趋势，但是只有T2 和T3期的TKT 水平有统计学差异（P=0.01385）。而N0和 N1 期，M0 和 M1 期的 TKT 水平没有统计学差异，详见图5E、F。

图5 TKT 表达与预后、肿瘤分期的关系

2.6 TKT 表达与免疫浸润水平

采用TIMER 数据库分析TKT 与免疫浸润水平的相关性。如图 6 所示，TKT 与 B 细胞、CD8+T 细胞、中性粒细胞和树突状细胞呈正相关的关系，但是没有统计学意义（P>0.05）。而同样编码TKT（转酮醇酶）同工酶的基因TKTL1 和TKTL2［12］与免疫浸润有显著相关性：TKTL1 与 B 细胞、CD8+T 细胞、CD4+T 细胞、巨噬细胞、中性粒细胞和树突状细胞水平均有显著相关性（P<0.05），且均为正相关（r>0）；TKTL2 表达只与巨噬细胞水平显著相关（r=-0.179，P=0.0191）。

图6 胰腺癌中TKT、TKTL1、TKTL2 表达水平与免疫浸润水平的相关性

3 讨 论

本研究通过对数据集GSE118556 和GSE18670表达谱的差异分析和KEGG 代谢通路富集分析发现，一碳代谢为排名较前的富集通路，由此推测一碳代谢可能在胰腺癌CTCs 形成和维持中起关键作用。而一碳代谢相关基因TKT 在两个数据集的DEGs 中重复出现，TKT 可能在整个调控过程中处于中心地位。生存分析发现，TKT 表达水平高的胰腺癌患者，总体生存期更短。此外，TKT 的表达水平与肿瘤分期、组织学分级、T 分期相关。而编码 TKT 的同工酶基因 TKT、TKTL1、TKTL2 中，TKTL1、TKTL2 均不同程度上调，与免疫浸润水平相关。

癌细胞通过重新调整自身的代谢来支持其增强的增殖水平。本研究发现，一碳代谢在PANC⁃1小鼠胰腺癌模型和胰腺癌患者CTCs 中，均位于富集的代谢通路前列。一碳代谢包括叶酸和蛋氨酸循环，通过此途径，细胞产生一碳单元（也称为甲基），并用于合成重要的代谢前体和甲基化反应。首先，嘌呤和嘧啶核苷酸的生物合成都需要一碳单元。癌细胞的增殖能力普遍增强，因此对于核苷酸的需求增加，通过增强一碳代谢以应对需求。其次，一碳代谢为癌细胞内的甲基化反应提供甲基供体。肿瘤细胞常出现DNA 甲基化模式的改变，一碳代谢产生的SAM（S⁃腺苷甲硫氨酸）可以为蛋白质翻译后修饰的甲基化、DNA 甲基化、调节RNA 甲基化的过程［13］提供甲基。此外，一碳单元还能影响NADH/NADPH 的产生，从而改变多种代谢和合成途径［14］。本研究中发现，CTCs 内一碳代谢通路相关基因表达升高，证明了CTCs 可能有着比一般癌细胞更强的增殖能力，并且部分蛋白质、DNA 的甲基化水平可能出现改变。

此外，我们发现在两个数据集的一碳代谢通路中，TKT 均显著性升高。TKT 基因编码的TKT（转酮醇酶）催化磷酸戊糖途径中几个非氧化分支的关键反应，这些反应在利用葡萄糖合成核糖的过程中起关键作用。TKT 已经被发现与阿尔兹海默病、糖尿病等疾病相关［15，16］。某些肿瘤细胞甚至有超过85%的核糖是通过该途径合成。TKT 的上调，使得肿瘤细胞能够在非氧化的情况下利用葡萄糖合成核糖［17］，而胰腺癌周围充满间质细胞而缺少血管，处于缺氧的微环境中，TKT 的上调大大增强癌细胞在缺氧情况下的生存能力。

此外，以往的研究发现CTCs 可以与免疫细胞相互作用。有研究报道，CTCs 与白细胞（WBC）相互作用形成的CTC⁃WBC 簇会加速乳腺癌小鼠的转移发生［18］。而在转移性乳腺癌研究中，发现CTCs 的数量和中性粒细胞与淋巴细胞的比率（NLR）相关，且CTCs 和 NLR<3 的患者，疾病复发的风险高8 倍［19］。本研究中，编码转酮醇酶（TKT）的同工酶的三个基因（TKT、TKTL1、TKTL2）中，TKTL1、TKTL2 均与免疫浸润相关，说明转酮醇酶的上调可能与CTCs 转移能力增强相关。

本研究创新地把胰腺癌CTCs 与代谢重编程联系在一起，并发现CTCs 中一碳代谢通路上调、通路中的TKT 显著上调且与预后相关。因此CTCs可能是由于胰腺癌细胞内一碳代谢的改变，增殖能力增强，并造成肿瘤细胞脱落并进入血液循环，从而使肿瘤转移到其他部位。然而本研究仅对GEO 数据库中的数据集进行生物信息学分析，仍需要进一步的体内和体外实验来验证一碳代谢对CTCs 的影响。另外，一碳代谢影响CTCs 的具体机制和作用仍有待深入探索。

综上所述，本研究为一碳代谢在CTCs 生成和维持的代谢重编程中的作用提出了新的观点。一碳代谢和TKT 的靶向治疗，可能会为减缓胰腺癌肿瘤转移提供新的可能。