|
摘要 胃癌(GC)是全球第四大癌症相关死亡原因。晚期胃癌患者预后较差，生存期缩短。寻找新的预测GC预后的生物标志物是迫切需要的。线粒体自噬是对受损线粒体的选择性降解，以维持细胞内稳态，已被证明具有促肿瘤和抗肿瘤作用。本研究结合单细胞测序数据和转录组学筛选与GC进展相关的线粒体相关基因(MRGs)，并分析其临床价值。RT-qPCR和免疫化学(IHC)进一步验证了基因表达谱。将单细胞测序数据与MRGs交叉后，共鉴定出18个DE-MRGs。MRG评分较高的细胞主要分布在上皮细胞簇中。上皮细胞与其他类型细胞之间的细胞间通讯显著上调。作者建立并验证了一个可靠的基于DE-MRGs (GABARAPL2和CDC37)和传统临床病理参数的nomogram模型。GABARAPL2和CDC37表现出不同的免疫浸润状态。鉴于hub基因与免疫检查点之间的显著相关性，靶向GC中的MRGs可能为接受免疫治疗的患者补充更多的益处。总之，GABARAPL2和CDC37可能是GC的预后生物标志物和候选治疗靶点。
图 1 研究的工作流程 本文思路可重复也可扩展，欢迎联系小编咨询分析 研究结果 胃癌中癌及癌旁组织的单细胞图谱 7个GC样本参与了单细胞测序分析。经过严格的质量控制，保留39475个高质量细胞用于进一步分析，其中26318个来自肿瘤组织，13157个来自癌旁组织。使用Seurat包对细胞进行聚类和降维。随后，作者采用t-SNE技术将高维scRNA-seq数据可视化，并成功地将细胞分为六个簇(图2A)，分别标注为上皮细胞、成纤维细胞、内皮细胞、淋巴细胞、髓系细胞和浆细胞。在探索所有6个细胞簇的分布时，我们注意到每个样本的细胞比例是不同的。上皮细胞在三个样本(5846-t、6342-t和6592-t)中占比最高，其余样本中主要细胞类型为淋巴细胞(图2B)。此外，这些细胞亚类用12个公认的细胞标记进行了验证，如图2C-D所示。每个聚类的前20个显著差异表达基因被绘制在热图上(图2E)。
图 2 GC中的细胞类型分类 线粒体相关差异表达基因(DE-MRGs)的鉴定与评分 为了阐明MRG与不同细胞群之间的相关性，作者将1588个差异表达细胞标记(log FC > 0.25)与83个MRGs取交集。作者获得了18个DE-MRGs(图3A)。作者分析了18个DE-MRGs之间的相关性(图3B)。TOMMF7与UBA52 (R = 0.18, P < 0.001)、RPS27 (R = 0.21, P < 0.001)呈正相关。JUN与UBB (R = 0.20, P < 0.001)、UBC (R = 0.24, P < 0.001)呈显著正相关。同样，UBA52与RPS27A呈正相关，具有统计学意义。同时，分析了18个DE-MRGs在各细胞类型中的表达情况。ATF4、JUN、TOMM7、UBA52、RPS27A、UBB、UBC和SQSTM1广泛分布于各种细胞簇中。MAP1LC3B、FIS1、VDAC1和CDC37主要在上皮细胞中表达，而FOXP3在成纤维细胞中相对高表达。这些基因可能在不同细胞类型之间协同作用，共同促进胃癌的发生。
图 3 基于单细胞RNA测序的线粒体相关基因差异表达鉴定 用AUCell R包对所有单细胞进行评分。一般来说，DE-MRGs高表达的细胞具有较高的AUC值。当AUC值设置为0.5时，4590个细胞的AUC值(MRG评分)高于阈值(图4A)。进一步分析发现，这些高MRG评分细胞主要集中在上皮细胞(图4B)。对MRG得分高的细胞群体中的差异表达基因( DEGs )进行GO和KEGG分析。GO结果显示，DE-MRGs主要参与细胞质转译、ATP代谢过程、氧化磷酸化、有氧呼吸、细胞呼吸、ATP合成偶联电子传递、线粒体ATP合成偶联电子传递、需氧电子传递链、有机化合物氧化产生能量等。前体代谢产物和能量的生成(图4C)。KEGG结果显示DE-MRGs主要与核糖体、冠状病毒病-COVID-19、氧化磷酸化、帕金森病、化学致癌-活性氧、糖尿病心肌病、朊病毒病、肌萎缩性侧索硬化症、亨廷顿病、产热有关(图4D)。
图 4 18个DE-MRGs评分 上皮细胞群的细分及轨迹分析 作者分离了9763个上皮细胞，并根据不同的免疫细胞标记将其分为8个亚组，包括以KRT18标记的内皮细胞1 ( Epi1 )，以MUC5AC标记的内皮细胞2 ( Epi2 )，以CEACAM5标记的内皮细胞3 ( Epi3 )，以ITLN1标记的杯状细胞，以FABP1标记的肠细胞，以MUC6标记的胃窦基底腺黏液细胞( GMC )，以PGA4标记的主细胞和以CHGA (图5A-B)标记的肠内分泌细胞。随后，我们分析了18个DE-MRGs在不同上皮细胞亚群(图5C-T)中的表达情况。结果表明，ATF4、JUN、TOMM7、UBA52、RPS27A、UBB、UBC和SQSTM1广泛分布于各种上皮细胞中。相反，HIF1A、BNIP3L、FOXO3、CITED2和RAB7A在任何亚组中均无明显表达。值得注意的是，GABARAPL2基因、MAP1LC3B、FOXO5、FIS1、VDAC1和CDC37在Epi1细胞中高表达。
图 5 上皮细胞的细分及轨迹分布 作者还对所有上皮细胞的细胞分化轨迹进行了模拟分析。如图5U所示，将七种不同的状态用不同的颜色标记出来。蓝色越深，分化越早，表明上皮细胞随时间由右向左分化(图5W)。Epi1细胞主要出现在分化初期，可向Epi2和Epi3细胞分化(图5V)。18个DE-MRGs的表达随分化状态的不同而变化(图6)，沿分化轨迹，JUN、UBB、UBC、SQSTM1、GABARAPL2、MAP1LC3B、HIF1A、BNIP3L、FOXO3、CITED2、RAB7A等DE-MRGs水平在过渡过程中逐渐升高。而UBA52、RPS27A、VDAC1和CDC37则呈现相反的趋势。同时，ATF4、TOMMF7、FIS1在轨迹分析中无明显趋势。
图 6 18个DE-MRGs在分化伪时间轨迹中的相对表达谱。 细胞通讯分析 CellChat用于描述复杂的细胞间通信网络。在GC样本中，上皮细胞与成纤维细胞、内皮细胞和免疫细胞（包括淋巴细胞、髓系免疫细胞和浆细胞）的相互作用增强；成纤维细胞与其他类型细胞的相互作用一般减弱(图7A)。作者分析了上皮细胞与其他类型细胞之间的配体-受体相互作用网络。结果显示MIF-(CD74+ CXCR4)和MIF-(CD74+ CD44)在胃癌细胞中明显增强(图7B)。进一步分析亚群间的细胞-细胞通讯也显示，在GC组织中有多对配体-受体对被激活，如MIF-(CD74)+ Cxcr4)， MIF -(cd74+ MDK-SDC4(图 7C)。总之，肿瘤组织中的细胞-细胞通讯模式与正常细胞不同，这可能与胃癌的发生发展有关。
图 7 细胞间通讯分析 正常和肿瘤组织中MRGs表达的差异 作者比较了450名患者的Bulk转录组测序数据。从TCGA-STAD数据集中，分析肿瘤组织与正常组织之间的差异表达基因。共筛选出4393个DEGs，其中上调的基因1949个，下调的基因2444个。通过火山图和热图(图8A、B)，在正常和肿瘤组织中共鉴定出10个MRGs。结果表明，E2F1在GC中表达量较高，而CITED、PINK1、PRKN、JUN、GABARAPL1、ATG9B、MAPK10、BNIP3、MAP1LC3C在GC中表达量较低。将单细胞测序数据与bulk-RNA测序数据交叉后，最终获得了233个基因(图8C)。根据GO功能富集分析，233个基因主要与脂多糖响应、细菌源分子响应、伤口愈合、金属离子响应、机械刺激响应、活性氧响应、肽酶活性调节相关(图8D)。
图 8 基于转录组学鉴定差异表达基因(DEGs)。 一致性聚类分析 为了探索18个DE-MRGs表达与GC亚型之间的联系，作者对来自TCGA队列的415例GC患者进行了一致性聚类分析。累积分布函数(CDF)曲线表示两两相似度得分小于或等于给定阈值的比例，用于确定最优的聚类数量。根据CDF, k = 2时得到最优除法(图9A-C)。415例胃癌患者可以很好地分为两个簇，簇1为205例，簇2为210例。我们比较了两个簇间18个DE-MRGs的表达(图9D)。11个基因(ATF4、JUN、TOMM7、UBA52、RPS27A、UBB、GABARAPL2、CITED2、FIS1、VDAC1、CDC37)在cluster-1中高表达，2个基因(HIF1A、FOXO3)在cluster-2中过表达。此外，两组间5个基因(UBC、BNIP3L、MAP1LC3B、RAB7A和SQSTM1)的表达差异无统计学意义。主成分分析显示了两个聚类之间的差异(图9E)。ESTIMATE包用于计算所有415个GC样本的免疫评分。结果表明，两组之间的免疫评分无统计学意义(图9F)。使用ssGSEA算法来分析两类免疫细胞的浸润情况。作者发现簇 2显著富集于中央记忆CD8+ T细胞，效应记忆CD8+ T细胞，记忆B细胞，中性粒细胞，调节性T细胞，滤泡辅助性T细胞，17型辅助性T细胞，效应记忆CD4+ T细胞，单核细胞，2型辅助性T细胞(图9G )。簇1中活化的CD8 + T细胞免疫浸润显著高于簇2。
图 9 一致性聚类分析 基于GABARAPL2和CDC37的预后风险模型 作者从TCGA数据集中提取了378例具有完整生存信息的GC患者的临床特征。单变量Cox回归分析筛选DE-MRGs, 5个基因(GABARAPL2、MAP1LC3B、HIF1A、CITED2、CDC37)与GC预后显著相关。作者进行LASSO回归来缩小候选基因的范围，最后保留两个基因(GABARAPL2和CDC37)，构建基于最优值(图10 A - B)的预后预测模型。结合临床资料，采用以下公式计算风险评分: risk score = (expression of GABARAPL2*(0.410) + expression of CDC37*(-0.455) + age*(0.033) + stage*(0.556)。在TCGA-STAD队列中，高危患者的生存期比低危患者短(P < 0.0001)(图10C)。ROC曲线在1年、3年和5年的AUC分别为0.7、0.66和0.68，表明作者的预后模型在预测OS方面表现良好(图10D)。在两个验证数据集(GSE62254和GSE15459)中也可以得到类似的结果。高危组和低危组的K-M曲线有显著性差异(图10E, G)，在GSE62254队列中，ROC曲线显示出良好的预测能力，1年的AUC为0.78,3年为0.75,5年为0.72(图- 10f)。在GSE15459队列中，1年、3年、5年的AUC值分别为0.76、0.81、0.83(图10H)。 此外，作者使用森林图来显示单变量和多变量Cox回归结果(图10I-J)。有意义的变量包括GABARAPL2、CDC37、年龄和分期。用加权风险评分和上述变量建立了一个预后nomogram模型(图10K)。总评分采用每个风险的定量评分计算。校准曲线用于评价该列线图预测3年生存率的准确性，显示了良好的预测能力(图10L)。
图 10 预后风险模型的构建与验证 GABARAPL2和CDC37的预后价值验证 图11A验证的两个预后生物标志物的RNA表达水平(GABARAPL2, CDC37)。RT-qPCR结果证实，与癌旁正常组织相比，肿瘤组织中GABARAPL2 (P = 0.0028)和CDC37 (P = 0.0003) mRNA表达明显上调。免疫组化结果(图11B)在蛋白质水平上支持了这一发现。在62个组织切片中，GABARAPL2蛋白在与中位OS较短相关的GC组织中高表达(高vs低= 15.97 vs 20.7个月)。GC组织中也检测到CDC37蛋白的上调，而其过表达则预示着更长的中位os和更好的预后(高vs低= 16.3个月vs 20.7个月)。图11C绘制了不同表达水平的GABARAPL2和CDC37的Kaplan-Meier曲线。这些发现进一步验证了TCGA-STAD队列转录组数据的结果。
图 11 免疫浸润分析 GABARAPL2和CDC37与免疫浸润的相关性 为了验证预后DE-MRGs ( GABARAPL2基因、CDC37 )与免疫的相关性，使用CIBERSORT探究了20个肿瘤浸润免疫细胞在高危组和低危组之间的差异。过滤Cibersort P ＜ 0.05，GABARAPL2基因与M2型巨噬细胞、静息树突状细胞、CD8 + T细胞呈正相关( Fig 11D )，CDC37与CD8+ T细胞、调节性T ( Treg )细胞、M2型巨噬细胞、静息NK细胞、M1型巨噬细胞呈正相关( Fig 11E )。免疫评分与CD8+ T细胞、记忆性B细胞、M1型巨噬细胞、活化的记忆性CD4+ T细胞、M2型巨噬细胞呈正相关(图11F )。随后，作者分析了高危组和低危组的免疫浸润差异(图11G)。高危组M2巨噬细胞浸润明显增加(P = 0.029)，而高危组M0巨噬细胞浸润相对较少(P = 0.015)。此外，作者还分析了GABARAPL2、CDC37与7个常见的免疫检查点基因(PDCD1、CD274、CTLA4、LAG3, CD276, ENTPD1, NT5E)(图11H)。GABARAPL2与CD247呈负相关，与NT5E呈正相关。CDC37与CD274、LAG3呈正相关，与ENTPD1、NT5E呈负相关。 基因突变分析 基于TCGA数据集的体细胞突变数据，作者使用' maftools '包对高危组和低危组患者进行单核苷酸突变分析。在高危组中，TTN ( 64 % )，MUC16 ( 37 % )，TP53 ( 37 % )，LRP1B ( 32 % )和SYNE1 ( 31 % )是前5位突变基因(图12A )，而在低危组中，前5位频繁突变的基因是TTN，TP53，MUC16，SYNE1，LRP1B (图12B )。两组肿瘤突变负荷( TMB )无统计学差异(图12C )。根据TMB的中位值，可将GC患者分为高TMB组和低TMB组。如图12D所示，高TMB组的存活率明显高于低TMB组( P = 0.018)。 基于TCGA的拷贝数变异(copy number variation, CNV)数据，作者分析了肿瘤的免疫浸润和功能障碍。与低危组相比，高危组Gistic得分更高，扩增和缺失突变更多(图12E, F)。同时，作者分析了TMB、微卫星不稳定性(MSI)、肿瘤免疫功能障碍和排除(TIDE)、免疫功能障碍、和免疫排斥，选择P < 0.05作为统计学意义阈值。作者观察到TMB与MSI (R = 0.723, P < 0.001，图12G)、TIDE与排除(R = 0.612, P < 0.001，图12H)之间呈正相关;功能障碍和TIDE (R = 0.219, P < 0.001，图12J)。相反，作者也观察到功能障碍与MSI呈负相关(R =， P < 0.001，图12I)。
图 12基因突变分析 结论 本研究结合了单细胞测序和转录组学，以识别有意义的DE-MRGs作为GC治疗的预后生物标志物和治疗靶点。开发了一个可靠的nomogram模型来预测胃癌患者1年、3年和5年的OS。RT-qPCR和IHC已经验证了基因表达谱。在其他GEO数据集中进一步证实了GABARAPL2和CDC37的预测价值。MRG评分较高的细胞主要分布在上皮细胞中。在GC中，上皮细胞与其他簇状细胞之间的通信显著上调。GABARAPL2和CDC37表现出不同的免疫浸润状态。鉴于hub基因和免疫检查点之间的显著相关性，靶向肿瘤细胞中的MRGs有望补充当前GC治疗中的免疫治疗。此外，GABARAPL2和CDC37的潜在机制还需要在体外和体内进行更多的研究。 |
|
|
