scvi-tools 深度学习单细胞分析：功能详解

import scvi

# 基础用法
scvi.model.SCVI.setup_anndata(adata, layer='counts', batch_key='batch')
model = scvi.model.SCVI(adata)

# 关键参数
model = scvi.model.SCVI(
    adata,
    latent_dim=30,        # 潜在维度：20-50（默认 10）
    n_layers=2,           # 神经网络层数：1-3
    n_hidden=128,         # 隐藏层大小：64-256
    gene_likelihood='nb',  # 分布假设：'nb'（负二项）或 'zinb'
    dispersion='gene-batch',  # 离散度参数：'gene', 'gene-batch'
)

# 渐进式训练（快速收敛）
model.train(max_epochs=50)  # 初始训练
model.train(max_epochs=100, plan_kwargs={'lr': 1e-3})  # 降低学习率继续训练

# 早停策略
model.train(
    max_epochs=400,
    early_stopping=True,
    early_stopping_patience=15,  # 15 轮无改善则停止
    early_stopping_min_delta=0.001
)

# 准备数据：参考数据有标签，查询数据标签设为 'Unknown'
adata.obs['cell_type'] = adata.obs['cell_type'].astype(str)
adata.obs.loc[adata.obs['batch'] == 'query', 'cell_type'] = 'Unknown'

# 设置模型
scvi.model.SCANVI.setup_anndata(
    adata,
    layer='counts',
    batch_key='batch',
    labels_key='cell_type',  # 标签列
    unlabeled_category='Unknown'  # 未标记类别
)

# 从 scVI 模型初始化（加速收敛）
scvi_model = scvi.model.SCVI(adata)
scvi_model.train(max_epochs=100)
scanvi_model = scvi.model.SCANVI.from_scvi_model(
    scvi_model,
    unlabeled_category='Unknown'
)
scanvi_model.train(max_epochs=50)

# 获取预测标签和概率
adata.obs['predicted_cell_type'] = scanvi_model.predict()
adata.obs['prediction_confidence'] = scanvi_model.predict(soft=True).max(axis=1)

# 如果有部分真实标签，可计算准确率
from sklearn.metrics import accuracy_score, classification_report

mask = adata.obs['cell_type'] != 'Unknown'
true_labels = adata.obs.loc[mask, 'cell_type']
pred_labels = adata.obs.loc[mask, 'predicted_cell_type']

print(f"准确率：{accuracy_score(true_labels, pred_labels):.2%}")
print(classification_report(true_labels, pred_labels))

# AnnData 结构要求：
# adata.X -> RNA 计数
# adata.obsm['protein_counts'] -> 蛋白质计数 (n_cells × n_proteins)

import scvi

# 设置模型
scvi.model.TOTALVI.setup_anndata(
    adata,
    batch_key='batch',
    protein_expression_obsm_key='protein_counts'  # 蛋白数据位置
)

# 训练
model = scvi.model.TOTALVI(adata)
model.train(max_epochs=400)

# 获取去噪的蛋白质表达（纠正技术噪声）
protein_expression = model.get_normalized_expression(
    adata,
    transform='log',  # 对数变换
    get_protein_background=True,
    include_n_background_proteins=10
)

adata.obsm['protein_denoised'] = protein_expression[1]

# 比较原始与去噪后的蛋白质表达
import matplotlib.pyplot as plt

protein_idx = 0  # 选择某个蛋白
plt.scatter(
    adata.obsm['protein_counts'][:, protein_idx],
    adata.obsm['protein_denoised'][:, protein_idx],
    alpha=0.5
)
plt.xlabel('Raw counts')
plt.ylabel('Denoised expression')
plt.title('TotalVI denoising')
plt.show()

# AnnData 结构要求：
# adata.X -> RNA 计数
# adata.obsm['mode2'] -> ATAC 峰计数 (n_cells × n_peaks)

scvi.model.MULTIVI.setup_anndata(
    adata,
    batch_key='batch',
    modalities_key='modality'  # 标识每个细胞是哪种模态
)

# 训练
model = scvi.model.MULTIVI(
    adata,
    n_genes=adata.shape[1],
    n_regions=adata.obsm['mode2'].shape[1]
)
model.train(max_epochs=400)

# 获取联合潜在表示
adata.obsm["X_multivi"] = model.get_latent_representation()

# 分析模态特异性
modality_specificity = model.get_modality_specificity()
# 值接近 0：RNA 和 ATAC 一致
# 值接近 1：模态特异性

# 1. 准备参考数据（scRNA-seq）
scvi.model.SCVI.setup_anndata(adata_ref, layer='counts', batch_key='batch')
ref_model = scvi.model.SCVI(adata_ref)
ref_model.train()

# 2. 准备空间数据
scvi.model.DESTVI.setup_anndata(
    adata_spatial,
    layer='counts',
    batch_key='batch'
)

# 3. 创建 DestVI 模型
spatial_model = scvi.model.DESTVI.from_scvi_model(ref_model, adata_spatial)
spatial_model.train(max_epochs=200)

# 4. 获取细胞类型比例
cell_type_abundance = spatial_model.get_proportions()
# shape: (n_spots, n_cell_types)

adata_spatial.obs['cell_type_composition'] = [
    {ct: prop for ct, prop in zip(cell_types, row)}
    for row in cell_type_abundance
]

# 1. 保存预训练的参考模型
ref_model = scvi.model.SCVI(adata_ref)
ref_model.train()
ref_model.save('reference_model/', overwrite=True)

# 2. 加载新数据并映射
adata_query = sc.read_h5ad('new_data.h5ad')

# 使用 scArches 方式训练（冻结部分网络）
query_model = scvi.model.SCVI.load_query_data(
    'reference_model/',
    adata_query,
    freeze_dropout=True,  # 冻结 dropout
    freeze_expression=True,  # 冻结编码器部分层
)

query_model.train(max_epochs=50)  # 快速训练

# 3. 获取映射后的表示
adata_query.obsm['X_mapped'] = query_model.get_latent_representation()

# 簇间差异表达
de_results = model.differential_expression(
    groupby='leiden',
    group1='0',  # 目标簇
    group2='1',  # 对照簇
    mode='vanilla',  # 'vanilla': 简单比较，'change': 最小检测限
    batch_correction=True  # 批次校正
)

# 结果列：
# - lfc_mean: 对数折叠变化（后验均值）
# - lfc_median: 对数折叠变化（后验中位数）
# - prob_de: 差异表达概率
# - bayes_factor: 贝叶斯因子（>3 为强证据）

# 筛选显著基因
significant_genes = de_results[
    (de_results['bayes_factor'] > 3) &
    (abs(de_results['lfc_mean']) > 0.5)
]
print(significant_genes.head(10))

# 找到每个簇的标志基因
from scvi.tools import differential_expression

cluster_markers = {}
for cluster in adata.obs['leiden'].unique():
    de_results = model.differential_expression(
        groupby='leiden',
        group1=cluster,
        group2='rest'  # 与其他所有簇比较
    )
    cluster_markers[cluster] = de_results[
        de_results['bayes_factor'] > 3
    ].head(20)

# 检查 GPU 可用性
import torch
print(f"CUDA available: {torch.cuda.is_available()}")
print(f"GPU count: {torch.cuda.device_count()}")

# 指定 GPU 设备
model.train(use_gpu=True, gpus=[0])  # 使用第一块 GPU

# 多 GPU 训练（大数据集）
model.train(use_gpu=True, gpus=[0, 1], accelerator='ddp')

# 自定义回调
from pytorch_lightning.callbacks import ModelCheckpoint

checkpoint = ModelCheckpoint(
    dirpath='checkpoints/',
    monitor='elbo_validation',
    save_top_k=1,
    mode='min'
)

model.train(
    max_epochs=400,
    callbacks=[checkpoint],
    check_val_every_n_epoch=5  # 每 5 轮验证一次
)

# 重新加载最佳模型
model = scvi.model.SCVI.load('checkpoints/last.ckpt', adata)

# 数据流式加载（内存不足时）
model.train(
    max_epochs=400,
    data_loader_kwargs={
        'batch_size': 1024,  # 增大批大小
        'num_workers': 4  # 多进程加载
    }
)

# 分层训练（先训练子集）
adata_subset = adata[:10000, :]  # 训练 1 万细胞
model.train(max_epochs=200)
# 再用全数据微调
model.train(max_epochs=100)

# 获取不确定性估计
latent, uncertainty = model.get_latent_representation(
    adata,
    give_mean=False,  # 返回采样而非均值
    return_dist=True  # 返回均值和标准差
)

adata.obsm['X_scvi_uncertainty'] = uncertainty
sc.pl.umap(adata, color='X_scvi_uncertainty', cmap='Reds')

# 获取每个细胞类型的标志基因重要性
importance = model.get_feature_correlation(
    adata,
    groupby='leiden',
    group1='0'
)

# 可视化 top 基因
import seaborn as sns
top_genes = importance.nlargest(20, importance.columns[0])
sns.barplot(x=top_genes.iloc[:, 0], y=top_genes.index)
plt.title('Top contributing genes for cluster 0')
plt.show()

# ============ 模型设置 ============
scvi.model.SCVI.setup_anndata(adata, layer='counts', batch_key='batch')
scvi.model.SCANVI.setup_anndata(adata, layer='counts', batch_key='batch', labels_key='cell_type', unlabeled_category='Unknown')
scvi.model.TOTALVI.setup_anndata(adata, batch_key='batch', protein_expression_obsm_key='protein_counts')

# ============ 模型训练 ============
model.train(max_epochs=400, batch_size=128, use_gpu=True, early_stopping=True)

# ============ 获取表示 =-----------
adata.obsm["X_scvi"] = model.get_latent_representation()
denoised = model.get_normalized_expression(adata, transform='log')

# ============ 差异表达 ============
de_results = model.differential_expression(groupby='leiden', group1='0', group2='1')

# ============ 模型保存 ============
model.save('model_dir/', overwrite=True)
loaded_model = scvi.model.SCVI.load('model_dir/', adata)