この記事のゴール
QC を通したきれいなデータから、「細胞のグループ(クラスタ)」を見つけるところまでを、4ステップで実行できるようになること。各ステップで「何のために」を押さえ、代表的な選択肢と典型的な誤りも整理します。
全体の流れ(まずここを押さえる)
4ステップは、次の順番でデータを変換していきます。ここで一番こんがらがるのは「PCA」と「UMAP」の役割の違いなので、先に全体像で整理します。
ポイントは2つです。
- PCA は解析の土台(近傍グラフ・クラスタリング・UMAP すべての入力)。
- UMAP は可視化専用。クラスタリングは UMAP の2D座標ではなく、近傍グラフ(=PCA空間)の上で行います。
1. 正規化(Normalization)
なぜ必要か
細胞ごとに、検出された総 mRNA 分子数(シーケンス深度)が大きく違います。この違いの多くは技術的要因(捕捉効率・PCR・シーケンス量)で、生物学的な発現差ではありません。揃えないと「深く読まれた細胞=高発現」という偽の差が生まれます。
基本のコード(shifted log)
最も標準的なのが shifted log(深度を揃えて log を取る)です。本シリーズもこれを既定にします。
pythonimport scanpy as sc
# adata: ダブレット除去済み。生カウントは adata.layers["counts"] に退避済み
sc.pp.normalize_total(adata, target_sum=1e4) # 各細胞の総カウントを 1e4 に揃える(CP10K)
sc.pp.log1p(adata) # log(1 + x) で分散を安定化
代表的な手法(迷ったら shifted log)
| 手法 | Scanpy | 入力 | 向く用途 |
|---|---|---|---|
| shifted log(★本記事の既定) | normalize_total + log1p |
生カウント | 分散安定化。次元削減・DEG に手堅い |
| analytic Pearson residuals | sc.experimental.pp.normalize_pearson_residuals |
生カウント | 希少細胞型の検出・HVG 選択に有利 |
| scran pooling | R/scran(別環境) | 生カウント | 深度のばらつきが大きい時に頑健 |
💡 どれを選ぶか:大規模ベンチマーク(Ahlmann-Eltze & Huber, 2023)と sc-best-practices は、shifted log が次元削減・差次的発現遺伝子(DEG、differentially expressed genes)に手堅い一方、analytic Pearson residuals は生物学的に変動する遺伝子の選択や希少細胞型の同定に向く、と整理しています(Lause et al., 2021)。まず shifted log、希少細胞型を狙うなら Pearson residuals が実用的です。
⚠️ 落とし穴
– target_sum=1e4(CP10K)は慣習で、絶対的な正解ではありません。None にすると中央値で正規化します。
– scVI 経路には正規化を渡さない:scVI(single-cell variational inference)は生カウントの統計モデルを内部で扱うため、counts レイヤー(生)を入力にします(→ 統合:scVI の記事)。
2. 特徴選択(Highly Variable Genes, HVG)
なぜ必要か
数万ある遺伝子の大半は、細胞間でほぼ一定です。一定の遺伝子は構造の手がかりにならず、ノイズと計算量を増やすだけです。変動の大きい遺伝子(HVG)だけに絞ると、細胞集団の構造が見えやすくなり、計算も軽くなります。
基本のコード(分散ベース)
python# 変動の大きい遺伝子を分散ベース(seurat flavor)で選ぶ
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata) # 選ばれた遺伝子を確認
adata.raw = adata # 正規化済みの全遺伝子を退避(後で DEG 等に使う)
adata = adata[:, adata.var.highly_variable].copy() # HVG だけにスライス
flavor(選択アルゴリズム)の比較
flavor と「入力データの種類」の対応を間違えないことが最重要です。
| flavor | Scanpy 指定 | 入力データ | 備考 |
|---|---|---|---|
| seurat(分散ベース・★本記事) | flavor="seurat"(既定) |
正規化・log 後 | min_mean / max_mean / min_disp で指定 |
| seurat_v3 | flavor="seurat_v3", layer="counts" |
生カウント | n_top_genes で本数指定。scVI と相性が良い |
| pearson_residuals | sc.experimental.pp.highly_variable_genes(flavor="pearson_residuals") |
生カウント | 生物学的に変動する遺伝子に強い |
⚠️ 落とし穴
– seurat(既定)は正規化・log 後のデータが前提。seurat_v3 と pearson_residuals は生カウントが前提(layer="counts")。入力と flavor がちぐはぐだと HVG が無意味になります。
– 本数の目安は 2000〜3000。複数サンプルがあるなら batch_key を指定して、特定サンプルだけの変動に引っ張られないようにします。
3. 次元削減(Dimensionality Reduction)
ここは役割の違う2段階を混同しないことが最重要です(冒頭の図を再確認)。
3-1. PCA(線形)— 解析の土台
PCA はノイズを落として主要な変動軸(主成分)に圧縮します。近傍グラフ・クラスタリング・UMAP の入力になる、解析の土台です。
python# (任意)技術ノイズの影響を回帰で除く
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])
sc.tl.pca(adata, n_comps=50, svd_solver="arpack") # 主成分に圧縮
sc.pl.pca_variance_ratio(adata, log=True) # 何 PC 使うかの目安
🔧 regress_out は任意:技術ノイズ(総カウント・ミトコンドリア割合)の影響を回帰で除く処理です。効果はデータ次第で、近年は省略する例も増えています(生物学的シグナルまで削るリスクがあるため)。元の解析に合わせて入れていますが、外して比較するのも有効です。標準化(sc.pp.scale)を挟むパイプラインもあります。
💡 PC 数の決め方:pca_variance_ratio の曲線が寝るあたり(elbow)を目安に 30〜50 を選びます。多すぎるとノイズを拾い、少なすぎると構造を取りこぼします。
📌 統合する場合:複数サンプルのバッチ効果をそろえるなら、PCA の後に Harmony や scVI を挟みます。その場合は近傍グラフを use_rep="X_pca_harmony" や use_rep="X_scVI" で作ります。
3-2. UMAP / t-SNE(非線形)— 可視化「専用」
UMAP(McInnes et al., 2018)は、高次元の構造を2Dに落として見るためのものです。
pythonsc.pp.neighbors(adata, n_neighbors=15, n_pcs=50) # PCA 空間で近傍グラフを作る
sc.tl.umap(adata)
sc.pl.umap(adata, color="leiden")
| 手法 | Scanpy | 用途・特徴 |
|---|---|---|
| UMAP(標準) | sc.tl.umap |
局所構造を保持。大域的な距離は保存しない |
| t-SNE | sc.tl.tsne |
局所構造重視。大規模で遅く、大域構造はさらに弱い |
⚠️ UMAP の3大誤用
1. UMAP 座標でクラスタリングしない。クラスタリングは近傍グラフ(=PCA空間)で行います。2D 座標は歪んでいます。
2. 点間の距離・クラスタの大きさ・配置を定量解釈しない。UMAP は局所構造を優先し、大域的な距離は保存しません。「A と B が近い/遠い」を生物学的に読むのは危険です(Chari & Pachter, 2023)。
3. パラメータ(n_neighbors, min_dist)で見た目が大きく変わる。1枚の UMAP を「真実」と思わないこと。
4. クラスタリング(Clustering)
Leiden を使う(Louvain ではなく)
現在の標準はグラフベースの Leiden 法です。近傍グラフ上で「密につながった細胞のかたまり」を見つけます。かつて広く使われた Louvain 法には、まれに内部がつながっていない(disconnected)クラスタを作る欠陥があり、Leiden はこれを解消して連結性を保証します(Traag et al., 2019)。
pythonsc.tl.leiden(
adata, key_added="leiden",
flavor="igraph", n_iterations=2, directed=False, # 現行の推奨指定
resolution=1.0,
)
sc.pl.umap(adata, color="leiden")
🔧 更新メモ:sc.tl.leiden(adata) をオプション無しで呼ぶと旧実装が使われ FutureWarning が出ます。上記のように flavor="igraph", n_iterations=2 を明示するのが現行の推奨です。Louvain(sc.tl.louvain)は別パッケージが必要で性能でも劣るため、新規解析では Leiden を使います。
resolution は「粒度のつまみ」
resolution を上げるとクラスタは細かく(多く)なり、下げると粗く(少なく)なります。唯一の正解はありません。
python# 複数の解像度で試し、安定する粒度を選ぶ
for res in [0.3, 0.5, 0.8, 1.0]:
sc.tl.leiden(adata, key_added=f"leiden_r{res}",
flavor="igraph", n_iterations=2, directed=False, resolution=res)
⚠️ 落とし穴:クラスタ数は「真実」ではない
– クラスタリングは探索的な仮説生成であって、確定した細胞型ではありません。
– 複数解像度を比較し、既知マーカーで意味が通る粒度を選びます(→ 次回のアノテーション)。
– 過剰クラスタリングは「ノイズの島」を、過少クラスタリングは「異なる細胞型の融合」を生みます。
全体を通したコード
ここまでを1つの関数にまとめると、次のようになります(統合が必要なら PCA の後に Harmony / scVI を挟みます)。
pythonimport scanpy as sc
def standard_preprocess(adata, n_pcs=50, n_neighbors=15, resolution=1.0, random_state=0):
"""正規化 → HVG → 回帰 → PCA → 近傍グラフ → UMAP → Leiden を通しで実行して返す。"""
# 1. 正規化(shifted log)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
# 2. 高変動遺伝子(HVG, 分散ベース)→ raw に退避 → HVG だけにスライス
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata.raw = adata
adata = adata[:, adata.var.highly_variable].copy()
# 3. 技術ノイズの回帰(任意)→ PCA
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])
sc.tl.pca(adata, n_comps=n_pcs, svd_solver="arpack")
# 統合する場合はここで:sc.external.pp.harmony_integrate(adata, "batch")
# 4. 近傍グラフ → UMAP → Leiden
sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs)
sc.tl.umap(adata, random_state=random_state)
sc.tl.leiden(adata, key_added="leiden", flavor="igraph",
n_iterations=2, directed=False, resolution=resolution)
return adata
adata = standard_preprocess(adata)
sc.pl.umap(adata, color="leiden")
この関数のなかで実際に起きているのは、AnnData に対する列の追加とスライスです。とくに adata.raw = adata と adata[:, adata.var.highly_variable] の2行は、何を退避して何を捨てているのかが分かれ目になります。ここが曖昧なまま進むと、後の差次的発現解析で「HVG に絞った後の遺伝子しか見えない」といった形で跳ね返ってきます。
本記事で繰り返した「flavor と入力データの対応を間違えない」という注意点も、突き詰めれば adata.X と adata.layers["counts"] のどちらを渡しているかという pandas と NumPy の話です。この層をきちんと踏みたい場合は、『改訂 独習Pythonバイオ情報解析』が、NumPy の配列と pandas の DataFrame を土台から扱ったうえで Scanpy に進む構成になっています。関数を呼ぶだけの状態から、引数の意味とデータへの副作用を判断できる状態に移りたい段階に向いています。
![]()
つまずきやすいポイント(まとめ)
- 正規化:既定は shifted log。希少細胞型狙いなら Pearson residuals。scVI には生カウントを渡す。
- HVG:flavor と入力データ(生 vs 正規化済み)の対応に注意。
seuratは正規化後、seurat_v3/pearson_residualsは生カウント。 - 次元削減:PCA は解析の土台、UMAP は可視化専用。UMAP の距離・配置を定量解釈しない。クラスタリングは近傍グラフ上で行う。
- クラスタリング:Leiden(
flavor="igraph")。resolution は粒度のつまみで、クラスタ数は真実ではない。
次の記事:細胞型アノテーション(マーカー遺伝子の同定) — 得られたクラスタに意味を与えます。マーカー遺伝子の同定から始める細胞型アノテーションを扱います。
関連記事
- 📖 scRNA-seqの前処理 入門 ─ 全体像と進め方 … このシリーズの全体像(ハブ)
- 📖 データの読み込みと品質管理(QC) / 📖 ダブレット検出(scrublet・SOLO)… 本記事の前段
- 📖 サンプル統合(Harmony・scVI)… 複数サンプルを扱う場合は PCA の後に
- 📖 細胞型アノテーション手法の全体像 / 📖 手動アノテーション(マーカー遺伝子の同定) … 次の段階
- 📖 AnnDataのデータ構造 / 📖 scvi-tools の全体像 … 土台・AI手法の地図
参考文献
- Ahlmann-Eltze, C., & Huber, W. (2023). Comparison of transformations for single-cell RNA-seq data. Nature Methods, 20, 665–672. doi:10.1038/s41592-023-01814-1
- Lause, J., Berens, P., & Kobak, D. (2021). Analytic Pearson residuals for normalization of single-cell RNA-seq UMI data. Genome Biology, 22, 258. doi:10.1186/s13059-021-02451-7
- Traag, V. A., Waltman, L., & van Eck, N. J. (2019). From Louvain to Leiden: guaranteeing well-connected communities. Scientific Reports, 9, 5233. doi:10.1038/s41598-019-41695-z
- McInnes, L., Healy, J., & Melville, J. (2018). UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv:1802.03426
- Chari, T., & Pachter, L. (2023). The specious art of single-cell genomics. PLOS Computational Biology, 19(8), e1011288. doi:10.1371/journal.pcbi.1011288
- 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


コメント