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') を実行します。Scrublet がまだサポートされていない GPU のみの環境では、doublets_method を省略します。
   • 返された 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 環境では、RAPIDS のグラフプリミティブを呼び出すために ov.pp.neighbors(..., method='cagra') を使用します。
   • ノートブックの種類に応じて、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 を実行する手順も紹介します。
   • IMPORTANT - 事前チェック：クラスタリング結果（例：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')）と、QC の概要や埋め込みを保持するための Matplotlib による図のエクスポート（plt.savefig(...)）を保存するよう推奨します。
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 のインポートエラーが発生した場合は、インストールされている CUDA ドライバーに対応する RAPIDS バージョンと conda 環境が一致していることを確認します（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」のようなエラーが発生した場合：
      • 依存する処理の前に、必要な前処理（neighbors、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. バッチ列がカテゴリ型でない：メモリ効率を高めるためにカテゴリ型に変換してください
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

別の方法 — バッチ対応 HVG には cell_ranger フレーバーを使用する：

# 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