记录下最近阅读的一篇于2024年6月在Cell Reports Medicine
上发表的文章"scRank infers drug-responsive cell types from untreated scRNA-seq data using a target-perturbed gene regulatory network"，记录下翻译过程，以备后期随时翻阅。
引言
细胞对药物的反应表现出巨大的异质性，这主要归因于细胞群体的多样性。在不同细胞的特定生物网络中，药物与其靶点（如受体或酶）相互作用会产生不同的效果。例如，激酶抑制剂可能会阻止肿瘤细胞增殖，而对正常细胞影响最小，因为不同的途径激活有所不同。因此，理解异质的药物反应需要对复杂组织中的细胞亚群进行表征。传统上，组织转录组分析仅提供聚合的表达数据，掩盖了不同细胞群体的独特分子特征。这些独特的分子特征在介导细胞对治疗的不同反应中至关重要。最近的单细胞 RNA 测序（scRNA-seq）技术进展提供了在单细胞水平上解析细胞异质性的机会。这一方法揭示了大量信息，包括细胞类型的识别、其空间分布、细胞间通信模式以及理解疾病进展和细胞对药物反应的基因调控机制。因此，药物开发策略已从针对特定因素转向针对关键细胞类型，从而为精准医学开辟了新途径。最近的研究，包括使用治疗集群识别可药物作用的细胞群体和深度学习方法在疾病背景下优先考虑细胞类型，进一步加深了我们对药物开发中细胞靶标的理解。此外，将整体表达数据与 scRNA-seq 分析相结合的工具为确定表型相关的细胞群体提供了补充见解。尽管 scRNA-seq 取得了进展，但识别对药物治疗有反应的关键细胞类型仍是一个重要挑战，因为这需要理解药效动力学，包括药物与其靶点的直接相互作用及其对细胞途径的后续级联效应。这些知识对于理解细胞水平的药物作用机制和开发副作用更少的精确治疗策略至关重要。虽然之前大多数分解药物反应的努力都集中在识别基因表达或蛋白质活性显著变化上，但哪个细胞类型对药物治疗表现出最高反应性的问题尚未得到充分解决。
结果
scRank 概述
我们提出了 scRank（https://github.com/ZJUFanLab/scRank），这是一种基于细胞类型对指定药物潜在基因网络活动进行排名的计算工具（图 1）。药物扰动不仅对单一靶点产生影响，还会影响下游基因。这些与靶点相关的下游基因往往在共同的路径中共同表达，通常协同作用以响应药物治疗。药物通过与靶点的相互作用，不仅通过直接靶点调节发挥作用，还通过这些共同表达的下游基因启动一系列信号事件，最终导致特定的生物功能。因此，对药物影响的综合评估必须超越仅仅检查靶基因表达，涵盖药物在更广泛的生物网络中的整体影响。我们使用 scRank 的方法利用网络方法来预测药物在不同细胞类型中的影响程度。我们假设化合物扰动通过抑制靶点，应该沿分子网络传播，最终影响网络的整体结构和强度。通过扰动基于共表达关系重构的基因表达谱中的分子网络中的药物直接靶基因节点，我们可以捕捉药物的系统性影响。我们将这种扰动效应建模为计算机模拟药物扰动所产生的全局网络差异和靶相关信号的局部网络扩散的组合。通过使用 tpGRN 实施这一假设，我们使用扰动评分来评估计算机模拟药物扰动对每个细胞类型的影响并对它们进行排序。由于这种计算机模拟药物扰动，我们的方法允许仅基于未处理的数据推断药物反应中的细胞类型优先级。
图 1：scRank 方法工作流程 > (A) 未处理的单细胞 RNA 测序数据和直接靶基因。
(B) 基于 tpGRN 对药物响应的细胞类型进行排名的概念。首先构建细胞类型特异性的 GRNs。在确定靶基因后，通过切断靶基因节点的所有边缘创建 tpGRN。通过比较未处理的 GRNs 和 tpGRNs，量化药物扰动以对细胞类型进行排名。
(C) 结合局部和全局扰动效应，以优先考虑每种细胞类型的药物反应
(D) 构建 GRN 的过程。抑制剂药物的靶基因从 DGIdb 收集，前 2,000 个高变异基因和转录因子被考虑用于构建 GRN。使用主成分回归获得共表达邻接矩阵。对于每种细胞类型，使用相同的基因节点创建特定的 GRN。
(E) 通过流形对齐计算全局扰动效应的示意图
(F) 使用网络扩散计算局部扩散扰动效应的示意图。Dd 表示药物靶点节点的距离，Dn 表示 1 跳基因节点 n 的距离，Wout 表示输出边的权重，Win 表示输入边的权重，Deg 表示基因节点的度。
使用合成真实数据集验证 scRank
我们首先在三种情况下对合成数据上的 scRank 进行测试，以评估其在识别药物反应细胞类型方面的性能。在我们的研究中，合成数据集设计用于模拟疾病状态下的细胞转录组。为实现这一点，我们使用 SERGIO 从预定义的 GRN 模拟单细胞基因表达数据集，其中疾病相关基因相互作用并形成疾病模块。我们假设如果药物的靶基因位于高活性疾病模块内，那么药物的扰动效应在细胞网络内将更为显著。因此，具有活性靶基因的细胞类型应对这种扰动表现出更显著的反应，从而产生细胞类型排名的真实情况（图 2A）。具体来说，我们开发了三个不同的场景，每个场景都有定义好的药物反应细胞类型的真实排名（图 2B）。这种方法使我们能够在不同情况下评估模型的预测准确性。
图 2：scRank 在模拟数据集中的基准测试 > (A) 生成模拟真实数据的示意图。我们使用预定义的 GRNs 通过 SERGIO 模拟三种情况下不同细胞类型中具有不同模块活动的疾病 scRNA 数据。我们假设具有高药物相关模块活性的细胞类型可能会受到较大的扰动。药物反应细胞类型排名的真实情况可以通过预定义的模块活动获得。
(B) 每种情况下的细胞类型真实排名。
(C) 四种细胞类型的 3,000 个模拟细胞，具有递增的模块 1 活性梯度。scRank 在不同细胞类型中的扰动评分按比例缩放后以点线图呈现。数据以中位数 ±SD 表示。每种细胞类型的数据点数量为 25。scRank 预测的具有递增细胞数量的细胞类型的平均排名以热图表示，其中星号表示每个数据集的排名最高的细胞类型。
(D) 六种细胞类型的 3,000 个模拟细胞，可根据它们的二进制模块 1 活性分为高低两组。scRank 计算的两组的扰动评分按比例缩放后以箱线图表示（最小值，第 25 百分位数，中位数，第 75 百分位数和最大值）。每组的数据点数量为 75。p 值通过 Wilcoxon 秩和检验计算，随细胞数量增加
(E) 四种细胞类型的 3,000 个模拟细胞，具有细胞类型特异性模块高活性。条形图为 scRank 计算的不同模块基因的扰动评分。数据以中位数 ±SD 表示。每个模块的数据点数量为 25。不同类型扰动下排名最高的细胞类型随细胞数量增加而变化的点线图。
scRank 与其他方法的性能比较
为了进一步评估 scRank 在识别药物靶向细胞类型方面的性能，我们将其与现有方法进行了比较。这些方法基于不同条件下（治疗前和治疗后）相同细胞类型之间的相对 DEGs 数量或可分离性来优先排序对药物反应最强的细胞类型。差异表达分析方法，如 Wilcoxon 秩和检验、MAST、二模态似然比检验和 DESeq2，专注于识别统计显著的 DEGs。
需要注意的是，这些方法都需要具有治疗前和治疗后数据集，而 scRank 只需要治疗前数据集和药物信息，不需要任何治疗后数据。因此，为了将 scRank 与这些方法进行比较，我们首先将它们应用于合成的配对条件数据集（场景 4，图 3A）。数据集包含两种条件：疾病和治疗。对于疾病条件，四种细胞类型在模块 1 的活性上表现为梯度增加，代表疾病中异常激活的基因模块。考虑到引入药物抑制模块 1 中的基因，模块 1 的活性应当降至各细胞类型的基础水平。我们推测，基因模块活性较高的细胞类型将受到药物更大的影响，从而生成细胞类型的真实排名。因此，采用 Spearman 秩检验来比较预测的细胞类型排名与真实排名。通过增加数据集中的细胞数量，各种方法的细胞类型排名与真实情况一致。然而，scRank 在使用四个基准数据集进行的推断中表现出稳健性能，在推断药物反应细胞类型方面优于现有方法（图 3A）。
图 3：scRank 优于现有方法的性能
(A) 模拟数据中具有治疗前和治疗后条件的四种细胞类型。跨细胞类型和条件的共表达网络以热图表示。scRank 计算的疾病条件下各细胞类型的扰动评分以箱线图表示（最小值，第 25 百分位数，中位数，第 75 百分位数和最大值）。每种细胞类型的数据点数量为 25。p 值使用 Wilcoxon 检验计算。∗p ≤ 0.05 和 ∗∗∗p ≤ 0.001。scRank 与现有细胞类型优先级方法（Augur，通过 Wilcoxon 秩和检验的差异表达分析）在模拟数据中的性能比较通过点线图表示，细胞数量递增。星号表示特定数据集中所有细胞类型的最佳优先级。
(B) 具有真实靶细胞类型的数据集的低维投影。靶细胞类型用虚线圈出。
(C) scRank 与其他方法（Augur 和差异表达分析方法，包括 Wilcox, MAST, bimod 和 DESeq2）在九种响应细胞类型和两种非响应细胞类型中的性能比较。星号表示与真实情况一致的响应和非响应细胞类型的最佳优先级。
(D) 各方法的平均性能。数据以箱线图表示（最小值，第 25 百分位数，中位数，第 75 百分位数和最大值）。每种方法的数据点数量为 9
(E) scRank 与其他药物反应预测方法（beyondcell 和 scDEAL）在体外样本中的性能比较，平均准确率分别为 scRank 的 77.8%，beyondcell 的 72.2% 和 scDEAL 的 64.8%。数据以箱线图表示（最小值，第 25 百分位数，中位数，第 75 百分位数和最大值）。每种方法的数据点数量为 18。beyondcell 方法包括三种策略（beyondcell, bc_SwitchBinary, bc_BCS）来排名细胞类型。scDEAL 方法包括两种策略（scDEAL, scDEAL_binary）来排名细胞类型。
(F) 箱线图表示 scRank 预测的响应和非响应细胞系的排名，涵盖了 53 种药物中 179 个细胞系的分析。每组的数据点数量为 159。p 值通过单侧配对 Wilcoxon 检验计算。
识别髓母细胞瘤中对 vismodegib 敏感和耐药的肿瘤亚型
作为在真实 scRNA-seq 数据集上的基准扩展，我们选择了一个基准数据集，该数据集是从携带定义的药物敏感和耐药细胞类型的小鼠产生的，以展示 scRank 如何识别药物靶向细胞类型及其相关的生物机制（图 4A）。使用预处理条件下的细胞和 vismodegib 直接靶基因 Smo 作为输入数据（图 4A），scRank 成功识别了药物敏感和耐药的肿瘤亚型（Node_B 和 Node_A），分别基于其最高和最低的扰动评分，这与原始论文的发现一致（图 4B）。
图 4：验证髓母细胞瘤中 vismodegib 药物敏感和耐药肿瘤亚型的预测
(A) 采用两种治疗条件下携带髓母细胞瘤的小鼠的 scRNA 数据。scRank 的输入为车辆 scRNA 数据和 vismodegib 的靶基因 Smo。t- 分布随机邻嵌（t-SNE）图按每种细胞类型中预测的最高排名百分比着色。
(B) 两种条件下原始样本的 t-SNE 图按细胞类型着色。Node_A 到 Node_D 代表肿瘤细胞。根据治疗前后相对差异，Node_A 和 Node_B 分别定义为耐药和敏感肿瘤亚型，如原始研究中所述。
(C) 预处理条件下 Node_A 和 Node_B 的 Smo 相关子网络。网络节点按基因模块类别着色。
(D) 两个肿瘤亚型的 GRN 热图及其差异（治疗后减去治疗前）
(E) 不同模块基因的基因本体（GO）分析。
(F) Node_A 和 Node_B 中差异表达基因的热图。条形代表每个样本中的平均表达。
(G) 基因集合富集分析（GSEA）确定的 Node_B（上部）和 Node_A（下部）中显著激活的通路。
(H) Node_A 和 Node_B 标志基因的生存概率分析。在耐药肿瘤亚型（Node_A）中，标志基因的高表达表现出较低的生存概率。在药物敏感肿瘤亚型（Node_B）中，标志基因的高表达表现出较高的生存概率。p 值通过双侧 log-rank 检验计算。
展示对心肌梗死中丹参酮 IIA 响应的巨噬细胞亚群及其另一潜在靶标
巨噬细胞在导致心肌梗死的病理生理过程中发挥多种作用，表现出广泛的促炎和抗炎效应，这由其异质性表型决定。这些不同的巨噬细胞亚群也可能对药物表现出不同的反应。因此，识别关键的药物响应巨噬细胞亚群对于加深对药理机制的理解和发现细胞水平的新治疗靶标至关重要。因此，我们将 scRank 应用于我们之前研究中的小鼠心脏巨噬细胞的 scRNA-seq 数据集，包含三种条件（假手术、心肌梗死和丹参酮 IIA 治疗）（图 6A）。我们使用心肌梗死的 scRNA-seq 数据和丹参酮 IIA 的直接靶基因 Ctsk 作为输入数据。如图 6B 所示，巨噬细胞 5（M∅-5）是对丹参酮 IIA 响应的最高排名细胞类型，这与我们的原始结论一致，即该化合物可以显著减少 M∅-5 的比例和炎症效应（图 6C）。然后，我们放大了基于治疗前条件数据构建的最高排名和最低排名巨噬细胞亚型的 Ctsk 为中心的子网络。我们发现子网络被模块化为三个基因模块，药物靶基因 Ctsk 位于模块 1 中（图 6D 和 6E）。与 M∅-11 相比，M∅-5 表现出更高的模块 1 活性，主要富集在炎症相关功能中，通过基因本体分析确定（图 6D 和 6F）。此外，我们观察到丹参酮 IIA 显著降低了 M∅-5 中模块 1 的活性，但未在 M∅-11 中观察到，表明 M∅-5 亚群是治疗靶标。
图 6：利用 scRank 识别对丹参酮 IIA 响应的巨噬细胞亚群及其潜在靶标
(A) 来自 C57BL/6 小鼠心脏左心室的免疫细胞 scRNA 数据，包括三种条件（假手术，左前降支结扎 3 天后导致心肌梗死（MI），以及结扎和丹参酮 IIA 治疗 3 天后）。
(B)MI 和丹参酮 IIA 靶基因 Ctsk 的 scRNA 数据是 scRank 的输入。t-SNE 可视化的 MI 数据中，细胞类型根据通过 scRank 预测的优先级着色，从红色（最高排名）到黑色（最低排名）。用虚线圈出最高排名的巨噬细胞亚型 M∅-5。
(C)t-SNE 可视化数据，结合组根据实验条件着色。先前验证的巨噬细胞亚型 M∅-5 作为丹参酮 IIA 的真实靶细胞类型，用虚线圈出。
(D)在不同条件下（治疗前和丹参酮 IIA 治疗后）最高排名和最低排名巨噬细胞亚型（M∅-5 和 M∅-11）的 GRN 热图，包括 131 个 Ctsk 相关基因。
(E)在治疗前和治疗后条件下，M∅-5（顶部）和 M∅-11（底部）之间的药物靶相关模块的比较，根据 47 个模块 1 基因的平均边权重通过配对双侧 Wilcoxon 检验评估。
(F)不同模块基因的 GO 分析。
(G)scRank 用于基因排名的示意图。scRank 计算的 M∅-5 的 GRN 中每个基因的扰动评分，输入是靶细胞类型的 GRN，输出是基因的排名。
(H)不同实验条件下预测靶基因在 M∅-5 中的表达变化。预测靶基因在 M∅-5 中的表达变化以箱线图表示（最小值，第 25 百分位数，中位数，第 75 百分位数和最大值）。每种条件下的数据点数量分别为 93, 896 和 140。星号表示表达值显著上调（MI 对比假手术）或下调（MI 对比丹参酮 IIA） (∗∗∗p < 0.0001)。
(I)显示 CTSB 活性位点与丹参酮 IIA 之间形成氢键的 3D 结构和结合模式。
(J)显示丹参酮 IIA 与 CTSB 相互作用的传感器图。
(K) 巨噬细胞标记物 CD68（绿色）、CTSB（红色）和 DAPI（4′,6- 二脒基 -2- 苯基吲哚；蓝色）的免疫荧光染色，显示丹参酮 IIA 治疗抑制巨噬细胞 CTSB 水平。
讨论
准确识别对药物治疗负责的关键细胞类型仍然是一个挑战。大多数推断治疗反应的方法集中于已知药物敏感性标志物的基因表达，可能错过药物靶标与下游信号基因之间的关键关联。这些关联对于理解异质药物反应非常重要，因为药物扰动不仅对单一靶标产生影响，还会影响信号转导和 GRNs。
在这项研究中，我们提出了 scRank，这是一种使用 tpGRN 对药物响应细胞类型进行排序和推断的新方法。其基本假设是，不同细胞类型的治疗反应取决于其与药物靶标相关的内在基因网络。scRank 的主要目标是模拟 GRN 中的计算机模拟药物扰动，并量化其在不同细胞类型中的效果。为此，scRank 创建了 tpGRNs 以模拟计算机模拟药物扰动，并将 tpGRN 与原始 GRN 进行比较，以测量药物引入的网络变化，从而对细胞类型进行排序。比较网络的常用方法是使用流形对齐算法将不同网络的节点投射到共享的流形空间，在那里计算节点相似性以评估网络结构的差异。然而，这种方法无法捕捉 tpGRN 中的动态变化，这在准确测量药物引起的网络变化方面可能是一个挑战。为了解决这个问题，我们将药物扰动在 tpGRN 中建模为具有全局和局部效应。具体而言，扰动的靶基因可以在每个基因节点上引起全局变形，结果变形可以沿着靶点为中心的子网络扩散。在实践中，距离与药物靶标及其相关的子网络信号相结合，以量化药物效应。虽然 tpGRN 的创建可能简化了响应药物扰动的复杂网络重组，但其设计旨在通过网络中的基因 - 基因相关性反映潜在药物靶标的遗传扰动。我们对 GRNs 和模拟的 tpGRNs 与来自处理细胞的 GRNs 进行比较，验证了这种方法可以准确模拟药物效应。我们所描述的方法，结合药物和生物网络以模拟和评估计算机模拟药物扰动，是基于未处理的 scRNA-seq 数据推断药物响应细胞类型的首例。
研究的局限性
scRank 的一个局限性是其主要集中于作为抑制剂作用的药物，这部分是由于可用数据的性质。我们目前访问的大多数数据都涉及抑制性药物，这指导了 scRank 的开发和测试主要在这种背景下。然而，认识到扩大 scRank 范围的需求，我们还对激动剂进行了初步测试，调整 scRank 的方法以模拟其效果。具体而言，为了在网络中反映激动剂的激活效应，我们将从靶基因节点出发的边缘权重设置为其最大可能值，从而生成 tpGRN。随着更多激动剂药物的单细胞数据变得可用，我们计划进一步完善和扩展 scRank 的能力。