Transl Oncol | 机器学习衍生的自然杀伤细胞特征可预测透明细胞肾细胞癌的预后和治疗反应

‍

导语

今天给同学们分享一篇文章“Machine learning-derived natural killer cell signature predicts prognosis and therapeutic response in clear cell renal cell carcinoma.”，这篇文章发表在Transl Oncol期刊上，影响因子为4.5。

结果：

鉴定 NK 相关标记基因

预处理后，共包括 37,962 个细胞用于后续分析。使用“FindClusters”功能，将细胞分为 17 个簇（Fig. 1 A）。根据已发表的论文和 CellMarker 网站鉴定每个簇的细胞类型，并成功注释了 10 种细胞类型（ Fig. 1 B）。这些亚群高表达和低表达的标记基因如 C 和图 S1A 所示 Fig. 1 ，其中 NK 细胞显示 GNLY、NKG7 和 FGFBP2 等基因的高表达。

透明细胞肾细胞癌中 NK 细胞相关基因的鉴定

（A） 17 个细胞簇的 t-SNE。（B）根据已发表的论文和 CellMarker 注释的细胞类型的 t-SNE。（C）每种细胞类型的高表达和低表达的标记基因。（D） NK 细胞相关基因的功能富集分析。

为了进一步探索 156 个已鉴定的 NK 相关基因的生物学功能，进行了功能富集分析。GO 分析揭示了标记基因与白细胞介导的细胞毒性、细胞杀伤和白细胞介导的免疫等生物过程之间的显着关联（Fig. 1 D）。此外，KEGG 分析表明，这些 NK 细胞标记基因参与与自然杀伤细胞介导的细胞毒性和抗原加工和呈递途径相关的免疫反应（ Fig. 1 D）。富集分析的结果有力地支持作者对 NK 细胞群的准确注释。

蜂窝通信网络的描述

肿瘤微环境中不同类型细胞之间的相互作用在癌症的发生和发展中起着关键作用。细胞通讯可以通过直接的物理接触或通过分子（如激素、生长因子、趋化因子、充当配体的细胞因子）发生 [16 ]。利用“cellchat”R 包来描述细胞通讯网络可以进一步了解细胞间相互作用。

首先，研究了各种细胞类型之间通信的概率，表示为它们单独发出的信号数量（Fig. 2 A， S1B）。此外，还说明了每个细胞组的主要输出和输入信号。在 NK 细胞的情况下，输出信号主要集中在 CCL 、 TGF-β 、 IFN-II. 和 CX3C 等信号通路，而输入信号主要涉及 CCL 、 IL2 、 CX3C 和 BAG 信号通路（ Fig. 2 B）。这表明 CCL 通路在 NK 细胞与其他细胞群之间的通讯中起着至关重要的地位（ Fig. 2 C），其中 NK 细胞作为影响因素，其中 CCL4-CCR5 和 CCL5-CCR5 配体-受体对的贡献最大（ Fig. 2 D， S1C）。 Fig. 2 E 以热图的形式呈现 CCL 通路内细胞亚群之间的相互作用模式。很明显，NK 细胞主要向 T 细胞发出信号，主要接收来自 T 细胞和巨噬细胞的信号，这意味着 NK 细胞可能通过与 T 细胞和巨噬细胞相互作用来改变肿瘤细胞的生物学行为。

蜂窝通信网络的描述

（A）不同细胞类型之间的相互作用次数。（B）显示每种细胞类型的主要输出和输入信号的热图。（C） CCL 信号通路中的通讯相互作用。（D） CCL 信号通路中每个配体-受体对的贡献。（E） NK 细胞主要与 CCL 信号通路中的 T 细胞和巨噬细胞通讯。

基于机器学习的综合程序构建 NK 相关预后模型

对来自 TCGA KIRC 数据集的肿瘤和正常样本进行差异分析。火山图（Fig. 3 A）说明了 5189 个差异表达基因（DEG）（P < 0.05，logFC > 1）。随后与 NK 相关基因的交叉确定了 95 个 NK 相关 DEGs。进行单变量 cox 回归分析，筛选出 36 个可能与预后相关的 NK 相关 DEGs（图 S2A）。在 TCGA-KIRC 和 E-MTAB-1980 数据集中，作者拟合了 101 个预测模型，并进一步计算了训练集和验证集中每个模型的 C 指数。有趣的是，作者发现最佳模型是 Coxboost 和 plsRcox 的组合，它们的平均 C 指数最高（0.747）（ Fig. 3 B），并确定了最终的 18 个基因集（表 S4）。根据“survminer”包确定的最佳临界值，将 TCGA 患者分为高危组和低危组。Kaplan-Meier （K-M）曲线显示，与高危组相比，低风险组的总生存期（OS）和无进展生存期（PFS）明显更好（P < 0.0001）（ Fig. 3 C， D）。作者在外部 E-MTAB-1980 验证集（P < 0.0001）（ Fig. 3 E）中进一步证实了这一结果。这些数据表明，作者通过基于机器学习的综合程序构建的模型可以有效区分不同风险人群患者的预后，表明其作为临床预后分层工具的潜力。

构建 NK 相关预后模型

（A）火山图显示了 ccRCC 中正常样本和肿瘤样本之间的 DEGs。下降，下调;Ns，不显著;上升，上调。（B）共有 101 种预测模型，并进一步计算了每个模型在训练和验证数据集中的 C 指数。（C、D、E）基于风险评分的训练（TCGA）和验证（E-MTAB）数据集中生存分析的 Kaplan-Meier 曲线。

NK 相关预后模型的验证

多变量 cox 回归验证了风险评分可以作为独立的预后因素（图 S2B）。为了进一步证实作者建立的模型的有效性和敏感性，基于训练集绘制了时间依赖性的 ROC 曲线，1 、 3 和 5 年的 AUC 均高于 0.72，值分别为 0.724 、 0.741 和 0.753 （ Fig. 4 A）。同样，外部验证集中 1 、 3 和 5 年的 AUC 也达到了良好的水平（ 0.777 、 0.851 和 0.846 ）（ Fig. 4 B），表明该模型具有相当的有效性和稳健性。

NK 相关预后模型的验证

（A， B）训练数据集（A）和验证数据集（B）的 1 年、3 年、5 年 ROC 曲线。（C， D）根据 NK 细胞相关预后风险评分和其他临床病理特征得出的预后列线图，用于训练数据集（C）和验证数据集（D）的 KIRC 患者 1 年、3 年和 5 年生存率的预测。（E， F）训练数据集（E）和验证数据集（F）的校准曲线。（G）训练数据集中预测模型、风险评分、无风险评分的预测模型和临床病理因素的 C 指数可视化。（H）风险评分、等级、分期和组合的 DCA 曲线。

接下来，使用风险评分、年龄、性别、肿瘤分期和组织学分级一起作为临床特征构建列线图。列线图显示对应于不同变量的不同分数。通过结合每个变量的个体评分，总分可用于预测 1 、 3 和 5 年的生存概率（Fig. 4 C， D）。值得注意的是，校准曲线显示列线图预测的生存概率与实际观察值之间具有令人满意的一致性，证明了该模型出色的预测准确性（ Fig. 4 E， F）。此外，作者计算了模型的 C 指数，结果表明，结合风险评分和其他临床特征的模型的 C 指数为 0.774，优于任何其他预测模型（ Fig. 4 G）。作者继续使用 DCA 来评估预测模型在临床决策中的效用和价值。增加风险评分的预测模型明显显示出更好的益处（ Fig. 4 H）。最后，这些结果在验证队列中也得到了证实（图 S2C、S2D）。综上所述，作者的 NK 相关预后模型可以有效预测 ccRCC 的预后，有望为临床决策提供有益的信息。

NK 细胞的单细胞伪时间分析

从单细胞数据中提取 NK 细胞的亚群数据，并使用 R 包“monocle3”进行伪时间分析。图 S3A 说明了 NK 细胞群的伪时间值，其中较高的伪时间值代表细胞分化的后期。图 S3B 描述了在 NK 细胞群发育过程中表现出最显着变异的六个基因。 Fig. 5 A 和 B 展示了作者的基因特征中涉及的基因的分布和趋势，其中 CD7 、 PLAC8 和 METRNL 在 NK 细胞发育过程中表现出表达差异。METRNL 表现出表达增加，而 CD7 和 PLAC8 表现出下降趋势。结果提供了对 NK 细胞细胞状态和分化途径的顺序的见解。提取 NK 细胞群发育过程中表达差异显著的基因（P < 0.05）用于后续的 GO 和 KEGG 分析。富集分析结果表明，这些基因可能参与关键的生物学过程，如 RNA 转录、蛋白质翻译、抗原加工和呈递、NF-κB 信号通路以及白细胞趋化性（ Fig. 5 C）。

NK 细胞的单细胞伪时间分析

（A）显示参与基因特征的 18 个基因分布的特征图。（B） CD7、PLAC8 和 METRNL 在不同伪时间的变化趋势。（C）对发育过程中表现出显著差异表达的基因的 GO 和 KEGG 分析。图中标记的基因是变化最显着的 20 个基因。

功能富集分析为通路提供了新的见解

揭示了不同风险组患者临床病理特征的关联，作者发现较高的组织学分级和晚期肿瘤分期与较高的风险评分相关。验证集验证了上述发现（图 S3C）。为了进一步研究可能与作者的基因特征相关的途径，进行了相关性分析以过滤掉与风险评分呈正相关的基因（Pearson 相关性，R > 0.2 & p < 0.05）。然后将基因集用于功能富集分析，以进一步探索与风险评分相关的生物过程和途径。在 GO 分析中，在细胞-基质粘附、循环系统中的血管过程等方面观察到富集（图 S4A），这表明风险评分相关基因与细胞间通讯之间存在密切关联。在 KEGG 分析中，在神经活性配体-受体相互作用和细胞因子-细胞因子受体相互作用等途径中观察到富集（图 S4B）。此外，在高危组和低危组之间进行了差异分析，然后是 GSEA。结果表明，差异表达基因与 EMT、G2-M 检查点、IL6-JAK-STAT3 信号传导、干扰素γ反应等通路呈正相关，与脂肪酸代谢和胆汁酸代谢等代谢通路呈负相关（图 S4C）。 EMT 、 IFN-γ 反应和 JAK/STAT3 通路的显着富集表明特征与癌症增殖、迁移、侵袭和代谢之间存在潜在的紧密联系，这可能至少部分解释了在高危人群患者中观察到的较差临床结果。

高低风险组突变状态比较

据报道，基因突变的逐渐积累与肿瘤的发生发展密切相关 [20 ]。体细胞突变数据的分析有助于更深入地了解导致不同临床结局的因素。作者从 KIRC 患者的 TCGA 中获得了体细胞突变数据，并阐明了高低风险组之间体细胞突变的差异。高低风险组突变基因的比例表现出显著差异。低风险组前 5 位突变基因为 VHL （48 %）、 PBRM1 （ 45 %）、 TTN （17 %）、 SETD2 （9 %）和 MUC16 （ 8༅ ）（ Fig. 6 A）。在高危组中，前五个突变基因是 VHL （46 %）、PBRM1 （39 %）、BAP1 （17 %）、TTN （17 %）和 SETD2 （16 %）（ Fig. 6 B）。尽管高危组和低危组之间的肿瘤突变负荷（TMB）水平没有统计学上的显着差异（图 S5A），但 TMB 较高的患者的预后明显低于 TMB 较低的患者（P = 0.036）（ Fig. 6 C）。将风险评分与 TMB 相结合后，作者观察到高 TMB 和高危评分患者的生存时间远短于低 TMB 和低风险评分患者（P < 0.0001）（ Fig. 6 D），表明风险评分加 TMB 的组合可能提高预测 ccRCC 患者临床结局的准确性。

高危组和低危组之间突变的差异。

（一、二）训练集（A）和验证集（B）的体细胞突变瀑布图。（C）比较高 TMB 水平和低水平组的 KM 生存曲线。（D）比较四组的 KM 生存曲线。

高危组别和低危组别的免疫浸润情况

为了进一步探讨高危组和低危组之间免疫浸润的差异，作者使用“CIBERSORT” 和 “ssGSEA” 进行了额外的分析，以检查不同组的免疫细胞浸润模式。与低风险组相比，高危组观察到 CD8+ T 细胞、调节性 T 细胞、MDSCs 和 CD56 亮 NK 细胞的浸润显着增加。相比之下，高危组静息 NK 细胞、M1 巨噬细胞、肥大细胞和中性粒细胞的浸润率显着降低。（ Fig. 7 A，S5B）。此外，对几种常见免疫检查点基因表达差异的分析表明，高危组 PDCD1、CTLA4、LAG3 和 TIGIT 显着上调。然而，发现由 CD274 编码的 PD-L1 在高危组中的表达较低。这一 Fig. 7 观察结果表明，高危组可能具有增强的抑制抗肿瘤免疫反应的能力，最终导致生存时间缩短。采用“ESTIMATE”算法计算肿瘤浸润评分和免疫评分。结果表明，与低风险组相比，高危组表现出免疫评分升高和肿瘤浸润评分降低（ Fig. 7 C）。随后的相关分析表明，风险评分与免疫评分呈正相关（R = 0.25，P < 0.05），而与肿瘤浸润呈负相关（R = −0.14，P < 0.05）（图 S5C）。此外，利用 cibersort、ssGSEA、EPIC、estimate 和 MCP_counter 等多种算法全面分析两组免疫浸润差异。（ Fig. 7 D）.

高危组别和低危组别的免疫浸润

（A）高危组和低危组之间 28 个免疫细胞组分的比较。（B）高危组和低危组之间免疫检查点的不同表达。（C）箱形图说明了不同风险组之间免疫评分和肿瘤浸润的差异。（D）显示由多种算法计算的免疫浸润景观的热图。*p < 0.05、**p < 0.01 和 ***p < 0.001。

基因特征作为预测治疗反应的潜在工具

首先，利用 R 包 “oncoPredict” 计算每位患者 545 种药物的 IC50 值，并选择 ccRCC 临床治疗中的几种常用药物（舒尼替尼、阿西替尼、替沃扎尼、帕唑帕尼、索拉非尼、替西罗莫司、多柔比星、吉西他滨）进行研究。结果表明，高危组舒尼替尼、阿西替尼、替沃扎尼、帕唑帕尼、索拉非尼、替西罗莫司、多柔比星和吉西他滨的 IC50 值显著降低，表明酪氨酸激酶抑制剂（TKI）药物和 mTOR 抑制剂可能在高危组中表现出增强的治疗反应。值得注意的是Fig. 8 ，发现阿昔替尼的 IC50 值表现出类似的趋势，尽管没有达到统计学意义。

治疗敏感性分析

（A）高危组和低危组之间常用的 8 种药物的估计 IC50 的差异。（B）显示高危组和低危组免疫功能障碍评分和排除评分的小提琴图。（C）显示高危组和低危组 IPS 评分的小提琴图。（D） IMmotion150 队列中高危组和低危组之间反应率的差异。（E） IMmotion150 和 IMmotion151 免疫治疗队列中高危组和低危组之间的 PFS 差异。

为了逃避免疫系统的抗肿瘤作用，肿瘤细胞一方面可以通过免疫抑制因子阻碍免疫细胞浸润，另一方面，与 ccRCC 一样，尽管细胞毒性免疫细胞浸润水平升高，但这些免疫细胞往往处于功能抑制状态 [21 ]。采用肿瘤免疫功能障碍和排除（TIDE）算法综合评估每位患者的免疫排除和免疫功能障碍评分，预测肿瘤免疫逃逸能力和免疫检查点治疗的疗效。结果显示，与低风险组（ Fig. 8 B）相比，高危组免疫功能障碍评分升高，排除评分较低。IPS 评分越高，免疫原性越高，对免疫治疗的反应概率越高，高危组的 IPS 评分越过低风险组，意味着风险评分较高的患者对免疫治疗的反应效果可能增强（ Fig. 8 C）。作者还从 IMmotion150 队列中选择了 77 名接受 atezolizumab 治疗的患者，并分析了高危组和低危组之间的反应率差异。结果表明，高危组对治疗无反应的患者比例（SD/PD）显着更高（ Fig. 8 D）。最后，作者检查了 IMmotion150 和 IMmotion151 队列中高危和低危患者之间 PFS 的差异。高危组患者在两个队列中始终表现出较差的预后（ Fig. 8 E）。

ARHGAP9 的功能验证实验

考虑到 ARHGAP9 在模型中具有最大的系数，可以合理地推断其在模型框架中的关键作用。因此，作者进行了一系列实验来阐明其在肾细胞癌中的功能意义。作者的初步研究集中在 RCC 组织中 ARHGAP9 与邻近正常组织相比的差异表达。qPCR 和 IHC 分析证实，ARHGAP9 在 mRNA 和蛋白质水平（ Fig. 9 A、B、C）在 RCC 组织中显著过表达。作者使用两种人 RCC 细胞系 786-O 和 769P 进行了体外实验。在集落形成测定中，转染 siARHGAP9 的 RCC 细胞表现出克隆形成能力显著降低（ Fig. 9 D）。在 CCK8 测定中，ARHGAP9 敲低细胞对舒尼替尼表现出增强的耐药性，与作者之前的生物信息学预测一致（ Fig. 9 E）。最后，作者探讨了 ARHGAP9 对 RCC 细胞中 PD-L1 表达的影响。流式细胞术显示，ARHGAP9 敲低 ARHGAP9 导致 RCC 细胞中 PD-L1 表达上调（ Fig. 9 F、G）。总的来说，这些结果强调了 ARHGAP9 在调节 RCC 增殖、耐药性和免疫微环境调节中的潜在作用。这些发现表明，ARHGAP9 是 RCC 病理生理学不可或缺的一部分，可能是一个关键的治疗靶点。

患者样本和肾癌细胞的实验验证

（A）肿瘤和正常组织之间 mRNA 表达水平的差异。（二、三）IHC 揭示了肿瘤组织和正常组织之间 ARHGAP9 的差异蛋白表达水平。图像右上角的比例尺代表 50 μm。（D）敲低 ARHGAP9 显著降低两种肾癌细胞系的增殖速率。（E）敲低 ARHGAP9 显著增加对舒尼替尼的耐药性。（女，克）敲低 ARHGAP9 诱导 PD-L1 表达增加。MFI，中位荧光强度。数据表示为平均 ± SD，重复三次。*，p < 0.05;**， p < 0.01;， p < 0.001;， p < 0.0001.

总结

通过利用单细胞数据和批量测序数据，作者探索了 NK 细胞与各种细胞亚群之间错综复杂的信号相互作用，研究了 NK 细胞群在不同阶段的基因变化。随后，通过基于机器学习的集成程序建立了与 NK 细胞相关的 18 个基因特征，这不仅在预后预测方面表现出出色的疗效和稳健性，而且还揭示了风险评分与肿瘤微环境的复杂性之间的密切相关性，以及对免疫疗法的反应。最后，作者进行了 qPCR、IHC 和一系列细胞功能测定，以研究 ARHGAP9 在作者的患者样本和肾癌细胞系中的表达模式和功能。作者的研究不仅为 NK 细胞在 ccRCC 微环境中的关键作用提供了新的见解，而且还建立了具有潜在临床应用的预测模型。对这篇文章感兴趣的老师，欢迎扫码咨询！

往期推荐

纯生信选刊

非肿瘤生信

预后模型

单基因生信

单细胞系列