5.基于python的scRNA-seq细胞状态分析-细胞扰动

参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

背景

前面学习了不同处理条件下的基因差异表达，和不同处理条件下的细胞组成变化。这延申出另一个问题，我们想在不同处理条件下，看看究竟是什么细胞受到的影响最大。

在这里，有人会想到，直接观察差异表达基因数量，或者哪一类细胞变多了，或者对差异表达基因进行通路富集分析去评估通路的变化，这三种思路都是间接体现细胞的状态，我们需要一个直接的评价指标，去评估条件施加对细胞的影响。

首先阐明扰动的定义：由外部影响所造成的细胞暂时性或永久性改变。最常见的是CRISPR-Cas9实验带来的基因组扰动现象，其次是由外界干预所带来的转录水平或者蛋白水平的改变。

尽管我们可以从实验的角度去衡量扰动的效果，但是对单细胞层面，不同细胞的扰动效果，仍有着较大的发展空间，具体来说有以下几点：

扰动响应：根据不同的外界干预方法，我们可以预测扰动后的组学特征。一般我们可以通过预测特征和真实值的相关性来评估，还可以通过表型的测量值进行评估，例如IC50，不同剂量反应下的AUC（曲线下面积），毒性等。
靶点和机理：我们通常会采用实验手段去分析药物的靶点与机理，但实际上，对于一种未知作用的化合物而言，我们分析其作用模式，也可以通过扰动响应的方式进行。
扰动相互作用：由于干预一般不是独立的，我们有时候会去预测扰动的组合效应，以了解基因组与药物组合之间的相互关联的影响。
化学性质：我们除了通过组学层特征的改变去分析药物的扰动响应外，我们也可以根据现有的化合物知识，包括R基团，药效基团等，通过已有的数据来直接预测药物的扰动效应。

下面将分析两个方面，包括"查找受扰动影响最大的细胞类型"以及"预测单细胞对扰动的转录反应"

这里依然使用kang数据集，基于10x液滴scRNA-seq的PBMC（外周血单核细胞），来自8名红斑狼疮患者在 INF-β 治疗前后 6 小时的测量。

补充：PBMC是一个细胞大类，所以一般机器学习的细胞类型注释其实是细胞亚型注释

将label重命名为condition以提高可读性，并且将ctrl替换成control，stim替换成stimulated：

import omicverse as ov
import scanpy as sc
import pertpy as pt

adata=ov.read('./data/kang.h5ad')
adata.obs.rename({"label": "condition"}, axis=1, inplace=True)
adata.obs["condition"].replace({"ctrl": "control", "stim": "stimulated"}, inplace=True)

识别受扰动影响最大的细胞类型

扰动很少对所有细胞产生相同的影响，Augur是一种量化响应程度的方法。Augur 的目标是根据细胞类型对给定单细胞基因表达数据的实验扰动的反应，对细胞类型进行排序或优先排序。基本思想是，在分子测量领域，对诱导扰动反应强烈的细胞比反应很少或没有反应的细胞类型更容易分为受扰动组和未受扰动组。通过测量每种细胞类型内实验标签（例如处理"treatment"和对照"control"）的预测程度来量化这种可分离性。Augur 训练一个机器学习模型，在多次交叉验证运行中预测每种细胞类型的"处理条件"标签，然后根据衡量模型准确性的指标分数对细胞类型响应进行优先级排序。

Augur有两种分类器：随机森林和逻辑回归。

在Kang中，处理条件有Control和stimulated两种类别，同时我们选择随机森林作为分类器：

ag_rfc = pt.tl.Augur("random_forest_classifier")
loaded_data = ag_rfc.load(adata, label_col="condition", cell_type_col="cell_type")

v_adata, v_results = ag_rfc.predict(
    loaded_data, subsample_size=20, n_threads=4, select_variance_features=True, span=1
)

print(v_results["summary_metrics"])

可视化排序结果和UMAP结果：

lollipop = pt.pl.ag.lollipop(v_results)

ov.utils.embedding(v_adata,
                   basis='X_umap',
                    frameon='small',
                   color=["augur_score", "cell_type", "label"],
                   cmap='Reds',
                   ncols=2)

fig1
fig2
据观察，CD14+单核细胞受IFN-β影响最大，而巨核细胞受影响最小。这与原始论文中各细胞类型差异表达基因的统计分析吻合。既然与差异表达基因统计一致，那为什么不用差异表达基因数量来衡量扰动强弱，这是因为，除了关心扰动强弱外，我们还关心扰动的主要干预基因，这些基因不一定是差异表达最高的基因，但他们对模型的贡献度是可以量化的：

important_features = pt.pl.ag.important_features(v_results)

fig3
这些是扰动相关的所有基因。

不同处理中受扰动影响的细胞类型优先级

在上述分析中，我们仅研究了Ctrl和Stim两个组，但在真实情况中，我们通常会进行多种药物处理，或者是不同时间段的药物干预分析。在Augur中，提供了排列测试来进行差异优先级划分。

这里使用另一个数据集，该数据给小鼠提供可卡因，并在戒断后48小时（withdraw_48h_Cocaine）和15天（withdraw_15d_Cocaine）采集前额皮质，测量技术为scRNA-seq：

adata=ov.read('./data/bhattacher.h5ad')
adata.obs.rename({"label": "condition"}, axis=1, inplace=True)
#ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=2000,)
ag_rfc = pt.tl.Augur("random_forest_classifier")
sc.pp.log1p(adata)

# 分析5天
# Default mode
bhattacherjee_15 = ag_rfc.load(
    adata,
    label_col='condition',
    condition_label="Maintenance_Cocaine",
    treatment_label="withdraw_15d_Cocaine",
    cell_type_col="cell_type",
)

bhattacherjee_adata_15, bhattacherjee_results_15 = ag_rfc.predict(
    bhattacherjee_15, random_state=None, 
    n_threads=4, 
)
bhattacherjee_results_15["summary_metrics"].loc["mean_augur_score"].sort_values(
    ascending=False
)

# 分析48小时
# Default mode
bhattacherjee_48 = ag_rfc.load(
    adata,
    label_col='condition',
    condition_label="Maintenance_Cocaine",
    treatment_label="withdraw_48h_Cocaine",
    cell_type_col="cell_type",
)

bhattacherjee_adata_48, bhattacherjee_results_48 = ag_rfc.predict(
    bhattacherjee_48, random_state=None, 
    n_threads=4,
)

bhattacherjee_results_48["summary_metrics"].loc["mean_augur_score"].sort_values(
    ascending=False
)

scatter = pt.pl.ag.scatterplot(bhattacherjee_results_15, bhattacherjee_results_48)

fig4
在戒断后，Oligo的响应从0.5上升到0.8，这为后续研究提供了一个潜在的启发。对于几乎没有变化的细胞类型，说明戒断对其影响没有提高，有可能是可卡因已经造成了永久性损伤。

预测细胞的扰动响应

除了对已知细胞评估扰动响应，我们还可以根据已有的细胞，预测数据集中未知细胞的扰动响应。对于分析未知的细胞，从实验角度需要重新测序，这带来了额外的代价，scGen是一种自编码器，学习数据的潜在表示，估计对照（未处理）和扰动（处理）细胞之间的差异向量，然后将估计的差异向量添加到对照组感兴趣的细胞类型中，以预测每个单细胞的基因表达反应。

构建模拟数据集

以Kang数据集为例，假设CD4 T cells在治疗后没有被测序仪捕获，构建数据如下，首先选择HVG加快计算速度：

import omicverse as ov
import scanpy as sc
import pertpy as pt

adata=ov.read('./data/kang.h5ad')
adata.obs.rename({"label": "condition"}, axis=1, inplace=True)
adata.obs["condition"].replace({"ctrl": "control", 
                                "stim": "stimulated"}, inplace=True)
adata=ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=2000,)
adata.raw = adata
adata = adata[:, adata.var.highly_variable_features]

遵循假设（治疗后的CD4 T cells没有被捕获），我们删除stimulated下的所有CD4 T cells，命名为adata_t，同时存放stimulated下的所有CD4 T cells，命名为cd4t_stim：

adata_t = adata[
    ~(
        (adata.obs["cell_type"] == "CD4 T cells")
        & (adata.obs["condition"] == "stimulated")
    )
].copy()

cd4t_stim = adata[
    (
        (adata.obs["cell_type"] == "CD4 T cells")
        & (adata.obs["condition"] == "stimulated")
    )
].copy()