单细胞多样本整合之scVI和scANVI
我也有用python做过BBKNN和Harmony：
scVI（Single-Cell Variational Inference）是基于深度学习的整合算法，多次被评价为最佳的整合算法之一。scVI发表于2018年的Nature Methods上，后续经过多次的更新改进，文章题目是《Deep generative modeling for single-cell transcriptomics》。
而scANVI（single-cell ANnotation using Variational Inference）算法和scVI类似，是同一个团队做的算法，发表在2021年的Mol Syst Biol上，文章题目是《Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models》。
注意，scVI和scANVI都是用于单细胞转录组数据分析的深度生成模型，它们都是scvi-tools包的一部分。scVI可以用于数据集成、降维、聚类、差异表达等任务，而scANVI可以用于半监督的细胞类型注释，即利用部分已知的细胞类型标签来推断其他细胞的状态（需要预先对每个批次的单细胞数据进行预注释）。scANVI可以看作是在scVI的基础上增加了一个分类器，因此可以从预训练的scVI模型初始化。
- 官网教程在https://docs.scvi-tools.org/en/stable/tutorials/notebooks/harmonization.html，
- github在https://github.com/scverse/scvi-tools
需要强调的是用scVI/scANVI整合分析的话比较吃资源，如果有GPU的话速度会快5-10倍；如果没有GPU，仅用CPU运行scVI，还需谨慎。因为CPU模式的scVI运行速度实在是慢，特别是单细胞大样本数据的运行。
一. 环境部署
mamba create -n scvi python=3.9
which python
#/home/data/fuli09/miniconda3/envs/scvi2/bin/python
/home/data/fuli09/miniconda3/envs/scvi2/bin/python -m pip install scvi-tools anndata numpy scanpy scib certifi scib-metrics pymde scvi-colab
#Please be sure to install a version of PyTorch that is compatible with your GPU (if applicable).
还需要在python里运行：
from scvi_colab import install
install()
如果需要衔接Jupyter notebook使用的话，需要在scvi这个环境里安装一下几个插件，然后就能在Jupyter notebook里选择scvi这个环境（ps：我用的是JupyterLab，和Jupyter notebook是差不多的）：
mamba install -y nb_conda_kernels ipykernel
python -m ipykernel install --user --name scvi --display-name scvi
JupyterLab
二. 运行示例数据
import scanpy as sc
import scvi
from rich import print
from scib_metrics.benchmark import Benchmarker
from scvi.model.utils import mde
sc.set_figure_params(figsize=(4, 4))
%config InlineBackend.print_figure_kwargs={'facecolor' : "w"}
%config InlineBackend.figure_format='retina'
Step1. 读入数据
adata = sc.read(
"data/lung_atlas.h5ad",
backup_url="https://figshare.com/ndownloader/files/24539942",
)
adata
AnnData object with n_obs × n_vars = 32472 × 15148
obs: 'dataset', 'location', 'nGene', 'nUMI', 'patientGroup', 'percent.mito', 'protocol', 'sanger_type', 'size_factors', 'sampling_method', 'batch', 'cell_type', 'donor'
layers: 'counts'
Step2. 数据集预处理
该数据集已经按照scIB手稿中的描述进行了处理。通常，scvi工具中的模型期望以与Scanpy/Seurat相同的方式进行标准分析。
值得注意的是高度可变的基因选择问题。虽然scVI和scANVI都可以在运行时使用所有基因，但作者通常建议过滤基因以获得最佳整合性能。同时，这也有助于消除由于批次特异性基因表达而产生的批次特异性变异。
一般来说，用户可以使用Scanpy的标准流程执行高变基因选择。
adata.raw = adata # keep full dimension safe
sc.pp.highly_variable_genes(
adata,
flavor="seurat_v3",
n_top_genes=2000,
layer="counts",
batch_key="batch",
subset=True,
)
Step3. 使用scVI进行整合分析
第一步，根据count数据和batch信息创建模型：
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="batch")
vae = scvi.model.SCVI(adata, n_layers=2, n_latent=30, gene_likelihood="nb")
#Now we train scVI. This should take a couple of minutes on a Colab session
vae.train()
将结果添加至obsm，命名为X_scVI：
adata.obsm["X_scVI"] = vae.get_latent_representation()
最后，对数据集进行聚类，并可视化：
sc.pp.neighbors(adata, use_rep="X_scVI")
sc.tl.leiden(adata)
adata.obsm["X_mde"] = mde(adata.obsm["X_scVI"])
sc.pl.embedding(
adata,
basis="X_mde",
color=["batch", "leiden"],
frameon=False,
ncols=1,
)
image-20230723214135712
因为该数据集已做注释，所以这里可以根据注释结果来评价scVI的效果：
sc.pl.embedding(adata, basis="X_mde", color=["cell_type"], frameon=False, ncols=1)
image-20230723214212112
Step4. 使用scANVI进行整合
作者说，如果输入的数据有预注释的信息，使用scANVI可以产生更准确的整合结果。
由于上述已经训练了一个scVI模型，我们将使用它来初始化scANVI。在初始化scANVI时，我们为它提供labels_key。由于scANVI也可以用于具有部分观察到的注释的数据集，因此我们需要为它指定与未标记细胞对应的类别名称。由于我们没有未标记的细胞，我们可以给它任意一个不是现有细胞类型名称的随机名称。
lvae = scvi.model.SCANVI.from_scvi_model(
vae,
adata=adata,
labels_key="cell_type",
unlabeled_category="Unknown",
)
lvae.train(max_epochs=20, n_samples_per_label=100)
将结果添加为"X_scANVI"
adata.obsm["X_scANVI"] = lvae.get_latent_representation(adata)
结果可视化：
adata.obsm["X_mde_scanvi"] = mde(adata.obsm["X_scANVI"])
sc.pl.embedding(
adata, basis="X_mde_scanvi", color=["cell_type"], ncols=1, frameon=False
)
image-20230723215034705
Step5. 计算整合效果
这里我们使用scIB-metrics包，它包含scIB基准测试套件中使用的度量的可伸缩实现。我们可以使用这些度量来评估集成的质量。
我们可以看到，使用标签信息和scANVI的额外训练改进了捕获生物保护的指标(cLISI, Silhouette标签)，而不会牺牲太多的批校正能力(iLISI, Silhouette批)。
bm = Benchmarker(
adata,
batch_key="batch",
label_key="cell_type",
embedding_obsm_keys=["X_pca", "X_scVI", "X_scANVI"],
n_jobs=-1,
)
bm.benchmark()
bm.plot_results_table(min_max_scale=False)
image-20230723215710408
df = bm.get_results(min_max_scale=False)
print(df)
Isolated labels Leiden NMI Leiden ARI \
Embedding
X_pca 0.50606 0.724169 0.474548
X_scVI 0.540021 0.750638 0.58713
X_scANVI 0.536124 0.771481 0.639793
Metric Type Bio conservation Bio conservation Bio conservation
Silhouette label cLISI Silhouette batch \
Embedding
X_pca 0.557418 1.0 0.878176
X_scVI 0.535189 0.994187 0.900951
X_scANVI 0.588762 0.998546 0.88585
Metric Type Bio conservation Bio conservation Batch correction
iLISI KBET Graph connectivity \
Embedding
X_pca 0.003306 0.223596 0.728504
X_scVI 0.111625 0.339272 0.899708
X_scANVI 0.105794 0.349628 0.922797
Metric Type Batch correction Batch correction Batch correction
PCR comparison Batch correction Bio conservation \
Embedding
X_pca 0.0 0.366716 0.652439
X_scVI 0.891112 0.628534 0.681433
X_scANVI 0.809051 0.614624 0.706941
Metric Type Batch correction Aggregate score Aggregate score
Total
Embedding
X_pca 0.53815
X_scVI 0.660273
X_scANVI 0.670014
Metric Type Aggregate score
三. 评价和总结
我用这个算法跑了一下示例数据，还用了自己的数据也试了一下，使用后的第一个想法是，等我以后有钱了，我一定要买一台高配服务器，加GPU的那种。目前贫穷的我，是用不起scvi-tools这种高贵的算法的（纯CPU跑真的很慢，非常焦虑）。
另外，目前scVI可以用R语言运行了，下期推文进行介绍。
- END -
腾讯云开发者
