pathway-enrichment

パスウェイエンリッチメント

概要

エンリッチメント解析では、「遺伝子の中で、どの生物学的機能が過剰に存在するか？」を調べます。差次的発現解析、スクリーニング、クラスタリングの後に行う標準的な最後のステップです。主な手法は2つあり、適切な手法を選ぶことが最も重要です。

• ORA（過剰代表解析） — 閾値を設定した遺伝子リスト（例：padj < 0.05）を取り、遺伝子セットとの重複が偶然より多いかを、Fisherの正確確率検定／超幾何検定で調べます。ツール：Enrichr、g:Profiler。
• GSEA（遺伝子セットエンリッチメント解析） — 遺伝子の全順位リスト（閾値なし）を取り、各遺伝子セットが上位または下位に集中しているかを調べます。Preranked GSEAでは、遺伝子ごとのスコア（例：DESeq2の stat）を使います。影響が広範囲に及ぶものの微妙な場合に適しています。

このスキルでは、これらの解析、その基盤となる遺伝子セットデータベース、そして結果を誤らせたり、論文に掲載できないものにしたりする解釈上の落とし穴を一元的に扱います。

このスキルを使う場面

次のような目的がある場合に、このスキルを使います。

• 遺伝子リストから、エンリッチされたGO用語／KEGG／Reactome／WikiPathways／MSigDB Hallmarkセットを見つける。
• DESeq2、edgeR、limma、またはScanpy rank_genes_groups の出力に対して、GSEA／preranked GSEAを実行する。
• サンプル／細胞ごとにパスウェイ活性をスコア化する（ssGSEA、GSVA）。
• エンリッチメント結果を解釈、重複除去、可視化する、または論文掲載用の表／図を作成する。
• ORAとGSEAのどちらを使うか、遺伝子セットライブラリを選ぶか、背景を選ぶか、遺伝子IDの問題を修正するかを判断する。

Enrichrを手早く単発で検索する場合は、gget スキル（gget enrichr）のほうが軽量です。生のパスウェイ／相互作用API（Reactome、KEGG、STRING）については、database-lookup スキルを参照してください。完全で信頼性の高いエンリッチメントワークフローには、このスキルを使います。

適切な手法の選択

状況	手法	ツール／エントリーポイント
個別のヒットリスト（DE遺伝子、スクリーニングヒット、クラスターマーカー）がある	ORA	gp.enrichr(...) またはg:Profiler
スコア付きの完全な順位リスト（検定した全遺伝子＋スコア）がある	Preranked GSEA	gp.prerank(...)
発現行列＋クラスラベルがある	GSEA	gp.gsea(...)
サンプル／細胞ごとのパスウェイスコアが必要	ssGSEA / GSVA	gp.ssgsea(...)、gp.gsva(...)
カスタム背景、または500以上の生物種が必要	カスタムドメインを使ったORA	g:Profiler（domain_scope='custom'）
TF／シグナル伝達の活性（PROGENy、DoRothEA）を調べたい	活性推定	references/databases-and-gene-sets.md（decoupler）を参照

迷った場合は、閾値を設定したリストならORA、スコア付きの順位テーブルならGSEAを使います。リストに閾値を設定してからGSEAに渡してはいけません。GSEAが必要とする順位情報が失われるためです。

セットアップ

uv pip install gseapy gprofiler-official
# gseapy pulls pandas, numpy, scipy, matplotlib. Network access is needed for
# Enrichr, g:Profiler, and MSigDB downloads. For fully offline ORA, use a local
# GMT file with gp.enrich() (see references/gseapy.md).

利用可能な遺伝子セットライブラリを確認して一覧表示します（名称は時間とともに変わるため、決してむやみにハードコードしないでください）。

import gseapy as gp
names = gp.get_library_name(organism="human")   # 200+ Enrichr libraries
print([n for n in names if "Reactome" in n or "KEGG" in n or "Hallmark" in n])

クイックスタート

ヒットリストに対するORA（gseapy + Enrichr）

import gseapy as gp

# Enrichr libraries expect HGNC gene SYMBOLS (human: UPPERCASE). Map IDs first if needed.
genes = [g.strip() for g in open("deg_symbols.txt") if g.strip()]

enr = gp.enrichr(
    gene_list=genes,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2023",
               "KEGG_2021_Human", "Reactome_2022"],
    organism="human",
    outdir=None,            # in-memory; set a path to also write tables/plots
)
res = enr.results
sig = res[res["Adjusted P-value"] < 0.05].sort_values("Adjusted P-value")
print(sig[["Gene_set", "Term", "Overlap", "Adjusted P-value", "Combined Score", "Genes"]].head(20))

DESeq2の結果からPreranked GSEAを実行する

import gseapy as gp
import pandas as pd

res = pd.read_csv("deseq2_results.csv", index_col=0)   # index = gene symbols
# Rank by the test statistic (sign = direction, magnitude = evidence). This is
# more stable than ranking by log2FoldChange, which is noisy for low-count genes.
rnk = res["stat"].dropna().sort_values(ascending=False)
rnk.index = rnk.index.str.upper()
rnk = rnk[~rnk.index.duplicated(keep="first")]

pre = gp.prerank(
    rnk=rnk,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2023"],
    min_size=15, max_size=500,        # drop tiny/huge sets (noisy or generic)
    permutation_num=1000, seed=123,   # seed = reproducible p-values
    threads=4, outdir=None,
)
out = pre.res2d.sort_values("FDR q-val")
print(out[["Term", "ES", "NES", "NOM p-val", "FDR q-val", "Lead_genes"]].head(20))

stat 列がない場合は、sign(log2FoldChange) * -log10(pvalue) から順位を作成します。

基本的なワークフロー

信頼性の高い解析を行うには、次の手順に沿って進めます。途中の手順（IDの種類、背景）で、結果が気づかないうちに誤ることが最もよくあります。

ステップ1 — 入力を明確にし、手法を選ぶ

対象の遺伝子、生物種、遺伝子ごとのスコアがあるか（→ GSEA）、それともリストだけか（→ ORA）、またどの比較を表すか（解釈では方向が重要）を確認します。

ステップ2 — 遺伝子IDを適切な名前空間に変換する

Enrichr／MSigDBライブラリは遺伝子シンボル（ヒトは大文字、マウスは先頭のみ大文字）をキーにしています。Ensembl／Entrez IDを持っている場合は、先に変換します。gp.Biomart、g:Profilerの g:Convert、および mygene については、references/databases-and-gene-sets.md を参照してください。IDの不一致に気づかないことが、「有意な結果が何もない」最大の原因です。

ステップ 3 — 質問に合った遺伝子セットライブラリを選択する

Hallmark（幅広いテーマ）→ GO:BP（メカニズム）→ KEGG/Reactome/WikiPathways（キュレーション済み経路）→ C7（免疫）などを選びます。50 個のライブラリを実行するのではなく、生物学的なテーマに合った 2～4 個を選びましょう。カタログと選択のガイダンスについては、references/databases-and-gene-sets.md を参照してください。

ステップ 4 — 背景のユニバースを設定する（ORA のみ）

背景には、ゲノム全体ではなく、アッセイで検出される可能性があった遺伝子（例：発現または検査されたすべての遺伝子）を指定する必要があります。誤った背景を使うと、有意性が過大評価されます。Enrichr は固定の背景を使用します。背景が重要な場合は、domain_scope='custom' と background を指定した g:Profiler、または明示的な背景を指定した gp.enrich() を使用してください。根拠については、references/interpretation.md を参照してください。

ステップ 5 — 解析を実行する

クイックスタートのパターン、または同梱の scripts/run_enrichment.py を使用します。GSEA では必ず seed を設定し、permutation_num を報告してください。

ステップ 6 — 調整済み p 値でフィルタリングする

生の p 値ではなく、Adjusted P-value（ORA、Benjamini–Hochberg）または FDR q-val（GSEA）を使用します。一般的なカットオフは 0.05 です。また、「ヒット」が 2000 遺伝子のセットに含まれる 1 遺伝子だけではないことを確認するため、重複数／遺伝子数もチェックしてください。

ステップ 7 — 可視化する

ドットプロット、棒グラフ、エンリッチメントマップ、GSEA のランニングスコアプロットは gseapy に組み込まれています（gp.dotplot、gp.barplot、gp.enrichment_map、gp.gseaplot）。references/gseapy.md を参照してください。

ステップ 8 — 冗長性を減らして解釈する

GO では特に、非常によく似た用語が多数返されます。エンリッチメントマップ（用語間の類似性）、リーディングエッジの重複、または親用語を使ってまとめ、代表的な用語を報告してください。解釈の枠組みと論文用テーブルの形式については、references/interpretation.md を参照してください。

ヘルパースクリプト

scripts/run_enrichment.py は ORA または GSEA をエンドツーエンドで実行し、定型処理（シンボルのクリーンアップ、重複除去、NA の削除、DESeq2 テーブルからの順位付けの作成、ライブラリごとの FDR フィルタリング）を行い、結果テーブルとドットプロットを書き出します。

# ORA from a hit list (one gene symbol per line)
python scripts/run_enrichment.py ora \
  --genes deg_symbols.txt \
  --libraries MSigDB_Hallmark_2020 GO_Biological_Process_2023 KEGG_2021_Human \
  --organism human --outdir results/

# Preranked GSEA from a DESeq2 results CSV (auto-builds the rank from `stat`)
python scripts/run_enrichment.py gsea \
  --deseq2 deseq2_results.csv \
  --libraries MSigDB_Hallmark_2020 GO_Biological_Process_2023 \
  --organism human --outdir results/ --seed 123

# Preranked GSEA from an explicit 2-column rank file (gene,score)
python scripts/run_enrichment.py gsea --rnk ranked_genes.csv --outdir results/

すべてのオプション（背景ファイル、FDR カットオフ、セットサイズの最小値／最大値、置換回数）を確認するには、python scripts/run_enrichment.py --help を実行してください。

よくある落とし穴

以下は、誤った結果や再現性のない結果の主な原因です。

1. 遺伝子 ID／生物種の不一致 — シンボルと Ensembl、人間とマウスの大文字・小文字の違いなどです。ID を対応付けて organism を正しく設定してください。設定を誤ると、一致する遺伝子が気付かないうちにほぼゼロまで減少します。
2. 誤った背景（ORA） — 検査済み／発現遺伝子セットではなくゲノム全体を使用すると、p 値が過大評価されます。背景が重要な場合は、カスタム背景を設定してください。
3. GSEA の前にしきい値を適用する — GSEA には完全な順位付きリストが必要です。カット済みリストを使うのは ORA の場合のみです。
4. GSEA の順位付けに log2FoldChange のみを使う — 低カウント遺伝子では不安定です。stat または sign(LFC) * -log10(p) を優先してください。
5. ライブラリ間の多重検定 — FDR はライブラリ内で計算されます。多くのライブラリを実行すると、検定数が増えます。ライブラリごとの FDR を報告し、慎重に解釈してください。
6. 冗長な GO 用語 — 同じ用語のバリエーションを 40 件報告せず、まとめて代表的な用語を示してください。
7. 有意性 ≠ 関連性 — 重複数と遺伝子セットのサイズを確認してください。小さなセットは簡単に有意になります。
8. ORA に対して短すぎる／長すぎるリスト — 10 遺伝子未満では検出力が不足し、2000 遺伝子を超えると特異性が失われます（代わりに GSEA を検討してください）。
9. 再現性に必要なメタデータがない — Enrichr／GO ライブラリにはバージョンがあり、時間とともに変化します。ライブラリ名と日付を記録し、GSEA の seed を設定してください。

他のスキルとの連携

• 上流（遺伝子の取得元）： pydeseq2（DE 遺伝子と GSEA 用の stat）、scanpy（rank_genes_groups のマーカー／スコア）、depmap／pytdc（スクリーニングのヒット）、プロテオミクスのスキル（pyopenms、matchms）。
• データベース／ID： database-lookup（Reactome、KEGG、STRING、Gene Ontology API）、gget（クイックパスの gget enrichr、ID 対応付け用の gget info）、bioservices。
• 下流： scientific-visualization（カスタム図）、networkx（エンリッチメントマップのグラフ）、scientific-writing／literature-review（解釈と引用）、statistical-analysis（多重検定の詳細）。

参考ファイル

詳しい情報が必要な場合は、該当するファイルをお読みください。

• references/gseapy.md — gseapy API の全容：enrichr、オフラインの enrich、prerank、gsea、ssgsea、gsva、Msigdb、Biomart、get_library_name／read_gmt、すべてのプロット、結果列の意味、GMT／オフラインでの使用方法、トラブルシューティング（レート制限、空の結果）。
• references/databases-and-gene-sets.md — GO、KEGG、Reactome、WikiPathways、MSigDB コレクション、Enrichr ライブラリの命名、g:Profiler のソース、生物種の扱い、遺伝子 ID 変換、質問に応じたライブラリの選択、Reactome／STRING API と decoupler による活性推論へのリンク。
• references/interpretation.md — ORA と GSEA の統計、背景ユニバースの選び方、多重検定法（BH と g:SCS、Bonferroni）の比較、先行遺伝子、冗長性の削減、効果と有意性、論文用テーブルのテンプレート、再現性チェックリスト。

リソース

• gseapy ドキュメント：https://gseapy.readthedocs.io/ · リポジトリ：https://github.com/zqfang/GSEApy
• g:Profiler：https://biit.cs.ut.ee/gprofiler/ · Python クライアント：https://pypi.org/project/gprofiler-official/
• Enrichr：https://maayanlab.cloud/Enrichr/ · MSigDB：https://www.gsea-msigdb.org/gsea/msigdb/
• GSEA メソッド：Subramanian ら（2005）PNAS、DOI: 10.1073/pnas.0506580102

関連スキル

• bulk-rnaseq：バルク RNA-seq のエンドツーエンドオーケストレーター — 生の FASTQ リードを受け取り、QC とトリミング（FastQC、fastp/Trim Galore）、アラインメント、定量化まで実行します…
• bioservices：40 以上のバイオインフォマティクスサービスに対応する統合 Python インターフェース。
• gget：gget CLI と Python パッケージを通じて、20 以上のバイオインフォマティクスデータベースおよび解析サービスを照会します。Ensembl の遺伝子検索、情報取得などに対応します…