scRNA-seq解析：機能エンリッチメント・パスウェイ解析

scRNA-seq
2026.06.28

 

📚 この記事について：「scRNA-seq解析 実践シリーズ」下流解析編。差次的発現（DE）で得た遺伝子リストを、既知の遺伝子セット（GO・経路・転写因子）に照らして生物学的に解釈する方法を、3つのアプローチで整理します。

🔙 前の記事：差次的発現（DE）　／　🔜 次の記事：細胞組成の差（differential abundance）

前提：DE 結果（log2FC・補正 p 値・統計量）があること（→差次的発現（DE））。Python の基本。

この記事のゴール

DE 遺伝子を生物学的な言葉（経路・機能・転写因子）に翻訳する3つのアプローチ――ORA／GSEA／活性推定――の違いと使い分け、Python での実行を理解すること。

なぜエンリッチメントが必要か

DE で「変化した遺伝子のリスト」が出ても、それだけでは「生物学的に何が起きたか」は分かりません。そこで、既知の遺伝子セット（GO、KEGG/Reactome 経路、転写因子の標的群など）に照らし、「どの機能・経路がまとまって動いたか」を読み取ります。これが機能エンリッチメントです。

3つのアプローチ ― 入力が違う

エンリッチメントには大きく3系統があり、必要な入力が異なります。

機能エンリッチメントの3アプローチ。(a) ORA は「有意DE遺伝子のリスト＋背景」から過剰代表を検定（GSEApy enrichr）。(b) GSEA は「全遺伝子を統計量で順位づけ」してセットの偏りを見る（GSEApy prerank）。(c) 活性推定は「発現/統計量＋重み付きシグネチャ（PROGENy/CollecTRI）」で経路・TF活性を出す（decoupler、細胞ごとも可）。
機能エンリッチメントの3アプローチ。(a) ORA は「有意DE遺伝子のリスト＋背景」から過剰代表を検定（GSEApy enrichr）。(b) GSEA は「全遺伝子を統計量で順位づけ」してセットの偏りを見る（GSEApy prerank）。(c) 活性推定は「発現/統計量＋重み付きシグネチャ（PROGENy/CollecTRI）」で経路・TF活性を出す（decoupler、細胞ごとも可）。
  • (a) ORA（過剰代表解析）：有意 DE 遺伝子のリストが、ある遺伝子セットに偏って含まれるかをFisher検定。しきい値で区切ったリストが必要。
  • (b) GSEA（遺伝子セット濃縮解析）：全遺伝子を統計量で順位づけし、セットがランキングの上位/下位に偏るかを見る。しきい値が不要。
  • (c) 活性推定（フットプリント）：重み付きの応答シグネチャ（経路は PROGENy、TF は CollecTRI/DoRothEA）を使い、経路・転写因子の活性を推定。単一細胞ごとにも出せる。

(a) ORA ― 遺伝子リストの過剰代表

有意 DE 遺伝子のリストを、GO/KEGG などに照らします。背景（universe）の選び方が結果を左右します。

pythonimport gseapy as gp

# 有意な DE 遺伝子（前記事の res から：上方制御の例）
sig_genes = res[(res["padj"] < 0.05) & (res["log2FoldChange"] > 1.0)].index.tolist()

# 背景＝その細胞型で「発現していた」全遺伝子（重要：全ゲノムではない）
background = res.index.tolist()

enr = gp.enrichr(
    gene_list=sig_genes,
    gene_sets=["GO_Biological_Process_2023", "KEGG_2021_Mouse"],
    background=background,
    outdir=None,
)
enr.results.head()      # Term, Adjusted P-value, Overlap, Genes ...

⚠️ 背景は「発現遺伝子」にする
背景に全ゲノムを使うと、「そもそもこの細胞で発現する遺伝子」の偏りを拾って偽の有意が出ます。背景はその文脈で発現していた遺伝子に揃えます。

(b) GSEA ― 順位づけ全体を使う

しきい値で区切らず、全遺伝子を統計量で順位づけして使います。弱い変化が多数集まる経路を拾えるのが利点です。

pythonimport gseapy as gp

# 全遺伝子を統計量で順位づけ（cutoff 不要）。例：DESeq2 の Wald 統計量 stat
rnk = res[["stat"]].dropna().sort_values("stat", ascending=False)

pre = gp.prerank(
    rnk=rnk,
    gene_sets="MSigDB_Hallmark_2020",
    min_size=15, max_size=500, permutation_num=1000,
)
pre.res2d.head()        # Term, NES, FDR q-val ...

💡 順位づけの量が肝
GSEA は「並び」に意味を持たせる必要があります。方向と強さを反映する量（Wald 統計量や符号付き −log10 p）で並べます。log2FC のみだと分散の大きい遺伝子に引っ張られることがあります。

(c) 活性推定 ― 経路・転写因子の活性

遺伝子セットの「重なり」ではなく、重み付きの応答シグネチャで活性を推定します。単一細胞ごとに経路・TF 活性を出せるため、scRNA-seq と相性が良い方法です。

pythonimport decoupler as dc

# 経路活性：PROGENy（応答遺伝子に重みづけしたシグネチャ）
progeny = dc.get_progeny(organism="mouse", top=500)
dc.run_mlm(mat=adata, net=progeny)
adata.obsm["mlm_estimate"]     # 細胞ごとの経路活性

# 転写因子活性：CollecTRI（高信頼の TF–標的）
collectri = dc.get_collectri(organism="mouse")
dc.run_ulm(mat=adata, net=collectri)
adata.obsm["ulm_estimate"]     # 細胞ごとの TF 活性

🔧 環境メモ：pip install gseapy decoupler。gseapy は Enrichr/MSigDB のライブラリ名を指定します（生物種に注意：上記は Mouse 版）。decoupler の関数名（get_progeny / get_collectri / run_mlm / run_ulm）はバージョンで変わることがあるため、公式チュートリアルに合わせてください。pseudobulk の統計量を入力にすればサンプルレベルの活性も出せます。転写因子活性をさらに掘り下げるなら 遺伝子制御ネットワーク・転写因子活性（GRN／TF activity） も参照。

遺伝子セット資源とツールの早見表

データリソース 内容 主な用途
MSigDB（Hallmark/C2/C5…） 厳選された遺伝子セット（GO/KEGG/Reactome 含む） ORA・GSEA
GO（BP/MF/CC） 遺伝子オントロジー ORA・GSEA
KEGG / Reactome 代謝・シグナル経路 ORA・GSEA
PROGENy 経路の「応答遺伝子」シグネチャ（重みつき） 活性推定（経路）
CollecTRI / DoRothEA 転写因子–標的（重みつき） 活性推定（TF）
ツール 言語 手法
GSEApy Python ORA（enrichr）・GSEA（prerank）
decoupler Python 活性推定（PROGENy/CollecTRI ほか）
g:Profiler Python/Web ORA
fgsea R GSEA

落とし穴

⚠️ つまずきやすい点
– ORA の背景（universe）：全ゲノムではなく発現遺伝子を背景に。ここを誤ると結論が変わる。
– GSEA の順位づけ：方向と強さを反映する量で並べる（Wald 統計量など）。
– 多重比較：多数の遺伝子セットを試すので FDR 補正必須。
– 遺伝子セットの偏り：よく研究された遺伝子・経路に偏る。GO は階層で冗長になりやすい。
– エンリッチメント ≠ 因果：有意でも、その経路が原因とは限らない。

まとめ

  • DE 遺伝子は、既知の遺伝子セットに照らして生物学的に解釈する。
  • 3アプローチ：(a) ORA（リスト＋背景・過剰代表）、(b) GSEA（順位づけ全体・しきい値不要）、(c) 活性推定（重み付きシグネチャで経路・TF 活性、単一細胞ごとも可）。
  • Python では GSEApy（ORA/GSEA）と decoupler（活性推定）が中心。
  • 注意：ORA の背景は発現遺伝子、GSEA は意味のある順位づけ、多重比較補正、マウス/ヒトのオーソログ、エンリッチメント≠因果。

次の記事：細胞組成の差（differential abundance） — 条件間で細胞型の割合が変わるかを、Milo / scCODA で検定する方法に進みます。

関連記事

下流解析の地図と前後

活性推定の発展先

参考文献

  • Subramanian, A., Tamayo, P., Mootha, V. K., et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. PNAS, 102, 15545–15550. doi:10.1073/pnas.0506580102
  • Liberzon, A., Birger, C., Thorvaldsdóttir, H., et al. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Systems, 1, 417–425. doi:10.1016/j.cels.2015.12.004
  • Fang, Z., Liu, X., & Peltz, G. (2023). GSEApy: a comprehensive package for performing gene set enrichment analysis in Python. Bioinformatics, 39, btac757. doi:10.1093/bioinformatics/btac757
  • Kuleshov, M. V., Jones, M. R., Rouillard, A. D., et al. (2016). Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Research, 44, W90–W97. doi:10.1093/nar/gkw377
  • 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
  • Schubert, M., Klinger, B., Klünemann, M., et al. (2018). Perturbation-response genes reveal signaling footprints in cancer gene expression. Nature Communications, 9, 20. doi:10.1038/s41467-017-02391-6
  • Müller-Dott, S., Tsirvouli, E., Vázquez, M., et al. (2023). Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities. Nucleic Acids Research, 51, 10934–10949. doi:10.1093/nar/gkad841

 

コメント

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