scVI 深度学习单细胞 — 变分自编码器驱动的单细胞分析框架¶

# 创建 conda 环境（需要 Python 3.11+）
conda create -n scvi_env python=3.11 -y
conda activate scvi_env

# 安装 scVI-tools（最新版 1.4.2）
pip install scvi-tools==1.4.2

# 安装 PyTorch（CPU 版本）
pip install torch --index-url https://download.pytorch.org/whl/cpu

# 若有 GPU（CUDA 12.x），安装 GPU 版 PyTorch
# pip install torch --index-url https://download.pytorch.org/whl/cu121

# 安装 Scanpy（数据处理依赖）
pip install scanpy

# 验证安装是否成功
python -c "import scvi; print(scvi.__version__)"

import scvi                            # 导入 scVI-tools 主包
import scanpy as sc                    # Scanpy 配合使用
import anndata as ad                   # AnnData 数据格式

# ── 读取数据 ──────────────────────────────────────────
adata = sc.read_h5ad('preprocessed.h5ad')

# 必须在 AnnData 中保存原始整数计数（scVI 需要原始计数作为输入）
# 如果 adata.X 已经是标准化后的数据，需要从 raw 或 layers 恢复
adata.X = adata.layers['counts'].copy()   # 使用原始计数层

# ── 设置 SCVI 模型 ────────────────────────────────────
# 注册 AnnData 数据，告诉 scVI 哪一列是批次信息
scvi.model.SCVI.setup_anndata(
    adata,
    layer='counts',                    # 原始计数层名称
    batch_key='batch',                 # 批次信息列名（用于批次校正）
    categorical_covariate_keys=['donor']  # 其他分类协变量（可选）
)

# 创建 SCVI 模型（VAE 结构）
model = scvi.model.SCVI(
    adata,
    n_layers=2,                        # 神经网络层数（默认2层够用）
    n_latent=30,                       # 潜在空间维度（类似 PCA 的 PC 数）
    gene_likelihood='nb'               # 负二项分布，适合 scRNA-seq 计数数据
)

# 训练模型（CPU 约需 10-30 分钟，GPU 更快）
model.train(
    max_epochs=400,                    # 最大训练轮数
    batch_size=128,                    # 每批次样本数
    early_stopping=True                # 自动提前停止，避免过拟合
)

import scvi
import scanpy as sc
import matplotlib.pyplot as plt

# ── 完整 scVI 分析流程 ────────────────────────────────

# 1. 读取多批次数据（假设有两个样本需要整合）
adata = sc.read_h5ad('multi_batch.h5ad')

# 2. 基础预处理（保留原始计数）
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata.layers['counts'] = adata.X.copy()   # 保存原始计数
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(
    adata, n_top_genes=3000,
    subset=True,                           # 只保留高变基因
    layer='counts',
    flavor='seurat_v3'                     # 基于 counts 层计算高变基因
)

# 3. 设置并训练 SCVI 模型
scvi.model.SCVI.setup_anndata(adata, layer='counts', batch_key='sample')
vae = scvi.model.SCVI(adata, n_latent=20, n_layers=2)
vae.train(max_epochs=300, early_stopping=True)

# 4. 获取批次校正后的细胞嵌入（替代 PCA 用于后续分析）
adata.obsm['X_scVI'] = vae.get_latent_representation()  # 潜在空间坐标

# 5. 用 scVI 嵌入做聚类
sc.pp.neighbors(adata, use_rep='X_scVI', n_neighbors=15)  # 用 scVI 嵌入建 KNN 图
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)
sc.pl.umap(adata, color=['leiden', 'sample'])              # 检查批次效应是否去除

# 6. scVI 差异表达分析（比 Wilcoxon 更准确，考虑了批次）
de_df = vae.differential_expression(
    adata,
    groupby='leiden',                      # 按聚类比较
    group1='0',                            # 第一组
    group2='1'                             # 第二组
)
print(de_df.head(20))                      # 查看前20个差异基因

# 7. 基因归因/插补（填充 dropout 造成的零值）
normalized_expr = vae.get_normalized_expression(
    adata, library_size=1e4              # 归一化到每细胞 10000 counts
)
adata.layers['scVI_normalized'] = normalized_expr.values

# 8. 保存模型和结果（方便后续加载）
vae.save('scvi_model_dir/', overwrite=True)
adata.write_h5ad('scvi_result.h5ad')

# 后续加载模型
# loaded_model = scvi.model.SCVI.load('scvi_model_dir/', adata=adata)

# 安装
pip install scvi-tools==1.4.2 torch scanpy

# 快速流程（3步走）
scvi.model.SCVI.setup_anndata(adata, layer='counts', batch_key='batch')
vae = scvi.model.SCVI(adata, n_latent=20)
vae.train()
adata.obsm['X_scVI'] = vae.get_latent_representation()

# 其他模型
# SCANVI  - 半监督，支持细胞类型标签迁移
# totalVI - 蛋白+RNA 多组学整合（CITE-seq）
# scANVI  - 跨数据集标签迁移