基于人类癌细胞的深度学习模型预测药物反应和协同作用

导语

大多数进入临床试验的药物基本会失败，这通常与对药物反应机制的不完全理解有关。机器学习技术在更好地预测药物反应方面有着巨大的希望，但由于缺乏可解释性和对单一疗法的关注，大多数技术尚未进入临床实践。本文通过开发DrugCell来应对这些挑战，DrugCell是一种可解释的人类癌细胞深度学习模型，接受1,235 种肿瘤细胞系对 684 种药物的反应训练。肿瘤基因型诱导细胞子系统的状态，这些状态与药物结构相结合，以预测对治疗的反应，同时了解药物反应的生物学机制。DrugCell 预测在细胞系中是准确的，并且还能对临床结果进行分层。体外和患者来源的异种移植物。DrugCell 为构建预测医学的可解释模型提供了蓝图。

背景介绍

今天小编为大家带来一篇基于深度学习预测药物反应和协同作用发表在38分+ Cancer Cell的思路。题目为   Predicting Drug Response and Synergy Using a Deep Learning Model of Human Cancer Cells。

数据介绍

研究设计

为了在计算上表示癌症基因型，本文根据癌细胞系百科全书(CCLE) 在基因本体(GO)术语注释的基因中选择了人类癌症中最常突变的前15%基因。这个过程产生了3008个基因，后来被称为“药物细胞基因”，用于模型构建。这些基因被组织成一个嵌套基因集的层次结构，代表不同尺度的细胞子系统，基于从GO生物过程层次结构中提取的术语。如果术语至少有10个DrugCell基因，并且与所有子术语不同，则从GO中保留术语，定义为比任何子术语至少有30个DrugCell基因(part_of和is_a层次术语关系都被考虑在内)。每隔一个术语就从层次结构中删除，而将其子术语直接分配给父术语，以保持层次结构的连接。为了进一步降低模型的复杂性，本文通过删除层次结构底层子系统(没有任何子子系统)之上的所有超过五个父子关系的子系统，将层次结构限制为五个子系统的最大深度。得到的层次结构由2086个子系统组成，定义了用于基因型嵌入的DrugCell分支，也称为VNN。

为了获得足够大的药物基因组学数据集用于模型训练，从癌症药物敏感性基因组学数据库(GDSC)和癌症治疗反应门户(CTRP)中检索原始药物敏感性数据。这些数据涵盖了509,294对(细胞系，药物)。在这些数据中，在两个存储库中冗余测量的24,923对数据完整地保留在训练数据集中，因为这样的重复有利于减少模型过拟合。为了在两个数据集之间进行标准化，本文计算了剂量反应曲线下的面积(AUC)，使AUC = 0表示完全细胞杀伤，AUC = 1表示没有效果。曲线是通过以分段线性方式连接单个响应点而不是使用s型曲线拟合而创建的。然后将此分段线性拟合的AUC归一化为跨越测试浓度范围的零曲线下的面积。计算的AUC值与之前对该数据集的分析高度一致(r2 = 0.87)。为了标准化跨数据集的药物表示，查询了CTRP或GDSC中使用的每种化合物的PubChem条目，以获得基于数据集中提供的药物名称或InChIKey的异构SMILES表记法。在初始搜索中没有匹配的化合物被手动注释。为了在计算上表示化学结构，本文使用RDKit (http://www.rdkit.org/)计算Morgan指纹(半径= 2)，该指纹通过迭代获得分子中每个原子的不同路径，将每种化学结构分解为分子片段。这些片段被散列成一个长度为2048的位向量，用于模型训练。数据集被过滤后，只代表前15%最常突变的基因(n = 3,008)。每个细胞系基因型被表示为一个位载体，横跨3008个DrugCell基因，表明该细胞系中每个基因的突变状态(0 =野生型;1 =变异)。

DrugCell VNN(基因型嵌入分支)是采纳了DCell protoco算法，并进行了轻微修改。每个子系统 在 DCell中，以及在DrugCell的子系统层次结构中，被分配一个数字k 来表示它的多维状态。 子系统的状态，用输出向量O(s)表示。 被定义为它的c 子子系统和 g个直接注释的基因，连接在输入向量I(s)中。

结果解析

01、设计与训练可解释神经网络的用于药物反应预测

细胞药物反应是一种复杂的现象，它取决于生物和化学因素。目前使用这两种因素的药物反应黑箱模型已经开始达到预测性能的极限。因此，本文的目标是设计一个模型，以保持这种高水平的预测能力，同时获得模型预测的机制可解释性。为了在可解释的模型中捕获药物反应的两个决定因素，本文将DrugCell设计为具有两个分支的神经网络(图1A)。第一个分支是VNN，它从基因本体(GO)数据库中记录的2086个生物过程中提取，对人类细胞中分子子系统的层次组织进行建模。从涉及小蛋白质复合物(如β-连环蛋白破坏复合物)到更大的信号通路(如MAPK信号通路)到总体细胞功能(如糖酵解)的每一个子系统，都被分配了一组人工神经元来代表该子系统的状态(图1B)。神经元的连通性被设置为反映生物层次结构，因此神经元只接受来自子子系统的输入，只向父系统发送输出，连接权在训练过程中确定。每个子系统使用多个神经元允许细胞子系统具有多功能，具有不同的状态，能够沿多个维度采用一系列值。层次结构的输入层映射到基因的突变状态。VNN输出的6个神经元对应于层次的根，代表了基于基因型的整个细胞的嵌入状态(图1B)。总的来说，VNN使用了12516个神经元，分层分布在六个不同的层上。DrugCell的第二个分支是嵌入药物摩根指纹的传统人工神经网络(ANN)，这是化学结构的规范向量表示(图1C)。模型的两个分支，即VNN嵌入细胞基因型和ANN嵌入药物结构的输出，被组合在单层神经元中，然后将其整合以生成给定基因型对特定治疗的响应(图1A)。

图1

为了训练模型，本文协调了两个大型癌症药物筛选资源的数据:癌症治疗反应门户(CTRP) v2和癌症药物敏感性基因组学(GDSC)数据库。合并后的数据集包括509294个细胞系-药物对，涵盖684种药物和1235个细胞系。所有主要组织类型均有体现，其中造血和肺谱系最为普遍。每个细胞系基因型用一个二元载体表示，记录癌症中前15%最常突变基因的突变状态(1=突变，0=未突变)(n = 3,008;每个细胞系中位突变基因=73个)。每种药物的化学结构在Morgan指纹载体中平均由81个激活位表示，每个位通常代表少于10个分子片段。drug-cell经过训练，将每个基因型药物对与相应的药物反应联系起来，通过剂量-反应曲线下的面积来测量。

02、可解释的药物反应模型没有性能损失  

本文首先试图在5倍交叉验证中使用预测和观察到的AUC值之间的Spearman相关性来评估drug-cell的预测准确性。所有细胞系-药物对的总准确度为rho=0.80(图2A)。通过单独计算每种药物的预测精度，进一步深入了解，揭示了具有非常高预测精度的药物亚群，总体分布范围更广(范围为-0.29至+0.83，中位数为0.37)。这些精度显著高于弹性网(中位数rho = 0.35)，弹性网是一种最先进的回归技术，用于许多先前的药物反应预测方法(图2B)。DrugCell的逐药预测性能与具有匹配神经元数量、层数和连接的传统黑箱神经网络没有显著差异(图2C)。它也可以与之前将药物的化学特征纳入反应预测(例如，结构和理化性质，如溶解度、亲脂性和分子量)的努力相媲美，并且它优于仅使用生物特征预测反应的模型。最后，由于即使在没有其他信息的情况下，组织类型的知识也可以预测药物反应，本文认为这些模型的一些性能可能是由于它们能够从输入数据(即其突变谱)中识别细胞系的组织类型。因此，本文将DrugCell与仅训练药物结构和组织标签的等效神经网络模型进行了比较。drug-cell的表现大大优于该组织模型(中位数rho = 0.18;图2D)，表明该模型已经从起源组织以外的体细胞突变中获得了信息。

drug-cell预测最准确的化合物来自不同的靶点类别，包括化疗药物和靶向治疗(例如，靶向PLK1的GSK461364，靶向Src的KX2-391;图2E)。DrugCell保持了训练数据的特异性，因为它的预测是针对个别类别的药物的(例如，MEK抑制剂的预测是高度特异性的)，而不是简单地反映一般药物的毒性。药物的预测性能与用于训练的细胞系药物对的数量没有很强的相关性，也与化合物的结构复杂性。本文确实发现，引发更大范围细胞系反应的化合物往往更可预测。同样，单个细胞系和组织类型，引起大范围的反应，通常是高度可预测的。

图2

03、DrugCell学习介导特异性药物反应的机制

在评估了预测能力之后，本文接下来转向机械解释。这一任务得到了两个模型分支的帮助，它们从细胞内化学结构对药物活性的影响(药物包埋，图1A)中剖析了基因型对细胞系统配置的影响(基因型包埋)。本文通过绘制顶部两个主成分(图3A-3E)来直观地检查这些嵌入。从VNN中嵌入的基因型揭示了根据已知赋予特定药物敏感性的突变的基因型分离，例如BRAF中的激活突变(图3A)促进对MEK抑制剂selumetinib的敏感性(图3B)。基因型嵌入还区分了导致耐药的突变，例如EGFR、LKB1或BRAF(图3C)的突变，这些突变赋予了对β-家族抑制剂JQ-1的耐药性(图3D)。本文类似地检查了VNN内单个子系统的drug-cell嵌入，发现许多子系统的活性与实验测量的子系统活性一致，通过使用反相蛋白质阵列(rppa)对蛋白质丰度和磷酸化状态进行独立分析;)。例如，DrugCell在MAPK级联调控的子系统嵌入中准确捕获MAPK通路活性，其与ERK1/2磷酸化显著相关。总体而言，大多数DrugCell子系统与这些子系统的RPPA测量值具有良好的相关性。其他精确捕获的子系统包括蛋白水解、PI3K信号传导调节和细胞周期阻滞。

从人工神经网络中对药物包埋的检查显示，在主要药物靶标类别中，药物根据其作用机制进行了分层(图3E)。化学结构包埋中每对药物之间的距离与它们的整体化学相似性无关，这与前人对药物活性和化学结构的研究一致(Breinig et al.， 2015)。由于训练数据仅由药物和类药物分子组成，因此化学结构包埋并没有根据膜透性、溶解度或药效学性质。总之，这些结果表明，drug-cell能够了解控制药物敏感性和耐药性的基因型的关键特征，以及控制药物生物活性的化学结构特征。

由于DrugCell的VNN是根据包含人类细胞的生物子系统的层次结构构建的，因此它的输出(基因型嵌入)是该层次结构中特定子系统状态变化的结果。为了确定这些子系统中最重要的子系统，本文使用预测能力指标(RLIPP)的相对局部改进，根据子系统状态比其子子系统状态更能预测药物反应的程度对子系统进行评分。作为概念的初步证明，本文使用RLIPP评分来识别细胞对紫杉醇(紫杉醇)反应的重要子系统，紫杉醇是一种稳定微管的药物(图2E和3F)。在紫杉醇得分最高的子系统中，许多子系统是代谢过程(图3G和3H)，包括对cAMP的反应(最高分)以及葡萄糖和葡萄糖反应时的胰岛素分泌。通过检查证实，这些子系统的状态有能力对紫杉醇敏感细胞系和耐药细胞系进行分层(例如，对cAMP子系统的反应，图3I)。鉴于这些潜在的代谢途径，假设紫杉醇的功效可能受到代谢扰动的调节。因此，将A427细胞暴露于三种不同的处理-紫杉醇，糖酵解抑制剂2-脱氧葡萄糖(2-DG)或两者的组合-并发现该组合比任何一种单独的化合物都有效(图3J)。

对下一个(第二重要的)子系统，泛素蛋白转移酶活性的调节进行了类似的分析(图3G)。通过蛋白酶体抑制剂硼替佐米将紫杉醇与干扰泛素依赖性蛋白降解结合使用。本文发现这些处理是拮抗的，这与最近的研究结果一致，表明糖酵解受细胞骨架上泛素连接酶的负物理调节。多西紫杉醇的姊妹化合物泛素和子系统也被鉴定出来。值得注意的是，这些药物细胞通路没有被早期的基因突变分析发现，并且与紫杉醇敏感系与耐药系mRNA表达差异分析发现的通路不同。与drug -cell分析中出现的糖酵解扰动不同，差异表达途径建议的联合治疗不能成功提高紫杉醇的疗效。

除了紫杉醇之外，我们还研究了其他药物的重要子系统，我们发现其中一些子系统与先前确定的药物敏感性机制相对应，而其他许多子系统是值得进一步研究的新途径。特别是，我们使用II型误差方差分析模型对GDSC数据集进行了早期分析，报告了60种泛癌症诊断基因突变的药物。对于许多药物，DrugCell在前子系统(4种药物)或前10个子系统(14种药物，子系统的上0.4百分位)中恢复了先前报道的诊断基因。然而，对于绝大多数药物(56种药物)，DrugCell通过咨询额外的或不同的标记物获得了比以前报道的更好的预测性能。

图3

考虑到新型药物反应途径的范围，本文试图系统地研究指明的机制(图4A)，重点研究MEK1抑制剂曲美替尼(trametinib);olaparib, PARP1抑制剂;坚果素-3，一种稳定和激活p53的MDM2拮抗剂。三种药物靶点(MEK1, PARP1, TP53)的CRISPR敲除与定制CRISPR/Cas9文库中的每个基因的敲除相结合，该文库具有广泛的癌症信号通路(MCF7细胞;图4b)。确定了对每种药物反应的前5个重要子系统(RLIPP分析;图4C-4E)，以及CRISPR文库覆盖的这些子系统中的基因。MAPK1与曲美替尼子系统中基因的组合破坏(图4C)导致的细胞杀伤明显高于随机不重要子系统的基因(图4F)。在奥拉帕尼子系统中，PARP1与基因的组合破坏也观察到类似的细胞杀伤效应(图4G)。相比之下，TP53与nutlin-3子系统基因的组合破坏(图4E)对细胞生长的影响与随机情况没有显著差异(图4H)。这一结果在意料之中，因为TP53敲除与nutlin-3相比具有相反的作用，从而导致p53激活。这些结果与紫杉醇的初步结果一起，为DrugCell确定的顶级反应途径的重要性提供了系统的支持。

图4

04、识别的子系统带来了协同药物组合的机会  

药物协同作用的平行通路抑制理论(Yeh et al.， 2009)认为，如果两种药物抑制调节共同基本功能的不同通路，则两种药物将具有协同作用(图5A)。drug -cell模型的分支结构(图1A)反映了这种平行通路结构，其中药物的生物活性由药物包埋分支学习，平行通路由基因型包埋分支学习(图5B)。因此，预测药物反应的重要子系统可能代表协同药物联合的机会。上述分析正是使用这种平行性来命名联合治疗(即2-DG与紫杉醇协同作用)。

为了进一步探索这一概念，研究使用RLIPP评分对DeepSynergy数据库中调节对25种药物敏感性的子系统进行排名，其中25种药物的所有对都在39个细胞系中进行了测试(图5C)。然后，我们分析了每个化合物的前5个和后5个drug - cell子系统，以指定协同和非协同药物组合。我们发现，与预测的非协同或随机组合相比，drug cell提名的药物组合具有很强且显著的增效细胞杀伤效果(图5D)。 其中一个例子是依托泊苷，一种导致DNA损伤的拓扑异构酶抑制剂。顶层依托泊苷子系统包括主要的激酶信号通路PI3K- akt(调控PI3K活性，PI3K;图5E)和RAF-MEK-ERK (ERK1/ERK2级联负调控，ERK;图5e)。事实上，在DeepSynergy测试的大多数细胞系中，依托泊苷与AKT和MEK具有很强的协同作用(图5F)。本文在A549细胞中使用CRISPR/Cas9基因编辑技术，通过单独或与PI3K-AKT信号通路(PIK3CA)或RAF-MEK-ERK信号通路(MAP2K1)核心基因联合，删除依托opo苷的靶点TOP2，进一步验证了观察到的协同作用。我们观察到，与单基因敲除相比，PIK3CA或MAP2K1缺失TOP2都表现出细胞活力的显著丧失(图5G)。APC的子系统(β-catenin破坏复合体)没有被RLIPP识别(表S5)，没有显示出相同的模式(图5G)。同样，依托泊苷与蛋白酶体抑制剂硼替佐米没有协同作用(图5F)，这与DrugCell没有识别蛋白酶体子系统(图5E)是一致的。 进一步的检查表明，DrugCell捕获的PI3K信号、ERK信号和etopo苷敏感性之间的关系可以通过整合6个基因突变状态的逻辑函数大致近似(图5H和5I;星方法)。其中，FLT1 (Das et al.， 2005)和PIN1 (Mathur et al.， 2011)先前已被证明可调节etopo苷反应，而DUSP1、PIK3R4、SRC和RPS6KA6则未被证实。单独考虑，这些基因中的任何一个在癌细胞系中很少发生突变，预测依托泊苷敏感性与耐药性的能力有限(突变频率0.9%-8.9%;优势比<2;图5j)。然而，作为一个集成电路，这些基因突变聚集在PI3K或ERK子系统上，创造了一个强大的基于网络的药物反应生物标志物(图5j)。我们还注意到，这两种途径仅代表完整drug-cell模型的一部分，该模型预测依托泊苷敏感性的比值比为14.3。

图5

05、药物细胞改善患者来源的异种移植模型的无进展生存  

接下来，本文希望超越细胞系来预测和解释患者来源的异种移植模型(PDX;图6A)。为此，本文查阅了PDX百科全书，其中对399种不同组织类型的PDX肿瘤进行了筛查，共针对40种不同的单一疗法和27种联合疗法。每种PDX的基因型也已确定，并将其提供给DrugCell，以预测对每种单药治疗的反应。如果预测的AUC低于所有PDX-药物对的预测中位数，我们认为PDX肿瘤对治疗敏感(drug- cell(+));否则，该肿瘤被标记为不敏感(DrugCell(−))。drug-cell(+)肿瘤的无进展生存期(PFS)明显长于drug-cell(-)肿瘤(2.19个月vs 1.58个月，p = 9.4 × 10−10,log rank检验)。然而，考虑到这些PDX肿瘤对单药治疗的总体不敏感，对应于两种drug - cell类别观察到的较短的PFS，本文希望评估drug-cell能够建议有效的药物组合的程度。我们使用RLIPP评分对六种主要药物的药物反应进行排序，并将此列表过滤到含有次要药物靶点的子系统。观察到的这些(主要的，次要的)组合的PFS被用来估计沿受试者工作特征曲线(ROC;图6A)。我们发现DrugCell能够准确地识别出PDX肿瘤中与有效药物组合相对应的子系统(auROC = 0.75;图6B)，假阳性和假阴性相对较少(图6C)。

图6

例如，drug - cell对PI3K抑制剂BKM-120的分析发现，ERK1 + ERK2级联的负调控是BKM-120应答的重要子系统，提示PI3K + MAPK途径抑制剂(BKM-120 + encorafenib)的联合作用。与单药治疗相比，该联合治疗显著提高了整个PDX组的PFS(图6D)。同样，DrugCell发现DNA损伤反应，p53类介质的信号转导导致细胞周期阻滞是abraxane反应的重要子系统，建议联合化疗与诱导DNA损伤和细胞周期阻滞的药物(abraxane +吉西他滨)。这种组合同样显著改善了PFS(图6D)。对于未被DrugCell优先考虑的组合(未在RLIPP的前20%的子系统中)，这些组合确实未能显著改善PFS(图6C和6E)。这些结果表明，药物细胞在指导肿瘤患者联合治疗的设计方面具有实用价值。

06、DrugCell预测雌激素受体阳性转移性乳腺癌患者对mTOR和CDK4/6抑制剂的反应
最后，我们试图评估药物细胞是否可以在临床上用于将癌症患者分为反应性和非反应性患者群体。 我们获得并分析了221例雌激素受体(ER)阳性转移性乳腺癌患者的综合临床试验数据，这些患者接受了多轮治疗，除了mTOR抑制剂(依维莫司)或CDK4/6抑制剂(ribociclib)治疗外，还包括ER拮抗剂(氟维司汀)。 对于这项分析(STAR方法)，我们使用预训练的drug - cell模型预测了患者对mTOR或CDK4/6抑制的反应。 如果预测患者对任何一种治疗都敏感，本文认为患者是drug-cell(+)型;如果预测患者对两种治疗都不敏感，我们认为患者是drug - cell(-)型。 DrugCell(+)患者的总生存期明显长于DrugCell(-)患者(48.2个月vs 33.6个月，p=0.018;图7a)。

接下来，我们通过对mTOR和CDK4/6抑制剂进行RLIPP分析，探究了drug -cell(+)和drug-cell(−)患者之间差异敏感性的机制。值得注意的是，我们发现这两种药物反应都受到er相关子系统的调节(图7B和7C)，这与它们在er阳性乳腺癌中的应用一致(Hare和Harvey, 2017;Pernas et al.， 2018)。我们还发现，这两种药物的主要作用机制都在最重要的途径中，PI3K信号对于mTOR抑制剂的反应特别重要，CDK活性对于CDK4/6抑制剂的活性很重要(图7B和7C)。有趣的是，在CDK4/6抑制剂中也发现了TOR信号传导(图7B)，并且在mTOR抑制剂中发现了CDK活性(图7C)，这表明这些药物可能是一种有效的联合治疗，这一发现得到了最近的临床前研究的支持(Michaloglou等人，2018;Occhipinti et al.， 2020)。 

关于特定的基因改变，我们发现DrugCell(+)患者比DrugCell(-)患者更有可能携带AKT1突变(图7D)。相比之下，drug - cell(−)患者具有先前与耐药相关的基因突变，包括ESR1 (Reinert等人，2017)，RB1 (Condorelli等人，2018)和PTEN (Costa等人，2020)(图7D)，这表明我们根据导致治疗耐药的复杂突变模式对患者进行了分层。引人注目的是，AKT1突变状态本身并不能预测治疗反应，AKT1突变患者实际上倾向于更短的总生存期(35.8个月对43.1个月)，尽管这种差异没有统计学意义(图7E)。这一分析说明了drug - cell如何能够比单基因标记研究更精确、更深入地有效指导临床治疗决策。

图7

小编总结

本文作者探索了一个可解释的深度学习模型的结构和功能的人类癌细胞对治疗的反应。这项工作将预测建模推向了药物反应背后的生物机制的系统表示，这是精准医学的一个关键方向。在模型预测之后，获得机械解释使实验学家或临床医生参与对生物功能的推理。例如，对drug - cell的依托泊苷模型的分析确定了一小部分对细胞反应很重要的子系统，并且靶向药物可用。该分析促使我们进行后续实验，以MAPK或PI3K途径靶向遗传和药理学上的拓扑异构酶II;两种组合均表现出显著的协同效应。这种人类推理和后续实验的参与有助于极大地提高对机器学习模型预测的问责制和信任。相比之下，传统的黑箱预测模型只产生一个模型输出——药物反应——没有进一步的信息来建立信任。

未来的工作可能还会选择将突变与其他水平的分子信息(如表观遗传状态、基因表达或微环境影响)整合起来。这种整合可以通过预处理多层信息来完成，以获得每个细胞系或肿瘤的基因评分概况，然后将其输入DrugCell。额外级别的信息也可以通过在现有的神经网络分支上添加新的可见的或传统的神经网络分支来集成。或者，特定突变对基因功能的影响可以通过组合注释依赖耗尽评分(Rentzsch等人，2019)等指标来纳入，或者通过将基因结构域作为层次结构的额外层来纳入。 另一个机会是根据组学数据而不是文献整理(GO)来构建drug - cell系统层次结构，就像之前在出芽酵母中所做的那样(Kramer等人，2014;Ma et al.， 2018)。数据驱动的，而不是文献策划的，层次结构有可能将新的基因子系统关联以及全新的子系统合并到模型中。它也有可能修改和定制GO中通用的子系统定义，以适应与癌症相关的特定环境。例如，我们发现当前形式的DrugCell包含许多子系统，这些子系统具有基于GO命名约定的误导性标签。例如，迷宫发育是曲美替尼的顶级子系统之一，最初令人困惑，但经过进一步检查，它对应于众所周知参与癌症增殖的MAPK级联基因(例如，MAP2K1, MAPK1, GRB2, FGFR2)。将数据驱动的层次结构整合到DrugCell中，提供了重新标记这些子系统并修改其特定基因内容的途径。最后，鉴于DrugCell输入了一个完整的药物结构，它有可能被用来设计新的化合物。利用强化学习在药物设计方面的进步，就有可能针对任何给定的基因组背景设计出具有最大功效的化合物。