我也有用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模型初始化。
需要强调的是用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'
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'
该数据集已经按照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,
)
第一步,根据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
作者说,如果输入的数据有预注释的信息,使用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
这里我们使用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 -