Bohrium
robot
新建

空间站广场

论文
Notebooks
比赛
课程
Apps
我的主页
我的Notebooks
我的论文库
我的足迹

我的工作空间

任务
节点
文件
数据集
镜像
项目
数据库
公开
使用深度生成模型scVI进行单细胞批次效应去除实践
python
notebook
深度学习
机器学习
生物信息学
单细胞分析
中文
pythonnotebook深度学习机器学习生物信息学单细胞分析中文
Gongyq
更新于 2024-09-09
推荐镜像 :scvi-tools:v1.0
推荐机型 :c4_m30_1 * NVIDIA P100
赞 1
使用深度生成模型scVI进行单细胞批次效应去除实践
目录
背景介绍
前置准备
加载必要的库
设置全局参数
加载和准备数据
创建并训练模型
保存和加载
获取模型输出
校正批次效应
没有批次校正的可视化
批次校正后的可视化 (scVI)
对scVI潜在空间聚类
批次校正后的差异表达分析

使用深度生成模型scVI进行单细胞批次效应去除实践

scVI论文：https://www.nature.com/articles/s41592-018-0229

📖 上手指南
本文档可在 Bohrium Notebook 上直接运行。如果是在案例广场，你可以直接运行；如果是保存到自己的空间中，你可以点击界面上方按钮 开始连接，选择 scvi-tools:v1.0 镜像及 c4_m30_1 * NVIDIA P100 节点配置，稍等片刻即可运行。

代码
文本

目录

代码
文本

背景介绍

scvi是Nature Methods发表的一个单细胞RNA测序(scRNA-seq)数据的集成分析平台。scvi是一个快速的多功能方法，整合了来自机器学习、统计模型、概率图模型等领域的思想和模型，应用到了单细胞数据的多种下游分析上来，相比其他针对单一下游分析的方法来说，从生物信息学的角度来讲是一个很大的进步。同期的Nature Methods期刊上也发表了一篇News来介绍这个工作，里面用一个瑞士军刀的比喻来说明scvi功能的全面。

scvi

单细胞测序数据将人们对生物特性的认知细化到了细胞层面，极大地促进了生物发现的进程。但是这种“高分辨率”的信息捕获技术也有着自身的局限性，比如测量灵敏度的不确定性，批次效应（batch effect，指即便是同一种细胞同一测序技术不同次的测序结果也存在着偏差），还有技术本身带来的噪声等等，进而影响了基于单细胞测序技术的下游分析的质量，比如细胞聚类（可以找到一组细胞中的常见子类，并发现新的子类或者亚群）、基因差异表达（找到不同细胞群之间表达量显著差异的基因）、轨迹分析（trajectory analysis，研究细胞发展变化的过程）等等。

scVI提出了基于深度神经网络的层次贝叶斯模型 (Hierarchical Bayesian Model) 来解决这些问题，其中贝叶斯模型中的条件概率由神经网络来实现。这是一个生成式模型，模型结构如下图所示，将scRNA-seq数据中每个细胞作为一个样本，其基因表达作为特征，通过encoder的神经网络与重参数化将高维的基因表示压缩到低维隐空间（比如说10维）；之后基于单细胞RNA测序数据基因表达量服从零膨胀负二项分布 (ZINB) 的假设，再利用decoder的神经网络将隐空间映射到基因表达分布参数的后验估计上。scvi的两大优势在于：

  • 模型考虑了scRNA-seq数据中的两大问题：library size 和 batch effect
  • 利用单一生成式模型为一系列的下游分析提供了数据基础
scvi-workflow

scvi-tools包含了一系列由Yosef团队开发的单细胞深度生成概率模型，scVI也被整合到了其中，本notebook将演示如何使用scvi-tools中的scVI来进行单细胞数据的批次效应去除。

代码
文本

前置准备

加载必要的库

代码
文本
[1]
import os
import tempfile
import scanpy as sc
import scvi
import seaborn as sns
import torch
代码
文本

设置全局参数

代码
文本
[2]
scvi.settings.seed = 0
print("Last run with scvi-tools version:", scvi.__version__)
Seed set to 0
Last run with scvi-tools version: 1.1.6
代码
文本

你可以更改下面的save_dir来决定本notebook产生的结果文件的存储路径。

代码
文本
[3]
sc.set_figure_params(figsize=(6, 6), frameon=False)
sns.set_theme()
torch.set_float32_matmul_precision("high")
save_dir = tempfile.TemporaryDirectory()

%config InlineBackend.print_figure_kwargs={"facecolor": "w"}
%config InlineBackend.figure_format="retina"
代码
文本

加载和准备数据

首先，让我们加载内置的PBMC数据集。scvi-tools 有许多“内置”数据集，并支持加载任意 .csv、.loom 和 .h5ad (AnnData) 文件。

代码
文本

所有的scvi-tools模块都需要用AnnData对象作为输入。

代码
文本
[5]
adata = scvi.data.purified_pbmc_dataset(save_path=save_dir.name)
adata
INFO     Downloading file at /tmp/tmpe7bl6sag/PurifiedPBMCDataset.h5ad                                             
Downloading...: 157054it [00:12, 12251.40it/s]                              

AnnData object with n_obs × n_vars = 105868 × 21932
,    obs: 'cell_types', 'barcodes', 'labels', 'batch'
代码
文本

现在，我们对数据进行预处理，以删除表达水平极低的基因和其他异常值。对于这些任务，我们使用Scanpy 预处理模块来进行。

代码
文本
[6]
sc.pp.filter_genes(adata, min_counts=3)
代码
文本

在scRNA-seq分析中，对数据进行标准化是很常见的做法。scvi-tools不使用这些值，但考虑到它们在其他任务以及可视化中的流行程度，我们将它们单独存储在anndata对象中（通过 .raw 属性）。

除非另有说明，否则scvi-tools模型需要原始计数（而不是标准化后的）。

代码
文本
[7]
adata.layers["counts"] = adata.X.copy() # preserve counts
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata # freeze the state in `.raw`
代码
文本

最后，我们进行特征选择，以减少用作scvi-tools模型输入的特征（在本例中为基因）的数量。对于scVI，建议使用1,000到10,000个HVG，视具体情况而定。

代码
文本
[11]
sc.pp.highly_variable_genes(
adata,
n_top_genes=1200,
subset=True,
layer="counts",
flavor="seurat_v3",
batch_key="batch",
)
代码
文本

现在让我们运行setup_anndata()，它会提醒scvi-tools各个矩阵在anndata中的位置。使用正确的参数运行此函数很重要，这样scvi-tools才会知道数据集有批次、注释等。例如，如果批次已在scvi-tools中注册，则后续模型将校正批次效应。

在此数据集中，我们只有一个分类协变量“batch”，使用batch_key参数指定，如果有多个协变量需要校正，则可以使用categorical_covariate_keys参数和continuous_covariate_keys来注册这些协变量（离散或连续）。

代码
文本
[18]
scvi.model.SCVI.setup_anndata(
adata,
layer="counts",
batch_key="batch"
)
代码
文本

如果在运行“setup_anndata”后修改了adata，请在创建模型实例之前再次运行“setup_anndata”。

代码
文本

创建并训练模型

代码
文本
[19]
model = scvi.model.SCVI(adata)
model
SCVI model with the following parameters: 
,n_hidden: 128, n_latent: 10, n_layers: 1, dropout_rate: 0.1, dispersion: gene, gene_likelihood: zinb, 
,latent_distribution: normal.
,Training status: Not Trained
,Model's adata is minified?: False
,
代码
文本


使用 GPU 时，scVI模型运行速度会更快。




代码
文本
[20]
model.train()
GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
/opt/mamba/lib/python3.10/site-packages/lightning/pytorch/trainer/connectors/data_connector.py:441: The 'train_dataloader' does not have many workers which may be a bottleneck. Consider increasing the value of the `num_workers` argument` to `num_workers=7` in the `DataLoader` to improve performance.
Epoch 76/76: 100%|██████████| 76/76 [09:26<00:00,  7.39s/it, v_num=1, train_loss_step=137, train_loss_epoch=128]`Trainer.fit` stopped: `max_epochs=76` reached.
Epoch 76/76: 100%|██████████| 76/76 [09:26<00:00,  7.45s/it, v_num=1, train_loss_step=137, train_loss_epoch=128]

代码
文本
保存和加载


保存包括保存模型神经网络权重，以及用于初始化模型的参数。


代码
文本
[21]
model_dir = os.path.join(save_dir.name, "scvi_model")
model.save(model_dir, overwrite=True)
代码
文本
[22]
model = scvi.model.SCVI.load(model_dir, adata=adata)
INFO     File /tmp/tmpe7bl6sag/scvi_model/model.pt already downloaded                                              
/opt/mamba/lib/python3.10/site-packages/scvi/model/base/_utils.py:66: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.
  model = torch.load(model_path, map_location=map_location)

代码
文本
获取模型输出


建议将scvi-tools的输出存储回原始anndata，这样就能用Scanpy进行进一步分析。


代码
文本
[23]
SCVI_LATENT_KEY = "X_scVI"
latent = model.get_latent_representation()
adata.obsm[SCVI_LATENT_KEY] = latent
latent.shape
(105868, 10)
代码
文本
model.get...()函数默认使用用于初始化模型的 AnnData。也可以查询anndata的子集，甚至可以使用完全独立的anndata对象，只要anndata有相同的数据结构即可。


代码
文本
[25]
adata_subset = adata[adata.obs.cell_types == "cd4_t_helper"]
latent_subset = model.get_latent_representation(adata_subset)
latent_subset.shape
INFO     Received view of anndata, making copy.                                                                    
INFO     Input AnnData not setup with scvi-tools. attempting to transfer AnnData setup                             

(22426, 10)
代码
文本
[26]
denoised = model.get_normalized_expression(adata_subset, library_size=1e4)
denoised.iloc[:5, :5]

,
,
,  
,    
,      
,      
,      
,      
,      
,      
,    
,  
,  
,    
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,    
,  
,
RP4-758J18.13VWA1RP5-892K4.1RP11-181G12.2PRDM16
00.0926750.0346990.0459960.0337640.010116
10.0919820.0174940.0255300.1318600.089071
20.0610150.0195500.0737920.0664600.010783
30.0104240.0023320.0168000.1435820.005972
40.0231810.0152230.0237270.0601230.011572

,
代码
文本
我们将标准化的值重新存储到anndata中。


代码
文本
[27]
SCVI_NORMALIZED_KEY = "scvi_normalized"
adata.layers[SCVI_NORMALIZED_KEY] = model.get_normalized_expression(library_size=10e4)
代码
文本
校正批次效应


Scanpy 是一个功能强大的Python库，用于对scRNA-seq数据进行可视化和下游分析。我们在此展示如何用scanpy可视化scVI去除数据批次效应前后的结果。


代码
文本
没有批次校正的可视化




我们使用UMAP定性评估我们的低维细胞embedding。使用scVI生成的embedding可以直接作为PCA获得的embedding的直接替代，如下所示。




首先，我们通过绘制原始计数数据前30个PCA主成分的UMAP结果来证明这里的批次是一个干扰变异（nuisance variation）。


代码
文本
[28]
# run PCA then generate UMAP plots
sc.tl.pca(adata)
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=20)
sc.tl.umap(adata, min_dist=0.3)
代码
文本
[33]
adata.obs.batch = adata.obs.batch.astype('category')
sc.pl.umap(
    adata,
    color=["cell_types", "batch"],
    ncols=2,
    frameon=False,
)
代码
文本
我们发现，虽然细胞类型通常分离得很好，但是干扰变异在数据的变化中起着很大的作用。


代码
文本
批次校正后的可视化 (scVI)


现在，让我们尝试使用scVI潜在空间生成相同的UMAP图，看看scVI是否成功缓解了数据中的批次效应。


代码
文本
[34]
# use scVI latent space for UMAP generation
sc.pp.neighbors(adata, use_rep=SCVI_LATENT_KEY)
sc.tl.umap(adata, min_dist=0.3)
代码
文本
[36]
sc.pl.umap(
    adata,
    color=["batch"],
    frameon=False,
)
代码
文本
我们可以看到，scVI成功校正了由于批次造成的干扰变化。


代码
文本
对scVI潜在空间聚类


scvi-tools与scanpy的接口使得使用scanpy从scVI的潜在空间对数据进行聚类，然后将其重新嵌入scVI（例如，用于差异表达）变得容易。


代码
文本
[37]
# neighbors were already computed using scVI
SCVI_CLUSTERS_KEY = "leiden_scVI"
sc.tl.leiden(adata, key_added=SCVI_CLUSTERS_KEY, resolution=0.5)
/tmp/ipykernel_7526/3580189653.py:3: FutureWarning: In the future, the default backend for leiden will be igraph instead of leidenalg.

 To achieve the future defaults please pass: flavor="igraph" and n_iterations=2.  directed must also be False to work with igraph's implementation.
  sc.tl.leiden(adata, key_added=SCVI_CLUSTERS_KEY, resolution=0.5)

代码
文本
[38]
sc.pl.umap(
    adata,
    color=[SCVI_CLUSTERS_KEY],
    frameon=False,
)
代码
文本
批次校正后的差异表达分析


例如，1-vs-1 DE测试非常简单：


代码
文本
[42]
de_df = model.differential_expression(
    groupby="cell_types", group1="cd4_t_helper", group2="naive_cytotoxic"
)
de_df.head()
DE...: 100%|██████████| 1/1 [00:01<00:00,  1.51s/it]


,
,
,  
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,  
,  
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,  
,
proba_deproba_not_debayes_factorscale1scale2pseudocountsdeltalfc_meanlfc_medianlfc_std...raw_mean1raw_mean2non_zeros_proportion1non_zeros_proportion2raw_normalized_mean1raw_normalized_mean2is_de_fdr_0.05comparisongroup1group2
RP5-892K4.10.99280.00724.9264470.0000040.0000330.00.25-3.179379-3.1232451.344204...0.0000890.0015060.0000890.0012550.0064160.079633Truecd4_t_helper vs naive_cytotoxiccd4_t_helpernaive_cytotoxic
CTD-2154B17.40.99220.00784.8458000.0000030.0000150.00.25-2.612882-2.6092930.988461...0.0000000.0003350.0000000.0002510.0000000.025562Truecd4_t_helper vs naive_cytotoxiccd4_t_helpernaive_cytotoxic
PCSK60.99160.00844.7710870.0000030.0000160.00.25-2.675164-2.6440901.080938...0.0000890.0015060.0000890.0014220.0074940.081850Truecd4_t_helper vs naive_cytotoxiccd4_t_helpernaive_cytotoxic
FAM171A20.99000.01004.5951190.0000020.0000140.00.25-2.820170-2.8166931.138890...0.0000000.0004180.0000000.0003350.0000000.023196Truecd4_t_helper vs naive_cytotoxiccd4_t_helpernaive_cytotoxic
S100B0.98880.01124.4805770.0001420.0033850.00.25-9.200255-9.32450514.401048...0.0099880.6471200.0072240.3207560.58408338.634869Truecd4_t_helper vs naive_cytotoxiccd4_t_helpernaive_cytotoxic

,
5 rows × 22 columns

,
代码
文本
我们还可以进行1-vs-all DE测试，将每种细胞类型与数据集的其余部分进行比较：


代码
文本
[43]
de_df = model.differential_expression(
    groupby="cell_types",
)
de_df.head()
DE...: 100%|██████████| 10/10 [00:38<00:00,  3.88s/it]


,
,
,  
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,  
,  
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,    
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,      
,    
,  
,
proba_deproba_not_debayes_factorscale1scale2pseudocountsdeltalfc_meanlfc_medianlfc_std...raw_mean1raw_mean2non_zeros_proportion1non_zeros_proportion2raw_normalized_mean1raw_normalized_mean2is_de_fdr_0.05comparisongroup1group2
AL928742.120.99980.00028.5169430.0002470.0000030.00.257.3336447.3381812.282522...0.0121960.0001570.0083290.0000210.6016310.002720Trueb_cells vs Restb_cellsRest
MS4A10.99960.00047.8236210.0069950.0000880.00.258.1293928.2136031.953751...1.3876550.0113280.6547350.00630668.8708570.402582Trueb_cells vs Restb_cellsRest
HLA-DOB0.99940.00067.4179640.0014910.0000250.00.257.1528907.2294681.737764...0.2832880.0039150.2306400.00279814.7465570.104858Trueb_cells vs Restb_cellsRest
CD79B0.99940.00067.4179640.0139110.0004380.00.255.3568305.3643570.985360...3.0693240.0918650.8981660.072591158.5021674.351209Trueb_cells vs Restb_cellsRest
RP11-16E12.20.99920.00087.1300860.0002150.0000080.00.256.0850506.0835021.904058...0.0386710.0005640.0358950.0004281.8700000.016863Trueb_cells vs Restb_cellsRest

,
5 rows × 22 columns

,
代码
文本
我们现在使用DE结果提取每个cluster的top marker：


代码
文本
[45]
markers = {}
cats = adata.obs.cell_types.cat.categories
for c in cats:
    cid = f"{c} vs Rest"
    cell_type_df = de_df.loc[de_df.comparison == cid]
    cell_type_df = cell_type_df[cell_type_df.lfc_mean > 0]
    cell_type_df = cell_type_df[cell_type_df["bayes_factor"] > 3]
    cell_type_df = cell_type_df[cell_type_df["non_zeros_proportion1"] > 0.1]
    markers[c] = cell_type_df.index.tolist()[:3]

sc.tl.dendrogram(adata, groupby="cell_types", use_rep="X_scVI")  
代码
文本
[47]
sc.pl.dotplot(
    adata,
    markers,
    groupby="cell_types",
    dendrogram=True,
    color_map="Blues",
    swap_axes=True,
    use_raw=True,
    standard_scale="var",
)
代码
文本
我们还可以使用“layer”选项可视化经过scVI标准化后的基因表达值。


代码
文本
[48]
sc.pl.heatmap(
    adata,
    markers,
    groupby="cell_types",
    layer="scvi_normalized",
    standard_scale="var",
    dendrogram=True,
    figsize=(8, 12),
)
代码
文本
python
notebook
深度学习
机器学习
生物信息学
单细胞分析
中文
pythonnotebook深度学习机器学习生物信息学单细胞分析中文
已赞1
推荐阅读
公开
多分辨率反卷积模型DestVI
pythonnotebookTutorial生物信息学单细胞分析空间转录组机器学习深度学习中文
pythonnotebookTutorial生物信息学单细胞分析空间转录组机器学习深度学习中文
Gongyq
更新于 2024-09-10
2 赞
公开
因果推断中的工具变量（DeepIV）
Deep LearningnotebookCausalityInstrumental Variable
Deep LearningnotebookCausalityInstrumental Variable
wenh@dp.tech
发布于 2023-09-26
1 转存文件