代谢重编程是肿瘤细胞的主要特征。由肿瘤微环境恶化驱动的葡萄糖、脂质和氨基酸代谢重编程与肿瘤的发生发展存在密切联系，近年来得到了越来越广泛的关注。今天小编为大家解析的这篇文章就是以肺腺癌（lung adenocarcinoma, LUAD）转移相关的代谢重编程为侧重点，基于生信信息学方法对公共数据库中的数据进行挖掘，进一步揭示其潜在的细胞间的变化和可能的致病机理。这篇文章于今年3月发表于Communications Biology(IF = 6.268)。

一、摘要

本研究的主要目的是探究与肺腺癌（lung adenocarcinoma, LUAD）转移相关的代谢重编程，揭示其潜在的细胞间的变化。基于TCGA数据库中394例LUAD患者的基因表达谱数据，作者鉴定出11个与转移相关的代谢基因，这些基因涉及在糖酵解和脂质代谢途径中。并且，研究通过无监督聚类算法定义了3个代谢重编程表型：MP-I、MP-II、MP-III。在TCGA-LUAD队列和6个GEO表达谱共计1235个样本中，MP-III型样本的糖解酶水平最高、脂质代谢水平最低，并表现出较高的转移能力和较低的生存率。分析结果表明，TP53和KEAP1基因的突变以及SETD2和PBRM1的缺失可能是导致MP-III型患者代谢重编程的主要因素。LUAD样本单细胞转录组数据进一步验证从正常到MP-II、MP-III，再到MP-I的代谢轨迹。细胞间通讯分析结果发现，在ANGPTL通路中MP-III型细胞群与内皮细胞和成纤维细胞具有独特的相互作用，而在VEGF通路中其与内皮细胞的相互作用更强。鉴于此，研究认为，糖酵解-脂质失衡模式能够导致与转移相关的代谢重编程表型。该研究为进一步了解致癌驱动因子和肿瘤微环境之间的相互作用提供参考，对于LUAD转移的治疗提供策略。

二、材料与方法

1. 数据来源

本研究中，作者从TCGA、GEO和CPTAC数据库共获取10个LUAD队列。两个转移性LUAD数据集（TCGA-LUAD和GSE11969）的患者类型包括远端转移、淋巴结转移和非转移性。其中，TCGA-LUAD队列共包括394例LUAD患者样本（25例远端转移、160例淋巴结转移以及209例非转移性样本）。此外，对于TCGA-LUAD和GSE11969这两个数据集，作者分别下载了390和385个突变和拷贝数变异（CNV）数据。GSE11969数据集共包括18个远端转移、15个淋巴结转移和35个非转移性LUAD样本的mRNA表达谱数据。用于生存分析的5个LUAD表达谱，纳入的标准为：手术切除后未接受辅助治疗的I期患者样本数要超过50例。

此外，作者选择了2个LUAD单细胞表达谱数据：GSE131907和GSE123902。在GSE131907数据集中，作者从11位tLung患者、10位脑转移患者和11个正常肺组织中共提取了208506个细胞。在GSE123902数据集中，作者从7位tLung患者的原发性组织中提取了18511个细胞。

2. 数据预处理

TCGA-LUAD数据集，测序平台为Illumina HiSeq 2000 platform，fpkm表达矩阵经log2(x+1)转换；对于单细胞表达谱（scRNA-seq）数据，测序平台为Illumina HiSeq 2500。将基因的Ensembl ID转换成symbol。芯片数据GSE31210、GSE50081和GSE68465，测序平台为Affymetrix；GSE42127数据集测序平台为Illumina、GSE13213测序平台为Agilent。体细胞突变数据从Illumina Genome Analyzer DNA Sequencing GAIIx平台获取，包含17821个非同义突变。拷贝数变异（CNV）数据使用GISTIC算法进行处理。

3. 差异分析、富集分析

差异基因的鉴定使用student’s t-test，基于超几何分布模型对差异基因进行富集分析，数据集来源于KEGG数据库（https://www.genome.jp/kegg/pathway.html）。

4. 代谢重编程表型的鉴定

首先，作者基于代谢相关基因的表达相关性使用层次聚类的方法将患者分成不同的簇（或组）。在聚类之前，去除掉表达量均为0的样本或细胞，并对表达矩阵进行归一化处理。接下来，作者通过公式计算MP-score以评估不同簇的患者之间糖酵解和脂代谢基因之间的表达差异。计算公式如下：

5. 相关性及生存分析

使用Wilcoxon秩和检验来评估不同簇患者之间缺氧评分、干性评分、增殖评分、免疫评分和肿瘤突变负荷（TMB）的差异。基因间的表达相似性、突变情况、拷贝数变异（CNV）和蛋白表达情况使用Pearson相关系数进行评估。生存时间被定义为从最初的手术切除到死亡或最近一次随访的时间。使用Kaplan–Meier法绘制生存曲线。使用单因素cox模型分析患者的生存状态与临床特性之间的相关性。

6. 单细胞转录组分析

基于前20个主成分，使用R package Seurat中的RunPCA和JackStraw函数来对上皮细胞进行无监督聚类。使用RunUMAP函数对聚类结果进行可视化。使用Featureplot函数对每个细胞中mRNA的平均表达水平进行可视化。使用SingleR对单细胞类型进行注释。使用R package CellChat分析细胞间通讯。配体-受体互作及其相关的信号通路从CellChat数据库（https://www.cellchat.org）下载。该方法是根据两个细胞群间的基因表达、信号因子和细胞百分比来推断配体-受体间的潜在相关作用强度。

三、结果

3.1 鉴定LUAD转移相关的代谢重编程表型

首先，基于TCGA数据库中的表达谱数据，使用Student’s t-test方法进行差异分析（筛选标准为：FDR < 0.05），在转移和非转移性LUAD样本间共鉴定出431个差异基因，其中，显著下调表达基因151个，显著上调表达基因280个。分析结果发现，上调表达基因显著富集到HIF-1信号通路，该通路是糖酵解途径中的重要通路。先前的研究发现，该通路能够通过激活糖酵解基因的转录来增强糖酵解作用，在癌症的发生发展过程中起关键作用。相比之下，下调表达基因显著富集到脂代谢相关的生物途径，包括甘油磷脂代谢（Glycerophospholipid metabolism）和醚脂代谢（Ether lipid metabolism）等。由此作者推测，糖代谢和脂代谢可能在诱导LUAD转移过程中扮演不同的角色（图1 a-c）。

因此，作者将显著富集到糖酵解/糖异生（Glycolysis/Gluconeogenesis）途径的8个下调基因（ALDOA, ENO1, GAPDH, GPI, LDHA, PGAM1, PGM2, TPI1）以及显著富集到甘油磷脂代谢（Glycerophospholipid metabolism）和醚脂代谢（Ether lipid metabolism）这两个脂代谢通路的3个上调基因（PLPP1, GPD1L, PLD3）提取出来，进行更进一步的研究。

首先，作者基于上述11个基因在全部LUAD患者样本中的表达模式进行无监督聚类。分析结果将全部LUAD患者划分成3个簇（图1 d），由于这些基因的表达模式与糖、脂代谢途径显著相关，因此作者将这些簇定义为代谢重编程表型（metabolic reprogramming phenotypes, MPs）。接下来，作者计算了每个簇的MP得分，并根据得分的由低到高将簇分别命名为：MP-I、MP-II、MP-III。通过构建混淆矩阵，作者发现，MP-I型样本的患者类型主要为转移性，MP-II型样本的患者类型主要为淋巴结转移性，而MP-III型样本的患者类型主要为远端转移性（图1 f）。此外，生存分析结果表明，MP-III患者的生存率显著低于MP-I和MP-II患者（图1 g）。

3.2 不同代谢表型的分子特性和临床性状差异

分析结果表明，MP-III型患者处于stage IV和stage III的比例显著高于MP-I和MP-II。此外研究发现，从MP-I ~ III，患者的缺氧评分、干性评分和增殖评分逐渐增加。MP-II型患者的免疫评分最高，而MP-III型患者的免疫评分最低，两者之间的免疫评分存在显著差异（p < 0.0001）,见图2。

3.3 LUAD代谢表型的验证

上述结果表明，MP-III型大多为未发生任何转移的早期患者。由此推断，这些患者可能会有较高的转移风险，进而导致较差的预后。因此，作者首先基于11个代谢基因的表达谱数据，对来自5个公共数据集的stage I患者进行层次聚类。分析结果发现，stage I患者也可聚类为3种表型（图3 a-e）。生存分析结果表明，在GSE31210、GSE50081、GSE13213、GSE42127数据集中，MP-III型的stage I患者生存率均显著低于MP-I和MP-II型。除此之外，作者发现，在全部五个数据集中，MP-I和MP-II型的生存差异并不显著，但是可以看出在GSE13213、GSE42127和GSE68465这3个数据集中，相较于MP-I型，MP-II型患者的生存率是偏低的，但差异并不显著。

此外，研究从CPTAC数据库中获取了110例LUAD患者的表达谱数据，从蛋白水平上评估11个代谢基因的表达模式。从分析结果中可以看到，基于11个代谢基因的蛋白表达谱，同样将110例LUAD患者聚成了3个表型。此外，研究观察到每个基因的mRNA表达水平和蛋白表达水平均呈显著正相关（图4）。

3.4 在单细胞水平上描述代谢表型的转录轨迹

基于单细胞数据集GSE131907，26436个上皮细胞被分成了14个簇，其中包括正常细胞簇（Cluster 1、11、13）、原代肿瘤细胞簇（Cluster 4、8、9、10）、脑转移瘤细胞簇（Cluster 0、2、3、5、6、7、12），见图5 a。图5 b表现为每个细胞中11个代谢基因的平均表达水平。接下来，作者基于11个基因的表达矩阵对样本进行聚类以确定上皮肿瘤细胞的代谢表型，分析结果表明，MP-III型患者的脑转移瘤细胞的比例显著高于MP-I和MP-II型（图5 e）。

接下来，作者根据11个代谢基因的mRNA表达变化对上皮细胞进行排序，以构建代谢轨迹。轨迹中的7个转录状态分别代表不同的细胞代谢轨迹。每个细胞状态中各表型所占比例见图5 f，可以很明显的看出，正常上皮细胞显著富集在state 1和state 4这两个状态，然而，处于state 4阶段的细胞11个代谢基因的表达水平普遍较低，这一点在TCGA数据集中并没有观察到。糖-脂质代谢失衡的代谢轨迹似乎主要是从部分正常上皮细胞开始的，主要在State 1末期分化成MP-I型肿瘤细胞，然后形成一个具有两个主要细胞命运的分支结构：Cell fate 1和Cell fate 2。通过追踪Cell fate 1的代谢轨迹，作者发现MP-I型肿瘤细胞逐渐分化成MP-III型（state 3）。此外，对Cell fate 2代谢轨迹的追踪结果表明，MP-II和MP-III型肿瘤细胞主要富集在State 6和State 7这两个独立的分支。

3.5 揭示代谢表型下潜在的细胞间通讯

在该部分分析中，作者从患者原发性肿瘤组织中提取6372个肿瘤上皮细胞、2373个干细胞和35506个免疫细胞。如图5 a所示，作者在10个细胞群两两间共鉴定出41个具显著意义的配体-受体对。正如所推断的那样，作者发现MP-III型细胞群表现出强烈的信号输出，其次为MP-II和MP-I型细胞群（图6 b）。

接下来，与3种代谢表型相关的重要配体-受体对被进一步划分为15个信号通路（图6 c）。各配体-受体对在信号通路中的相对贡献率以及在各细胞群中的作用如图6 b、c所示。此外，分析结果发现，在3个MP细胞群中，MP-II是EGF信号通路的主要受体，可能单独发挥作用，也可能是与髓细胞、MAST、NK和内皮细胞共同作用（图7 g-i）。

四、总结

在本研究中，分析结果表明在转移性LUAD样本中显著上调和下调表达的差异基因分别富集在糖酵解和脂代谢途径中，由此推测糖脂代谢的失衡可能会促进LUAD的转移。不足之处在于，（1）首先，本研究是基于11个代谢基因的mRNA表达谱来定义代谢表型的，因此，糖酵解和脂质代谢活性需要通过估计其代谢物来验证；（2）其次，在本研究中发现的细胞间相互作用应使用更多的方法来进行验证，因为它们是通过配体及其受体的mRNA表达预测的；（3）尽管MP-III群体与其他群体的相互作用在另一个scRNA seq数据集中得到了验证，但想要获得更可靠的结果仍需进一步的实验验证；（4）此外，其他改变肿瘤微环境的代谢重编程模式，如代谢物等，仍待进一步探究。总的来说，研究鉴定了可能与LUAD转移相关的三种代谢表型，发现这一代谢重编程背后的四种致癌事件及其与其他细胞间的相互作用（尤其是ANGPTL信号通路）为防止LUAD转移提供了候选治疗策略。