scRNA-seq解析：サンプル統合とバッチ補正 ― Harmony（GPU不要の古典的手法）

📚 この記事について：「scRNA-seq解析 実践シリーズ（前処理編）」の記事です。全体像と手法比較は「前処理の全体像とツール選択」を参照。本記事は統合（バッチ補正）の古典的手法 Harmony を扱います。コードは初めてでも追えるよう、各ステップで「何のために」を先に説明します。

🖥 動作環境：CPU で完結（GPU 不要）。

🔀 対になる手法：AI 版（GPU 必須）の「サンプル統合：scVI」。統合はどちらか一方を選びます。

🔙 前の記事：「データの読み込みとQC」　／　🔜 次の記事：「ダブレットの検出：scrublet」

前提：QC で得たクリーンな AnnData（annotated data、scanpy のデータ形式）／Python の基本。

この記事のゴール

複数サンプルをまたいで、「測定のクセ（バッチ効果）」だけを消し、本物の生物学的な違いは残したデータを、軽量・普通のPCで動く Harmony で作ること。

📷 たとえると：違うカメラで撮った写真の「色味」をそろえる作業です。同じ被写体（同じ細胞型）が、カメラ（サンプル）ごとに違う色に写ってしまう。そのカメラ差だけを補正し、被写体そのもの（本物の細胞型の違い）は変えません。

0. バッチ効果ってなに？ なぜ補正するの？

同じ細胞型でも、測定した日・試薬のロット・10x のチップなどが違うと、発現のパターンが系統的にずれます。これがバッチ効果です。

補正しないとどうなるか？ 次元削減やクラスタリングで、「細胞型」ではなく「サンプル」でグループが割れてしまい、解釈が壊れます。

💡 統合のゴールは「混ぜる」と「残す」の両立 良い統合とは、次の2つを同時に満たすことです。

🔵 同じ細胞型は、サンプルをまたいで混ざる（バッチ補正）

🔴 違う細胞型・違う状態は、分かれたまま（生物差の保存）

この2つは綱引きの関係です。手法ごとに「どちらをどれだけ重視するか」で性格が分かれます。

⚠️ 最大の落とし穴：消しすぎ（過補正） バッチ補正を強くかけすぎると、本物の細胞型の違いや処理の効果まで消えます。 確かに「そろった」けれど、サンプル間の違い（生物学的意味のある差）も消えています。「WT と KO がきれいに重なった！」は統合の成功ではなく、差を消してしまった失敗かもしれません。だから統合後は必ず評価します（最終節）。

1. 準備：正規化・特徴選択（HVG）

統合の前に、共通の下ごしらえをします。ここでは「個性の出る遺伝子」だけに絞る高変動遺伝子（HVG、highly variable genes）選択を行います。seurat_v3 という方法は生（なま）のカウントを必要とするので、先に生カウントを別の場所へ退避してから選びます。

import numpy as np
import scanpy as sc
import anndata as ad

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

# 生カウントを保存（seurat_v3 は生カウントを要求するため）
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,
)

🔧 環境メモ：flavor="seurat_v3" は scanpy の HVG 選択アルゴリズムの名前で、内部で scikit-misc を使います（R や Seurat 本体は不要）。未導入なら pip install scikit-misc。

💡 なぜ遺伝子を絞るの？ 遺伝子の大半は、どの細胞でもほぼ一定で「個性」が出ません（シグナルよりノイズを足すだけ）。変動の大きい遺伝子だけに絞ると、細胞集団の構造が見えやすくなり、計算も軽くなります。「人を見分けるのに、答えが分かれる質問だけ使う」のと同じ発想です。

2. Harmony で統合する

仕組み（やさしく）

Harmony は 主成分分析（PCA、principal component analysis：情報を圧縮して扱いやすくする手法）の空間で動きます。やっていることはシンプルで、

  1. 似た細胞をゆるくグループ分けし、
  2. 各グループの中で「カメラ（バッチ）ごとの偏り」を打ち消し、
  3. これを何度か繰り返してカメラ差をそろえる、

というものです（Korsunsky et al., 2019）。深層学習を使わない軽量・高速で、普通のPCで十分な方法です。

コード（正規化 → PCA → Harmony → 近傍 → クラスタリング）

まず正規化（公平に比較できるよう数値をそろえる）と対数変換をして、HVG に絞って PCA をかけます。

import scanpy as sc

# 正規化と対数変換
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata                       # 全遺伝子の正規化済みデータを退避

# HVG に絞ってスケーリング → PCA
adata_h = adata[:, adata.var.highly_variable].copy()
sc.pp.scale(adata_h, max_value=10)
sc.tl.pca(adata_h, n_comps=30)

次に Harmony を適用します。

import harmonypy

def harmony_integrate_robust(adata, batch_key, basis="X_pca", adjusted="X_pca_harmony"):
    ho = harmonypy.run_harmony(adata.obsm[basis], adata.obs, [batch_key])
    Z = ho.Z_corr
    if Z.shape[0] != adata.n_obs:        # 古い向き (PCs, cells) の場合は転置
        Z = Z.T
    adata.obsm[adjusted] = np.asarray(Z) # → (細胞数, PC数) で格納
    return adata

harmony_integrate_robust(adata_h, "batch")   # obsm["X_pca_harmony"] が作られる

# ★ 補正後の埋め込みで近傍グラフを作る（ここが最重要）
sc.pp.neighbors(adata_h, use_rep="X_pca_harmony", n_neighbors=15)
sc.tl.umap(adata_h)
sc.tl.leiden(adata_h, key_added="leiden_harmony",
             flavor="igraph", n_iterations=2, directed=False, resolution=1.0)

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

⚠️ 落とし穴①：Harmony の結果を「使い忘れる」（最重要） sc.pp.neighbors() は、何も指定しないと X_pca（補正していないPCA）を使います。Harmony を回しても、use_rep="X_pca_harmony" を指定しなければ補正は一切反映されません。「Harmony したのにサンプルで割れる」原因のNo.1がこれです。補正した結果を指定して使う——ここだけは必ず守ってください。

⚠️ 落とし穴②：harmonyはカウントデータそのものを補正しているわけではなく、UMAP上のデータ点を補正していることに注意が必要です。この後に、遺伝子発現量解析などを行う場合この点を考慮する必要があります。

💡 なぜ regress_out を入れないの？ 古い手順では、PCA の前に sc.pp.regress_out(["total_counts","pct_counts_mt"]) を挟むことがよくありました。しかし今のベストプラクティスでは、既定では行いません（Heumos et al., 2023）。技術指標の回帰除去は、本物の生物シグナルまで削るリスクがあり、効果も限定的だからです。

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

「UMAP（uniform manifold approximation and projection：高次元データを2次元の地図にする手法）できれいに重なった」だけでは不十分です。

# ① バッチが混ざっているか（混ざっていれば技術差は除けている）
sc.pl.umap(adata_h, color="batch", frameon=False)

# ② 既知マーカーで細胞型が保たれているか（消えていないか）
known_markers = ["GeneA", "GeneB", "GeneC"]   # 自分の系のマーカーに置換
sc.pl.umap(adata_h, color=known_markers, color_map="RdBu_r", frameon=False)

💡 数値で評価したいなら scib パッケージ（📖 詳しくは「scibとは？」記事。Python 完結なら scib-metrics） 目で見るのに加え、scib で数値化できます。代表的な指標は、

  • バッチ混合（よく混ざったか）：iLISI、kBET（高いほど良い）
  • 生物保存（細胞型が残ったか）：cLISI、ARI、NMI（既知ラベルとの一致度）

4. もう一方の選択肢：AI 手法の scVI（GPU 必須）

統合にはもう一つ、深層学習ベースの scVI という選択肢があります（📖 scVI 自体の詳細は「scVIとは？」記事）。強いバッチ効果や大規模データに強く、この後のダブレット検出（SOLO）・アノテーション・マルチモーダル解析へ発展させやすい一方、実用上 GPU が必須です。

📌 どちらを選ぶ？（詳しくは全体像の記事）

  • 小〜中規模／生物シグナルが明瞭／GPU が無い・素早く回したい → 本記事の Harmony
  • 大規模／強いバッチ効果／GPU が使え、発展解析に繋げたい → scVI（別記事）

まとめ

  • 統合のゴールは「カメラ差（バッチ）を消し、被写体の違い（生物差）は残す」の両立。消しすぎ（過補正）に注意し、必ず評価する。
  • Harmony は軽量・普通のPCで完結。ポイントは3つ：neighbors で use_rep="X_pca_harmony" を必ず指定／harmonypy 2.0 は直接呼びで回避／regress_out は使わない。
  • GPU が使え、強いバッチ・大規模・発展性を求めるなら scVI（別記事） を検討。

次の記事：ダブレットの検出 ― scrublet（GPU不要） — 2細胞が1つに化けた人工物を、普通のPCで完結する古典的手法で検出・除去します。

参考文献

  • Korsunsky, I., Millard, N., Fan, J., et al. (2019). Fast, sensitive and accurate integration of single-cell data with Harmony. Nature Methods, 16, 1289–1296. doi:10.1038/s41592-019-0619-0
  • 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
  • Single-cell best practices（統合の章）: https://www.sc-best-practices.org/cellular_structure/integration.html

コメント

タイトルとURLをコピーしました