- 著者: Tianzhi Wu, Erqiang Hu, Shuangbin Xu, Meijun Chen, Pingfan Guo, Zehan Dai, Tingze Feng, Lang Zhou, Wenli Tang, Li Zhan, Xiaocong Fu, Shanshan Liu, Xiaochen Bo, Guangchuang Yu
- Corresponding author: Xiaochen Bo (Beijing Institute of Radiation Medicine); Guangchuang Yu (Southern Medical University)
- 雑誌: Innovation (Cambridge (Mass.))
- 発行年: 2021
- Epub日: 2021-07-01
- Article種別: Original Article
- PMID: 34557778
背景
大規模オミクスデータ（ゲノム・トランスクリプトーム・プロテオーム・エピゲノム）の解析では、差次発現遺伝子や関心ゲノム領域である ROI (Regions of Interest) から生物学的意味を導くため、遺伝子セットエンリッチメント解析が不可欠のステップである。しかし、従来のツールの多くは、(1) ヒト・マウスなど限られたモデル生物にしか対応せず、(2) GO (Gene Ontology) や KEGG (Kyoto Encyclopedia of Genes and Genomes) など特定のアノテーションデータベースのみをサポートし、(3) 最新データベースへのアップデートが不定期であるという課題を抱えていた。Wadi et al. (2016) の報告では、ツールの約 42% が 5 年以上の古いアノテーションを使用しており、最新アノテーションと比較すると 26% のパスウェイしか検出できていなかったことが示されている。このような古いアノテーションの使用は、生物医学研究において機能的意義を過小評価し、誤解を招く結果につながる可能性がある。
初代の clusterProfiler について報告した Yu et al. (2012) の先行研究や、GSEA の基本アルゴリズムを提唱した Subramanian et al. ProcNatlAcadSciUSA 2005、さらにゲノム領域アノテーションの重要性を示した McLean et al. (2010) などの研究が存在するが、依然として非モデル生物への対応や最新アノテーションの自動同期という点において、既存ツールには大きな「課題」が残されており、多様なオミクスデータを統合的に解釈するための汎用的なプラットフォームが「不足」しているという「ギャップ」が存在した。特に、非モデル生物の機能注釈を迅速に行うためのシステムや、エピゲノムデータとトランスクリプトームデータをシームレスに紐付けるワークフローは「未開拓」であり、研究者が独自のアノテーションデータをインポートして解析する自由度も「不十分」であった。これらの背景から、より広範な生物種と多様なデータタイプに対応し、最新のアノテーションを常に利用できる汎用的なエンリッチメント解析ツールの開発が強く求められていた。
目的
本研究の目的は、Bioconductor パッケージである clusterProfiler をバージョン 4.0 として全面刷新し、以下の主要な機能強化を実現することである。第一に、5,000 種超の生物種への対応を拡大し、モデル生物に限定されない広範な生物学的研究を支援する。第二に、GO・KEGG・WikiPathways・Reactome・MSigDB (Molecular Signatures Database) ・ユーザー定義アノテーションなど、多様な生物学的知識データベースを統一インターフェースで利用可能にする。第三に、ChIP-seq (Chromatin Immunoprecipitation sequencing) や ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) などの非コード領域データから得られるゲノム ROI に対する機能エンリッチメント解析を実装し、エピゲノムデータの解釈能力を向上させる。第四に、複数条件や時系列実験設計における多重比較に対応する機能を提供し、複雑な実験デザインからの生物学的洞察を効率的に抽出する。第五に、dplyr/ggplot2 との完全統合による tidy インターフェースを提供し、データ操作と可視化の効率を大幅に向上させる。
結果
GO解析の強化と多生物種対応:
enrichGO() および gseGO() 関数は、Bioconductor が半年ごとに更新する OrgDb パッケージを通じて、ヒト、マウス、ショウジョウバエ、酵母、線虫など 20 種の主要モデル生物の GO 解析を提供する。さらに、AnnotationHub パッケージを介して UCSC、Ensembl、NCBI、STRING などの各種プロバイダのゲノムワイドアノテーションにオンラインアクセスすることで、事実上 5,000 種超の生物種に GO 解析を拡張した。GO 解析の主要課題であった冗長 GO タームの問題には simplify() 関数で対処する。乳癌差次発現遺伝子 (fold change > 2) の BP-ORA 実施例では、オリジナル結果では細胞分裂・免疫関連の冗長タームが密集してネットワーク全体が不明瞭であったが、simplify() 適用後 (Figure 1B) は複数の機能モジュール (細胞分裂・免疫応答・走化性シグナル等) が明確に分離し、解釈が大幅に向上した。この解析では n=194 の差次発現遺伝子が用いられ、そのうち 34 個が「mitotic nuclear division」に、36 個が「nuclear division」にアノテーションされた (p.adjust 値はそれぞれ 6.70e-23, 1.70e-19) (Figure 1)。
KEGGパスウェイ解析と最新データベース維持: 多くのツールが 2011 年 7 月の KEGG FTP の有料化以降に更新を停止しているのに対し、clusterProfiler は KEGG データをパッケージ内に保持せず、KEGG web API (Application Programming Interface) 経由でリアルタイムに最新データを取得する設計とした。これにより 6,000 種以上の KEGG 収録生物種に対応し、常に最新の経路情報を使用した解析が可能である。GSEA 解析例では乳癌データの上位 5 パスウェイを可視化し (Figure 2A)、UpSet plot による core enriched 遺伝子の重複構造解析では、細胞周期・DNA 複製経路の交差にある MCM (Minichromosome Maintenance) 遺伝子群が腫瘍診断バイオマーカー候補として浮上した (Figure 2B)。MCM 遺伝子群は細胞周期と DNA 複製パスウェイの重複部分に位置し、その発現値は単独のパスウェイに属する遺伝子よりも高いことが示された。IL-17 シグナル伝達経路とプロテアソーム経路の交差は、インターフェロンγ (IFNγ) という 1 つの遺伝子のみに関連しており、IL-17 シグナル伝達経路が炎症反応を誘導し、IFNγ がプロテアソーム形成を調節することが示唆された。
汎用インターフェースによる多データベース・非モデル生物対応:
enricher() および GSEA() の 2 つの汎用関数が、ユーザー提供の GMT (Gene Matrix Transposed) 形式またはデータフレーム形式のアノテーションを受け入れ、WikiPathways・MSigDB・Disease Signatures・CCLE (Cancer Cell Line Encyclopedia) ・InterPro・Clusters of Orthologous Groups・Mouse Phenotype Ontology・ヒト細胞マーカー・COVID-19 関連遺伝子セットなど任意のデータベースに対して ORA および GSEA を実行できる。WikiPathways GSEA の実施例では、卵巣癌細胞 (n=8 replicates) のシスプラチン・オキサリプラチン時系列処理 (0/2/6/24 時間) データの 8 グループ一括解析が compareCluster(formula=Gene~Time+Treatment, fun="GSEA") の 1 行で実行でき、DNA 損傷・細胞周期停止経路が両薬剤共通の後期パスウェイとして同定された (Figure 4)。
ゲノム領域エンリッチメント解析と複数条件比較:
ChIPseeker パッケージの seq2gene() 関数を使用し、ChIP-seq・ATAC-seq 等で得られたゲノムリージョンをホスト遺伝子、近位遺伝子、フランキング遺伝子の 3 カテゴリに対して多対多マッピングし、コーディング遺伝子IDに変換した後に clusterProfiler で機能エンリッチメント解析を実行する。CBX6 の ChIP-seq データ (GSM1295076) 適用例では、ENCODE/ChEA 転写因子遺伝子セットを背景にした ORA により、CBX6 結合領域遺伝子が POU5F1 (OCT4) ・TRIM28 (KAP1) ・SUZ12・EZH2 の標的遺伝子と有意に重複することを確認した (p<0.001) (Figure 3)。この解析では、CBX6 結合領域の遺伝子と POU5F1、TRIM28、SUZ12、EZH2 の標的遺伝子との間に有意な重複 (p<0.001) が認められた。さらに、compareCluster() 関数は R の標準統計モデル式による指定に対応し、時系列×処理条件の複雑な実験デザインを一括処理できる。卵巣癌細胞のシスプラチン・オキサリプラチン時系列処理データ (8 グループ、n=8 replicates) の WikiPathways GSEA 結果は dot plot に X 軸 (時間) ・ファセットパネル (薬剤) で同時表示され、DNA 損傷・細胞周期停止経路が両薬剤共通の後期パスウェイとして同定された (Figure 4)。
tidyインターフェースとggplot2統合:
clusterProfiler は dplyr の動詞を enrichResult/gseaResult/compareClusterResult の S4 オブジェクトに拡張し、データフレームと同様に操作可能にした。as.data.frame() で CSV エクスポートも容易である。ggplot2 の全文法が enrichment result オブジェクトに直接適用可能であり、lollipop chart・bar chart 等のカスタム可視化が数行のコードで実現する。enrichplot パッケージは DOT plot・bar plot・enrichment map network・UpSet plot・gene-concept network・ridgeline plot 等の豊富なパブリケーション品質グラフを提供する (Figure 5)。例えば、GO エンリッチメント結果 (ORA) をリッチファクターで可視化した lollipop chart (Figure 5A) や、GSEA の正規化エンリッチメントスコア (NES) を可視化した bar chart (Figure 5B) が示された。
考察/結論
clusterProfiler 4.0 は、従来ツールが抱えていた「限られた生物種」「古いアノテーション」「特定データベースへの依存」「複雑な実験デザイン非対応」という 4 つの本質的課題を一括解決した。常時最新データベースへのオンラインアクセス設計 (KEGG の web API リアルタイム取得) は、旧アノテーション問題を根本から回避する革新的アプローチであり、今後の解析再現性・比較可能性の向上に貢献する。
先行研究との違い: 本研究は、従来の GREAT (Genomic Regions Enrichment of Annotations Tool) が限られた参照ゲノムのみに対応しエピゲノム ROI 解析が不可能であった点や、Blast2GO が可視化・統合解析機能が限定的であった点と異なり、これらの制約を解消した。また、これまでの clusterProfiler v3 以前は GO・KEGG に限定され、R の tidyverse 哲学との整合性も不十分であったが、バージョン 4.0 ではこれらすべての制約が解除されたことで、単一の解析環境でゲノムから機能注釈・可視化・報告書作成まで完結できるようになった。特に、Subramanian et al. ProcNatlAcadSciUSA 2005 で提唱された GSEA の概念を、多様なデータベースと生物種に適用可能とした点で、既存ツールとは一線を画している。
新規性: 本研究で初めて、ORA と GSEA の統一インターフェースを多様なデータベースと生物種に適用可能とし、dplyr/ggplot2 と完全統合された tidy インターフェースを新規に提供した。これにより、ユーザーが解析目的に応じて手法を自由に選択しながら一貫した作業フローを維持することを可能にし、コードの再利用性を高める。本研究で初めて、これらの機能が統合された形で提供され、多様なオミクスデータ解析における汎用性と効率性を大幅に向上させた。また、非モデル生物に対するゲノム ROI の機能解釈を可能にした点も新規性がある。
臨床応用: EV (Extracellular Vesicles) 関連の生物医学研究との関連では、EV プロテオミクス・EV トランスクリプトーム・EV-miRNA シーケンシングにおける差次発現タンパク質や候補 miRNA のパスウェイ解析・機能注釈として広く活用されることが期待される。CRISPR/Cas9 スクリーンによる EV 生合成機構解析、ChIP-seq・ATAC-seq による EV 産生細胞のエピゲノム解析においても本ツールが標準的解析基盤となり得る。これらの応用は、基礎研究から臨床応用への橋渡しを促進する臨床的意義を持つ。例えば、腫瘍診断バイオマーカー候補として浮上した MCM 遺伝子群の解析は、乳癌の早期診断や治療標的の特定に繋がる可能性がある。
残された課題: バージョン 4.0 以降は API の安定性を優先しており、将来の変更時には少なくとも 1 年間の後方互換性を保証することを表明している。2020 年の年間 2,500 件超の引用実績が示すように、clusterProfiler は生物医学オミクス解析の標準ツールとして今後もさらに広い利用が見込まれる。残された課題としては、非常に大規模なゲノムデータセット (植物の多倍体ゲノム等) での処理速度の最適化、およびネットワーク可視化の大規模遺伝子セットへの対応改善が挙げられる。また、特定の生物学的コンテキストにおけるアノテーションの網羅性や、新しいオミクスデータタイプ（例：シングルセルオミクス）への適用性についても、今後の検討課題である。
方法
clusterProfiler 4.0 は、R パッケージとして Bioconductor エコシステム内で実装された。主要な機能と検証方法は以下の通りである。
エンリッチメント解析: ORA (Over-Representation Analysis) には enrichGO()、enrichKEGG()、enrichMKEGG()、enrichWP()、および汎用インターフェースである enricher() などの関数群を提供する。GSEA (Gene Set Enrichment Analysis) には gseGO()、gseKEGG()、gseMKEGG()、gseWP()、および汎用インターフェースである GSEA() などの関数群を提供する。
比較解析: compareCluster() 関数は、複数遺伝子リストの一括解析を可能にする。この関数は、R の標準統計モデル式 (Gene ~ Time + Treatment など) をサポートし、複雑な実験デザインにおける複数条件間の機能プロファイルの比較を容易にする。卵巣癌細胞のシスプラチン・オキサリプラチン処理時系列データ (GEO: GSE8057) の 8 グループ一括解析に適用された。
データ操作: clusterProfiler は、dplyr パッケージの動詞 (filter, arrange, group_by, mutate, select, summarise など) を enrichResult、gseaResult、compareClusterResult の S4 オブジェクトに拡張適用した。as.data.frame() メソッドにより、結果を CSV (Comma Separated Values) ファイルとして容易にエクスポートできる。また、[[ 演算子を再定義し、ORA 結果からは各タームに注釈付けされた遺伝子リストを、GSEA 結果からは leading edge 遺伝子リストを直接抽出できるようにした。
ゲノム領域エンリッチメント解析: ChIPseeker パッケージの seq2gene() 関数を使用し、ChIP-seq や ATAC-seq などで得られたゲノムリージョンをホスト遺伝子、近位遺伝子、フランキング遺伝子に多対多マッピングし、コーディング遺伝子IDに変換する。その後、clusterProfiler を用いて機能エンリッチメント解析を実行する。
検証データの解析と統計解析: パッケージの動作検証のため、ヒト乳癌由来の MCF-7 細胞株および A549 細胞株から得られた公開トランスクリプトームデータ、および CBX6 の ChIP-seq データ (GSM1295076) を使用した。統計的有意性は、主に p 値および調整済み p 値 (p.adjust) を用いて評価された。GO タームの冗長性除去には、GOSemSim パッケージの意味的類似度計算が用いられ、類似度>0.7 のタームが統合された。乳癌差次発現遺伝子の同定には fold change > 2 を基準とし、統計的検定には Fisher’s exact test、GSEA のスコアリング、および複数群比較における調整 p 値算出のために FDR (False Discovery Rate) 制御法が適用された。また、2群間の発現変動の有意差検定には Student t-test および Wilcoxon rank-sum test が用いられた。