分享一篇2022年2月18日发表在frontiers in Cell and Developmental Biology（影响因子 6.684）《Single-Cell Transcriptome Reveals the Metabolic and Clinical Features of a Highly Malignant Cell Subpopulation in Pancreatic Ductal Adenocarcinoma》。胰腺导管腺癌（PDAC）是一种恶性肿瘤，死亡率很高。PDAC表现出明显的异质性以及代谢途径的改变，与它的恶性进展有关。在这项研究中，我们基于单细胞转录组技术，探索了PDAC高度恶性亚组的代谢和临床特征。
根据恶性基因的表达，在单细胞分辨率下确定了一个高度恶性的细胞亚群。根据代谢途径基因组，分析了不同细胞类型的代谢情况。进行了体外实验来验证标记基因的生物学功能。根据五个糖酵解标记基因，区分了具有高度恶性细胞亚群的PDAC患者亚群。五个糖酵解高度恶性相关基因特征被用来构建糖酵解高度恶性相关基因特征（GHS）评分。
介绍
代谢重编程已被认为是癌症的一个共同特征。肿瘤细胞利用重编程采用几种代谢途径，以满足细胞生长的大量能量需求。在PDAC中，肿瘤在高度纤维化和结缔组织增生的微环境中发展，导致血管受压和血液灌注不足。因此，PDAC细胞的葡萄糖摄取率大多是中等的。然而，它在营养不足的环境中表现出较高的生长速度。这些特征反映了新陈代谢在PDAC进展中的关键作用。因此，近年来开发了更多关注PDAC代谢靶点的治疗方式。然而，目前还缺乏合适的生物标志物来区分代谢亚型。单细胞RNA测序（scRNA-seq）具有在单个细胞水平分析转录组信息的特点，与传统的批量RNA-seq相比，scRNA-seq的特点将深入了解恶性肿瘤细胞的代谢异质性。
结果
单细胞RNA测序发现了一个高度恶性的导管细胞亚群
scRNA-seq数据来自24名PDAC患者，通过标准化的单细胞处理步骤，最终确定了10个不同的细胞簇，包括针状细胞、B细胞、导管1细胞、导管2细胞、内皮细胞、内分泌细胞、成纤维细胞、巨噬细胞、星状细胞和T细胞。每个细胞簇的标记基因都与已知的细胞标记一致（图1A）。为了进一步验证细胞亚群的分配，我们使用InferCNV R软件包进行了推断基因拷贝数分析。使用HMM方法预测了导管和尖锐湿疣细胞的CNV分数。导管2细胞的CNV分数明显高于导管1细胞（图1B），这表明导管2细胞是一个恶性细胞亚群。为了研究恶性细胞之间在恶性程度上的区别的可能性，我们整理了136个参与PDAC恶性进展的基因。scRNA-seq数据归一化后，我们根据收集到的PDAC恶性基因对恶性导管2细胞进行评分，发现导管2细胞的恶性评分呈双峰分布。这一发现意味着PDAC中的一些细胞表现出高度恶性倾向的特征。对数据进行高斯拟合后，我们根据恶性评分将恶性的导管细胞分为高恶性和低恶性组（图1C）。
图1. 识别PDAC中高度恶性的细胞亚群。(A）40,986个细胞的基因表达的t-SNE降维，显示24个PDAC患者的主要细胞类型。(B) 通过t-SNE分析对导管1、导管2和尖锐湿疣细胞的拷贝数变异（CNV）评分。(C)导管2细胞的恶性评分的密度图。(D）比较导管2细胞中高度恶性和低度恶性细胞的UMAP图。(E) 相应的恶性评分水平的UMAP图[从低（灰色）到高（红色）的比例] 。(F）高度恶性细胞生物过程（BP）GO术语（最多五个）的柱状图，按FDR（P<0.05）排序。
高度恶性细胞亚群的代谢重编程
我们参考了一种公共方法来量化不同类型细胞的代谢途径。这使得我们能够进一步调查不同细胞类型之间代谢途径的变化和整体特征，特别是高度和低度恶性细胞亚群之间的代谢途径。不同的细胞类型表现出不同的代谢概况。在免疫细胞中，巨噬细胞表现出最高的代谢活性。核黄素的代谢在巨噬细胞中明显上调。在成纤维细胞中，糖胺聚糖的生物合成被明显上调。引人注目的是，高度恶性细胞亚群在所有细胞类型中显示出最高的代谢活性。此外，低度恶性细胞的代谢途径比导管一细胞的代谢途径更活跃。这表明，代谢活动与肿瘤细胞的恶性进展相关（图2A，B）。在所有79条代谢途径中，有44条在高度恶性的细胞亚群中高度表达。过度表达的代谢途径主要集中在碳水化合物代谢、辅助因子和维生素的代谢、糖类的生物合成和代谢以及脂质代谢。为了进一步揭示高度恶性细胞亚群的重要代谢途径，我们对所有代谢途径的得分与细胞恶性程度的得分进行了皮尔逊相关分析。嘧啶代谢、糖酵解/葡萄糖生成、嘌呤代谢、酪氨酸代谢和药物代谢-其他酶与肿瘤细胞的恶性程度表现出明显的正相关（图2C和补充图S2C），表明这五个代谢途径可能与PDAC的恶性进展有关。随后，我们比较了这五个代谢途径在导管1细胞和低度和高度恶性细胞中的表达。发现在高度恶性的肿瘤细胞中，糖酵解的表达最为明显（图2D）。这些分析表明，加强糖酵解是高度恶性细胞亚群的主要代谢特征。
图2. 高度恶性细胞亚群的代谢重编程。(A) 表达热图显示了胰腺导管腺癌的不同细胞类型中京都基因和基因组百科全书（KEGG）代谢途径的丰富程度。(B) 每种细胞类型的代谢得分的维林图。Wilcoxon排名测试，p = 0.00084。(C) 所有癌细胞中与恶性分数相关性最高的代谢途径（前五位）的相关热图（Pearson相关）。红色，正相关；蓝色，负相关；白色，无相关。(D) 条形图显示了三种细胞类型中五条代谢途径的得分，并使用方差分析计算了P值。管状1细胞，灰色柱状图；低度恶性，蓝色柱状图；高度恶性，橙色柱状图。****p < 0.0001.
高度恶性细胞亚群的糖酵解标记基因
为了进一步揭示PDAC肿瘤细胞恶性发展过程中糖酵解的关键基因，我们用|logfc|>1和P值<0.01来比较高度和低度恶性细胞，发现共有364个基因在高度恶性细胞中明显上调（图3A和补充表S3）。在这364个基因中，确定了糖酵解基因，包括EN O 1、LDHA、PKM、PGK1和PGM1（图3B，C）。此外，在这364个基因中，有16个基因属于以前确定的5个恶性肿瘤相关的代谢途径。蛋白质-蛋白质相互作用（PPI）网络证实了糖酵解标记基因在恶性肿瘤相关代谢基因中的中心地位。此外，我们试图分析五个糖酵解基因与患者预后的关系。EN O 1（p < 0.0001）、LDHA（p = 0.00013）、PKM（p = 0.00021）、PGK1（p = 0.00082）和PGM1（p = 0.015）的高表达与PDAC患者的生存期明显较差有关（图3D），表明这五个基因可作为高度恶性细胞亚群的糖酵解标志物。
图3
PDAC中五个糖酵解标志基因的mRNA和蛋白质表达概况
为了研究PDAC中五个糖酵解标志基因的mRNA表达水平，我们使用了基因表达总库（GEO）的数据集进行进一步分析。GSE62452数据集共包含130个样本，包括肿瘤和配对的正常活检样本。将表达谱归一化后，与配对的正常活检组织相比，PDAC中五个糖酵解标志基因的mRNA表达水平明显升高（图4A）。此外，这五个基因的mRNA表达在TCGA队列中也进行了类似的研究。同样，在肿瘤组织中也观察到这五个糖酵解标志基因的mRNA水平升高。在使用TCGA程序和GTEx程序进行的分析中也观察到了这一发现。随后的步骤是利用HPA的免疫组织化学（IHC）分析来检查PDAC中五个糖酵解标记基因的阳性表达。这五个基因的mRNA和蛋白表达在PDAC中都是过度表达的（图4B）。这些结果与之前scRNA-seq数据的分析结果一致。
糖酵解标志基因参与PDAC细胞的恶性进展
为了研究糖酵解标志基因在PDAC细胞中的影响，我们在PANC-1和CFPAC-1细胞系中分别沉默了五个糖酵解标志基因。每个基因都使用了两个独特的沉默位点。西方印迹表明，与对照组细胞相比，EN O 1、LDHA、PKM2、PGK1和PGM1的表达明显下降（图4C）。用CCK8增殖测定法测量各组在不同时间点的吸光度。如图所示，在PANC-1和CFPAC-1细胞系中，与shCont组相比，糖酵解标记基因的沉默降低了增殖能力（图4D，E）。在正常的胰腺导管细胞系HPDE中，这五个糖酵解基因的敲除对增殖率的影响没有肿瘤细胞那么明显。这些数据证实了这五个糖酵解标志基因在PDAC的恶性进展中的贡献。
图4. 糖酵解标志基因参与PDAC的恶性进展。(A) GSE62452队列中肿瘤和配对的正常组织之间五个糖酵解标志基因的表达水平。蓝色，正常；橙色，肿瘤。(B) HPA数据库中10名PDAC患者的EN O 1、LDHA、PKM、PGK1和PGM1蛋白的表达水平。比例尺：20 µm。(C) PANC-1细胞和CFPAC-1细胞用sh-对照、sh-ENO1、sh-LDHA、sh-PKM2、sh-PGK1和sh-PGM1转染。使用Western blot验证敲除效率。(D) 分别敲除5个糖酵解标志基因后，用CCK8测定PANC-1细胞的生存能力。(E) 分别敲除5个糖酵解标志基因后，用CCK8法测定CFPAC-1细胞的活力。*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001。
为了确定高度恶性细胞亚群的临床特征，我们从TCGA（n = 146）获得了五个标记基因的表达谱，并使用UMAP来降低维度（图5A）。五个标记基因表现出线性上升的特征，根据四分位数间的分界线，将患者分为C1、C2和C3组。(图5B)。五个糖酵解标志基因在C3亚型中表达量最高，并显示出明显的预后差异，C3亚型的生存率最低（P = 0.043）（图5C）。在GSE62452（n = 65）队列中，类似的增加趋势和生存差异令人惊讶地保持着（p = 0.011）（图5D-F）。糖酵解标志基因的表达仅在PDAC患者中与病理等级相关。在GSE62452队列中也观察到类似的结果。
此外，为了研究C3亚型的突变和免疫特征，我们分析了TCGA数据库中146名患者的突变数据。从突变的数量来看，C3亚型表现出最高的突变数量（图5G）。具体来说，C3亚型突变率最高的前10个基因是KRAS、TP53、CKDN2A、CKDN2A-DT、CKDN2B、SMAD4、MTAP、RN7SL151P、KLHL9和IFNE，这反映了C3亚型的突变负担过重。此外，我们根据细胞的已知功能，将免疫细胞分为抗肿瘤和促肿瘤。使用免疫细胞的基因特征来计算免疫浸润分数。我们观察到多种抗肿瘤免疫细胞的C3亚型明显减少，包括活化的CD8 T细胞、效应记忆CD4 T细胞、效应记忆CD8 T细胞、1型T辅助细胞和17型T辅助细胞。值得注意的是，我们观察到C3亚型中促肿瘤免疫细胞的浸润水平很低，包括未成熟的树突状细胞、巨噬细胞、骨髓源性抑制细胞（MDSCs）和浆细胞树突状细胞（图5H）。这一发现表明，几乎所有的免疫细胞在C3亚型中都受到抑制
图5. 高度恶性细胞亚群患者的临床特征。(A) 基于TCGA队列的五个糖酵解标志物基因表达值的UMAP图，并分为C1、C2和C3组。(B）TCGA队列中C1组（低表达）、C2组（中位表达）和C3组（高表达）的五个糖酵解标志基因的表达水平。(C) TCGA队列中C1、C2和C3组的总生存率。(D) 基于GSE62452队列的五个糖酵解标记基因表达值的UMAP图，并分为C1、C2和C3群。(E) GSE62452队列中C1组（低表达）、C2组（中位表达）和C3组（高表达）的五个糖酵解标志基因的表达水平。(F) GSE62452队列中C1、C2和C3组的总生存率。(G) C1、C2和C3组中突变的数量比较。(H）C1、C2和C3组中抗肿瘤和促肿瘤免疫细胞的浸润分数。*P＜0.05；**P＜0.01；***P＜0.001；****p＜0.0001。
由糖酵解标记基因构建的GHS评分是一个独立的预后因素
为了评估5个糖酵解标志基因作为PDAC患者预后标志物的价值，通过使用上述公式，整合5个糖酵解标志基因，建立GHS评分。进行单变量和多变量的Cox回归分析，进一步检验GHS评分对其他临床病理变量的独立作用。我们的结果显示，GHS评分是PDAC患者的独立预后因素（HR = 13.67；95%置信区间（CI）1.46-127.66；P = 0.02；图6A，B）。
图6
此外，我们通过生存分析和ROC分析测试了GHS的可靠性。在TCGA队列中，Kaplan-Meier生存分析显示，GHS分数较高组的预后明显较差（P = 0.00041）。GHS评分在1年（曲线下面积[AUC]=0.68）、3年（AUC=0.57）和5年（AUC=0.66；图6C）后有明显的预后价值。同样，在GSE62452队列中也获得了同样可靠的预后价值。较高的GHS分数组表现出较低的生存率（P = 0.0063）。1年、3年和5年的时间依赖性ROC分析的GHS评分分别为0.56、0.74和0.94（图6D），表明GHS评分可用于预测PDAC患者的预后。
小结
在这项研究中，我们结合了 scRNA-seq、bulk RNA-seq 和临床数据进行生物信息学分析。在 PDAC 患者中发现了以高度糖酵解为主要特征的高度恶性细胞亚群。以前的研究大多使用批量 RNA-seq 进行全面的基因组分析。与经典的bulk RNA-seq 相比，scRNA-seq 可以准确区分bulk 组织中的不同细胞类型，拓宽了对肿瘤生物学的理解。我们揭示了PDAC中一个以高糖酵解速率为主要特征的高度恶性细胞亚群，并获得了5个与PDAC发育密切相关的生物学功能的糖酵解标记基因。这些基因可以作为标志物来确定PDAC患者的代谢亚型和预后。
