シーケンス型 空間トランスクリプトーム:cell2location 実践 ─ 割合ではなく、細胞数を推定する

Spatial transcriptome
📚 この記事について
cell2location を実際に動かします。2段階のワークフロー、自分で決めなければならない2つのハイパーパラメータ、そして結果をどう検証するかを、コードとともに追います。
🔙 前の記事:シーケンス型 空間トランスクリプトーム:デコンボリューションの原理と手法分類
🔜 次の記事:シーケンス型 空間トランスクリプトーム:手法比較(CARD・RCTD・Tangram)
📌 前提scRNA-seq解析:細胞型アノテーション手法の全体像scvi-tools の全体像

空間トランスクリプトーム解析(spatial transcriptomics)のうちシーケンス型(Visium・Stereo-seq など NGS で読み出す方式)は、デコンボリューションが要ります。その中で cell2location(Kleshchevnikov et al., Nature Biotechnology, 2022)は、割合ではなく細胞数そのものを推定するという、他にない性質を持っています。階層ベイズモデルなので計算コストは高いのですが、16 手法のベンチマーク(Li et al., Nature Methods, 2022)でも Spotless(Sang-aram et al., eLife, 2023)でも総合上位に来ており、最初に試す価値があります。

1. cell2location が持っている3つの性質


  • 絶対量を返す:細胞型ごとの「細胞数」を推定する。割合ではないので、細胞密度の違いが残る。
  • 稀な細胞型を拾える:フィーチャ間で統計的な情報を借り合う設計なので、1つのフィーチャだけでは弱いシグナルも拾える。
  • 技術差をモデル化する:参照と空間データで遺伝子の捕まりやすさが違うことを、モデルの中で吸収する。

代償は計算コストです。GPU がほぼ必須で、GPU を使っても RCTD などより時間がかかります。その代わり、ベイズモデルなので推定の不確かさが分布として返ります

2. 2段階のワークフロー


図:cell2location は2段階で動く
図:cell2location は2段階で動く

段階1は参照 scRNA-seq だけを使います。負の二項回帰をあてて、細胞型ごとの発現シグネチャを推定します。このとき batch_key にドナーや実験日を指定すると、ドナー間の差を吸収した、きれいなシグネチャが得られます。

段階2で、そのシグネチャを空間データに当てます。ここで初めて空間データが登場します。

⚠️ どちらの段階でも、渡すのは生のカウント
負の二項モデルは、カウントの分布そのものを扱います。正規化・対数変換した行列を渡すと、モデルの前提が崩れます。adata.X に生のカウントが入っていることを、毎回確認してください。scvi-tools は正規化済みのデータを検知すると警告を出しますが、警告を無視しても処理自体は走ってしまいます。

3. 環境構築


bash
# GPU があることを前提にする。CPU でも動くが、実用的な速度にはならない
conda create -y -n c2l python=3.10
conda activate c2l
pip install cell2location

# GPU が見えているかを確認する
python -c "import torch; print(torch.cuda.is_available())"

4. 段階1 ─ 参照からシグネチャを作る


python
import scanpy as sc
from cell2location.models import RegressionModel
from cell2location.utils.filtering import filter_genes

# 参照 scRNA-seq。生のカウントを渡す
adata_ref.X = adata_ref.layers["counts"].copy()

# 情報量の少ない遺伝子を落とす
selected = filter_genes(
    adata_ref,
    cell_count_cutoff=5,
    cell_percentage_cutoff2=0.03,
    nonz_mean_cutoff=1.12,
)
adata_ref = adata_ref[:, selected].copy()

# batch_key にドナーや実験日を入れると、その差を吸収したシグネチャができる
RegressionModel.setup_anndata(
    adata=adata_ref,
    batch_key="donor",
    labels_key="cell_type",
)

mod = RegressionModel(adata_ref)
mod.train(max_epochs=250, batch_size=2500, train_size=1, lr=0.002)

# 事後分布を要約して、シグネチャを取り出す
adata_ref = mod.export_posterior(
    adata_ref, sample_kwargs={"num_samples": 1000, "batch_size": 2500}
)
inf_aver = adata_ref.varm["means_per_cluster_mu_fg"]
inf_aver.columns = adata_ref.uns["mod"]["factor_names"]
print(inf_aver.shape)   # (遺伝子, 細胞型)
引数 何を指定するか 注意点
batch_key ドナー、実験日、10x のバージョンなど ここを指定しないと、ドナー間の差がシグネチャに混ざる
labels_key 細胞型のラベル この解像度が、そのまま結果の解像度になる
categorical_covariate_keys その他の技術共変量 任意
💡 遺伝子フィルタの意味
filter_genes は、ほとんど検出されない遺伝子を落とします。cell_count_cutoff は「最低これだけの細胞で検出されていること」、cell_percentage_cutoff2 は「細胞の何%以上で検出されていること」、nonz_mean_cutoff は「検出された細胞での平均発現がこれ以上あること」。3つの条件を組み合わせて、稀だが特異的な遺伝子(稀な細胞型のマーカー)を残しつつ、ノイズだけの遺伝子を落とします。

5. 2つのハイパーパラメータ ─ ここが実質のすべて


cell2location で自分が決めるのは、たった2つです。そしてこの2つを間違えると、結果は意味を失います

図:自分で決めるハイパーパラメータは、2つだけ
図:自分で決めるハイパーパラメータは、2つだけ

6. N_cells_per_location ─ 組織像から数える


「1つのフィーチャに、平均で何個の細胞がいるか」を指定します。この値が、推定される細胞数の全体のスケールを決めます。大きすぎれば水増しになり、小さすぎれば細胞を取りこぼします。

公式ドキュメントは、組織像から推定することを明示的に求めています。チュートリアルの値(30)はリンパ節の値であって、あなたの組織の値ではありません。

python
import numpy as np
from stardist.models import StarDist2D

# H&E 画像から核を検出して、1スポットあたりの細胞数を見積もる
model = StarDist2D.from_pretrained("2D_versatile_he")
labels, _ = model.predict_instances(he_norm)
n_nuclei = int(labels.max())

# 組織が載っている面積を µm² に直す
area_um2 = tissue_area_px * (um_per_px ** 2)
density = n_nuclei / area_um2                # 核の数 / µm²

# Visium のスポットは直径 55 µm の円。面積は およそ 2,400 µm²
spot_area = np.pi * (55 / 2) ** 2
print("N_cells_per_location の目安:", round(density * spot_area, 1))
⚠️ チュートリアルの 30 を、そのまま使わない
リンパ節は細胞が非常に密に詰まった組織です。脳、肝、筋、脂肪では、まったく違う値になります。核を数えるのに 10 分かければ、その後の解析すべてが信用できるものになります。StarDist(Schmidt et al., 2018)を使えば、H&E 画像から核を自動で数えられます。この手順は、Visium HD のビン→細胞再構成でも使うので、覚えておくと得をします。

7. detection_alpha ─ QC の診断と同じ判断


スライド内で、RNA の検出感度がどれだけばらつくかを、モデルにどう伝えるかのパラメータです。値が大きいほど強く正則化し、「感度は場所によらずほぼ一定」と仮定します。

意味 使う場面
200 強く正則化する(感度はほぼ一定と仮定) スライド内の技術変動が小さいとき
20 感度のばらつきを許す スライド内の技術変動が大きいとき(現在の既定値)

公式ドキュメントは、判断基準をこう述べています。総 RNA カウントの、バッチ内での変動が、組織像で見て「そこに細胞が多いから」では説明できないなら、detection_alpha=20 を使う。

💡 これは、QC でやった診断そのもの
「総カウントの変動は、細胞密度で説明できるのか、それとも技術由来なのか」──これは QC と組織検出の記事で、組織像と並べて見分けた問いと、まったく同じです。総カウントには細胞密度という生物学が入っている(Bhuva et al., Genome Biology, 2024)。その生物学で説明できない変動があるなら、それは技術ノイズであり、detection_alpha を下げてモデルに吸収させるべきだ、という理屈です。前処理でやった診断が、そのままハイパーパラメータの選択になります。
⚠️ 公式は「両方を試すこと」を推奨している
cell2location の開発元は、当初 200 を既定値にしていました。論文で使ったマウス脳とヒトリンパ節のデータは技術変動が小さく、強い正則化のほうが、推定した総細胞数と組織像の核数の一致がよかったからです。しかしその後、ヒト組織の Visium では技術変動が大きいことが多いと分かり、既定値は 20 に変わりました。20 と 200 の両方で走らせて、結果を比べてください。

8. 段階2 ─ 空間データに当てる


python
from cell2location.models import Cell2location

# 共通遺伝子だけに絞る
common = [g for g in inf_aver.index if g in adata_vis.var_names]
print("共通遺伝子:", len(common))
adata_vis = adata_vis[:, common].copy()
inf_aver = inf_aver.loc[common, :]

# ここでも生のカウントを渡す
adata_vis.X = adata_vis.layers["counts"].copy()
Cell2location.setup_anndata(adata=adata_vis, batch_key="sample")

mod = Cell2location(
    adata_vis,
    cell_state_df=inf_aver,
    N_cells_per_location=8,   # 組織像から数えた値に置き換える
    detection_alpha=20,       # 20 と 200 の両方を試す
)
mod.train(max_epochs=30000, batch_size=None, train_size=1)

adata_vis = mod.export_posterior(
    adata_vis,
    sample_kwargs={"num_samples": 1000, "batch_size": mod.adata.n_obs},
)

9. 結果を読む ─ なぜ q05 なのか


cell2location はベイズモデルなので、返ってくるのは1つの数字ではなく事後分布です。obsm には、平均・標準偏差・5% 分位点(q05)・95% 分位点(q95)が入ります。

図:返るのは1つの数字ではなく、分布
図:返るのは1つの数字ではなく、分布

実務では q05 を使うのが標準です。「少なくともこれだけの細胞はいる」という、保守的な推定だからです。

平均を使うと、推定に自信がない場所でも高い値が出ます。事後分布が広いということは「よく分からない」ということなのに、平均を取ると、それなりの値が返ってしまう。その結果、稀な細胞型を過剰に検出しがちになります。

python
import squidpy as sq

# q05 = 5% 分位点。「少なくともこれだけはいる」という保守的な推定
abund = adata_vis.obsm["q05_cell_abundance_w_sf"]
abund.columns = adata_vis.uns["mod"]["factor_names"]
print(abund.iloc[:3].round(2))

# 細胞型ごとに、組織の上に描く
for ct in abund.columns[:4]:
    adata_vis.obs[ct] = abund[ct].values
sq.pl.spatial_scatter(adata_vis, color=list(abund.columns[:4]), ncols=2)

# 推定された総細胞数。これが検証の要になる
adata_vis.obs["n_cells_est"] = abund.sum(axis=1).values
sq.pl.spatial_scatter(adata_vis, color="n_cells_est")

10. モデルの QC ─ 結果を信じてよいか


python
# 学習曲線。下がりきっていなければ、max_epochs を増やす
mod.plot_history(1000)

# モデルの QC。総 RNA カウント と 推定した総細胞数 の関係を見る
mod.plot_QC()
# 直線的に並ぶのが望ましい。大きく外れる点があれば、
# その場所は技術的な問題を抱えているか、参照が合っていない

# 推定した総細胞数と、組織像の核数を突き合わせる
# ここがずれるなら、N_cells_per_location が合っていない
  • 学習曲線:ELBO が下がりきっていなければ、max_epochs を増やす。
  • plot_QC():総 RNA カウントと、推定された総細胞数の関係を見る。直線的に並ぶのが望ましい。大きく外れる点は、技術的な問題か、参照が合っていない場所。
  • 推定された総細胞数 vs 組織像の核数:これが最も重要な検証。系統的にずれるなら、N_cells_per_location が合っていない。
⚠️ そして、マーカーを組織の上に投影する
モデルの QC が通っても、それは「モデルがデータに当てはまった」というだけです。生物学的に正しいかは別の話です。割り当てられた細胞型ごとに、そのマーカー遺伝子の発現を組織の上に描き、分布が重なるかを必ず確認してください。組織学的にありえない場所に細胞型が出ているなら、参照の欠落・spillover・QC の甘さの3つを疑い、前処理に戻ります。

11. 発展 ─ 細胞型の共局在を見つける


cell2location には、推定した細胞数の行列を NMF で分解して、どの細胞型が一緒に出てくるかを見つける機能があります。

python
from cell2location import run_colocation

# 細胞型が一緒に出てくるパターン(微小環境)を NMF で見つける
res, adata_vis = run_colocation(
    adata_vis,
    model_name="CoLocatedGroupsSklearnNMF",
    train_args={
        "n_fact": [8],           # 微小環境の数
        "sample_name_col": "sample",
        "n_restarts": 3,
    },
)
# 「腫瘍細胞 + 疲弊 T細胞 + マクロファージ」のような組み合わせが見つかる
# これがニッチ解析の入口になる

「腫瘍細胞 + 疲弊 T細胞 + マクロファージ」のような組み合わせが、1つの因子として出てきます。これは組織の中の微小環境(ニッチ)にほかなりません。ニッチ解析への、そのままの入口になります。

12. つまずきやすい点


症状 原因 対処
全フィーチャで細胞数が同じくらいになる detection_alpha が大きすぎる 20 にして再実行する
推定した総細胞数が、核数と大きく違う N_cells_per_location が合っていない 組織像から数え直す
稀な細胞型が、いたるところに出る 平均を使っている q05 を使う
学習が終わらない max_epochs が大きすぎる、または GPU がない plot_history で収束を確認し、必要なだけに減らす
共通遺伝子が極端に少ない 参照と空間のプラットフォームが違いすぎる 参照を見直す。FFPE には FFPE 由来の参照が望ましい
結果が実行ごとに変わる 収束していない max_epochs を増やす。乱数シードを固定する

まとめ


  • cell2location は2段階。参照からシグネチャを作り、それを空間に当てる。どちらの段階でも生のカウントを渡す。
  • 自分で決めるのは2つのハイパーパラメータだけ。そしてこの2つがすべてを決める。
  • N_cells_per_location は組織像から数える。チュートリアルの 30 はリンパ節の値であって、あなたの組織の値ではない。
  • detection_alpha は、QC でやった診断そのもの。総カウントの変動が細胞密度で説明できるなら 200、できないなら 20。両方試す。
  • 返るのは事後分布q05(保守的な推定)を使う。平均を使うと、稀な細胞型を過剰に検出する。
  • 検証の要は、推定した総細胞数が、組織像の核数と合っているか
  • NMF による共局在解析が、そのままニッチ解析の入口になる。

関連記事


参考文献


  • Kleshchevnikov, V., Shmatko, A., Dann, E., et al. (2022). Cell2location maps fine-grained cell types in spatial transcriptomics. Nature Biotechnology, 40(5), 661–671. doi:10.1038/s41587-021-01139-4
  • Cable, D. M., Murray, E., Zou, L. S., et al. (2022). Robust decomposition of cell type mixtures in spatial transcriptomics. Nature Biotechnology, 40(4), 517–526. doi:10.1038/s41587-021-00830-w
  • Li, B., Zhang, W., Guo, C., et al. (2022). Benchmarking spatial and single-cell transcriptomics integration methods for transcript distribution prediction and cell type deconvolution. Nature Methods, 19(6), 662–670. doi:10.1038/s41592-022-01480-9
  • Sang-aram, C., Browaeys, R., Seurinck, R., & Saeys, Y. (2023). Spotless, a reproducible pipeline for benchmarking cell type deconvolution in spatial transcriptomics. eLife, 12, RP88431. doi:10.7554/eLife.88431
  • Bhuva, D. D., Tan, C. W., Salim, A., et al. (2024). Library size confounds biology in spatial transcriptomics data. Genome Biology, 25, 99. doi:10.1186/s13059-024-03241-7
  • Schmidt, U., Weigert, M., Broaddus, C., & Myers, G. (2018). Cell detection with star-convex polygons. In Medical Image Computing and Computer Assisted Intervention (MICCAI) 2018, 265–273. doi:10.1007/978-3-030-00934-2_30
  • Polański, K., Bartolomé-Casado, R., Sarropoulos, I., et al. (2024). Bin2cell reconstructs cells from high resolution Visium HD data. Bioinformatics, 40(9), btae546. doi:10.1093/bioinformatics/btae546
  • Wolf, F. A., Angerer, P., & Theis, F. J. (2018). SCANPY: large-scale single-cell gene expression data analysis. Genome Biology, 19, 15. doi:10.1186/s13059-017-1382-0
  • Palla, G., Spitzer, H., Klein, M., et al. (2022). Squidpy: a scalable framework for spatial omics analysis. Nature Methods, 19(2), 171–178. doi:10.1038/s41592-021-01358-2

コメント

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