这是一篇在2024年8月16日发表在Journal of Translational Medicine[IF:6.1/Q2]杂志上的文章《Machine learning‑based identification of biomarkers and drugs in immunologically cold and hot pancreatic adenocarcinomas》
研究背景：胰腺癌(PAAD)常表现为“冷”或免疫抑制的肿瘤环境，这与免疫检查点阻断治疗的抵抗有关；然而，其潜在的机制尚不完全清楚。在这里，研究团队旨在提高研究团队对发生在肿瘤微环境中的分子机制的理解，并确定生物标记物、治疗靶点和潜在的药物来改善PAAD的治疗。根据免疫热型和冷热型PAAD的不同亚型将患者分为不同的疾病结局。通过COX回归和加权相关网络分析，构建了一个新的基因特征，称为“热点肿瘤下调、预后和免疫相关基因”(DPIRGs)，用于通过机器学习(ML)建立PAAD的预后模型。综合分析DPIRGs在PAAD中的作用，并用ML鉴定能够区分PAAD免疫亚型和预测预后的生物标志物基因。利用公共单细胞转录和蛋白质组资源验证了生物标记物的表达。通过分子对接研究确定了使冷肿瘤变热的候选药物和相应的靶蛋白。以DPIRG信号为输入数据，从137个ML组合中选择生存随机森林和偏最小二乘回归的组合COX，构建了PAAD的最优预后模型。从遗传/表观遗传改变、免疫渗透、途径丰富和miRNA调控等方面探讨了DPIRGs的作用及其分子机制。确定了生物标志物和潜在的治疗靶点，包括PLEC、TRPV1和ITGB4等，并验证了生物标志物的细胞类型特异性表达。候选药物包括沙利度胺、SB-431542和博莱霉素A2，基于它们对DPIRG表达的有利调节能力。
Figure1 冷热肿瘤分型，并揭示了其显著的分子特征：(A) 聚类分析表明热肿瘤与冷肿瘤的基因表达差异明显；(B) 生存分析显示热肿瘤患者预后优于冷肿瘤患者；(C) 热肿瘤在免疫评分和基质评分上显著更高；(D)研究通过差异表达分析发现热肿瘤相较于冷肿瘤共有2055个基因上调、2565个基因下调 . （图1E、1F）GO分析表明，上调基因主要富集于免疫相关通路（如免疫反应激活、白细胞迁移等），而下调基因则与表皮发育和细胞骨架相关；（图1G、1H）KEGG分析显示，上调基因富集于免疫配体-受体相互作用和细胞因子信号通路，而下调基因与代谢通路相关；(I) GSVA分析进一步揭示热肿瘤在脂肪生成、移植物排斥和血管生成等通路中显著富集，反映了热肿瘤的免疫活性增强和代谢特性变化。
Figure2 预后和免疫相关基因特征的鉴定 使用 WGCNA 程序识别了用不同颜色表示的十个模块。 作者团队计算了每个模块与免疫簇或免疫细胞之间的相关性，发现粉色和青绿色模块与冷肿瘤相关，而黑色模块与热肿瘤相关，两者相关系数值均> 0.3（图2A）。 在粉色模块中，基因显着性（GS）和模块成员资格（MM）之间的相关系数为0.46（图2B），而青绿色和黑色模块的相关系数分别为0.81和0.23（图2C和D） ）。 选择GS > 0.3和MM > 0.5的基因进行进一步分析（图2B-D），其中包括165个热相关基因（HRG）和4183个冷相关基因（CRG）（图2E，F）。 DEG和WGCNA结果的交叉显示，上调DEG和HRG之间有118个重叠基因，下调DEG和CRG之间有375个重叠基因； 提取这些基因用于后续分析（图2E，F）。此外，相关性分析表明，粉色和青绿色模块中的基因可能负向调节CD8+ T细胞、M2型巨噬细胞和静息肥大细胞（R < − 0.2），正向调节静息NK细胞、M0型巨噬细胞、活化的 肥大细胞和嗜酸性粒细胞 (R > 0.2)。 黑色模块中的基因预计会上调幼稚 B 细胞、CD8+ T 细胞和记忆静息 CD4+ T 细胞 (R > 0.2)（图 2G），但会负向调节 M0 型巨噬细胞和激活的肥大细胞 (R < − 0.2)（ 图2G）。
Figure3 构建共识的机器学习驱动的预后模型为了开发共识模型，使用留一交叉验证 (LOOCV) 框架对 82 个 UPIRG、96 个 DPIRG 和所有 178 个预后免疫相关基因进行基于 ML 的整合程序。 在混合 TCGA + ICGC 训练队列中，作者团队基于单一和组合 ML 算法构建模型，并计算每个模型的 C 指数值。 此外，作者团队使用 TCGA-PAAD、ICGC-CA 和 ICGC-AU 数据集计算了训练模型的 C 指数值，然后计算了四个 GSE 数据集的平均 C 指数，以评估所有模型的预测能力 （图3A-C）。 结果表明，UPIRGs的最佳模型是Survival RF和Enet的组合（a = 0.1）（图3A）； 对于 DPIRG，它是 Survival RF 和 PlsRcox 的组合（图 3B）； 对于 PIRG，它是生存 RF 和 ANN 的组合（隐藏 = 7）（图 3C）。 基于这三个模型，结合 Survival RF 和 PlsRcox for DPIRGs 的模型被发现是具有最佳 C 指数的最佳预后模型，并被指定为混合模型。 为了评估 DPIRG 的预后意义，作者团队根据风险评分的中值将 PAAD 患者分为高风险组和低风险组。 根据TCGA + ICGC训练数据库，高危组和低危组的生存率存在显着差异（p < 0.001），其中低危组的生存概率更高（图3D）。 一致的是，TCGA + ICGC 队列中 1 年、2 年和 3 年生存率的时间依赖性 ROC 曲线的曲线下面积 (AUC) 值分别为 0.979、0.983 和 0.986（图 3E）。 作者团队还比较了两组之间的风险评分和临床状态，并证明高风险组的死亡率高于低风险组（补充图4A、B）。 高风险组的存活率较低（图3F）。
Figure4 使用四个 GSE 数据集来验证使用 DPIRG 预测 PAAD 患者预后的可行性。 作者团队还组合了四个验证数据集以形成组合测试。 K-M生存曲线显示，低危组患者的总生存期（OS）优于高危组（图4A-E）。 ROC曲线的AUC值为GSE28735为0.735，GSE62452为0.710，GSE78229为0.789，GSE85916为0.647，组合测试为0.760（图4F-J）。 作者团队还比较了四个 GSE 数据集中两组之间的风险评分和临床状态（补充图 4C-J），作者团队的模型显示高风险组中的患者死亡人数多于相应低风险组中的患者（图 4C-J）。 4K-O）。 总体而言，一个训练队列和四个验证队列的 KM 生存分析、时间依赖性 ROC 曲线分析和 C 指数值计算一致表明，DPIRG 特征可以潜在地预测外部验证队列中 PAAD 患者的结果。
Figure5PAAD 中 DPIRG 特征的预后价值在基于 ML 构建 PAAD 患者的预后模型后，作者团队还研究了是否可以将不同的预后模型分别应用于热肿瘤和冷肿瘤患者。 因此，使用单一和组合的ML算法来构建基于DPIRG的热肿瘤或冷肿瘤患者的共识模型。 根据所有模型的四个 GSE 验证数据集的平均 C 指数值，热肿瘤患者的最佳模型是生存 RF 和脊的组合，作者团队将其指定为热模型（图 5A），并且 对于冷肿瘤患者来说，最佳模型是 plsRcox 和 XGBoost 的组合，作者团队将其命名为冷模型（图 5B）。 Cox分析表明，热模型比冷模型更能将患者分为高危组和低危组； 根据热门模型以及 TCGA + ICGC 和四个 GSE 数据库的数据，风险组之间的生存差异更为显着（图 5C-G）。此外，作者团队将混合模型与热模型进行了比较。 作者团队发现，只有 GSE78229 数据集根据热模型（图 5F，左）比根据混合模型（图 4C）的生存率差异更大。 相比之下，使用混合模型的其他数据库产生了更显着的生存差异（图 4A、B、D）。 接下来，作者团队重点使用训练数据库和四个 GSE 数据库的数据来验证混合模型、热模型和冷模型预测 PAAD 预后的能力。 作者团队发现混合模型（图 4F-I）和热模型（图 5I-L，左）比冷模型（图 5I-L，右）表现更好，AUC 值更高。 总之，混合模型和热模型在不同数据集上都有优势，而冷模型在预后预测方面表现较差。
Figure6DPIRGs 的遗传改变和 DNA 甲基化分析 作者团队的研究结果表明，DPIRG 特征与 PAAD 患者的预后显着相关。 为了探索调节 DPIRG 表达的机制，作者团队使用 cBioPortal (http://www.cbiop ortal.org/) 来阐明 PAAD 中 DPIRG 的遗传改变。 一些基因突变被鉴定出来； 然而，突变频率相对较低（例如，ASPM，4％；AHNAK2，3％；DCST1，4％），其中PLEC的突变频率最高，为9％（图6A）。 作者团队还评估了基因拷贝数变异（CNV）和表达水平之间的关系。 检测到CNV与ASPM、TRPV1、PLEC、POLQ、DCST1、ITGB4、AHNAK2和GPR52表达水平之间呈正相关（图6B）。 此外，对具有这些 DPIRGs 的 CNV 的患者进行了生存分析，结果表明 TRPV1、SDHAP1、SCARNA9、SCARNA7、IGF2BP2-AS1、FAM86HP 和 FAM157A 的 CNV 与疾病特异性生存 (DSS)、进展显着相关。 -无生存期(PFS)和OS(图6C)。 接下来，探讨了DNA甲基化对这些DPIRG表达的影响，结果表明ITGB4、AHNAK2和XDH的DNA甲基化水平与其表达呈负相关（图6D）。 作者团队还发现，在 PAAD 患者中，TRPV1 DNA 甲基化与无病间隔 (DFI)、DSS、PFS 和 OS 相关（图 6E）； SNORA12 DNA 甲基化与 PFS、PFS 和 OS 相关； ITGB4 DNA 甲基化与 DFI、PFS 和 OS 相关（图 6E）。
Figure7预测 DPIRG 对于 PAAD 预后和通过 ML 进行聚类很有价值 为了了解 DPIRGs 影响 PAAD 患者肿瘤免疫景观和预后的潜在生物学功能和机制，作者团队首先使用五种 ML 算法（即 SVM、ANN、Boruta、RF 和 XGBOOST）在混合、 热肿瘤和冷肿瘤。 在混合肿瘤中，针对每种 ML 算法对 DPIRG 值进行归一化与患者生存最显着相关的前 5 个基因被确定为 AL591135.1、AL158201.1、AHNAK2、AK3P5 和 CEP295NL （图7A）。 热肿瘤中排名前5位的基因是AL627402.1、RGPD4、ASPM、SMARCE1P1和GJA1P1（图7B），冷肿瘤中排名前5位的基因是GLIPR1L1、AL158201.1、WASHC5AS1、AK3P5、 和 AL591135.1（图7C）
Figure8DPIRGs 在塑造 PAAD 亚型免疫景观中的作用 根据免疫亚型和风险水平，PAAD 患者被分为四组：热高组、热低组、冷高组和冷低组。 四组中热低组的预后最好，而热高组和冷高组的预后最差（图8A）。 根据混合数据库，分配到热低组的患者多于分配到冷低组的患者（33% vs. 17%），这表明低风险患者往往患有热肿瘤。 相比之下，高危患者均匀地分布在热组和冷组中（分别为 23% 和 27%）（图 8B）。 此外，作者团队进行ROC分析来预测热低肿瘤和冷高肿瘤患者的生存率，发现能够有效区分热低肿瘤和冷高肿瘤的前10个基因是AHNAK2、ITGB4、ACTBP7、PLEC、ANXA2P1 、FTH1P4、FTH1P12、KRT8P33、KRT18P7 和 AC087257.2，AUC 范围为 0.816 至 0.791（图 8C，补充图 9）为了检查各个 DPIRG 之间的潜在相互作用，作者团队估计了各个 DPIRG 之间的相关性，发现热低组中 DPIRG 表达水平的总体相关性大于冷高组 组（图8D，补充表3、4）。 接下来，为了检查个体 DPIRG 与不同 PAAD 风险组的关联，作者团队分析了个体 DPIRG 表达水平与使用 ssGSVA 算法计算的 DPIRG 评分之间的相关性（补充表 5-6）。 总体而言，冷高组中的 DPIRG 表达与 DPIRG 评分之间的相关性大于热低组中的相关性。 同样，在高冷组中，个体 DPIRG 表达水平与风险评分之间的相关性更大（图 8D，补充表 7、8），表明 DPIRG 可能更显着地影响高冷组中的 PAAD 风险水平。
Figure9PAAD亚型的miRNA和GSVA分析 为了探究 miRNA 是否参与 DPIRGs 的功能，作者团队搜索了其水平与 DPIRGs 表达显着相关的 miRNA（相关系数 > 0.7），并鉴定了热低组中的 16 个 miRNA 和冷高组中的 20 个 miRNA。 组（图9A）。 然后作者团队评估了这 36 个 miRNA 表达的差异，发现只有 hsa-mir-139 在 hot-low 组中表达量较高，而 hsa-mir-193a、hsa-mir-1248、hsamir 的表达水平较高。 -365a、hsa-mir-365b 和 hsa-mir-93 在冷高组中更大。 为了探讨 DPIRGs 在低热 PAAD 肿瘤和高冷 PAAD 肿瘤中是否发挥不同的生物学功能，作者团队进行了 GSVA； 相对于冷高组，鉴定为在热低组中富集的标志通路如图9C所示。 在两组中，幼稚 B 细胞与前 10 个富集标志通路呈正相关，而静息 NK 细胞、活化的肥大细胞和嗜酸性粒细胞与这些通路呈负相关（图 9E）。 一般来说，CD8+ T 细胞、M1 和 M2 巨噬细胞以及静息肥大细胞与冷高组的通路更显着相关，而 Tregs、滤泡辅助 T 细胞和活化的树突状细胞与热高组的通路更显着相关。 -低组（图9E）。 这些生物途径和各种免疫细胞之间的相互作用可能参与确定 PAAD 亚型的结果。
Figure10免疫治疗反应分析和药物预测 TIDE 评分与肿瘤免疫逃逸和免疫治疗耐药的可能性密切相关，TIDE 评分越高，预测免疫治疗反应率越低。 如图10A所示，热低肿瘤中的TIDE评分显着低于冷高肿瘤中的TIDE评分。 此外，对 PD-1 阻断有反应的患者的 DPIRG 评分显着低于无反应者（图 10B）。 作者团队还使用五种 ML 算法预测了影响 TIDE 分数的显着 DPIRG，并发现排名前 5 的基因是 RPL1P42、KRT18P7、MIR554、RBMXP2 和 FAM157A（图 10C）。 接下来，作者团队使用 OncoPredict 包计算 DPIRG 表达与现有药物反应之间的相关性。 23种药物与DPIRG表达水平呈负相关（补充图11，补充表15、16），这23种药物的IC50值在热低组和冷高组之间存在差异，热低组中的IC50值较高。 组，表明这些药物在高冷组中可能具有更有效的治疗效果（图10D）。 在热低组中，23 种药物对浆细胞、M2 型巨噬细胞和中性粒细胞产生正向调节，同时对滤泡辅助 T 细胞和嗜酸性粒细胞产生负向调节。 在高冷组中，药物与 Tregs 和静息 NK 细胞呈强烈负相关（图 10E）。 此外，这些药物与标志通路呈正相关，热低组的相关系数高于冷高组（图10F）。 接下来，使用五种ML算法来预测药物对预后的重要性，与患者生存相关的前8名药物是CAY10594、沙利度胺、AT13387、SB-431542、达沙替尼、博莱霉素A2、ML258和TGX-221（图1）。 10G）。 作者团队还使用了五种ML算法来预测可以显着调节DPIRG评分的药物，排名前8的药物是BRD-K02251932、WAY-362450、西咪替丁、semagacestat、SB-431542、博莱霉素A2、ML258和沙利度胺（图10H）。
Figure11药物分子与 DPIRG 的结合作者团队随后从ZINC15数据库下载了8个活性化合物的化学结构，即AT13387、卡铂、CAY10594、沙利度胺、TGX-221、达沙替尼、赛马西特和SB-431542。然后，作者团队选择了四个具有相应完整蛋白质结构的基因来探索基因与药物之间的结合模式(图11A-D)。根据对接分数，与GLIPR1L1结合最强的两种药物是SB-43154和赛马西坦(补充表17)。在GLIPR1L1蛋白分子表面发现了一个口袋，这使得SB-43154能够与其相互作用形成相对稳定的复合体(图11A)，而GLIPR1L1表面的另一个口袋与半乳糖相互作用形成复合体(图11B)。根据二维和三维分子相互作用可视化显示，Semagacestat与GLIPR1L1Asp84和Gln225氨基酸的相互作用比SB-43154更强。与TRPV1、PLEC和CEP295NL结合最强的两种药物也是SB-43154和赛马西特(补充表17)。虽然在TRPV1和Plec蛋白质的表面没有发现口袋，但由于药物和蛋白质之间有大量的氢键，这些蛋白质的结合保持相对稳定(图11C-F)。SB-43154和Semagacestat与CEP295NL蛋白弱结合，只与四个氨基酸相互作用(图11G，H，补充表17)。
Figure12PAAD 样品中细胞类型特异性 DPIRG 表达的验证 作者团队下载了 8 个 PAAD scRNA-seq 数据集，鉴定了 TME 中的 15 种细胞类型，并检查了几种 DPIRG 的细胞类型特异性表达，这些 DPIRG 在低热肿瘤和高冷肿瘤之间差异表达（图 8C），特别是那些表达为 与遗传/表观遗传调控相关（图6），包括AHNAK2、PLEC、ITGB4和XDH。 作者团队发现这些基因在恶性细胞中表达更强烈（图12A）。 在 HPA 网站上，提供了与 9 个 DPIRG 对应的蛋白质的 IHC 结果，包括 AHNAK2、ANKRD61、CEP295NL、DCST1、ITGB4、PLEC、RGPD3 和 ZC3H11B。 根据 IHC 分析，这些蛋白质在代表性 PAAD 样品中的分布和染色强度（从低表达到高表达）如图 12B 所示。
本文亮点：
- 免疫热型与冷型胰腺癌的分子特征揭示：通过对免疫热型（热肿瘤）和冷型（冷肿瘤）胰腺癌（PAAD）样本的基因表达差异进行深入分析，明确了两种肿瘤类型在免疫活性、代谢特性等方面的显著不同，并进一步确认了热肿瘤具有更高的免疫活性。
- DPIRGs基因特征的预后价值：通过加权相关网络分析（WGCNA）和机器学习（ML）算法，构建了新的基因特征DPIRGs（热点肿瘤下调、预后和免疫相关基因），并成功构建了基于DPIRGs的PAAD患者预后模型，展示了其在区分高风险与低风险患者中的潜力。
- 药物预测与免疫治疗策略：文章结合机器学习和分子对接，预测了能够调节DPIRGs表达的潜在药物，包括沙利度胺、SB-431542和博莱霉素A2等，这些药物可能对提高冷型肿瘤的免疫反应性起到重要作用。
- 细胞类型特异性验证：通过多种单细胞RNA测序数据集，验证了DPIRGs在PAAD不同细胞类型中的特异性表达，揭示了这些基因在肿瘤免疫微环境中的作用，特别是在与遗传和表观遗传调控相关的基因表达。
- 分子机制探讨与遗传分析：通过cBioPortal对DPIRGs的遗传改变进行分析，发现多个DPIRGs基因（如PLEC、TRPV1和ITGB4等）存在显著的基因突变和拷贝数变异，进一步探讨了它们在胰腺癌免疫逃逸和预后中的潜在作用。
这些亮点为PAAD的免疫机制研究和个性化治疗提供了新的思路，尤其是在优化免疫治疗和精准药物筛选方面具有重要应用价值。
当然了，本文也是提供了原始代码的，各位看官老爷自取（https:// github. com/ sangm m12/ ML_ Hot- cold.）
library(data.table)
library(GSVA)
library(ggplot2)
library(ggcor)
dir.create(paste0("./12.蝴蝶图"))
afdir <- paste0(getwd(),"/12.蝴蝶图")
ABC <- read.table("merged_data.txt", sep = "\t", header = TRUE, check.names = TRUE, row.names = 1)
immPath.score <- read.table("coldtcga.uniSigExp.txt", sep = "\t", header = TRUE, check.names = TRUE, row.names = 1)
# 加载 dplyr 包
library(dplyr)
# 根据列名合并两个数据框
immCorABC2 <- bind_rows(immPath.score, ABC)
# 去除含有NA的列
immCorABC2 <- immCorABC2[, colSums(is.na(immCorABC2)) == 0]
#immCorABC <- read.table("immCorABC.txt", sep = "\t", header = TRUE, check.names = TRUE, row.names = 1)
#immCorABC <- read.table("immCorABC.txt", sep = "\t", header = TRUE, check.names = TRUE, row.names = 1)
genei="riskscore"
# 循环计算相关性并绘制左下角
immCorABC <- NULL
for (i in rownames(immPath.score)) {
cr <- cor.test(as.numeric(immPath.score[i,]),
as.numeric(ABC),
method = "pearson")
immCorABC <- rbind.data.frame(immCorABC,
data.frame(gene = genei, ####？
path = i,
r = cr$estimate,
p = cr$p.value,
stringsAsFactors = F),
stringsAsFactors = F)
}
immCorABC$sign <- ifelse(immCorABC$r > 0,"Positive Relation","Negative Relation")
immCorABC$absR <- abs(immCorABC$r)
immCorABC$rSeg <- as.character(cut(immCorABC$absR,c(0,0.25,0.5,0.75,1),labels = c("0.25","0.50","0.75","1.00"),include.lowest = T))
immCorABC$pSeg <- as.character(cut(immCorABC$p,c(0,0.001,0.01,0.05,1),labels = c("<0.001","<0.01","<0.05","ns"),include.lowest = T))
immCorABC[nrow(immCorABC),"pSeg"] <- "Not Applicable"
immCorABC$rSeg <- factor(immCorABC$rSeg, levels = c("0.25","0.50","0.75","1.00"))
immCorABC$pSeg <- factor(immCorABC$pSeg, levels = c("<0.001","<0.01","<0.05","Not Applicable","ns"))
immCorABC$sign <- factor(immCorABC$sign, levels = c("Positive Relation","Negative Relation"))
#pdf(paste0(afdir,"/2",paste0(genei,"1.pdf"),width=10,height = 8)
p1 <- quickcor(t(immPath.score),
type = "upper",#绘制左下角的就是lower
show.diag = TRUE) +
geom_colour() +
anno_link(data = immCorABC,
mapping = aes(colour = pSeg, size = rSeg, linetype = sign),
spec.key = "gene",
env.key = "path",
diag.label = FALSE) +
scale_size_manual(values = c(0.5, 1, 1.5, 2)) +
scale_color_manual(values = c("#19A078","#DA6003","#7570B4","#E8288E","#65A818")) +
scale_fill_gradient2(low = "#9483E1",mid = "white",high = "#E11953",midpoint=0) +
remove_axis("x")
p1
#pdf(file =paste0(genei,"左下.pdf"),width=12,height = 9)
pdf(file =paste0(afdir,"/",genei,"右上.pdf"),width=12,height = 12)
print(p1)
#ggsave(filename = "左下.pdf", width = 10,height = 8)
dev.off()
写在最后：101机器模型我们也写过一次，基于此发出的文章也有很多，但是这篇的思路以及分析方法，真的让人眼前一亮，小编团队也实现了大部分图的复现，各位有条件的赶紧复现起来吧。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟