single-cell-clustering-and-batch-correction-with-omicverse

使用 omicverse 进行单细胞聚类和批次校正

概述

本技能提炼自单细胞教程 t_cluster.ipynb 和 t_single_batch.ipynb。当用户希望预处理 AnnData 对象、探索不同的聚类方法（Leiden、Louvain、scICE、GMM、主题/cNMF 模型），并使用 omicverse 工具评估或协调批次时，请使用本技能。

说明

1. 导入库并设置绘图默认值
   • 加载 omicverse as ov、scanpy as sc 和绘图辅助工具（使用齿状回演示数据时加载 scvelo as scv）。
   • 使用 ov.plot_set() 或 ov.utils.ov_plot_set()，以便在绘制嵌入图之前应用 omicverse 样式。
2. 加载数据并标注批次
   • 对于演示聚类，获取 scv.datasets.dentategyrus()；对于整合，使用 ov.read() 读取提供的 .h5ad 文件，并为每个队列设置 adata.obs['batch'] 标识。
   • 确认输入是稀疏数值矩阵；QC 步骤需要时，使用 adata.X = adata.X.astype(np.int64) 进行转换。
3. 运行质量控制
   • 执行 ov.pp.qc(adata, tresh={'mito_perc': 0.2, 'nUMIs': 500, 'detected_genes': 250}, batch_key='batch')，以剔除低质量细胞并检查各批次的汇总统计数据。
   • 保存中间过滤后的对象（adata.write_h5ad(...)），以便用户从干净的检查点继续操作。
4. 预处理并选择特征
   • 调用 ov.pp.preprocess(adata, mode='shiftlog|pearson', n_HVGs=3000, batch_key=None)，以进行归一化、对数转换并标记高变基因；设置 adata.raw = adata，并筛选 adata.var.highly_variable_features，用于后续建模。
   • 缩放表达量（ov.pp.scale(adata)），并使用 ov.pp.pca(adata, layer='scaled', n_pcs=50) 计算 PCA 得分。建议查看方差解释情况：ov.utils.plot_pca_variance_ratio(adata)。
5. 构建邻域图并进行基线聚类
   • 使用 sc.pp.neighbors(adata, n_neighbors=15, n_pcs=50, use_rep='scaled|original|X_pca') 或 ov.pp.neighbors(...) 构建邻域图。
   • 使用 ov.utils.cluster(adata, method='leiden'|'louvain', resolution=1)、ov.single.leiden(adata, resolution=1.0) 或 ov.pp.leiden(adata, resolution=1) 生成 Leiden 或 Louvain 标签；提醒用户，分辨率参数用于调整聚类粒度。
   • 重要 — 依赖项检查：聚类或绘图之前，务必检查前置条件：

     # Before clustering: check neighbors graph exists
     if 'neighbors' not in adata.uns:
         if 'X_pca' in adata.obsm:
             ov.pp.neighbors(adata, n_neighbors=15, use_rep='X_pca')
         else:
             raise ValueError("PCA must be computed before neighbors graph")
     # 按聚类绘图前：检查是否已执行聚类
     if 'leiden' not in adata.obs:
         ov.single.leiden(adata, resolution=1.0)
     

   • 使用 ov.pl.embedding(adata, basis='X_umap', color=['clusters','leiden'], frameon='small', wspace=0.5) 可视化嵌入，并确认簇之间已分离。绘图前，务必检查 color= 参数中的列是否存在于 adata.obs 中。
6. 探索高级聚类策略
   • scICE 共识聚类：实例化 model = ov.utils.cluster(adata, method='scICE', use_rep='scaled|original|X_pca', resolution_range=(4,20), n_boot=50, n_steps=11)，并通过 model.plot_ic(figsize=(6,4)) 检查稳定性，然后选择 model.best_k 个组。
   • 高斯混合模型：运行 ov.utils.cluster(..., method='GMM', n_components=21, covariance_type='full', tol=1e-9, max_iter=1000)，进行基于模型的分配。
   • 主题建模：拟合 LDA_obj = ov.utils.LDA_topic(...)，查看 LDA_obj.plot_topic_contributions(6)，通过 LDA_obj.predicted(k) 得出聚类结果，并可使用 LDA_obj.get_results_rfc(...) 进一步优化。
   • cNMF 程序：初始化 cnmf_obj = ov.single.cNMF(... components=np.arange(5,11), n_iter=20, num_highvar_genes=2000, output_dir=...)，执行因子分解（factorize、combine），通过 k_selection_plot 选择 K，并使用 cnmf_obj.get_results(...) 和 cnmf_obj.get_results_rfc(...) 将使用度分数回传。
7. 评估聚类质量
   • 使用 adjusted_rand_score(adata.obs['clusters'], adata.obs['leiden']) 将预测标签与已知参考进行比较，并报告每种方法（Leiden、Louvain、GMM、LDA 变体、cNMF 模型）的指标，以论证所选参数的合理性。
8. 使用多种布局进行嵌入
   • 使用 ov.utils.mde(...)，从不同的潜在空间（adata.obsm["scaled|original|X_pca"]、校正后的嵌入、主题组成）生成 MDE 投影，并通过 ov.pl.embedding(..., color=['batch','cell_type']) 或 ov.pl.embedding 绘图，以一致地检查簇/批次混合情况。
9. 执行批次校正与整合
   • 依次应用 ov.single.batch_correction(adata, batch_key='batch', methods='harmony'|'combat'|'scanorama'|'scVI'|'CellANOVA', n_pcs=50, ...)，生成存储在 adata.obsm 中的校正后嵌入（X_harmony、X_combat、X_scanorama、X_scVI、X_cellanova）。对于 scVI，注明潜在维度大小（n_latent=30）和 gene_likelihood="nb"；对于 CellANOVA，通过 control_dict 定义对照池。
   • 每次校正后，使用 ov.utils.mde 投影到 2D，并可视化 batch 与 cell_type，以检查混合情况和生物学特征保留情况。
10. 评测整合性能
    • 持久化最终对象（adata.write_h5ad('neurips2021_batch_all.h5ad', compression='gzip')），并在需要时重新加载。
    • 使用 scib_metrics.benchmark.Benchmarker 和嵌入列表（["X_pca", "X_combat", "X_harmony", "X_cellanova", "X_scanorama", "X_mira_topic", "X_mira_feature", "X_scVI"]），通过 bm.benchmark() 计算批次与生物学之间的权衡，并使用 bm.plot_results_table(min_max_scale=False) 汇总结果。
11. 常见故障排查
    • 确保在筛选 HVG 之前，adata.raw 保存的是未经缩放的对数归一化矩阵。
    • 聚类前确认 use_rep='scaled|original|X_pca' 字符串存在于 adata.obsm 中；若缺失，请重新运行预处理。
    • 运行 cNMF 或 scVI 时监控内存；对于较小的数据集，可调整 n_iter、components 或潜在维度。
    • 流水线依赖项错误：遇到诸如 "Could not find 'leiden' in adata.obs" 的错误时，务必检查并补齐前置条件：
      • Leiden/Louvain 聚类前 → 确保 'neighbors' in adata.uns
      • 按聚类绘图前 → 确保聚类列存在于 adata.obs 中
      • UMAP/嵌入前 → 确保已完成 PCA 或其他降维
    • 代码生成模式：生成多步骤代码时，应使用防御性检查，而不是假定先前步骤已成功完成。这样可以避免用户因步骤顺序不当或在不同会话中运行代码而导致错误级联。

示例

• “对齿状回细胞进行标准化，比较 Leiden、scICE 和 GMM 聚类，并报告与提供的 clusters 相比的 ARI 分数。”
• “使用 Harmony 和 scVI 对三个 NeurIPS 数据集进行批次校正，生成按 batch 和 cell_type 着色的 MDE 嵌入，并对这些嵌入进行基准测试。”
• “在预处理后的 AnnData 对象上拟合主题模型和 cNMF 模型，获取分类器优化后的聚类结果，并在 UMAP 上可视化生成的程序。”

参考资料

• 聚类演练：t_cluster.ipynb
• 批次整合演练：t_single_batch.ipynb
• 快速复制/粘贴命令：reference.md