大家晚上好！今天小编和大家分享一篇23年5月发表在(IF:7.300)杂志的文章《Machine learning-based glycolysis-associated molecular classification reveals differences in prognosis, TME, and immunotherapy for colorectal cancer patients》。作者通过共识聚类以及机器学习等方法，系统性地分析了糖酵解相关分子在结直肠癌患者中的突变、生存以及肿瘤微环境之间的关联。对于糖酵解相关方向的分析思路有兴趣的老师可以与我们联系，方法适用于其他分子。
糖酵解是能量转化的关键形式，将葡萄糖转化为丙酮酸，最终产生乳酸并为肿瘤细胞提供氧非依赖性能量。即使在有氧环境中，肿瘤细胞仍然可以限制糖酵解的能量代谢，通常称为需氧糖酵解。通过减少细胞抑制和凋亡，糖酵解为肿瘤细胞生长创造条件。越来越多的证据表明，糖酵解与结直肠癌的发展有关。许多糖酵解相关基因，包括乳酸、GLUT1、丙酮酸激酶M2、甘油醛-3-磷酸脱氢酶、烯醇化酶-1、乳酸脱氢酶5和己糖激酶2，目前被发现在结直肠癌中上调。同时也有证据表明，糖酵解对于控制肿瘤细胞与肿瘤微环境（TME）之间的相互作用至关重要。在这项研究中，作者使用无监督聚类方法使用糖酵解相关的基因对每个CRC患者进行分类。然后讨论了每个糖酵解相关簇（GAC）的临床、基因组学、TME和富集途径的特征。同时还将每种GAC类型映射到单细胞数据，分析了每个细胞中每种GAC样亚型的分布和功能。基于上述发现，作者整合了单细胞和普通RNA-seq的结果，以开发一个模型来预测每个CRC患者的GAC并对其进行验证，并探索了其在治疗中的作用。1，k-means方法和ConsensusClusterPlus包进行了共识聚类；2，survival and survminer包进行Kaplan-Meier生存分析；4，CMScaller包确定CRC患者的CMS归类；5，Seurat包对scRNA-seq数据进行处理；6，metafor包中的rma功能确定重要的细胞亚型；8，glmnet，randomForest，Boruta，XGBoost，e1071和caret包筛选重要基因；9，pROC和caret包计算ROC曲线下的AUC面积；10，pRRophetic包计算药物的半抑制浓度（IC50）值。
研究技术路线如图1。从MSigDB数据库（HALLMARK_GLCOLYSIS）收集198个糖酵解相关基因进行进一步研究。采用GO富集分析发现，这些基因主要富集到碳水化合物和衍生物（图2A)。PCA降维可视化了TCGA COAD/READ数据集中正常组织和肿瘤组织之间糖酵解相关基因的分布，结果表明糖酵解相关基因可有效将肿瘤与正常组织分离（图2B)。为了研究这些基因在肿瘤发生中发挥的潜在作用，我们使用共识聚类对TCGA和GEO数据集中CRC患者进行分类，并确定了三个糖酵解相关簇（GAC）如图2C。T-SNE分析显示每种GAC的分布呈现显著异质性（图2D)。
在GAC聚类的基础上，我们讨论了每种GAC相应的临床特征。三种GAC的总生存期（OS）和无病生存期（DFS）的生存曲线表明，GAC2的预后较差（TCGA-OS: p = 0.05; TCGA-DFS: p = 0.03; GEO-OS: p< 0.001; GEO-DFS: p = 0.011, log-rank test）（图2E–H)。临床分期与GAC的关系表明，GAC2的比例随分期的进展而增加，而GAC3的比例呈现相反趋势，这与GAC2的不良预后一致（图3C）。根据图3A，GAC1的Kras突变比例较大，而GAC3的Braf突变和MSI-H状态比例较高。此外，GAC2的年龄分布最年轻。因此以上结果表明GACs的临床特征具有异质性：1.GAC2预后较差，分期较晚，呈年轻化趋势；2.GAC3预后较好，分期较早，突变概率较高。3. GAC1在各种临床特征中处于中间位置。
利用“GSVA”算法，我们选择了四种表型（免疫相关、代谢相关、核苷和RNA相关和基质相关）进行热图展示（图2I)。这些途径包含“C2.cp.kegg”, “C2.cp.reactome”和“C2.go.bp”。通过“limma”包的分析，确定了差异途径。结果表明GAC1富集在与核苷相连的途径中。GAC2与基质和免疫功能有关。免疫和代谢性状存在于GAC3中。在此基础上，我们讨论了GAC与CMS的关系并创建了热图来显示GAC和CMS之间的关系（图3A)。令人惊讶的是，我们发现GAC1和GAC2分别与CMS2和CMS4重叠。而在GAC3患者中，CMS1和CMS3占据相等且最大的比例（图3A、B)。随后我们使用“CMScaller”包的经典途径进一步验证了这些发现。与CMS4相关的“EMT”和“TGF-Beta”途径在GAC2中富集（图3A）。我们用ssGSEA对糖酵解途径进行了量化，发现在GAC3中存在高表达（图3A、D）。此外，糖酵解与WNT、MSS、LGR5干细胞、HNF4A、CMS2和CMS4负相关，而与MSI、脂肪酸、DNA修复、CMS1、CMS3和细胞周期正相关（图3F)。GSE39582数据集的冲积图进一步证明GAC1指向CIN免疫下降和CIN Wnt-up，GAC3指向dMMR和KRASm，GAC2指向其余（图3E)。最后，GAC1共享了CMS2的大部分功能。GAC2与CMS4最相似，GAC3具有CMS1和CMS3的双重特征。
4.在scRNA-seq水平上将糖酵解相关簇映射到CRC细胞我们通过对GSE132465数据集进行预处理，整合了23个肿瘤样本。在bulk-seq的基础上进行单个细胞GAC分类的构建。基于Lee等人的研究进行细胞类型注释（图4A)。通过使用“AddModuleScore”，将三个类似GAC的子集群映射到单个细胞（图4B)。此外，糖酵解相关基因标记（glygene）在图4C，表明其主要存在于骨髓、基质和上皮细胞中。随后我们将不同GAC簇添加到每种细胞类型中（图4E)。有趣的是，我们发现上皮细胞中GAC1簇（69%）和GAC3簇（31%）占比较多。虽然T细胞和B细胞与上皮细胞具有共同的GAC特征，但GAC3簇在这些细胞中的浓度高于GAC1样，表明其与上皮细胞不同。基质细胞几乎是GAC2簇（95%），髓样细胞中GAC2簇占65%，GAC3簇占33%。每种细胞亚型的标记基因展示于火山图中（图4D）。显然，髓系-G2中富含药物反应、DNA甲基化和谷氨酸代谢途径，而髓系-G3则相反。T细胞-G2对药物有良好的反应。B细胞-G1、2的mRNA甲基化程度较高。重要的是，我们筛选出了与预后（OS和DFS）最相关的GAC簇临床亚群，其中的标记基因用于后续GAC预测因子的分析。Meta-pool分析显示T细胞-G2、T细胞-G3、基质-G2、髓样细胞-G2和上皮细胞-G1具有重要意义，因此我们确定这五个GAC样亚群为独立的预后因素（图4G，H)。在上皮细胞中，类GAC亚簇和CMS亚型的分布如图5A、B。研究发现，CMS和GAC类亚群之间的对应关系与bulk-seq的相当，表明了这种关系的稳健性。通过GO富集分析，我们发现上皮细胞-G3富集在转移特征中，而上皮细胞-G1则主要与RNA生物过程相关（图5C)。同样，通过GSVA分析，我们在两种类型的上皮细胞中显示了50个标志性特征，证明了GAC3样的厌氧性状和GAC1样上调氧化磷酸化的能力（图5E)。两种富集分析的结果显示，GAC1类上皮细胞会增殖，而GAC3类上皮细胞会入侵。此外，对两种上皮细胞亚型临床特征的研究也显示出了不同之处。在CMS亚型中，CMS2在scRNA和bulk-RNA水平上都是GAC1样，而在CMS1和CM3中发现了相当一部分GAC3样，这与之前的结果一致。就分期而言，CRC中晚期患者的GAC1样上皮细胞比例增加。在突变方面，GAC1样上皮细胞中含有大量TP53和APC突变，其他部位的突变不明显。对于MSS，大多数GAC3样细胞具有MSI-H特征（图5D)。根据拟时序分析，上皮细胞将从GAC3样转变为GAC1样（图5F)。同时，HLA-A和HLA-B主要存在于GAC3样细胞中，表明这些细胞具有更强的免疫反应倾向（图5G，H)。最后，采用细胞通讯分析研究非肿瘤细胞对肿瘤上皮的影响。Sromal-G2和骨髓-G2与上皮细胞之间的连接频率较高，但Epthelial-G3受B细胞G2的影响更大（图5I)。气泡图特别显示了配体-受体网络，从中我们可以判断SPP1-CD44在骨髓-上皮相互作用中值得注意。CD44被认为可促进结直肠癌的干性和侵袭性。MDK配体与多种受体相关，以调节基质细胞和上皮细胞之间的相互作用。GZMA-F2RL1在T细胞-G3和上皮-G3（图5J)。最后，增殖型上皮-G1和侵袭型上皮-G3均受TME调控。
我们收集了两组TME细胞特征（Lee和Charoentong）进行ssGSEA分析。与scRNA水平一样，我们发现基质细胞在bulk数据中主要集中在GAC2。B细胞、骨髓细胞和T细胞聚集在GAC2和3中。上皮细胞的表达与scRNA的结果一致。值得注意的是，GAC2拥有更多的调节性T细胞（Tregs）（图6A）。在Charoentong的特征中，巨噬细胞、T细胞和B细胞的表达趋势与Lee观察到的结果相似（图6B）。Cibersort算法计算了每个细胞的表达比例。在Lee的细胞特征中，增殖性EC在TME中的比例最高，增殖性巨噬细胞次之，其中GAC2的比例分别最低和最高（图6C）。在LM细胞特征中，同样的现象也出现在巨噬细胞中，而T细胞在所有细胞中所占比例次之，在GAC2中最低。值得注意的是，GAC2中M2巨噬细胞的比例更高（图6D）。根据估算方法，GAC2和GAC3的基质和免疫评分高于GAC1（图6E）。为了预测不同GAC对免疫疗法的反应，我们收集了MSS、dMMR、TMB和PD-1信息。结果表明GAC3的MSI-High和dMMR百分比更高（图6F），TMB得分更高（图6G）。而GAC2的PD-1、PDL-1和CTLA-4表达相对较高（图6H-J）。Lee氏细胞的相关性见图6K，其中我们显示了糖酵解与每种细胞之间的关系。GAC1的免疫细胞浸润最少，GAC2有较多的基质细胞和免疫细胞，但也有很大比例的肿瘤促进细胞（Tregs、M2巨噬细胞），而GAC3有较高的免疫细胞浸润。总之，GAC1相当于免疫凋亡型，GAC2相当于免疫炎症/排斥型，而GAC3则相当于免疫激活型。
我们确定了基因簇，以研究GAC之间DEGs存在的内在原因。为了得到最终的DEGs，我们对每对GAC之间DEGs求取交集（图7A）。然后将108个基因进行共识聚类，得到了三个聚类，即基因簇（图7B）。从临床特征来看，我们发现基因簇B有显著的GAC2患病率和不良预后（图7C）。图7D显示了这108个基因的聚类表达模式和其他临床特征。从富集分析来看，基因簇A的MSS和MYC上调。在基因簇B中，分化、MSI、糖酵解和脂肪酸被高度富集。基因簇C中的TGF-Beta和EMT被激活（图7E）。根据GO分析，基因簇B和C与细胞外基质合成和上皮细胞增殖有关，而基因簇A则与T细胞活化和O-糖过程有关（图7F）。我们利用两个ssGSEA特征发现，基因簇C组在大多数TME细胞中有强表达。基因簇B的免疫细胞相对活化（图7G、H）。综上所述，GACs的临床、生化和免疫学特征都得到了基因聚类的支持。
8.GAC中TMB、CNV、甲基化和CSC指数的探索本节重点介绍肿瘤干性和基因组学的研究。我们使用瀑布图在COAD和READ中显示了每个GAC的前10个突变基因。在COAD中，GAC1和GAC2的突变相似，而在GAC3中，我们没有发现TP53的高突变率。在READ中，我们检测到FATA4的高突变率，特别是在GAC2中。APC、TP53、TTN和KRAS是突变频率最高的基因（图8A、B)。此外，我们注意到TMB评分和糖酵解评分有良好的联系，GAC3具有高水平的TMB和糖酵解属性（图8C)。同时，GAC3具有较低的拷贝数变异，表明与促肿瘤作用相反（图8D、E)。GAC1的总体甲基化低于其他两种（图8F)。通过CNV计算发现CSC指数与糖酵解评分呈正相关，GAC2占CSC指数最低（图8G，H)。综上所述，我们认为GAC1为低甲基化类型，GAC2为低CSC类型，GAC3为低CNV和高TMB类型。
9.基于机器学习的GAC预测因子构建和GAC预测因子相关基因（GP基因）的描述我们使用四种机器学习方法为每个患者的GAC状态构建了一个预测模型。我们将TCGA称为训练集，将GEO称为验证集。整个过程显示在图9A。为了确保结果的可靠性，我们在bulk和单细胞水平上使用具有临床意义的基因进行后续分析。因此，GAC DEGs与重要细胞亚型（图4G，H）获得。为了保留对预测结果影响最大的基因，我们使用Lasso，SVM，Randomforest（Boruta）和XGBoost进行筛选。雷达图表明，4种方法在测试集和训练集中的AUC值较高。我们再次交叉了通过四种机器学习方法过滤的基因，以获得最终的基因集，称为GP基因。随后，使用GP基因进行多项式逻辑回归，构建最终的GAC预测因子。将GAC预测因子应用于验证后，在训练集和验证集中，我们获得了高水平的AUC和准确度值（训练集：AUC = 0.9984，准确率 = 0.9658；验证集：AUC = 0.9380，准确率 = 0.8114）（图9B)。我们进一步对GP基因进行了几项分析。热图显示了GP基因的表达顺序（图9C)。对GP基因进行了基因组学相关研究，包括TMB和CNV分析。我们发现EIF2S2，DPM1和ISG20作为CNV扩增，GMDS，BNIP3L和VCAN作为CNV缺失（图9D)。所有GP基因均显示在染色体环图中（图9F)。通过TMB分析，我们将VCAN，PXDN，FN1和NOTCH3鉴定为COAD和READ（图9E)。
10.用于GAC预测因子验证的每个GAC和免疫治疗队列的治疗策略化疗和靶向药物治疗可提高晚期结直肠癌的存活率。然而，一些患者仍然几乎没有获得任何益处或产生治疗耐药性。因此，有必要为每位CRC患者制定个性化的药物治疗方案。基于三种GAC在临床和生物学上的显著差异，我们假定每种GAC都有自己合适的治疗方法，然后研究了典型化疗药物和可能的小分子药物在三种GAC中的敏感性和拮抗作用。我们使用"pRRophetic "软件包和Cmap分析对每种GAC进行了治疗预测。如图10A所示，以IC50的形式得出了药物敏感性。对于GAC1，喜树碱和紫杉醇更敏感。吉非替尼和吉西他滨对GAC2的敏感性较高。Shikonin对GAC3具有更大的治疗潜力。此外，我们还进行了Cmap分析，以获取抵抗每种GAC的潜在化合物。较强的化合物功效源于较低的分数。30种化合物在图10B提供治疗策略。为了测试GAC聚类在免疫治疗队列中的作用，我们获得了两个数据集进行验证（IMVigor210和GSE78220)。在IMVigor210中应用GAC预测因子后，我们发现GAC2预后不良（p = 0.004，对数秩检验），以及免疫治疗反应（图10C，D)。在GSE78220，得出的结论相同，即GAC2的OS率不利（p = 0.128，对数排名检验）。即使GSE78220的p值不显著，每个GAC的生存模式清晰可见（图10E，F)。总体而言，我们探索了GAC的潜在化合物，并验证了GAC在免疫治疗领域的预测因子。
