题目：线粒体能量代谢与食管鳞状细胞癌的免疫抑制性肿瘤微环境和不良预后有关
杂志：Computational and Structural Biotechnology Journal
影响因子：IF=6.0
发表时间：2023年8月
研究背景
线粒体能量代谢（MEM）的重编程是肿瘤发生和癌症进展的一个重要标志。目前，还没有研究对食管鳞状细胞癌（ESCC）肿瘤微环境（TME）中的线粒体能量代谢进行检测，相关的药物靶点也尚未确定。
数据来源
研究思路
文章首先对单细胞数据进行质量控制、降维聚类和注释。通过KEGG数据库中MEM相关途径，提取MEM相关基因。使用AUCell R软件包分析单细胞数据集，以筛选富集分数较高的高分细胞群中与MEM相关基因。对高分细胞群进行通路富集分析和推断细胞分化，并分析了细胞间通讯网络。又对关键的MEM相关基因进行差异分析和富集分析。基于关键基因对患者分型，并构建了预后模型。使用互补DNA微阵列，通过qPCR验证了15例ESCC患者的预后基因转录水平。接着研究MEM相关基因与免疫微环境的生物学关系。利用体细胞突变数据，评估MEM相关风险评分对遗传变异变化的影响。最后评估了基于MEM相关预后基因的ESCC患者的药物敏感性。
主要结果
1 单细胞数据揭示ESCC的细胞异质性
利用GSE145370对14例病例（7例癌组织和7例癌旁组织）的单细胞RNA测序数据进行了分析。按照质量控制程序，共获得102268个过滤细胞。采用SNN（尖峰神经网络）算法进行聚类分析，得到了25个最佳细胞聚类（图1A）。利用气泡图确定了25个最佳细胞群中的排名靠前的标记基因（图1B）。这些标记基因在相应的细胞群中高表达，表明了细胞群特异性表达模式和细胞类型。
使用SingleR软件包识别了每个细胞群中的细胞类型（图1C）。细胞亚群0、3、4、5、6、12、15和20中共有50113个细胞被注释为T细胞，占细胞总数的49.002%。亚群1、10、16、17和22中共有16280个细胞被注释为B细胞，占细胞总数的15.919%。亚群2、8、9和19被注释为单核细胞（21075个；20.608%）。亚群7、13、14、23和24共有10773个细胞被注释为自然杀伤（NK）细胞，占细胞总数的10.534%。此外，亚群11被注释为中性粒细胞（3426个；3.350%），亚群18被注释为树突状细胞（DCs）（413个；0.404%），亚群21被注释为平滑肌细胞（188个；0.184%）。这些结果显示了数据集中细胞类型的多样性，并揭示了细胞异质性的存在。用小提琴图来显示各种细胞类型中差异最大的2个基因（图1D）。结果显示，与在其他细胞中的表达相比，CD3D和CD3E在T细胞中高表达，IGHG1和IGKC在B细胞中高表达，LYZ和CST3在单核细胞中高表达，FGFBP2和GNLY在NK细胞中高表达，PI3和CXCL8在中性粒细胞中高表达，CCL17和CCL22在DCs中高表达，DCN和CFD在平滑肌细胞中高表达。此外，每个细胞的前2个差异基因的表达量都显示在热图上。根据小提琴图预测，结果相似（图1E）。堆叠直方图用于显示每个样本中细胞类型的比例（图1F）。肿瘤组织中NK细胞的数量低于非肿瘤组织。
2 MEM相关基因在ESCC TME不同细胞群中的表达
将鉴定出的599个MEM相关基因与不同细胞群的DEGs相交，得到了121个标记的MEM相关基因（图2A）。根据平均对折变化（avg_log2FC）对这121个MEM相关基因进行排序，以确定在不同细胞状态或细胞群中表现出显著差异表达的基因。与MEM相关的前20个基因包括以下17个基因：SPP1、IL1B、APOE、SOD2、SERPINE1、GNG11、PTGS2、APP、CYBB、PPIF、FBP1、EGR1、COX7A1、MEF2C、LRP1、TUBB6和GNG5。热图显示，这17个基因的表达水平在不同细胞类型中存在显著差异（图2B）。SPP1、IL1B、APOE、PTGS2、CYBB、PPIF和FBP1在单核细胞中高表达，而在T、B和NK细胞中低表达。SOD2和GNG5在单核细胞、中性粒细胞和DC中高表达。此外，相关热图显示了17个MEM相关基因之间的相关性（图2C）。发现PTGS2与IL1B、IL1B与SOD2、IL1B与PPIF、PTGS2与SOD2以及SPP1与APOE之间存在显著相关性。
3 富集表达MEM相关标记基因的细胞亚群分析
使用AUCell来探索MEM相关标记基因的富集分数，并将富集分数高于0.086的细胞定义为“高分细胞”。结果得到了940个高分细胞，并用t-SNE观察了高分细胞群在所有细胞中的分布（图3A）。然后，使用SNN算法进行聚类分析，并根据PCA的降维结果使用t-SNE对高分细胞亚群进行聚类（图3B）。SingleR软件包用于识别不同的细胞类型并找到它们的DEGs。高分细胞群包括巨噬细胞（613个，65.213%）、T细胞（161个，17.128%）、上皮细胞（143个，15.213%）和B细胞（23个，2.447%）（图3B）。通过绘制小提琴图（图3C），确定了不同细胞类型中前2个差异标记基因。使用GO（图3D）和KEGG（图3E）分析软件包，利用clusterProfiler R软件包对不同细胞间的DEGs进行了富集。GO分析结果显示，高分细胞群主要富集的生物学过程与细胞活化有关，如细胞活化的正向调控、T细胞活化的调控和白细胞活化的正向调控。主要富集的细胞成分与核糖体有关，如细胞膜核糖体、核糖体亚基和核糖体。最后，主要富集的分子功能与免疫有关，如外源蛋白结合、主要组织相容性复合体（MHC）II类蛋白复合体和肽抗原结合（图3D）。此外，KEGG通路分析还显示，高分细胞群的DEGs主要富集于炎症、免疫和其他相关信号通路，如Th17细胞分化、Th1和Th2细胞分化以及T细胞受体信号通路（图3E）。此外，热图还显示了高分细胞群中17个MEM相关标记基因avg_log2FC top20的表达情况（图3F）。如图3F所示，巨噬细胞中17个MEM相关基因的表达量明显高于其他细胞类型。
4 高分细胞群的伪时间分析
使用Monocle R对高分细胞群进行伪时间分析，构建伪时间轨迹图。伪时间图从细胞群过程（图4A）和伪时间过程（图4B）两个方面着色，不同的分支沿伪时间序列方向传递。使用BEAM函数检验了不同支系之间与MEM相关的标记基因的差异，突出显示了两个支系之间存在差异的65个基因。这些结果在动态热图上可视化（图4C）。从第一个节点到根部是前分支，即与巨噬细胞相对应的状态；细胞命运1是与T细胞相对应的状态，细胞命运2包含与上皮细胞和巨噬细胞相对应的状态。考虑到其表达水平，根据其在不同细胞状态下表达模式的相似性，将差异标志物MEM相关基因分为六个不同的类别。这不仅显示了细胞群的过程，也显示了伪时的过程。通过伪时学的方向，观察到细胞群经历了不同的分支。这些结果提供了准时间过程中细胞状态变化和MEM动态的信息。
图4B描述了INV-H表型特有的前30条明显富集的通路和相关的MR。这些通路包括免疫系统（R-HSA-168256）、Toll样受体信号通路调控（WP1449）、II型干扰素信号传导（IFNG）（WP619）、纤维蛋白补体受体3信号通路（WP4136）、免疫系统中的细胞因子信号转导（R-HSA-1280215）、肿瘤微环境中免疫细胞与microRNAs的相互作用（WP4559）、结直肠癌中上皮到间质的转化（WP4239）、TGF-β信号通路（WP366）等。这些通路中的每一条都包含至少三种与INV-H表型（较差的OS）特异的不同MRs，从而表明这些MRs的较高活性和这些通路的富集不利于十种癌症中被归类为INV-H的患者的生存。
然后，对前30条通路进行了聚类，方法是利用2条富集通路所涉及的MR之间的重叠程度来估计富集通路集的相似性。得到相似性矩阵后，使用光谱聚类对通路进行聚类，将通路区分为有内聚力的组别（INV-H表型为5组）。通路按其所属的聚类进行颜色编码，与特定通路相关的MRs集被描述为邻接矩阵。观察到，显著富集的top通路大多是炎症（载脂蛋白E和miR-146在炎症中的作用、细胞因子和炎症反应）、免疫抑制（TGF-β信号通路、T细胞极化）、先天性免疫信号转导（Toll样受体信号转导、II型干扰素信号转导、白细胞介素信号转导）和转移前体（上皮细胞向间质转化）的标志性通路，如图4B所示。证明了INV-H表型及其在十种相关癌症中更差的生存预后。
针对INV-L表型的65个MRs也进行了类似的INV-H表型分析。在过度表达分析中，检测到70个GO术语和6条通路被显著富集。排名靠前的GO术语包括核酸代谢过程、杂环代谢过程、含核碱基化合物代谢过程、细胞芳香化合物代谢过程、基因表达等。它们主要与细胞内的代谢过程有关。
与INV-L表型相关的六条富集途径包括基因表达（转录）（R-HSA-74160）、RNA聚合酶II转录（R-HSA-73857）、通用转录途径（R-HSA-212436）、rRNA表达的负表观遗传调控（R-HSA-5250941）、B-WICH复合物正向调控rRNA表达（R-HSA-5250924）和染色质修饰酶（R-HSA-3247509）。如图4C所示，富集的通路分为两组，主要与转录调控有关。从图4C中，观察到x轴上的最大比值达到了约0.125，表明一条通路中最多只有八分之一的基因被过表达。在10种癌症中，这些通路的富集和相关MRs的较高活性有利于INV低组患者的生存。
综上所述，这些结果突显了TGF-β、Toll样受体信号通路、上皮细胞向间质转化通路等候选通路在多种癌症类型的高侵袭性癌症（INV-H表型）中的显著富集，可作为靶向通路来提高癌症侵袭性的生存率。
5 高分细胞群之间的细胞通讯
首先，使用热图（图5A）直观显示所有细胞群之间的相互作用强度。纵坐标上的每种细胞都是信号发送者。它们与横轴上其他细胞类型的交流强度显示在顶部和右侧两列。每个直方图代表相应细胞类型的交流强度之和。T细胞作为配体细胞或受体细胞与其他类型细胞的通讯强度最大。
此外，还分析了高分细胞群中不同类型细胞之间的相互作用。通讯网络图显示了不同类型细胞之间的相互作用数量（图5B）和相互作用信号的强度（图5C）。巨噬细胞之间有68对相互作用，作为配体细胞的巨噬细胞与T细胞之间有71对相互作用，作为受体细胞的巨噬细胞与T细胞之间有49对相互作用。相互作用信号强度分析结果表明，巨噬细胞本身之间的相互作用强度最强，作为配体细胞的巨噬细胞与T细胞之间的相互作用强度最强，作为受体细胞的巨噬细胞与T细胞之间的相互作用强度最强。此外，高分细胞群中有319对配体-受体。本研究采用热图显示了沟通概率最高的10对受体和配体（图5D），包括“SPP1-CD44”、“MIF-(CD74+CD44)”、“HLA-DRA-CD4”、“MIF-(CD74+CXCR4)”、“HLA-DPA1-CD4”、“HLA-DPB1-CD4”、“HLA-B-CD8A”、“SPP1-(ITGA5+ITGB1)”、“SPP1-(ITGAV+ITGB1)”和“HLA-DRB1-CD4”。
6 基于转录组数据中MEM相关标记基因的差异分析
首先，利用GSE111011数据集探讨了121个MEM相关标记基因在ESCC和癌旁组织中的表达，其中36个基因在ESCC和邻近组织中存在差异。如火山图所示，在36个DEGs中，18个基因上调，18个基因下调（图6A）。随后，利用GSE164158数据集探讨了121个MEM相关基因在ESCC和癌旁组织中的表达情况，其中28个基因在ESCC和癌旁组织中存在差异。火山图显示，在28个DEGs中，9个基因上调，19个基因下调（图6B）。7个基因在两个数据集中上调（图6C），包括TUBB、PSMB1、PSMA7、PSMA6、BID、APOE和SPP1；8个基因在两个数据集中下调（图6D），包括NDUFB2、VDAC2、NDUFA13、NDUFA11、GNAQ、MAP1LC3A、APPL1和KLF2。总体而言，热图显示，在GSE111011（P＜0.05，图6E）和GSE164158（图6F）数据集中，有15个基因在ESCC和癌旁组织之间存在显著差异。
为了探索其中的分子机制，使用clusterProfiler R软件包对15个关键的MEM相关基因进行了GO富集分析。在GO富集结果中，每个项根据P值从高到低排列，即富集越显著（P值越小），项越高（图6G）。因此，15个关键的MEM相关基因主要富集在与细胞有氧呼吸相关的生物过程中，如“氧化磷酸化”、“NADH脱氢酶复合物组装”和“线粒体呼吸链复合物I组装”。其主要富集的细胞成分与有氧呼吸有关，如“线粒体呼吸链复合物I”、“NADH脱氢酶复合物”和“呼吸链复合物I”。KEGG分析结果显示，15个MEM相关基因的主要信号通路与新陈代谢有关，如“氧化磷酸化”、“铁变态反应”和“胆固醇代谢”（图6H）。
7 基于关键MEM基因对TCGA-ESCC肿瘤样本进行分子分型
基于15个关键MEM相关基因和TCGA-ESCC样本普通转录组测序的FPKM表达数据，ConsensusClusterPlus被用来对ESCC样本进行分类。综合考虑一致性矩阵热图（图S2A）、累积分布曲线（图S2B）和delta面积曲线（图S2C），确定聚类数目为2。在两个亚型中，SPP1、APOE、BID、NDUFB2和NEUFA13的表达水平均有显著差异（图S2D）。根据这15个MEM相关基因进行了t-SNE分析，结果显示两种ESCC亚型聚成了两个不同的群（图S2E）。
此外，ssGSEA显示了16种免疫细胞的差异（图S2F）和13种免疫细胞在不同分子亚型中的富集途径（图S2G）。巨噬细胞在群集1中的含量高于群集2。在其他免疫细胞和免疫相关通路中，两个亚型的富集活性没有显著差异。
8 构建预后基因模型
为了评估15个关键MEM相关基因与ESCC患者预后之间的相关性，我们从92例TCGA-ESCC患者中随机挑选了56例组成训练集，并进行了单变量COX回归分析。结果发现，7个MEM相关基因（BID、MAP1LC3A、KLF2、APOE、APPL1、NDUFA13和TUBB）对ESCC预后有较大影响。然后，利用回归分析建立了拟合模型（图S3A），并通过交叉验证分析确定拟合模型的最佳lambda值为5（图S3B）。由此确定了5个与MEM相关的基因。最后，利用这5个基因进行了Cox回归分析，并构建了一个风险模型。
采用中位风险评分将患者分为高危和低危两组。Kaplan-Meier曲线显示，高风险患者的总生存率低于低风险患者（图S3C）。为了验证预测的准确性，还进行了ROC分析，结果表明风险评分能很好地预测ESCC患者的总生存率（OS）。在1年和2年时，曲线下面积（AUC）分别为0.78和0.86（图S3D）。最后，剩下的36例TCGA-ESCC患者被纳入验证集。同样，根据Kaplan-Meier曲线，高危组患者的总生存率低于低危组患者（图S3E）。根据时间ROC，风险评分对验证集患者OS的预测与训练集相似，1年和2年OS的AUC分别为0.77和0.88（图S3F）。此外，在92例TCGA-ESCC患者中，高危患者的OS低于低危患者（图S3G），风险评分预测OS的AUC也相对较高（图S3H）。
利用15例ESCC患者的cDNA芯片，通过RT-qPCR进一步验证了这四个预后基因的转录水平。ESCC组织中APOE和APPL1的表达水平高于正常食管组织，而ESCC组织中MAP1LC3A的表达水平较低（图7）。两种组织中NDUFA13的表达水平没有明显差异。
9 免疫特征的相关性分析
为了研究MEM基因与免疫微环境之间的生物学关系，使用CIBERSORT脚本分析了病变样本中的免疫细胞浸润情况。结果得出了22种免疫细胞的丰度值。值得注意的是，预后基因与一些免疫细胞密切相关（图8A）。例如，浆细胞的丰度与MAP1LC3A呈正相关，而T细胞CD8的丰度与APOE和NDUFA13呈正相关。此外，M2巨噬细胞与MAP1LC3A呈正相关。静息记忆CD4+T细胞的丰度与APOE呈负相关，与APPL1呈正相关。活化的DC与APOE呈负相关。此外，还研究了风险评分与免疫细胞丰度之间的相关性（图8B）。结果发现，风险评分与直流电活化和静息记忆CD4+T细胞的呈负相关。风险评分还与M2巨噬细胞和γ-δT细胞呈正相关。利用方框图还可以观察到高风险组和低风险组之间免疫细胞的差异（图8C）。低风险样本中静息记忆CD4+T细胞和活化DC的免疫浸润评分明显高于高风险样本，高风险样本中M2巨噬细胞的免疫浸润评分明显高于低风险样本（图8C）。
10 风险评分对TCGA-ESCC患者基因组变化的影响
为了评估MEM相关风险评分对遗传变异（包括单核苷酸多态性（SNP）和拷贝数变异（CNV））变化的影响，下载了体细胞突变数据。首先，使用maftools软件包显示TCGA-ESCC高危组和低危组中前20个基因的突变情况、包括TP53、TTN、NFE2L2、CSMD3、KMT2D、FLG、MUC16、NOTCH1、PIK3CA、ZNF750、DNAH5、HYDIN、NOTCH3、OBSCN、PCLO、SYNE2、DST、FAT3、PKHD1L1、RYR2（图9A）。在高危组和低危组中，PCLO突变的存在有显著差异。此外，还发现了高风险组和低风险组之间存在突变差异的基因（图9B）。结果显示，除PCLO外，DYNC2H1和KIF26B在高危组和低危组之间的突变样本数也有显著差异。
此外，还使用maftools软件包对TCGA-ESCC样本的突变特征进行了可视化分析（图9C）。结果显示，突变特征1、16、13、3、6、2、18和24在肿瘤样本中的出现频率高于正常样本。
11 基于MEM相关预后基因评估ESCC患者的药物敏感性
为了探讨MEM相关预后基因是否可用于评估ESCC对不同治疗药物的敏感性，从CellMiner数据库下载了60个细胞系的基因表达数据和24360种药物的IC50值。剔除信息不全的基因和药物后，分析了四个与MEM相关的预后基因与62种药物的IC50值之间的相关性。其中31种药物的IC50值与预后基因之间存在明显的相关性。APOE的表达与低霉素的IC50值呈正相关，而与7种药物的IC50值呈负相关，包括巴曲林和培福星。APPL1的表达与PX-316的IC50值呈正相关。MAP1LC3A的表达与冈田酸的IC50值呈正相关，与6-巯基嘌呤和6-硫鸟嘌呤等6种药物的IC50值呈负相关。NDUFA13的表达与15种药物的IC50呈正相关，包括羟基脲、美法兰、哌布溴铵、卡铂和氯霉素。预后基因与某些药物的相关性见图S6。
考虑到免疫疗法在治疗肿瘤中的重要性，使用TIDE算法评估了高风险和低风险患者对免疫疗法的敏感性。高风险组的TIDE评分高于低风险组（图10A）。此外，TIDE评分与风险评分也呈正相关（图10B），这表明免疫疗法对高危患者的效果较差。尽管排他性评分通常反映了免疫逃逸反应的强度，但高风险组和低风险组的排他性评分没有显著差异（图10C），而且与风险评分没有线性相关（图10D）。
为了进一步评估癌症患者免疫治疗的效果，估算了常见免疫检查点PD1、PD-L1、PD-L2、CTLA4与风险评分之间的相关性（图10E-L）。结果显示，高风险组PD-1的表达高于低风险组（图10E），且与风险评分呈正相关（图10F）。高危组与低危组的PD-L1表达量差异无统计学意义（图10G），但与风险评分呈正相关（图10H）。高危组PD-L2的表达高于低危组（图10I），且与风险评分呈正相关（图10J）。高风险组CTLA4的表达高于低风险组（图10K），且与风险评分呈正相关（图10L）。这些结果表明，针对这四个免疫检查点的免疫疗法可能对高危患者更敏感。
文章小结