scRNA-seq解析:差次的発現(DE)― マーカー検出と条件間比較(pseudobulk)を使い分ける

scRNA-seq

 

📚 この記事について:「scRNA-seq解析 実践シリーズ」下流解析編。差次的発現(differential expression, DE)には性質の異なる2つのタスクがあり、混同すると誤った結論につながります。両者の使い分けと、条件間比較で必須の pseudobulk を解説します。

🔙 前の記事下流解析の全体像 / 🔜 次の記事機能エンリッチメント・パスウェイ解析

前提:細胞型アノテーション済みのデータ(→アノテーション手法の全体像)。生カウントを adata.layers["counts"] に保持。Python の基本。


この記事のゴール

DE の2つのタスク――「マーカー検出」と「条件間比較」――を区別し、それぞれ正しい方法で実行できること。特に条件間比較では pseudobulk が必須である理由を理解します。


DE には2つの異なるタスクがある

同じ「差次的発現」でも、目的によってやり方がまったく違います

  • (a) マーカー検出:1つのクラスタ(細胞型)を残りの細胞と比べ、そのクラスタで高い遺伝子を見つける。目的は細胞型の特徴づけ・命名
  • (b) 条件間比較同じ細胞型の中で、条件(WT vs KO、健常 vs 疾患)間で変わる遺伝子を見つける。目的は条件(処理・遺伝型)の効果
DEの2つのタスク。(a) マーカー検出は「1クラスタ vs 残りの細胞」で、細胞を単位にしてよく、rank_genes_groups で行う(特徴づけ用)。(b) 条件間比較は「同じ細胞型内で条件群どうし」を比べ、個体(生物学的反復)を単位にする pseudobulk + DESeq2 で行う(条件の効果)。
DEの2つのタスク。(a) マーカー検出は「1クラスタ vs 残りの細胞」で、細胞を単位にしてよく、rank_genes_groups で行う(特徴づけ用)。(b) 条件間比較は「同じ細胞型内で条件群どうし」を比べ、個体(生物学的反復)を単位にする pseudobulk + DESeq2 で行う(条件の効果)。

⚠️ 混同が事故のもと
マーカー検出の手法(rank_genes_groups)を条件間比較に流用しないでください。後述のとおり、条件比較で細胞を独立サンプル扱いすると偽陽性が大量に出ます


(a) マーカー検出 ― 1クラスタ vs 残り

細胞型を特徴づけ・命名するための探索です。rank_genes_groups(既定は Wilcoxon)でクラスタごとに上位遺伝子を出します。ここは細胞を単位にしてよい(記述的な探索だからです)。

pythonimport scanpy as sc

# 各細胞型 vs 残りの細胞で、特徴的な遺伝子を出す(マーカー検出)
sc.tl.rank_genes_groups(adata, groupby="cell_type", method="wilcoxon")
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)

# 結果を表として取り出す
markers = sc.get.rank_genes_groups_df(adata, group=None)

💡 これは手動による細胞型アノテーションで使ったものと同じです。あくまで「どの遺伝子がそのクラスタを特徴づけるか」を見るためで、条件の効果を検定するものではありません


(b) 条件間比較は pseudobulk で

なぜ pseudobulk が必要か

同じ個体(マウス・患者)から取れた細胞は、互いに独立ではありません。それなのに細胞を1サンプルとして検定すると、「個体間のばらつき」ではなく「1個体内のばらつき」を見てしまい、p値が過小評価され偽陽性が爆発します(pseudoreplication)。

解決は pseudobulk個体×細胞型ごとに細胞の生カウントを合計して「サンプル(生物学的反復)」を作り、サンプル単位で bulk RNA-seq と同じ DESeq2 などにかけます。これで個体間のばらつきを正しく評価できます(bulk の発現変動解析はRNA-seq解析入門⑥:発現変動解析編も参照)。

pseudobulk の流れ。各個体×細胞型ごとに細胞の生カウントを合計し、行=個体サンプル(WT_1, WT_2, KO_1, KO_2…)の擬似バルク行列を作る。これを DESeq2 にかけ、log2 fold-change と補正 p 値を得て volcano で可視化する。細胞を直接サンプル扱いすると偽陽性が増える。
pseudobulk の流れ。各個体×細胞型ごとに細胞の生カウントを合計し、行=個体サンプル(WT_1, WT_2, KO_1, KO_2…)の擬似バルク行列を作る。これを DESeq2 にかけ、log2 fold-change と補正 p 値を得て volcano で可視化する。細胞を直接サンプル扱いすると偽陽性が増える。

⚠️ 生物学的反復が必須
pseudobulk は条件ごとに複数個体(最低でも各2〜3、多いほど良い)が必要です。1個体ずつしかないと、条件間比較は統計的に成立しません。

コード(decoupler + PyDESeq2)

pythonimport scanpy as sc
import decoupler as dc
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds  import DeseqStats

# 前提:layers["counts"]=生カウント、obs に sample(個体ID)・condition(WT/KO)・cell_type

# 1) 個体×細胞型ごとに生カウントを合計して pseudobulk を作る
pdata = dc.get_pseudobulk(
    adata, sample_col="sample", groups_col="cell_type",
    layer="counts", mode="sum", min_cells=10, min_counts=1000,
)

# 2) 対象の細胞型を1つ選ぶ(例:興奮性ニューロン ExN)
sub = pdata[pdata.obs["cell_type"] == "ExN"].copy()

# 3) 低発現の遺伝子を除く(推奨)
keep = dc.filter_by_expr(sub, group="condition", min_count=10, min_total_count=15)
sub = sub[:, keep].copy()

# 4) DESeq2(design = ~condition)で条件間比較
dds = DeseqDataSet(adata=sub, design_factors="condition", refit_cooks=True)
dds.deseq2()
stat = DeseqStats(dds, contrast=["condition", "KO", "WT"])
stat.summary()
res = stat.results_df                      # log2FoldChange, pvalue, padj
sig = res[(res["padj"] < 0.05) & (res["log2FoldChange"].abs() > 1.0)]

🔧 環境メモpip install decoupler pydeseq2。API はバージョンで変わることがあるため、decouplerget_pseudobulk / filter_by_expr)と pydeseq2DeseqDataSet / DeseqStats)は公式チュートリアルに合わせて使ってください。交絡(性別・処理日など)は design_factors=["batch", "condition"] のように追加します。

反復が少ないとき(2 vs 2 など)

「各条件2サンプルしかないので、検出力を上げたくて細胞単位(rank_genes_groups)で WT vs KO を比べた」――よくある選択ですが、これは擬似反復のため適切ではありません。要点は2つです。

  • 細胞単位の「検出力」は見かけ(偽物):有効な検定単位は個体(ここでは4)です。細胞を増やしても「個体4」という情報量は増えないのに、細胞単位の検定は数千細胞を有効サンプルサイズとみなし、p値が過度に小さくなります。結果、条件の効果ではなく1サンプルの偏りを「有意」と拾いがちです。
  • 2 vs 2 でも正解は pseudobulk:DESeq2/edgeR は各群2反復が分散推定の最小要件で、ちゃんと動きます。遺伝子間で分散情報を共有して安定化するため、素朴な遺伝子ごとの検定より頑健です。つまり方法は pseudobulk のままで、ただし検出力が限られ、強く一貫した効果しか拾えないという設計上の制約を受け入れることになります。

💬 すでに細胞単位で出した結果はどう扱うか
fold change(平均発現の差)は記述統計なので参考にできます。一方、p値・有意性・ランキングは有効 N が過大なため信頼できません。正式な統計としては報告しないのが安全です。
– 実務的には、いまのリストを探索的(仮説生成)と位置づけ、pseudobulk で検定し直すのが堅実です。
– 再検定後は頑健性を確認:有意な遺伝子が両方の WT・両方の KO で一貫しているか(サンプルレベルの PCA・ヒートマップで1サンプル駆動でないか)を見ます。
– 細胞単位の解像度を保ちたい場合は、サンプルを変量効果に入れた混合モデル(例:NEBULA)で擬似反復を補正できますが、まずは pseudobulk が簡便で実績豊富です。
– 本質的な検出力は反復数で決まります。可能なら各条件3個体以上が望ましく、2 vs 2 は誠実には「探索的」と明記します。


DE ツールの早見表

ツール 言語 主な用途 ひとこと
scanpy rank_genes_groups Python マーカー検出(1 vs 残り) 探索・特徴づけ向き(本記事 (a)
decoupler + PyDESeq2 Python 条件間比較(pseudobulk) 現在の推奨ワークフロー(本記事 (b)
edgeR / limma-voom R 条件間比較(pseudobulk) bulk の定番。pseudobulk で利用
muscat R 条件間比較(複数手法) sc 用 DS 解析の枠組み
MAST R 単一細胞レベル DE 細胞単位(混合モデルで反復を考慮可)

落とし穴

⚠️ つまずきやすい点
pseudoreplication:条件間比較で細胞を独立サンプル扱いしない。必ず pseudobulk(個体単位)で。
反復数:条件ごとに複数個体が要る。技術反復ではなく生物学的反復
低発現遺伝子の除去:検出力とFDRのために、ほとんど発現しない遺伝子は外す。
交絡:条件と技術要因(処理日・バッチ)が混ざっていないか。混ざるなら design に入れる。
効果量も見る:補正 p 値だけでなく log2FC を併用。有意 ≠ 生物学的に重要
細胞数の少ない細胞型:pseudobulk の合計が小さいと不安定。min_cells で足切り。


まとめ

  • DE は2タスク(a) マーカー検出(1クラスタ vs 残り・rank_genes_groups・細胞単位)と(b) 条件間比較(同一細胞型内・条件群どうし)。
  • 条件間比較は pseudobulk が必須。細胞を独立扱いすると偽陽性が爆発するため、個体×細胞型で合計してサンプル単位で DESeq2 にかける(decoupler + PyDESeq2)。
  • 生物学的反復低発現除去交絡の処理効果量の併用に注意。
  • 得られた DE 遺伝子は、次の機能エンリッチメントで生物学的に解釈する。

次の記事:機能エンリッチメント・パスウェイ解析 — DE 遺伝子の生物学的な意味(GO/KEGG・GSEA・経路/TF活性)を読み解く方法に進みます。


関連記事

下流解析の地図と前後

関連する前段・基盤


参考文献

  • Squair, J. W., Gautier, M., Kathe, C., et al. (2021). Confronting false discoveries in single-cell differential expression. Nature Communications, 12, 5692. doi:10.1038/s41467-021-25960-2
  • Zimmerman, K. D., Espeland, M. A., & Langefeld, C. D. (2021). A practical solution to pseudoreplication bias in single-cell studies. Nature Communications, 12, 738. doi:10.1038/s41467-021-21038-1
  • Crowell, H. L., Soneson, C., Germain, P.-L., et al. (2020). muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nature Communications, 11, 6077. doi:10.1038/s41467-020-19894-4
  • Murphy, A. E., & Skene, N. G. (2022). A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis. Nature Communications, 13, 7851. doi:10.1038/s41467-022-35519-4
  • Muzellec, B., Teleńczuk, M., Cabeli, V., & Andreux, M. (2023). PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics, 39, btad547. doi:10.1093/bioinformatics/btad547
  • Badia-i-Mompel, P., Vélez Santiago, J., Braunger, J., et al. (2022). decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinformatics Advances, 2, vbac016. doi:10.1093/bioadv/vbac016
  • Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. doi:10.1186/s13059-014-0550-8

 

コメント

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