2.基于Python的单细胞转录组数据处理-批次整合

scRNA-seq数据天然存在批次效应，批次效应是测量的基因表达水平的变化，这是测量不同组或“批次”的结果。例如，如果两个实验室从同一队列中采集样本，但这些样本的解离方式不同，则可能会出现批次效应。如果实验室 A 优化其解离方案以解离样品中的细胞，同时最大限度地减少对细胞的压力，而实验室 B 没有这样做，那么 B 组数据中的细胞很可能会表达更多与压力相关的基因。

注意，批次的划分情况很多，并不是说彻底消除批次效应就是好的，例如我们在探究肿瘤与性别的关系的时候，那么性别就应当保留，而不是当成批次效应去除。

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

整合模型的分类

批次整合模型一般包含三个步骤：1.降维，2.消除批次效应，3.生成优质嵌入。降维一般使用PCA，PCA可以降低数据的信噪比，提取细胞的差异信息。消除批次效应是模型的关键，这些方法经过长期发展被归纳为4类：全局模型，线性嵌入模型，基于图的模型，深度学习模型。

• 全局模型：把批次效应建模为加法或者乘法，ComBat是全局模型
• 线性嵌入模型：利用奇异值分解（SVD）以及其变体来获得嵌入数据，然后在嵌入中跨批次查找相似细胞的局部邻域，并以局部自适应方式（非线性）校正批次效应。常见方法包括MNN，Seurat，Harmony，Scanorama等
• 基于图的模型：这类方法通常运行速度最快。使用最近邻域图来表示每个批次的数据。通过强制连接不同批次的细胞，然后修剪细胞类型组成的边，可以校正批次效应，常见方法为BBKNN，这类方法之所以称为基于图的方法并且速度快，是因为它们不会改变细胞的embedding，而是通过修剪邻域图的边获得新的拓扑结构（新的KNN），新的KNN用于后续的聚类和可视化
• 深度学习模型：大多数深度学习模型基于VAE架构，在latent space中校正批次

下面我们比较不同的批次整合方法，我们始终要明白：不存在适用于所有场景的最佳方法。对于一些简单的批次效应校正任务来说，Seurat及Harmony的表现始终不错，而对于一些复杂的批次效应校正任务来说，scVI、scGen、scANVI和Scanorama则表现地更为出色。我们可以使用scib包来评估数据整合的效果，scib包中包含了10个指标。

查看批次效应

我们首先导入数据：

import omicverse as ov
print(f"omiverse version: {ov.__version__}")
import scanpy as sc
print(f"scanpy version: {sc.__version__}")
import scib
print(f"scib version: {scib.__version__}")
import scvi
print(f"scvi version: {scvi.__version__}")
ov.ov_plot_set()

#读取数据
adata1=ov.read('./data/neurips2021_s1d3.h5ad')
adata1.obs['batch']='s1d3'
adata2=ov.read('./data/neurips2021_s2d1.h5ad')
adata2.obs['batch']='s2d1'
adata3=ov.read('./data/neurips2021_s3d7.h5ad')
adata3.obs['batch']='s3d7'

#合并数据集
adata=sc.concat([adata1,adata2,adata3], merge='same')
adata.obs['batch'].unique()
# 删除多余的metadata
del adata.obsm
adata.obs = adata.obs[['cell_type', 'batch']]
adata.write_h5ad("./data/neurips2021_scib.h5ad", compression='gzip')

然后进行预处理：

adata = sc.read("./data/neurips2021_scib.h5ad")
#质控
adata=ov.pp.qc(adata,
              tresh={'mito_perc': 0.2, 'nUMIs': 500, 'detected_genes': 250},
              batch_key='batch')
#存放原始数据
ov.utils.store_layers(adata,layers='counts')
#归一化/高可变基因筛选
adata=ov.pp.preprocess(adata,mode='shiftlog|pearson',
                       n_HVGs=3000,batch_key='batch')

adata.raw = adata
adata = adata[:, adata.var.highly_variable_features]

adata.write_h5ad("./data/neurips2021_scib_1.h5ad", compression='gzip')

print(adata)
"""
AnnData object with n_obs × n_vars = 26950 × 3000
    obs: 'cell_type', 'batch', 'nUMIs', 'mito_perc', 'detected_genes', 'cell_complexity', 'n_genes', 'doublet_score', 'predicted_doublet', 'passing_mt', 'passing_nUMIs', 'passing_ngenes'
    var: 'feature_types', 'gene_id', 'mt', 'n_cells', 'percent_cells', 'robust', 'mean', 'var', 'residual_variances', 'highly_variable_rank', 'highly_variable_nbatches', 'highly_variable_intersection', 'highly_variable_features'
    uns: 'scrublet', 'layers_counts', 'log1p', 'hvg'
    layers: 'counts'
"""

通过在UMAP中按照batch和cell_type着色来直观感受批次效应：

import omicverse as ov
import scanpy as sc
import scib
import scvi
ov.ov_plot_set()
import matplotlib.pyplot as plt

adata = sc.read("./data/neurips2021_scib_1.h5ad")
ov.pp.scale(adata)
ov.pp.pca(adata, layer='scaled',n_pcs=50)

# GPU可视化UMAP, 自动包含计算邻域图
adata.obsm["X_mde_pca"] = ov.utils.mde(adata.obsm["scaled|original|X_pca"])
print(adata)

ov.utils.embedding(adata,
                basis='X_mde_pca',frameon='small',
                color=['batch','cell_type'],
                show=False)

对于CD14+ Mono，细胞被分为了不同的区域，这对应着分别属于不同的样本来源（批次），但实际上这些批次应该是对齐的。

全局模型整合批次-ComBat

这里使用全局模型的ComBat，需要提供批次标签：

adata_combat=ov.single.batch_correction(adata,batch_key='batch',
                                        methods='combat',n_pcs=50)
print(adata)
"""
AnnData object with n_obs × n_vars = 26950 × 3000
    obs: 'cell_type', 'batch', 'nUMIs', 'mito_perc', 'detected_genes', 'cell_complexity', 'n_genes', 'doublet_score', 'predicted_doublet', 'passing_mt', 'passing_nUMIs', 'passing_ngenes'
    var: 'feature_types', 'gene_id', 'mt', 'n_cells', 'percent_cells', 'robust', 'mean', 'var', 'residual_variances', 'highly_variable_rank', 'highly_variable_nbatches', 'highly_variable_intersection', 'highly_variable_features'
    uns: 'hvg', 'layers_counts', 'log1p', 'scrublet', 'scaled|original|pca_var_ratios', 'scaled|original|cum_sum_eigenvalues', 'batch_colors', 'cell_type_colors'
    obsm: 'scaled|original|X_pca', 'X_mde_pca', 'X_combat'
    varm: 'scaled|original|pca_loadings'
    layers: 'counts', 'scaled', 'lognorm'
"""

adata.obsm["X_mde_combat"] = ov.utils.mde(adata.obsm["X_combat"])
ov.utils.embedding(adata,
                basis='X_mde_combat',frameon='small',
                color=['batch','cell_type'],show=False)

上述过程的运行时间为：3m3.4s，整合结果如下所示：

可以发现ComBat的整合效果比较微弱，批次效应依然明显存在。

线性嵌入模型

Harmony整合

Harmony，是一种稳健灵活的多数据集整合算法，可满足无监督 scRNAseq 联合嵌入的关键挑战：1.扩展到大型数据集、2.识别微小簇、3.跨模态整合。Harmony的算法如下：

Harmony需要提供批次标签：

adata_harmony=ov.single.batch_correction(adata,batch_key='batch',
                                        methods='harmony',n_pcs=50)

adata.obsm["X_mde_harmony"] = ov.utils.mde(adata.obsm["X_harmony"])
ov.utils.embedding(adata,
                basis='X_mde_harmony',frameon='small',
                color=['batch','cell_type'],show=False)

上述过程的运行时间为：6m21.7s，整合结果如下所示：

看起来批次已经被混合在一起，批次整合表现良好。

scanorama整合

Scanorama可以识别和合并所有数据集之间的共享细胞类型，并准确地集成scRNA-seq数据的多个批次。算法如下：

Scanorama需要提供批次标签：

adata_scanorama=ov.single.batch_correction(adata,batch_key='batch',
                                        methods='scanorama',n_pcs=50)

adata.obsm["X_mde_scanorama"] = ov.utils.mde(adata.obsm["X_scanorama"])
ov.utils.embedding(adata,
                basis='X_mde_scanorama',frameon='small',
                color=['batch','cell_type'],show=False)

上述过程的运行时间为：5m3.0s，整合结果如下所示：

深度学习模型

基于变分自编码器-scVI

scVI是一种基于条件VAE的方法，VAE是一种试图降低数据维数的神经网络。CVAE中的条件部分是指在特定协变量（在该场景中协变量为batch）上调节降维过程，以便协变量不会影响低维表示。在基准研究中，scVI已被证明在一系列数据集上表现良好，同时保留了生物变异性。scVI直接对原始计数进行建模，因此为它提供counts矩阵而不是归一化表达矩阵非常重要。

CVAE和对抗训练都能帮助学习协变量无关的表示。

scVI整合如下，需要提供batch标签作为协变量：

# 使用scVI整合
adata1=adata.copy()
scvi.model.SCVI.setup_anndata(adata1, layer="counts", batch_key="batch")
vae = scvi.model.SCVI(adata1, n_layers=2, n_latent=30, gene_likelihood="nb")
vae.train()
adata.obsm["X_scVI"] = vae.get_latent_representation()

# 可视化
adata.obsm["X_mde_scVI"] = ov.utils.mde(adata.obsm["X_scVI"])
ov.utils.embedding(adata,
                basis='X_mde_scVI',frameon='small',
                color=['batch','cell_type'],show=False)

上述过程的运行时间为：32m0.6s，整合结果如下所示：

使用细胞标签作为先验知识–scANVI

scVI在整合时，没有利用细胞标签。在某些情况下，数据集中本身是携带细胞标签的，最常见的应用是将多个研究目标相关的公开数据集联合为一个整体来研究。

当获得细胞标签时，可以使用scANVI进行整合，scANVI是scVI的扩展。目的是保留细胞标签之间的差异，同时消除批次效应。与scVI相比，scANVI往往能更好地保留生物信号，但有时它在消除批次效应方面并不那么有效。虽然该案例为所有细胞提供了标签，但scANVI也可以用半监督方式运行，例如我们仅为某些细胞提供标签。

---

中间过程：标签统一----如果我们使用scANVI来整合已有标签的多个数据集，那么首先执行标签统一非常重要。这是指检查标签在多个数据集中是否一致的过程。例如，一个细胞在一个数据集中可能被注释为“T cell”，但同一类型的细胞在另一个数据集中可能被赋予标签“CD8+ T cell”。如何最好地统一标签是一个悬而未决的问题，一般来说我们可以根据我们的任务来确定，统一为主要的细胞类型是一种比较常见的策略。

---

如果我们还没有训练scVI模型，我们需要首先训练scVI。然后提供标注的细胞信息，对于无标注的细胞，通常来说我们会命名为Unknown或者Unlabelled：

lvae = scvi.model.SCANVI.from_scvi_model(
    vae,
    adata=adata1,
    labels_key="cell_type",
    unlabeled_category="Unknown",
)
lvae.train(max_epochs=20, n_samples_per_label=100)
adata.obsm["X_scANVI"] = lvae.get_latent_representation()

adata.obsm["X_mde_scANVI"] = ov.utils.mde(adata.obsm["X_scANVI"])
ov.utils.embedding(adata,
                basis='X_mde_scANVI',frameon='small',
                color=['batch','cell_type'],show=False)

上述过程的运行时间为：5m43.3s，整合结果如下所示：

基于细胞和特征的联合嵌入

当前大多数单细胞分析流程仅限于细胞嵌入，并严重依赖聚类，同时缺乏对不同特征类型之间的交互进行显式建模的能力。 此外，这些方法是针对特定任务量身定制的，因为不同的单细胞问题的表述方式不同。 为了解决这些缺点，SIMBA是一种图嵌入方法，它将单个细胞及其定义特征（例如基因、染色质可及区域和 DNA 序列）联合嵌入到一个共同的潜在空间中。通过利用细胞和特征的联合嵌入，SIMBA 可以研究细胞异质性、clustering-free的marker特征发现、基因调控推断、批次效应消除和多组学数据整合。SIMBA 提供了一个单一框架，允许以统一的方式表述不同的单细胞问题。

• SIMBA 将实验期间测量的各种特征和细胞表示共同嵌入到共享的潜在空间中，以完成单细胞数据分析中涉及的常见任务和单细胞基因组学中仍然存在的未解决问题的任务。左图是 SIMBA 编码的可能生物实体的示例，包括细胞、基因表达测量、染色质可及区域、TF 基序和reads中发现的 k-mer 序列。中间，SIMBA 将多种类型的实体嵌入到低维空间中。 所有以形状表示的实体（细胞-圆形；峰-三角形；基因-正方形；TF 基序-星形；k-mer-六边形）均按相关细胞类型着色（本例中为绿色、橙色和蓝色）。非信息性特征呈深灰色。在图中，每个实体都是一个节点，边表示实体之间的关系（例如，基因在细胞中表达，染色质区域在细胞中可访问，或者 TF 基序或 k-mer 存在于细胞的开放的染色质区域）。一旦在图中连接起来，这些实体就可以嵌入到共享的低维空间中，特定于细胞类型的实体嵌入到同一邻域中，非信息特征嵌入到其他地方。右边，可以使用 SIMBA 完成的常见单细胞分析任务。 不同的不透明度水平表示不同实验批次或模态的细胞。实线表示实验测量的边。虚线表示计算推断的边。

我们更换一个新的环境单独运行SIMBA，在此之前，保存之前的adata：

adata.write_h5ad("./data/neurips2021_scib_2.h5ad", compression='gzip')

然而，由于环境配置遇到困难，这里就不展示SIMBA的整合结果了。

整合结果评估

我们希望得到定量评估结果，比较不同算法的过程我们会称之为基准测试，在这里，我们将使用scib来执行这一过程：

adata.obsm['X_pca']=adata1.obsm['scaled|original|X_pca'].copy()
from scib_metrics.benchmark import Benchmarker
bm = Benchmarker(
    adata,
    batch_key="batch",
    label_key="cell_type",
    embedding_obsm_keys=["X_pca", "X_combat", "X_harmony",
                         'X_scanorama','X_scVI','X_scANVI'],
    n_jobs=8,
)
bm.benchmark()

bm.plot_results_table(min_max_scale=False)