scvi-tools 深度学习单细胞分析：应用场景

从真实科研场景中学习最佳实践

本章通过三个完整的科研案例，展示 scvi-tools 在实际研究中的应用。每个案例都来自已发表的高水平论文，我们简化了步骤但保留了核心思路。

案例 1：跨平台数据整合与细胞 atlas 构建

背景：某研究团队希望构建人类心脏发育的单细胞 atlas，数据来源包括：

5 个不同实验室（采用 10X v2、v3、Smart-seq2）
3 个发育阶段（胚胎期、胎儿期、成年期）
总计 50 万细胞

挑战：

批次效应极强（不同平台差异明显）
发育阶段差异与批次效应混杂
需要保留真实的生物学差异

解决方案：scVI + scANVI 联合流程

Step 1：分层数据整合

import scanpy as sc
import scvi

# 1.1 加载数据并标记来源
adata_list = []
for lab in ['lab1', 'lab2', 'lab3', 'lab4', 'lab5']:
    for tech in ['10xv2', '10xv3', 'smartseq2']:
        adata = sc.read_h5ad(f'{lab}_{tech}.h5ad')
        adata.obs['batch'] = f'{lab}_{tech}'  # 细粒度批次标签
        adata.obs['lab'] = lab  # 实验室标签
        adata.obs['tech'] = tech  # 技术标签
        adata_list.append(adata)

adata = adata_list[0].concatenate(adata_list, index_unique=None)

# 1.2 严格质量控制
sc.pp.filter_cells(adata, min_counts=1000)
sc.pp.filter_genes(adata, min_cells=20)

# 1.3 分层高变基因选择
# 使用 Seurat v3 方法，考虑批次信息
adata.layers['counts'] = adata.X.copy()
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=3000,
    batch_key='batch',
    flavor='seurat_v3',
    layer='counts',
    subset=True
)

Step 2：训练 scVI 模型

# 2.1 设置模型
scvi.model.SCVI.setup_anndata(
    adata,
    layer='counts',
    batch_key='batch'  # 使用细粒度批次标签
)

# 2.2 初始化模型
model = scvi.model.SCVI(
    adata,
    latent_dim=40,  # 较大的潜在维度以保留发育信息
    n_layers=3,  # 增加网络深度以学习复杂模式
)

# 2.3 训练
model.train(
    max_epochs=500,
    batch_size=256,
    early_stopping=True,
    early_stopping_patience=20,
    use_gpu=True
)

# 2.4 获取校正后的表示
adata.obsm['X_scvi'] = model.get_latent_representation()

Step 3：验证整合效果

# 3.1 聚类和可视化
sc.pp.neighbors(adata, use_rep='X_scvi')
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=1.0)

# 3.2 评估批次混合度
import scvi.tools

# 计算 LISI score（Local Inverse Simpson Index）
lisi = scvi.tools.compute_lisi(
    adata,
    label_keys=['batch', 'lab', 'tech'],
    embed_key='X_scvi'
)

print(f"批次 LISI: {lisi['batch'].mean():.2f}")  # >1.5 表示良好混合
print(f"实验室 LISI: {lisi['lab'].mean():.2f}")
print(f"技术 LISI: {lisi['tech'].mean():.2f}")

# 3.3 可视化
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
sc.pl.umap(adata, color='lab', ax=axes[0], show=False, title='By Lab')
sc.pl.umap(adata, color='tech', ax=axes[1], show=False, title='By Technology')
sc.pl.umap(adata, color='leiden', ax=axes[2], show=False, title='Clusters')
plt.tight_layout()
plt.show()

Step 4：使用 scANVI 进行细胞类型注释

# 4.1 利用已发表的心脏 atlas 作为参考
ref_adata = sc.read_h5ad('heart_atlas_reference.h5ad')

# 4.2 合并参考和查询数据
ref_adata.obs['cell_type'] = ref_adata.obs['cell_type'].astype(str)
adata.obs['cell_type'] = 'Unknown'

combined = ref_adata.concatenate(adata, index_unique='_')

# 4.3 训练 scANVI
scvi.model.SCANVI.setup_anndata(
    combined,
    layer='counts',
    batch_key='batch',
    labels_key='cell_type',
    unlabeled_category='Unknown'
)

scanvi_model = scvi.model.SCANVI(
    combined,
    unlabeled_category='Unknown',
    n_labels=len(ref_adata.obs['cell_type'].unique())
)

scanvi_model.train(max_epochs=200, use_gpu=True)

# 4.4 获取预测标签
adata.obs['predicted_cell_type'] = scanvi_model.predict()[
    adata.obs_names
]

# 4.5 评估置信度
pred_probs = scanvi_model.predict(soft=True)
adata.obs['prediction_confidence'] = pred_probs.max(axis=1)[
    adata.obs_names
]

# 过滤低置信度预测
high_conf_mask = adata.obs['prediction_confidence'] > 0.8
print(f"高置信度预测：{high_conf_mask.sum()}/{len(adata)} ({high_conf_mask.mean():.1%})")

结果与收获

整合质量：LISI score > 2.0，批次混合良好
生物学一致性：发育阶段在 UMAP 上形成连续轨迹
注释准确率：与专家手动注释一致率达 92%
效率：从手动注释 2 周缩短至自动注释 2 小时

案例 2：CITE-seq 数据揭示免疫治疗响应机制

背景：某肿瘤免疫研究团队分析了黑色素瘤患者治疗前后的 CITE-seq 数据（RNA + 27 个蛋白标志物），寻找预测免疫治疗响应的生物标志物。

挑战：

传统 RNA 分析无法捕捉表面蛋白变化
需要联合分析两种模态
蛋白数据噪声大，需要去噪

解决方案：totalVI 多模态分析

Step 1：准备 CITE-seq 数据

import scvi
import scanpy as sc

# 1.1 加载数据
adata = sc.read_h5ad('melanoma_citeseq.h5ad')

# 检查数据结构
print(f"RNA shape: {adata.shape}")
print(f"Protein shape: {adata.obsm['protein_counts'].shape}")

# 1.2 添加元数据
adata.obs['response'] = ['Responder' if x == 'R' else 'Non-responder'
                         for x in adata.obs['response']]
adata.obs['timepoint'] = ['Pre' if 'pre' in x else 'Post'
                          for x in adata.obs['sample_id']]

# 1.3 质量控制
sc.pp.filter_cells(adata, min_counts=1000)
sc.pp.filter_genes(adata, min_cells=20)

Step 2：训练 totalVI 模型

# 2.1 设置模型
scvi.model.TOTALVI.setup_anndata(
    adata,
    batch_key='batch',
    protein_expression_obsm_key='protein_counts'
)

# 2.2 初始化
model = scvi.model.TOTALVI(
    adata,
    latent_dim=30,
)

# 2.3 训练
model.train(max_epochs=400, use_gpu=True)

# 2.4 获取去噪的表达
rna_denoised, protein_denoised = model.get_normalized_expression(
    adata,
    transform='log',
    return_mean=True,
    return_protein_background=True
)

adata.layers['rna_denoised'] = rna_denoised
adata.obsm['protein_denoised'] = protein_denoised

Step 3：发现响应者特征

# 3.1 聚类
adata.obsm['X_totalvi'] = model.get_latent_representation()
sc.pp.neighbors(adata, use_rep='X_totalvi')
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

# 3.2 寻找差异表达特征
# 比较响应者与非响应者的 CD8+ T 细胞
cd8_t_cells = adata[adata.obs['leiden'].isin(['2', '5'])]

responder_mask = cd8_t_cells.obs['response'] == 'Responder'

# RNA 差异表达
de_rna = model.differential_expression(
    cd8_t_cells,
    groupby='response',
    group1='Responder',
    group2='Non-responder'
)

# 蛋白差异表达
de_protein = model.differential_expression(
    cd8_t_cells,
    groupby='response',
    group1='Responder',
    group2='Non-responder',
    modality='protein'  # 关键：指定分析蛋白
)

print("Top RNA markers:")
print(de_rna.head(10))

print("\nTop protein markers:")
print(de_protein.head(10))

Step 4：验证蛋白去噪效果

import matplotlib.pyplot as plt
import seaborn as sns

# 选择关键蛋白：PD-1 (PDCD1)
protein_idx = list(adata.uns['protein_names']).index('PDCD1')

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 原始 vs 去噪
axes[0].scatter(
    adata.obsm['protein_counts'][:, protein_idx],
    adata.obsm['protein_denoised'][:, protein_idx],
    alpha=0.5, s=10
)
axes[0].set_xlabel('Raw counts')
axes[0].set_ylabel('Denoised expression')
axes[0].set_title('PDCD1: Raw vs Denoised')

# 响应者 vs 非响应者
df = pd.DataFrame({
    'PDCD1': adata.obsm['protein_denoised'][:, protein_idx],
    'Response': adata.obs['response']
})

sns.boxplot(data=df, x='Response', y='PDCD1', ax=axes[1])
axes[1].set_title('PDCD1 expression by response')

plt.tight_layout()
plt.show()

结果与发现

新标志物：发现 CD8+ T 细胞上的 TIM-3 蛋白高表达与治疗响应相关
RNA-蛋白不一致：部分标志物在 RNA 水平无差异，但蛋白水平显著不同
预测模型：基于 5 个蛋白标志物构建的预测模型准确率达 85%

案例 3：空间转录组反卷解析析肿瘤微环境

背景：结直肠癌研究团队使用 Visium 空间转录组技术，希望解析肿瘤组织中不同区域的细胞类型组成。

挑战：

每个 spot 包含多个细胞（平均 5-10 个）
需要结合 scRNA-seq 参考进行反卷积
空间位置信息至关重要

解决方案：DestVI 空间反卷积

Step 1：准备参考和空间数据

import scvi
import scanpy as sc

# 1.1 加载 scRNA-seq 参考
ref_adata = sc.read_h5ad('colorectal_scrna.h5ad')

# 训练 scVI 模型
scvi.model.SCVI.setup_anndata(ref_adata, layer='counts', batch_key='batch')
ref_model = scvi.model.SCVI(ref_adata)
ref_model.train(max_epochs=400)

# 1.2 加载空间数据
spatial_adata = sc.read_visium('spatial_data/')
spatial_adata.layers['counts'] = spatial_adata.X.copy()

Step 2：训练 DestVI 模型

# 2.1 设置 DestVI
scvi.model.DESTVI.setup_anndata(
    spatial_adata,
    layer='counts',
    batch_key='batch'
)

# 2.2 从参考模型创建 DestVI
destvi_model = scvi.model.DESTVI.from_scvi_model(
    ref_model,
    spatial_adata
)

# 2.3 训练
destvi_model.train(max_epochs=200, use_gpu=True)

# 2.4 获取细胞类型比例
cell_type_proportions = destvi_model.get_proportions()

# 添加到空间数据
for i, ct in enumerate(ref_adata.obs['cell_type'].cat.categories):
    spatial_adata.obs[f'prop_{ct}'] = cell_type_proportions[:, i]

Step 3：可视化空间分布

import matplotlib.pyplot as plt

# 3.1 可观细胞类型分布
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
axes = axes.flatten()

cell_types = ['Tumor', 'T_cell', 'B_cell', 'Macrophage', 'Fibroblast', 'Endothelial']

for ax, ct in zip(axes, cell_types):
    sc.pl.spatial(
        spatial_adata,
        color=f'prop_{ct}',
        ax=ax,
        show=False,
        cmap='viridis',
        title=f'{ct} proportion'
    )

plt.tight_layout()
plt.show()

# 3.2 识别空间共存模式
import seaborn as sns

# 计算细胞类型相关性
prop_cols = [f'prop_{ct}' for ct in cell_types]
correlation = spatial_adata.obs[prop_cols].corr()

sns.clustermap(correlation, cmap='coolwarm', center=0)
plt.suptitle('Cell type co-occurrence patterns')
plt.show()

Step 4：识别肿瘤相关生态位

# 4.1 定义肿瘤富集区域
tumor_spots = spatial_adata.obs['prop_Tumor'] > 0.5

# 4.2 比较肿瘤与正常区域
tumor_area = spatial_adata[tumor_spots]
normal_area = spatial_adata[~tumor_spots]

# 计算差异
tumor_composition = tumor_area.obs[prop_cols].mean()
normal_composition = normal_area.obs[prop_cols].mean()

diff = tumor_composition - normal_composition
print("Enriched in tumor area:")
print(diff.sort_values(ascending=False).head(5))

print("\nDepleted in tumor area:")
print(diff.sort_values(ascending=True).head(5))

结果与发现

空间结构：发现肿瘤区域富集 M2 巨噬细胞和癌相关成纤维细胞
免疫排斥：T 细胞主要位于肿瘤边界，形成"免疫排斥"表型
治疗靶点：发现 CAF-T-mac 三细胞互作网络，可作为联合治疗靶点

案例总结：通用工作流程

通过以上三个案例，我们可以提炼出 scvi-tools 的通用工作流程：

1. 数据准备
   ├─ 加载数据（合并批次、添加元数据）
   ├─ 质量控制（过滤低质量细胞/基因）
   └─ 高变基因选择（考虑批次信息）

2. 模型选择
   ├─ scRNA-seq 整合 → scVI
   ├─ 有标签参考 → scANVI
   ├─ CITE-seq → totalVI
   ├─ 空间数据 → DestVI
   └─ 多模态 → MultiVI

3. 模型训练
   ├─ 设置 AnnData（指定 batch_key、labels_key）
   ├─ 初始化模型（选择合适的超参数）
   └─ 训练（启用 GPU、早停策略）

4. 下游分析
   ├─ 获取潜在表示
   ├─ 聚类和可视化
   ├─ 差异表达分析
   ├─ 标签转移
   └─ 特定模型分析（反卷积、去噪等）

5. 结果验证
   ├─ 评估批次混合度（LISI）
   ├─ 生物学一致性检查
   └─ 与独立数据集比较

最佳实践建议

1. 数据整合

✅ 使用细粒度的批次标签（实验室_技术_日期）
✅ 选择合适的高变基因数量（2000-4000）
✅ 验证整合效果（LISI score、可视化）

2. 模型训练

✅ 使用 GPU 加速（训练速度提升 10-50x）
✅ 启用早停策略（避免过拟合）
✅ 监控训练曲线（检查收敛性）

3. 结果解释

✅ 结合生物学知识验证聚类结果
✅ 使用多种指标评估整合质量
✅ 保存中间结果以便重现

4. 性能优化

✅ 大数据集：分批训练、流式加载
✅ 小数据集：降低 latent_dim
✅ 多模态数据：确保各模态质量

通过本章的案例学习，你已经掌握了 scvi-tools 在实际研究中的应用。在第五章中，我们将深入探讨这些模型背后的数学原理，帮助你理解为什么它们如此有效。

→ 继续阅读：第五章 - 工作原理