scDrug：从scRNA-seq到药物反应预测

scRNA-seq技术允许在转录组水平上对数千个细胞进行测量。scRNA-seq正在成为研究肿瘤微环境中细胞成分及其相互作用的重要工具。scRNA-seq也被用于揭示肿瘤微环境模式与临床结果之间的关联，并在复杂组织中剖析药物治疗的细胞特异性效应。scRNA-seq的最新进展推动了疾病和治疗靶点生物标志物的发现。虽然已经提出了利用scRNA-seq数据的基因表达来预测药物反应的方法，但还需要一个从scRNA-seq分析到药物发现的集成工具。scDrug作为一个完整pipeline，包括生成scRNA-seq聚类和预测药物治疗方法。scDrug管道由三个主要模块组成：用于鉴定肿瘤细胞亚群的scRNA-seq分析、细胞亚群的功能注释和药物反应预测。

来自：scDrug: From single-cell RNA-seq to drug response prediction
项目：https://github.com/ailabstw/scDrug

背景概述

scRNA-seq被用于分析高分辨率细胞组成，从而发现肿瘤异质性，并为定制化生物学任务提供了前所未有的机会。恶性肿瘤细胞表达特征的细节为药物治疗提供了靶点依据。药物重利用是基于已批准或正在研究的药物，开发针对不同疾病的治疗策略。为了连接药物发现和scRNA-seq分析两个领域，开发了scDrug。

scDrug可以从scRNA-seq分析到药物反应预测。在scDrug中，首先构建了scRNA-seq分析管道，用于对scRNA-seq数据进行全面分析。实现了在Python环境下对肿瘤细胞进行亚群聚类。

接下来，整合了两种不同的方法来预测针对癌细胞亚群的药物治疗，使用包括LINCS， GDSC和PRISM在内的公共数据集来全面表征癌细胞系的分子特征。具体来说，一种方法预测药物对特定肿瘤簇的敏感性，另一种方法预测药物对肿瘤簇的综合作用。scDrug提供预测结果，供领域专家对所选药物进行评价。实验结果表明，scDrug可以成功捕获细胞对药物治疗的反应。总之，scDrug允许研究人员探索肿瘤细胞的异质性，并找到有效治疗的候选药物。

数据和方法

scRNA-seq数据预处理

第一步是scRNA-seq数据分析，包括Scanpy的数据预处理、MAGIC对输入进行插补、Harmony进行批次校正、Louvain进行聚类、Scanpy差异表达基因DEG鉴定、GSEAPY功能富集注释，scMatch细胞类型注释。

在数据预处理中，过滤掉表达少于200个基因的细胞和表达少于3个细胞的基因，将线粒体基因占比的细胞保持在30%以下。剩余的数据进行归一化到每细胞10000个总数，再进行自然对数变换，高变基因搜索，缩放到单位方差和零均值。一旦需要数据输入，scDrug还集成了MAGIC来输入缺失值。接下来，应用主成分分析(PCA)，并根据需要使用Harmony消除批次效应。然后，计算前20个主成分的邻居图，并使用Louvain算法将细胞聚类成组。

自动分辨率聚类

为了确定聚类的分辨率，用户可以选择手动或自动分配。在自动模式下，对间隔为0.2的[0.4,1.4]区间内的分辨率值计算了chooseR中描述的基于次抽样的鲁棒性评分。对于给定的分辨率，使用定义为1.0的距离矩阵减去5次聚类的共聚类频率来计算平均轮廓分数，每次聚类在数据集的80%的随机子集上执行，不进行替换。将得分最高的分辨率作为最优聚类分辨率。

差异基因分析，细胞注释，功能富集

聚类后，scDrug使用默认参数的scanpy函数rank_genes_groups对每个聚类的基因进行排序，以识别DEGs。然后，scDrug用GSEAPY进行功能富集。此外，使用human GO_Biological_Process_2021库对log2倍变化大于2且p值和调整后的p值均低于0.01的DEG执行enrichment。对于细胞类型注释，使用所有基因表达，并计算其细胞的平均表达量作为每个簇的基因表达谱（GEP，gene expression profile）。接下来，应用scMatch，根据参考数据集的GEP对簇进行注释。

基于scRNA-seq数据分析的输出，包括一个AnnData对象，一个基因表达谱GEP，批校正、聚类和细胞类型注释结果的UMAPs，以及DEGs和GSEA文件。

生存分析

为了预测每个聚类对患者生存的影响，首先选择每个簇的前20个差异表达基因作为簇的特异性基因签名。然后，从TCGA数据库中下载不同癌症患者的bluk RNA profiles和相应的临床信息。为了评估每个患者的肿瘤簇活性，为每个患者构建了一个expression table，每个列代表一个簇的基因特征。对于每一簇及其所选的20个基因中的某一个，如果该患者的该基因表达高于所有患者的中位数表达，则赋值为1；否则，该值设置为0。按列求和（以下称为“activity score”）表示患者中每个簇的激活水平。对于每一簇，如果患者的活动得分在最高或最低四分位数，则将其分为“高表达”簇和“低表达”簇。最后，用Kaplan-Meier曲线和log-rank分析的p值比较两组的生存率（图S1）。下面示例中以A类型细胞为例，对比了A簇激活高低的两组患者，然后结合临床信息得到生存曲线。这里的生存分析是为了找到疾病相关的簇。
figs1

图S1：从 TCGA 中评估了 scRNA-seq 数据识别的每种细胞类型特异性特征的表达水平。对于每个簇（细胞类型），将患者分为两组，一组为高表达，另一组为低表达。最后，用 Kaplan-Meier 曲线和 log rank p 值分析比较了这两组的生存情况。

药物反应预测

在scDrug管道中，使用第一步生成的AnnData对象，并应用CaDRReS-Sc（Predicting heterogeneity in clone-specific therapeutic vulnerabilities using single-cell transcriptomic signatures）进行药物反应预测。CaDRReS-Sc是一个基于scRNA-seq数据的强大的癌症药物反应预测的机器学习框架，它估计细胞簇的halfmaximal inhibitory concentration（IC50）。基于CaDRReSSc框架，提供了两种预训练的预测模型GDSC和PRISM，用于预测细胞簇的药物反应。

这两个模型是使用GDSC和PRISM数据集的基因表达和药物反应数据，通过无样本偏差的目标函数进行训练的。通过计算实际药物反应值和预测药物反应值的Spearman相关系数来评估预测性能。按照升序排列，scDrug剔除了drug-wise系数低于第一个四分位数系数的药物。

“IC50”（半数最大抑制浓度）：表示在一定浓度下，化合物能够抑制蛋白质活性的程度

药物反应训练数据

对于GDSC模型，scDrug使用了226种药物在1074种癌细胞系中的反应数据（测量的IC50），数据来源于CaDRReS-Sc的GDSC数据集，GDSC数据集作为训练数据集。对于PRISM模型，scDrug使用PRISM Repurposing数据集（19Q4版本）作为训练数据，该数据集包含1448种药物对480种细胞系的反应。PRISM数据集以剂量-反应曲线下的面积（AUC）提供药物反应（不是来自IC50）。

GEP作为特征

对于GDSC模型，我们使用GDSC数据库中CaDRReS-Sc提供的1018个癌细胞系的基因表达数据，选择所有细胞系中共有的17419个基因作为特征基因进行模型训练。对于PRISM模型，从DepMap Portal（https://depmap.org/portal/）下载CCLE（Cancer Cell Line Encyclopedia）表达数据（21Q3版本），包含1,379个细胞系和19,177个基因。选择表达与PRISM AUC相关且绝对Pearson相关系数至少为0.2的8087个基因作为特征基因。scDrug计算了每个特征基因在细胞系间平均表达量的log2表达倍数变化。

预测模型的框架

为了预测细胞簇的IC50，scDrug计算了相对于AnnData的平均基因表达值的log2倍变化，并预测了每个细胞的IC50值。然后，平均IC50预测值确定每个簇的IC50。或者利用簇和其他簇之间的log2倍变化，直接预测簇的IC50。

模型从转录组学和药物反应中学习了潜在的药物-基因组学关系。CaDRReS-Sc中提出的模型定义为： $\widehat{s}_{iu}=\mu+b_{i}^{Q}+b_{u}^{P}+q_{i}\cdot p_{u}=\mu+b_{i}^{Q}+b_{u}^{P}+q_{i}(x_u W_{P})^{T}$ 其中， $s_{iu}$ 是药物 $i$ 对细胞系 $u$ 的观测药物反应（IC50）， $\widehat{s}_{iu}$ 表示预测的药物反应， $\mu$ 为总体平均药物反应， $b_{i}^{Q}$ 和 $b_{u}^{P}$ 分别是药物 $i$ 和细胞系 $u$ 的偏置项， $q_{i},p_{u}\in R^{f}$ 表示药物 $i$ 和细胞系 $u$ 在latent space下的f-dim表征。 $W_{P}\in R^{d\times f}$ 是将基因表达水平 $x_u\in R^{d}$ 投影到latent space的变换矩阵， $d$ 为基因数。也有简化的： $\widehat{s}_{iu}=b_{i}^{Q}+q_{i}\cdot p_{u}$ ，目标函数定义为： $L(\theta)=\frac{1}{2K}[\sum_{i}\sum_{u}(s_{iu}-\widehat{s}_{iu})^{2}+\lambda\sum_{d}||w_{d}||^{2}+\lambda\sum_{i}||q_{i}||^{2}]$ 其中， $K$ 是drug-cell pairs的总数， $\lambda$ 是L2正则化系数， $w_d$ 是 $W_{P}$ 中的向量。模型预测流程见图S2。
figs2