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() を取得します。統合の場合は、提供された .h5ad ファイルを ov.read() で読み込み、各コホートの 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 のラベルを生成します。解像度によってクラスタの粒度が調整されることをユーザーに伝えます。
   • 重要 - 依存関係の確認：クラスタリングやプロットの前に、必ず前提条件を確認します。

     # クラスタリングの前に、近傍グラフが存在することを確認
     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 スコアを報告してください。」
• 「3 つの NeurIPS データセットを Harmony と scVI でバッチ補正し、batch と cell_type で色分けした MDE 埋め込みを作成して、埋め込みをベンチマークしてください。」
• 「前処理済みの AnnData オブジェクトでトピックモデルと cNMF モデルを学習し、分類器で洗練されたクラスター割り当てを取得して、得られたプログラムを UMAP 上に可視化してください。」

参考資料

• クラスタリングの手順：t_cluster.ipynb
• バッチ統合の手順：t_single_batch.ipynb
• コピーしてすぐ使えるコマンド：reference.md