single-cell-preprocessing-with-omicverse

使用 omicverse 进行单细胞预处理

概述

当用户需要复现 omicverse 笔记本 t_preprocess.ipynb、t_preprocess_cpu.ipynb 和 t_preprocess_gpu.ipynb 中的预处理工作流时，请遵循本技能。这些教程使用 10x PBMC3k 数据集，涵盖 QC 过滤、归一化、高变基因（HVG）检测、降维和下游嵌入。

说明

1. 设置环境
   • 导入 omicverse as ov 和 scanpy as sc，然后调用 ov.plot_set(font_path='Arial')（或在旧版笔记本中调用 ov.ov_plot_set()），以统一图形样式。
   • 在笔记本中迭代时，建议使用 %load_ext autoreload 和 %autoreload 2，这样代码修改会生效，无需重启内核。
2. 准备输入数据
   • 从 10x Genomics 下载 PBMC3k 过滤后的矩阵（pbmc3k_filtered_gene_bc_matrices.tar.gz），并将其解压到 data/filtered_gene_bc_matrices/hg19/。
   • 使用 ov.io.read_10x_mtx(..., var_names='gene_symbols') 加载矩阵，并保留一个可写目录（如 write/）用于导出。
3. 执行质量控制（QC）
   • 在 CPU/CPU–GPU 流水线中运行 ov.pp.qc(adata, tresh={'mito_perc': 0.2, 'nUMIs': 500, 'detected_genes': 250}, doublets_method='scrublet')；在纯 GPU 流水线中省略 doublets_method，因为目前尚不支持 Scrublet。
   • 检查返回的 AnnData 摘要，确认双细胞率和 QC 阈值；并建议针对不同物种或测序深度调整截断值。
4. 转换前存储原始计数
   • 在 QC 后立即调用 ov.utils.store_layers(adata, layers='counts')，以便后续恢复和比较时仍可访问原始计数。
5. 标准化并选择 HVG
   • 使用 ov.pp.preprocess(adata, mode='shiftlog|pearson', n_HVGs=2000, target_sum=5e5) 执行 shift-log 标准化，然后通过 Pearson 残差检测 HVG（在 GPU 上设置 target_sum=None，以保留默认值）。
   • 对于 CPU–GPU 混合运行，演示如何使用 ov.pp.recover_counts(...) 反转标准化，并将重建的计数存储在 adata.layers['recover_counts'] 中。
6. 管理 .raw 和图层恢复
   • 使用 adata.raw = adata（或 adata.raw = adata.copy()）将标准化后的数据快照保存到 .raw，并展示如何使用 ov.utils.retrieve_layers(adata_counts, layers='counts') 比较标准化强度与原始强度。
7. 缩放、降维和嵌入
   • 使用 ov.pp.scale(adata) 缩放特征（缩放后的矩阵存放在图层中），然后调用 ov.pp.pca(adata, layer='scaled', n_pcs=50)。
   • 构建邻域图：
     • 基线笔记本使用 sc.pp.neighbors(adata, n_neighbors=15, n_pcs=50, use_rep='scaled|original|X_pca')。
     • CPU–GPU 流水线使用 ov.pp.neighbors(..., use_rep='scaled|original|X_pca')，以利用加速例程。
     • GPU 流水线使用 ov.pp.neighbors(..., method='cagra') 调用 RAPIDS 图原语。
   • 根据笔记本变体，通过 ov.utils.mde(...)、ov.pp.umap(adata)、ov.pp.mde(...)、ov.pp.tsne(...) 或 ov.pp.sude(...) 生成嵌入。
8. 聚类和注释
   • 构建邻域图后运行 ov.pp.leiden(adata, resolution=1) 或 ov.single.leiden(adata, resolution=1.0)；CPU–GPU 流水线还会展示如何在聚类前使用 ov.pp.score_genes_cell_cycle。
   • 重要 - 防御性检查：生成按聚类结果绘图的代码（例如 color='leiden'）时，务必先检查是否已执行聚类：

     # Check if leiden clustering exists, if not, run it
     if 'leiden' not in adata.obs:
         if 'neighbors' not in adata.uns:
             ov.pp.neighbors(adata, n_neighbors=15, use_rep='X_pca')
         ov.single.leiden(adata, resolution=1.0)
     

   • 使用 ov.pl.embedding(...) 或 ov.pl.umap(...) 绘制嵌入，按 leiden 聚类和标记基因着色。绘图前务必确认 color= 参数指定的列存在于 adata.obs 中。
9. 记录输出
   • 建议保存中间 AnnData 对象（adata.write('write/pbmc3k_preprocessed.h5ad')），并使用 Matplotlib 的 plt.savefig(...) 导出图形，以保留 QC 摘要和嵌入结果。
10. 笔记本专属说明
    • 基线（t_preprocess.ipynb）：专注于使用 Scanpy 邻域图的 CPU 执行；强调在 retrieve_layers 演示前后存储计数。
    • CPU–GPU 混合（t_preprocess_cpu.ipynb）：展示 Omicverse ≥1.7.0 的混合加速。加入计时魔法命令（%%time）以展示加速效果，并指出支持 doublets_method='scrublet'。
    • GPU（t_preprocess_gpu.ipynb）：需要支持 CUDA 的 GPU、RAPIDS 24.04 技术栈和 rapids-singlecell。说明 ov.pp.anndata_to_GPU/ov.pp.anndata_to_CPU 数据传输以及 method='cagra' 邻域图。指出当前的警告：纯 GPU 流水线依赖 RAPIDS 更新。
11. 故障排除提示
    • 如果 ov.io.read_10x_mtx 失败，请检查解压后的文件夹结构，并确保通过 var_names='gene_symbols' 获取基因符号。
    • 若 GPU 导入出错，请确认 conda 环境与已安装 CUDA 驱动对应的 RAPIDS 版本一致（nvidia-smi）。
    • 如果 ov.pp.preprocess 出现维度不匹配，请确保 QC 已过滤掉空条形码，以免 HVG 选择遇到零方差特征。
    • 当嵌入缺少预期字段（例如，缺少 scaled|original|X_pca）时，请重新运行 ov.pp.scale 和 ov.pp.pca，以重建缓存图层。
    • 流水线依赖错误：遇到类似 "Could not find 'leiden' in adata.obs or adata.var_names" 的错误时：
      • 始终检查依赖操作所需的预处理步骤（邻域图、PCA）是否已完成
      • 尝试按聚类结果为图着色前，检查 adata.obs 中是否存在聚类结果
      • 在生成的代码中使用防御性检查，以妥善处理不完整的流水线
    • 代码生成最佳实践：生成包含先决条件条件检查的健壮代码，而不要假设步骤总是按顺序完整执行。用户可能会在不同会话中运行各个步骤，或跳过中间步骤。

关键 API 参考——批次列处理

批次列验证——批次操作前必做

重要：在进行任何支持批次的操作（批次校正、整合等）之前，务必验证并准备好批次列。缺失值或 NaN 值会导致错误。

正确用法：

# Step 1: Check if batch column exists, create default if not
if 'batch' not in adata.obs.columns:
    adata.obs['batch'] = 'batch_1'  # Default single batch

# Step 2: Handle NaN/missing values - CRITICAL!
adata.obs['batch'] = adata.obs['batch'].fillna('unknown')

# Step 3: Convert to categorical for efficient memory usage
adata.obs['batch'] = adata.obs['batch'].astype('category')

# Now safe to use in batch-aware operations
ov.pp.combat(adata, batch='batch')  # or other batch correction methods

错误——请勿使用：

# WRONG! Using batch column without validation can cause NaN errors
# ov.pp.combat(adata, batch='batch')  # May fail if batch has NaN values!

# WRONG! Assuming batch column exists
# adata.obs['batch'].unique()  # KeyError if column doesn't exist!

常见的批次相关问题

1. 批次列中有 NaN 值：批次操作前务必使用 fillna()
2. 缺少批次列：使用前务必检查该列是否存在
3. 批次列不是分类类型：转换为 category 以提高内存效率
4. 数据类型混杂：分类前确保字符串类型一致

# Complete defensive batch preparation pattern:
def prepare_batch_column(adata, batch_key='batch', default_batch='batch_1'):
    """Prepare batch column for batch-aware operations."""
    if batch_key not in adata.obs.columns:
        adata.obs[batch_key] = default_batch
    adata.obs[batch_key] = adata.obs[batch_key].fillna('unknown')
    adata.obs[batch_key] = adata.obs[batch_key].astype(str).astype('category')
    return adata

高变基因（HVG）——小型数据集处理

小批次中的 LOESS 故障

重要：seurat_v3 HVG 风格使用 LOESS 回归，在小型数据集或每个批次的子集较小（每批次少于 500 个细胞）时会失败。表现如下：

ValueError: Extrapolation not allowed with blending

正确做法——使用 try/except 回退模式：

# Robust HVG selection for any dataset size
try:
    sc.pp.highly_variable_genes(
        adata,
        flavor='seurat_v3',
        n_top_genes=2000,
        batch_key='batch'  # if batch correction is needed
    )
except ValueError as e:
    if 'Extrapolation' in str(e) or 'LOESS' in str(e):
        # Fallback to simpler method for small datasets
        sc.pp.highly_variable_genes(
            adata,
            flavor='seurat',  # Works with any size
            n_top_genes=2000
        )
    else:
        raise

替代方案——使用 cell_ranger 风格进行批次感知 HVG：

# cell_ranger flavor is more robust for batched data
sc.pp.highly_variable_genes(
    adata,
    flavor='cell_ranger',  # No LOESS, works with batches
    n_top_genes=2000,
    batch_key='batch'
)

批次感知 HVG 的最佳实践

1. 计算 HVG 前检查批次大小：小批次（少于 500 个细胞）会导致 LOESS 失败
2. 批次大小差异较大时，优先选择 seurat 或 cell_ranger
3. 仅当所有批次都超过 500 个细胞时才使用 seurat_v3
4. 数据集大小未知时，务必使用 try/except 包装

# Safe batch-aware HVG pattern
def safe_highly_variable_genes(adata, batch_key='batch', n_top_genes=2000):
    """Select HVGs with automatic fallback for small batches."""
    try:
        sc.pp.highly_variable_genes(
            adata, flavor='seurat_v3', n_top_genes=n_top_genes, batch_key=batch_key
        )
    except ValueError:
        # Fallback for small batches
        sc.pp.highly_variable_genes(
            adata, flavor='seurat', n_top_genes=n_top_genes
        )

示例

• “下载 PBMC3k 计数数据，使用 Scrublet 运行 QC，使用 shiftlog|pearson 进行标准化，并在 CPU 上计算 MDE + UMAP 嵌入。”
• “在全新的 conda 环境中设置混合 CPU–GPU 工作流，在标准化后恢复原始计数，并在 Leiden 聚类前对细胞周期阶段进行评分。”
• “配置 RAPIDS 环境，将 AnnData 传输到 GPU，运行 method='cagra' 邻居搜索，并将嵌入结果返回 CPU 进行绘图。”

参考资料

• 详细演练笔记本：t_preprocess.ipynb、t_preprocess_cpu.ipynb、t_preprocess_gpu.ipynb
• 可快速复制粘贴的命令：reference.md