kegg-pathway-analysis

KEGG パスウェイエンリッチメント解析ガイド

概要

KEGG（Kyoto Encyclopedia of Genes and Genomes）パスウェイエンリッチメント解析では、差次的発現遺伝子の中で統計的に過剰に存在する生物学的パスウェイを特定します。このガイドでは、主な2つのエンリッチメント手法（ORA と GSEA）、遺伝子を発現方向別に分けるなどの重要なワークフロー上の判断、R の clusterProfiler と Python の gseapy の選択、そして信頼性に定評のある KEGG REST API への対処法を取り上げます。また、誤ったパスウェイ数を出したり、解析が停止したりする、繰り返し発生する失敗パターンについても説明します。

KEGG パスウェイ解析で最もよくあるエラーは次の3つです。（1）上方制御遺伝子と下方制御遺伝子を1回のエンリッチメント実行にまとめてしまい、実際のパスウェイシグナルが隠れる。（2）KEGG REST API のタイムアウトが発生し、代替手段がないために解析が失敗する。（3）完了しない可能性のある、見た目を整えるためだけのパスウェイ名の検索を試み、結果の報告が遅れる。このガイドでは、それぞれに対する具体的な解決策を紹介します。

主要な概念

ORA と GSEA

過剰代表解析（ORA）と遺伝子セットエンリッチメント解析（GSEA）は、パスウェイエンリッチメントの主な手法です。入力と統計的アプローチの両方が異なります。

ORA では、事前にフィルタリングした遺伝子リスト（例：padj < 0.05 かつ |log2FC| > 1.5 の遺伝子）を使い、KEGG パスウェイの構成遺伝子がバックグラウンド全体と比べてそのリストに過剰に存在するかを検定します。ORA では超幾何検定（Fisher の正確確率検定）を使用します。分かりやすい手法ですが、変化量の情報は失われ、選択した有意性のカットオフに大きく左右されます。

GSEA では、カットオフを設けず、すべての遺伝子を順位付けしたリスト（通常は log2 fold change または符号付き有意性統計量で順位付け）を使用します。順位付けされたリストを順にたどって累積エンリッチメントスコアを計算し、構成遺伝子が順位の上位または下位に集中するパスウェイを特定します。GSEA は、ORA では見逃される可能性のある、微妙な協調的変化を捉えます。

実際には、enrichKEGG()（clusterProfiler）または gp.enrichr()（gseapy）を使った ORA が、より一般的な出発点です。任意のカットオフを避けたい場合や、効果量が小さい場合は、gseKEGG() または gp.prerank() を使った GSEA が適しています。

エンリッチメントにおける発現方向

ORA を行う場合、遺伝子の発現方向（上方制御か下方制御か）は重要です。1つのパスウェイに、逆方向に制御される遺伝子が含まれることがあります。上方制御遺伝子と下方制御遺伝子を1つのリストにまとめると、それぞれの逆方向のシグナルが相殺され、エンリッチメントシグナルが弱まり、実際にエンリッチされているパスウェイが隠れてしまいます。上方制御遺伝子と下方制御遺伝子のセットを分けてエンリッチメントを実行することで、より正確で解釈しやすい結果が得られます。この分割は ORA では必須です。GSEA では符号付きの順位付けによって発現方向が本質的に考慮されますが、先頭のエッジ遺伝子の発現方向を解釈することも重要です。

KEGG の生物種コード

KEGG では、生物種ごとのパスウェイデータベースを識別するために、3文字（または4文字）の生物種コードを使用します。誤ったコードを使うと、エラーが表示されずに空の結果が返されます。一般的なコードは次のとおりです。

生物種	コード
ヒト	hsa
マウス	mmu
ラット	rno
ゼブラフィッシュ	dre
ショウジョウバエ	dme
C. elegans	cel
E. coli K-12	eco
P. aeruginosa PA14	pau
P. aeruginosa PAO1	pae
S. cerevisiae	sce
A. thaliana	ath

遺伝子 ID の形式も生物種によって異なります。真核生物では通常 Entrez gene ID が必要ですが、細菌では locus tag を使用します。ID の種類が一致しないことも、エラーが表示されない失敗の原因になります。

KEGG API の信頼性

KEGG REST API（rest.kegg.jp）にはレート制限があり、頻繁に応答が遅く、タイムアウトも発生しやすい傾向があります。clusterProfiler::enrichKEGG() と、KEGG への直接 HTTP リクエストは、どちらも予測できないタイミングで失敗することがあります。API の失敗に備えることは任意ではなく、KEGG を使うあらゆるワークフローに欠かせません。対策には、パスウェイデータを事前に取得してキャッシュすること、gseapy に同梱されたオフラインの遺伝子セットデータベースを使うこと、タイムアウトを設定した再試行ロジックを実装することなどがあります。

実行前の確認

解析コードを書く前に、ユーザーと以下の点を確認してください。

decisions:
  - id: D1
    param: organismCode
    kind: required
    source: upstream
    ask: "Which species should pathways be looked up for?"
    default: "carried from the dataset"

  - id: D2
    param: identifierMapping
    kind: required
    source: data
    depends_on: [D1]
    ask: "Which identifier type do the genes carry, and how should unmapped ones be handled?"
    default: null

  - id: D3
    param: testingApproach
    kind: required
    source: upstream
    ask: "Test a significant-gene list against pathways (over-representation), or the whole ranked list without a cut (enrichment)?"
    default: "over-representation when the input is a gene list"

  - id: D4
    param: backgroundSet
    kind: required
    source: upstream
    depends_on: [D3]
    ask: "Which genes form the background - everything measured in this experiment, or the whole annotated genome?"
    default: "the genes tested in this experiment"
    skip_if: "ranked-list enrichment, which uses the ranking rather than a background"

  - id: D5
    param: pathwayScope
    kind: optional
    source: user
    ask: "Should all pathway categories be tested, or only metabolism, signalling, or disease maps?"
    default: "all categories"

  - id: D6
    param: significanceThreshold
    kind: required
    source: user
    ask: "How strong must a pathway's evidence be, after correcting for the number of pathways tested?"
    default: "adjusted p-value 0.05"

D4は、最も見落とされやすく、回答に最も大きな影響を与える判断です。 実験で測定したのが発現遺伝子だけなのに、全ゲノムを背景として使うと、あらゆる組織特異的な経路が濃縮されているように見えてしまいます。濃縮の判定対象が、実験で測定できなかった遺伝子になり、実験の対象ではなくなってしまうためです。

判断の枠組み

Question: What enrichment analysis do you need?
|
+-- Have a pre-filtered DEG list (with cutoffs applied)?
|   +-- Yes --> ORA
|   |   +-- Using R? --> clusterProfiler::enrichKEGG()
|   |   +-- Using Python? --> gseapy.enrichr()
|   |   +-- KEGG API failing? --> gseapy with offline gene sets
|   +-- No, want cutoff-free analysis --> GSEA
|       +-- Using R? --> clusterProfiler::gseKEGG()
|       +-- Using Python? --> gseapy.prerank()
|
+-- Need to split by direction?
|   +-- ORA --> YES, always split up/down (mandatory)
|   +-- GSEA --> No split needed (direction encoded in ranking)
|
+-- KEGG API unreliable?
    +-- Try cached/pre-fetched data first
    +-- Fall back to gseapy offline databases
    +-- Use retry logic with short timeouts

シナリオ	推奨されるアプローチ	理由
Rで標準的なORAを実行する	clusterProfiler::enrichKEGG()を方向別に実行する	最も広く使われており、Bioconductorエコシステムと統合されています
Pythonで標準的なORAを実行する	KEGG_2021_Humanを指定してgseapy.enrichr()を使う	オフラインの遺伝子セットにより、APIへの依存を避けられます
カットオフなしで濃縮解析を行う	gseKEGG()またはgp.prerank()によるGSEA	微妙な協調的変化を捉え、恣意的な閾値を必要としません
KEGG APIが停止している	gseapyのオフラインデータベースに切り替える	gseapyにはKEGGの遺伝子セットがローカルにバンドルされています
条件間を比較する	条件ごとに上方・下方の濃縮解析を個別に実行する	条件間で方向を考慮した集合演算が可能になります
モデル生物以外を扱う	生物種コードを確認し、KEGGRESTで利用可能か調べる	コードが間違っていても、結果が空のまま返されることがあります

ベストプラクティス

1. ORAは必ず遺伝子の方向別に実行してください。 上方制御遺伝子と下方制御遺伝子に分け、それぞれに対してenrichKEGG()またはgp.enrichr()を実行してください。両者をまとめると遺伝子リストが膨らみ、濃縮シグナルが弱まり、経路数も不正確になります。両方向で有意となった経路の和集合を報告してください。

2. 背景となる遺伝子集合を明示的に指定してください。 有意な遺伝子だけではなく、検定したすべての遺伝子（差次的発現解析で得られた全遺伝子セット）をユニバースに設定してください。ユニバースを省略すると、KEGGデータベース内の全遺伝子がデフォルトで使われ、よく研究されている経路の有意性が過大評価されます。

3. 濃縮解析を実行する前に、KEGGデータを取得してキャッシュしてください。 解析の開始時に経路と遺伝子の対応データをダウンロードし、ローカルに保存してください。こうすることで、解析途中にKEGG APIが応答しなくなった場合の失敗を避け、解析の再現性も高められます。

4. 経路名を解決する前に、数値の結果を報告してください。 有意な経路IDの数またはリストを算出したら、その結果をすぐに出力してください。追加のKEGG API呼び出しでIDを人が読める名前に変換するのは表記上の処理にすぎず、タイムアウトして主要な結果が失われる可能性があります。

5. 多重検定補正を一貫して適用してください。 未調整のp値ではなく、調整済みp値（通常はBH法で補正されたp.adjust < 0.05）を使用してください。clusterProfilerとgseapyはどちらもデフォルトで補正を適用しますが、カットオフが調整済みの値に適用されていることを必ず確認してください。

6. 遺伝子IDの形式が生物種に合っていることを確認してください。 真核生物のKEGG経路ではEntrez遺伝子IDが必要であり、細菌ではlocus tagが必要です。形式が一致しないと、濃縮された経路が0件でも気づかれないことがあります。入力に遺伝子シンボルを使っている場合は、clusterProfilerのbitr()または同等のID変換を使用してください。

7. clusterProfilerで問題が発生した場合は、代替手段としてgseapyを使用してください。 KEGG APIの問題でenrichKEGG()が失敗した場合は、オフラインの遺伝子セットをバンドルしたgseapyのenrichr()関数（例：KEGG_2021_Human）を使えば、ネットワークに依存せず同等のORA結果を得られます。

よくある落とし穴

1. アップレギュレーション遺伝子とダウンレギュレーション遺伝子を1回のエンリッチメント解析にまとめること。 逆方向に制御された遺伝子を含む経路では効果が相殺され、実際よりも有意な経路が少なくなります。統合したリストから得られる結果は信頼できません。アンチパターンの例：

   # WRONG: combining up and down genes into one list
   all_sig_genes <- rownames(subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5))
   ekegg <- enrichKEGG(gene = all_sig_genes, ...)  # Will miss pathways
   

   • 回避方法: ORA の前に、必ず DEG を制御方向ごとに分けてください。エンリッチメント解析を2回（アップレギュレーション用とダウンレギュレーション用に1回ずつ）実行し、有意な経路を統合します。

2. 誤った KEGG 生物種コードを使用すること。 KEGG は無効なコードや不一致のある生物種コードに対して、エラーを表示せずに空の結果を返します。これは、株ごとに異なるコードが複数ある細菌種（例：PAO1 の pae と PA14 の pau）で特によく起こります。

   • 回避方法: エンリッチメント解析を実行する前に、KEGG の生物種リストでコードを確認してください。細菌の場合は、特定の株のコードを確認します。

3. 遺伝子 ID の種類が一致していないこと。 KEGG が Entrez ID（細菌の場合は遺伝子座タグ）を想定しているのに遺伝子シンボルを渡すと、エラーメッセージなしでエンリッチメントされた経路が0件になります。

   • 回避方法: 対象生物種で想定される ID の種類を確認してください。エンリッチメント解析の前に、clusterProfiler::bitr() または同等の方法を使って遺伝子シンボルを Entrez ID に変換します。

4. KEGG API のタイムアウトに対処しないこと。 KEGG REST API は頻繁にタイムアウトし、enrichKEGG() が解析の途中で失敗します。エラー処理がないと、解析全体が失われます。

   • 回避方法: KEGG API の呼び出しを、短いタイムアウト（30 seconds）を設定した再試行ロジックで囲んでください。解析の開始時に経路データを事前取得します。代替手段として gseapy を用意してください。

5. 経路名を解決するために回答を遅らせること。 結果の計算後に keggGet() を呼び出して経路 ID を人間が読みやすい名前に変換すると、タイムアウトによって数値の結果全体が失われる可能性があります。アンチパターンの例：

   # BAD: answer delayed by name lookup that may hang
   result_ids <- setdiff(pathways_condA, pathways_condB)
   names <- keggGet(result_ids)  # Can timeout -- answer never emitted
   count <- length(result_ids)
   

   • 回避方法: 経路 ID と件数をすぐに報告してください。主要な結果を確保した後で、tryCatch() または try/except ブロック内でのみ名前を解決します。

6. バックグラウンドユニバースを省略すること。 ユニバースパラメーターを指定しない場合、対象生物種の KEGG 遺伝子データベース全体がデフォルトで使用され、十分に注釈されたハウスキーピング遺伝子を含む経路の統計的有意性が過大評価されます。

   • 回避方法: enrichKEGG() に universe = rownames(res)（DE 解析で検定したすべての遺伝子）を必ず渡してください。

7. フィルタリングに調整済み p 値ではなく、生の p 値を使用すること。 多重検定補正なしで p < 0.05 の経路を報告すると、偽陽性が大幅に増加します。

   • 回避方法: 必ず p.adjust < 0.05 でフィルタリングしてください。フィルタリング対象の列が、生の p 値ではなく調整済みの値であることを確認します。

ワークフロー

1. ステップ 1：遺伝子リストを準備する

   • 差次的発現解析の結果から、有意な DEG をフィルタリングします（例：padj < 0.05、|log2FC| > 1.5）。
   • 上方制御遺伝子セットと下方制御遺伝子セットに分けます。
   • 遺伝子 ID を KEGG が想定する形式（Entrez ID または遺伝子座タグ）に変換します。

   library(clusterProfiler)
   
   # Filter significant DEGs
   sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5)
   
   # MANDATORY: Split by direction BEFORE running enrichment
   up_genes <- rownames(subset(sig_genes, log2FoldChange > 0))
   dn_genes <- rownames(subset(sig_genes, log2FoldChange < 0))
   
   cat("Up-regulated genes:", length(up_genes), "\n")
   cat("Down-regulated genes:", length(dn_genes), "\n")
   

2. ステップ 2：KEGG データを事前に取得する（推奨）

   • エンリッチメントを実行する前に、経路と遺伝子の対応関係をローカルにキャッシュします。
   • 判断ポイント：KEGG API にアクセスできる場合は、ライブクエリを実行します。アクセスできない場合は、オフラインの遺伝子セットに切り替えます（ステップ 2b）。

   R（KEGGREST パッケージ）：

   library(KEGGREST)
   pathway_list <- tryCatch(
     keggList("pathway", organism_code),
     error = function(e) NULL
   )
   

   Python（キャッシュを使用した requests）：

   import requests
   import json
   import os
   
   cache_file = f"kegg_{organism_code}_pathways.json"
   if os.path.exists(cache_file):
       with open(cache_file) as f:
           pathway_map = json.load(f)
   else:
       resp = requests.get(
           f"https://rest.kegg.jp/list/pathway/{organism_code}", timeout=30
       )
       if resp.ok:
           pathway_map = dict(
               line.split("\t") for line in resp.text.strip().split("\n")
           )
           with open(cache_file, 'w') as f:
               json.dump(pathway_map, f)
   

3. ステップ 3：それぞれの方向について個別にエンリッチメントを実行する

   • ORA を 2 回実行します。1 回は上方制御遺伝子、もう 1 回は下方制御遺伝子に対して実行します。
   • 常にユニバース（検定対象となったすべての遺伝子）を指定します。

   R（clusterProfiler）：

   ekegg_up <- enrichKEGG(gene = up_genes, organism = organism_code,
                          universe = rownames(res), pvalueCutoff = 0.05)
   ekegg_dn <- enrichKEGG(gene = dn_genes, organism = organism_code,
                          universe = rownames(res), pvalueCutoff = 0.05)
   
   up_pathways <- subset(as.data.frame(ekegg_up), p.adjust < 0.05)$ID
   dn_pathways <- subset(as.data.frame(ekegg_dn), p.adjust < 0.05)$ID
   
   all_sig_pathways <- union(up_pathways, dn_pathways)
   cat("Significant pathways (up):", length(up_pathways), "\n")
   cat("Significant pathways (down):", length(dn_pathways), "\n")
   cat("Total unique significant pathways:", length(all_sig_pathways), "\n")
   

   Python（gseapy）：

   import gseapy as gp
   
   up_genes = sig_genes[sig_genes['log2FoldChange'] > 0].index.tolist()
   dn_genes = sig_genes[sig_genes['log2FoldChange'] < 0].index.tolist()
   
   enr_up = gp.enrichr(gene_list=up_genes, gene_sets='KEGG_2021_Human',
                        organism='human', outdir=None)
   enr_dn = gp.enrichr(gene_list=dn_genes, gene_sets='KEGG_2021_Human',
                        organism='human', outdir=None)
   

4. ステップ 4：結果をすぐに報告する

   • 名前の解決を試みる前に、経路 ID と件数を出力します。
   • 必要に応じて、処理をブロックしない方法で経路名を解決します。

   result_ids <- setdiff(pathways_condA, pathways_condB)
   count <- length(result_ids)
   cat("Answer:", count, "pathways unique to condition A\n")
   cat("Pathway IDs:", paste(result_ids, collapse = ", "), "\n")
   
   # Optional: resolve names (non-blocking)
   names <- tryCatch(keggGet(result_ids), error = function(e) NULL)
   

5. ステップ 5：条件間の比較（該当する場合）

   • 条件ごとに上方制御と下方制御のエンリッチメントを個別に実行します（enrichKEGG の呼び出しは合計 4 回）。
   • 同じ方向の経路セットを比較します。
   • 方向ごとに固有の経路を統合した結果を報告します。

   # Condition A
   up_pathways_A <- subset(as.data.frame(ekegg_up_A), p.adjust < 0.05)$ID
   dn_pathways_A <- subset(as.data.frame(ekegg_dn_A), p.adjust < 0.05)$ID
   
   # Condition B
   up_pathways_B <- subset(as.data.frame(ekegg_up_B), p.adjust < 0.05)$ID
   dn_pathways_B <- subset(as.data.frame(ekegg_dn_B), p.adjust < 0.05)$ID
   
   # Find pathways unique to condition A in each direction
   unique_up <- setdiff(up_pathways_A, up_pathways_B)
   unique_dn <- setdiff(dn_pathways_A, dn_pathways_B)
   
   # Union
   unique_to_A <- union(unique_up, unique_dn)
   cat("Pathways in A but not B:", length(unique_to_A), "\n")
   cat("  From up-regulated:", length(unique_up), "-",
       paste(unique_up, collapse = ", "), "\n")
   cat("  From down-regulated:", length(unique_dn), "-",
       paste(unique_dn, collapse = ", "), "\n")
   

6. ステップ 6：再試行ロジックで API の失敗に対処する

   • KEGG API の呼び出しに失敗した場合は、オフラインデータに切り替える前に、タイムアウトを設定して再試行します。

   fetch_kegg_with_retry <- function(organism_code, max_retries = 3,
                                      timeout_sec = 30) {
     for (i in seq_len(max_retries)) {
       result <- tryCatch({
         R.utils::withTimeout(
           keggList("pathway", organism_code),
           timeout = timeout_sec
         )
       }, error = function(e) NULL)
       if (!is.null(result)) return(result)
       Sys.sleep(2)
     }
     warning("KEGG API unreachable after retries. Proceeding without pathway names.")
     return(NULL)
   }
   

参考資料

• KEGG：京都遺伝子ゲノム百科事典 -- 経路マップ、対象生物の一覧、REST API のドキュメントを掲載した公式 KEGG データベース
• clusterProfiler：遺伝子クラスター間の生物学的テーマを比較するための R パッケージ -- enrichKEGG 関数と gseKEGG 関数について解説した Guangchuang Yu による原著論文
• gseapy ドキュメント -- オフラインの KEGG 遺伝子セットに対応した、遺伝子セットエンリッチメント解析用の Python パッケージ
• KEGG REST API リファレンス -- KEGG の経路データにプログラムからアクセスするための技術ドキュメント

関連スキル

• gseapy-gene-enrichment -- Python ベースの遺伝子セットエンリッチメント解析です。clusterProfiler の KEGG API 呼び出しに失敗した場合の代替手段として、または Python ベースのワークフローの主要ツールとして使用します。
• deseq2-differential-expression / pydeseq2-differential-expression -- KEGG 経路エンリッチメントの入力に使用する DEG リストを生成する、上流の差次的発現解析です。