scRNA-seq解析：サンプル統合とバッチ補正 ― scVI（GPU必須のAIを使った手法）

この記事のゴール

複数サンプルをまたいで、バッチ効果（測定の技術的なズレ）だけを取り除き、生物学的な違いは保ったままのデータを、深層学習手法 scVI で作ること。強いバッチ効果や大規模データに強く、後続の発展解析（SOLO・アノテーション・マルチモーダル）に繋げられます。

0. バッチ効果とは／なぜ補正するのか

同じ細胞型でも、測定した日・試薬のロット・10x のチップなどが違うと、発現のパターンが系統的にずれます。これがバッチ効果です。補正しないと、クラスタリングで「細胞型」ではなく「サンプル」でグループが割れてしまい、解釈が壊れます。

1. 準備：正規化・特徴選択（HVG）・生カウントの保存

scVI は生（なま）のカウントを入力にします。正規化で X が上書きされる前に、生カウントを別の場所（レイヤー）へ退避しておきます。発現の変動が大きい遺伝子に絞る高変動遺伝子（HVG、highly variable genes）選択も、seurat_v3 で生カウントから行います。

pythonコピーimport numpy as np
import scanpy as sc
import anndata as ad

# combined: 前の記事（QC）で得たクリーンな AnnData（obs["batch"] にサンプル/群ラベル）
adata = combined.copy()

# ① 生カウントを保存（scVI が生カウントを必要とするため）
adata.layers["counts"] = adata.X.copy()

# ② HVG を「生カウント」から seurat_v3 で選択（バッチ考慮）
sc.pp.highly_variable_genes(
    adata, n_top_genes=2000, flavor="seurat_v3",
    layer="counts", batch_key="batch", subset=False,
)

2. scVI で統合する

仕組み

scVI は、各細胞の発現を低次元の潜在表現（数十次元のベクトル）に符号化する深層学習モデルで、中身は変分オートエンコーダ（VAE、variational autoencoder）です。動作は次の3点に整理できます。

• エンコーダが、生カウントを潜在表現 z に圧縮する
• デコーダが、z とバッチ情報から元の発現を再構成する
• カウントは負の二項分布でモデル化し、PCA を経由せず生データから直接学習する

要点は、バッチを「既知の共変量」としてモデルに与えることです。再構成のときにバッチ情報を別途渡すため、潜在表現 z 自体はバッチに依存しない方向へ学習されます。結果として、バッチ効果を分離した潜在表現 X_scVI が得られ、これを近傍グラフ・UMAP・クラスタリングに使います（Lopez et al., 2018）。

コード

まず HVG に絞り、生カウントとバッチを登録してから、モデルを学習させます。

pythonコピーimport scvi

# HVG に絞って学習するのが推奨（バッチ固有の発現変動を減らせる）
adata_scvi = adata[:, adata.var.highly_variable].copy()

# 生カウント（counts レイヤー）とバッチを登録
scvi.model.SCVI.setup_anndata(adata_scvi, layer="counts", batch_key="batch")

# モデル定義と学習（GPU を自動利用。大規模データほど epochs が自動調整される）
model = scvi.model.SCVI(adata_scvi, n_layers=2, n_latent=30, gene_likelihood="nb")
model.train()

# 潜在表現（バッチ補正済み）と denoise した発現量を取り出す
adata.obsm["X_scVI"] = model.get_latent_representation()
adata.layers["scvi_normalized"] = model.get_normalized_expression(library_size=1e4)

# scVI 潜在表現で近傍グラフ → UMAP → クラスタリング
sc.pp.neighbors(adata, use_rep="X_scVI", n_neighbors=15)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="leiden_scVI",
             flavor="igraph", n_iterations=2, directed=False, resolution=1.0)

sc.pl.umap(adata, color=["batch", "leiden_scVI"], ncols=2, frameon=False)

3. GPU 環境について

scVI は深層学習モデルなので、実用上 GPU が前提です。手元に GPU が無い場合は、Google Colab（無料の GPU ランタイム）やクラウド GPU でこのステップだけ実行し、結果（X_scVI など）を保存して持ち帰るのが現実的です。CPU でも動きますが、現実的なデータサイズでは学習に極端に時間がかかります。

4. 統合できたかを評価する（消しすぎていないか）

統合の評価では、「サンプルが混ざったか（バッチ効果が消えたか）」と「細胞型が保たれたか（消しすぎていないか）」の両面を確認します。

pythonコピー# ① バッチが混ざっているか
sc.pl.umap(adata, color="batch", frameon=False)
# ② 既知マーカーで細胞型が保たれているか
known_markers = ["GeneA", "GeneB", "GeneC"]   # 自分の系のマーカーに置換
sc.pl.umap(adata, color=known_markers, color_map="RdBu_r", frameon=False)

5. もう一方の選択肢：CPU で完結する Harmony

GPU が無い、または小〜中規模で素早く回したい場合は、古典手法の Harmony（PCA 空間で反復的にバッチ効果をそろえる）が軽量で実用的です（📖 詳しくは サンプル統合：Harmony 記事）。

まとめ

• scVI は生カウントを VAE で学習し、潜在表現を作る過程でバッチ効果を分離する。強いバッチ・大規模に強く、GPU 必須。
• 流れ：生カウントを counts に退避 → HVG（seurat_v3）→ setup_anndata(layer="counts") → train() → X_scVI で近傍グラフ。
• 確率的なので seed 固定・十分な細胞数が前提。GPU が無いなら Harmony（別記事） を。

次の記事：ダブレットの検出 ― SOLO（GPU必須） — 本記事で学習した scVI モデルを土台に、2細胞が1つに化けた人工物を検出・除去します。

関連記事

• 📖 scRNA-seqの前処理 入門 ─ 全体像と進め方 … このシリーズの全体像（ハブ）
• 📖 データの読み込みと品質管理（QC） … 本記事の前段（入力データ）
• 📖 サンプル統合：Harmony（古典・CPU） … 統合のもう一方の選択肢
• 📖 ダブレット検出：SOLO（AI・GPU） … 次の記事（scVI を土台にする）
• 📖 正規化・特徴選択・次元削減・クラスタリング … 統合後に続く段階
• 📖 AnnDataのデータ構造 ／ 📖 scvi-tools の全体像 … 土台・AI手法の地図

参考文献

• Lopez, R., Regier, J., Cole, M. B., Jordan, M. I., & Yosef, N. (2018). Deep generative modeling for single-cell transcriptomics. Nature Methods, 15, 1053–1058. doi:10.1038/s41592-018-0229-2
• Luecken, M. D., Büttner, M., Chaichoompu, K., et al. (2022). Benchmarking atlas-level data integration in single-cell genomics. Nature Methods, 19, 41–50. doi:10.1038/s41592-021-01336-8
• Heumos, L., Schaar, A. C., Lance, C., et al. (2023). Best practices for single-cell analysis across modalities. Nature Reviews Genetics, 24, 550–572. doi:10.1038/s41576-023-00586-w
• scvi-tools ドキュメント: https://docs.scvi-tools.org