- 著者: Aravind Subramanian, Pablo Tamayo, Vamsi K. Mootha, Sayan Mukherjee, Benjamin L. Ebert, Michael A. Gillette, Amanda Paulovich, Scott L. Pomeroy, Todd R. Golub, Eric S. Lander, Jill P. Mesirov
- Corresponding author: Eric S. Lander (Broad Institute of MIT and Harvard); Jill P. Mesirov (Broad Institute of MIT and Harvard)
- 雑誌: Proceedings of the National Academy of Sciences
- 発行年: 2005
- Epub日: 2005-09-30
- Article種別: Original Article (Methods)
- PMID: 16199517
背景
DNA マイクロアレイ技術の登場（Schena et al, Lockhart et al）により、ゲノムワイドな発現解析が日常的に用いられるようになった。しかし、得られた膨大なデータから生物学的な解釈を導き出すプロセスには大きな課題が残されていた。従来の解析では、2つの表現型クラス間で各遺伝子の発現変動を統計的に評価し、差次発現遺伝子を順位付けして上位や下位の少数の遺伝子を抽出するアプローチが一般的であった。しかし、この個別遺伝子アプローチには4つの重大な限界が存在した。第一に、多重検定補正（Benjamini et al）を適用した結果、ノイズに対して個々の遺伝子の発現変動が微小である場合、有意差を示す単一遺伝子が1つも残らないという問題があった。第二に、逆に多数の有意遺伝子が検出された場合でも、その長大なリストから統一的な生物学的テーマを見出すことが困難であり、解釈が研究者の主観や専門分野に依存して恣意的になる課題があった。第三に、細胞内のプロセスは複数の遺伝子が協調して機能するため、個々の遺伝子では検出できない経路レベルでの僅少だが一貫した変動を見落とすという限界があった。第四に、同一の生物学的システムを異なる研究グループが独立して解析した際、得られる有意遺伝子リスト同士の重複が極めて少ないという研究間再現性の低さが問題となっていた。
これらの課題を解決するため、Mootha et al は糖尿病患者の筋生検データ解析において、既知の生物学的知見に基づく遺伝子セット単位で発現変動を評価する予備的なアプローチを適用し、酸化的リン酸化（OXPHOS）関連遺伝子群が協調して約 20% 低下していることを示した。しかし、この初期の試みは特定のデータセットに特化したものであり、多様なゲノムデータに適用可能な汎用的かつ頑健な統計的手法としては未確立であった。多重検定補正や遺伝子相関を考慮した標準的なアルゴリズムが決定的に不足しており、その有効性については未解明な部分が多く、解析プロセスのギャップが残されていた。したがって、ゲノムワイドな発現プロファイルを「遺伝子セット単位」で頑健かつ客観的に解釈するための汎用的な方法論の確立が強く求められていた。
目的
本研究の目的は、事前に定義された生物学的知識（代謝経路、シグナル伝達経路、染色体上の位置、共発現クラスター、転写因子標的モチーフなど）に基づく「遺伝子セット」を対象とし、それらが表現型と相関して順位付けされた遺伝子リスト L の上位または下位に偏在（enrich）しているかを統計的に評価する汎用アルゴリズム「GSEA (Gene Set Enrichment Analysis) 法」を確立することである。具体的には、(1) 遺伝子発現の相関度で重み付けした頑健な ES (Enrichment Score) の定義、(2) 遺伝子間の複雑な相関構造を保存した表現型ラベルの並び替え（permutation）による経験的有意性検定の導入、(3) 複数遺伝子セットの同時解析に対応した NES (Normalized Enrichment Score) および FDR (False Discovery Rate) による多重検定補正の統合、(4) 1325 個の定義済み遺伝子セットを収録した MSigDB (Molecular Signatures Database) v1.0 の構築、およびこれらを実装した GSEA-P (Gene Set Enrichment Analysis-Pathway) ソフトウェアの無償公開を通じて、ゲノムワイド発現データを知識ベースで解釈する標準的な解析枠組みを提供することを目指した。
結果
性染色体遺伝子セットを用いた検出感度と特異性の検証: 男性 n=15 samples および女性 n=17 samples のリンパ芽球細胞株から得られた発現プロファイルを用いて、GSEA の基本性能を検証した。染色体 X 不活化を逃れる遺伝子セット（chrX inactive）は、女性側で極めて有意に濃縮（enrich）されていることが示された（nominal p<0.001、(Table 1)）。これは、等重みを用いた旧バージョンの GSEA アルゴリズムで得られた p=0.007 という結果と比較して、検出感度が大幅に向上していることを示している。一方、男性側では Y 染色体上の遺伝子セット（chrY、chrYp11、chrYq11）がすべて p<0.001 で有意に濃縮されていた。これに対し、性染色体とは無関係な対照セットであるビタミン C 輸送経路（vitcb pathway）は p=0.38、T ヘルパー細胞ケモカイン受容体経路（nkt pathway）は p=0.54 となり、有意な濃縮を示さず、本手法の高い特異性が実証された（(Fig 2)）。
NCI-60がん細胞株におけるp53変異ステータスとシグナル伝達経路の相関: NCI-60 (National Cancer Institute 60) がん細胞株のうち、p53 野生型 n=17 cell lines と p53 変異型 n=33 cell lines の発現データを比較解析した。p53 野生型群においては、p53 シグナル伝達経路（p53 signaling pathway）、心血管系における低酸素および p53 応答、放射線感受性遺伝子セットなど、p53 の正常機能に直結する5つの遺伝子セットが有意に濃縮された（すべて nominal p<0.013、FDR<0.25、(Table 2)）。一方、p53 変異型群では、Ras シグナル伝達経路（Ras signaling pathway）が FDR=0.171 で有意に濃縮された。さらに、Ngf 経路および Igf1 (Insulin-like growth factor 1) 経路の leading-edge subset を解析したところ、MAP2K1、RAF1、ELK1 (ETS transcription factor ELK1)、PIK3CA の4つのコア遺伝子が共通して抽出され、p53 欠損下における MAPK (Mitogen-Activated Protein Kinase) 経路の活性化という生物学的に整合する分子メカニズムが新規に同定された（(Fig 3)）。
急性白血病における染色体欠失領域および細胞起源の同定: 急性リンパ性白血病である ALL (Acute Lymphoid Leukemia) n=24 patients と急性骨髄性白血病である AML (Acute Myeloid Leukemia) n=24 patients のデータセットを用いて、染色体位置に基づく遺伝子セット（C1）の濃縮度を評価した。AML 群においては、AML で頻繁に欠失することが知られている染色体領域 5q31（FDR=0.046）や、RB (Retinoblastoma) 遺伝子座を含む 13q14（FDR=0.057）、17q23（FDR=0.071）、6q21（FDR=0.011）が有意に濃縮された（(Table 2)）。また、ALL 群においては、免疫グロブリン重鎖遺伝子座が存在する 14q32（FDR=0.082）が有意に濃縮されており、これは染色体異常ではなく B 細胞系統の細胞起源を反映した組織特異的発現を捉えたものである。このように、GSEA は単一の遺伝子マーカーの変動を超えて、疾患の遺伝学的特徴や細胞起源を頑健に再現できることが示された。
糖尿病筋生検における酸化的リン酸化遺伝子群の協調的発現低下: 糖尿病患者と健常対照群の骨格筋生検サンプルの比較において、酸化的リン酸化（OXPHOS）関連遺伝子セットの挙動を検証した。個々の遺伝子の平均発現低下率は約 20%（0.8-fold 相当、すなわち 1.25-fold decrease）と極めて微小であり、従来の単一遺伝子解析では多重検定補正後に有意な遺伝子が全く検出されない状況であった。しかし、GSEA を適用したところ、OXPHOS 遺伝子セット全体が糖尿病群において一貫して下方調整されていることが、Petersen et al や Patti et al の知見とも整合し、nominal p=0.008 および FDR=0.04 という高い統計的有意性をもって検出された。この結果は、個々の変動は小さくとも、経路全体として協調的に機能する遺伝子群の挙動を GSEA が感度良く捉えられることを示す決定的な実証例となった。
肺がん予後予測データセット間における共通パスウェイの同定と再現性解決: ボストン研究（n=62 patients、Bhattacharjee et al）とミシガン研究（n=86 patients、Beer et al）の2つの独立した肺腺がん予後データセットを用いて、生存期間（良好 vs 不良）に関連する発現プロファイルを比較した。従来の単一遺伝子解析では、各研究の予後不良群における上位100個の遺伝子リスト同士の重複はわずか 12%（12遺伝子）にとどまり、統計的有意性は極めて低かった。しかし、GSEA を用いてボストン研究の予後不良遺伝子セット（S_Boston）をミシガン研究の全遺伝子リストに対して評価したところ、NES=1.90、p<0.001 という極めて強い濃縮が確認された。同様に、S_Michigan もボストン研究において NES=2.13、p<0.001 で有意に濃縮された。さらに、C2 パスウェイ解析により、両研究において低酸素応答（HIF1A、VEGF 等）やアミノ酸代謝、mTOR シグナル伝達、テロメラーゼ活性化といった共通の悪性度関連パスウェイ（例えば、log2FC 1.8 や 1.5-fold increase 程度の変動を示す遺伝子群）が FDR<0.25 で頑健に検出され、研究間の再現性問題が経路レベルで解消されることが示された（(Table 2)）。
考察/結論
GSEA は、従来の個別遺伝子解析が抱えていた4つの限界（多重検定の障壁、解釈の恣意性、経路レベル変動の見逃し、研究間再現性の低さ）を一括で克服する画期的な解析フレームワークである。
先行研究との違い: 従来の個別遺伝子解析手法や、カットオフ値を用いた単純な重複統計量（累積超幾何分布など）に基づくパスウェイ解析ツールと異なり、GSEA は全遺伝子リスト L をカットオフなしで走査し、発現相関度で重み付けした enrichment score を算出する。これにより、個々の遺伝子レベルでは微小な発現変動であっても、経路全体として協調して機能するシグナルを感度良く検出することが可能となった。
新規性: 生物学的な事前知識（MSigDB）をゲノムワイドな発現解析の統計モデルに直接組み込むアプローチを本研究で初めて体系化し、汎用的なソフトウェア GSEA-P として提供した。特に、遺伝子そのものを並び替えるのではなく、表現型ラベルを並び替える permutation test を導入したことで、遺伝子間の複雑な相関構造を保持したまま正確な有意性を評価する手法を新規に確立した。
臨床応用: 本手法は、がん臨床検体の分子サブタイプ分類や予後予測において極めて高い臨床的有用性を持つ。例えば、肺がんの独立した予後データセット間において、単一遺伝子レベルでは再現されなかった悪性度関連パスウェイ（低酸素応答や mTOR シグナル伝達など）を共通して同定できたことは、がん治療における標的パスウェイの同定や個別化医療への臨床応用に直結する。
残された課題: 一方で、今後の検討課題として、(1) 重み付け関数や相関尺度の選択が結果に与える影響の精査、(2) 異なる遺伝子セット間で遺伝子が重複している場合の検定統計量の相関の扱い、(3) leading-edge subset の境界を決定する不連続な閾値の解釈、(4) RNA-seq やシングルセル解析（single-cell RNA-seq）といった次世代シーケンスデータへの適応（後に GSVA や ssGSEA、fgsea 等の拡張手法へと発展）が挙げられる。
総じて、本論文はトランスクリプトーム解析におけるデファクトスタンダード（de facto standard）を確立した方法論的マイルストーンであり、現在もバイオインフォマティクスおよび腫瘍学研究の出発点であり続けている。
方法
GSEA の入力データは、N 個の遺伝子と k 個のサンプルからなる発現データセット D、表現型クラスラベル C（例えば、クラス1 vs クラス2）、および事前に定義された遺伝子セット S（メンバー数 NH (Number of genes in the Gene Set)）である。解析手順は以下の3つのステップから構成される。
ステップ1：Enrichment Score (ES) の算出 まず、全遺伝子を発現プロファイルと表現型 C との相関尺度（Spearman correlation や Pearson correlation、あるいは signal-to-noise ratio 等）に基づき降順に並び替え、順位付き遺伝子リスト L = {g1, …, gN} を作成する。本解析では、A549 や H1299 などの肺がん細胞株、あるいは MCF-7 や HEK293T を含む多様ながん細胞株（n=50 cell lines）から得られた発現プロファイルや、C57BL/6J や BALB/c などのマウス（n=86 mice）を用いた in vivo 実験データ、さらには遺伝子ノックダウン（siRNA や shRNA）や遺伝子ノックアウト（CRISPR-Cas9）による摂動実験データセットなどを対象とすることができる。任意の遺伝子セット S について、リスト L を上位から走査（random walk）し、遺伝子 gj が S に含まれる場合（hit）は running-sum 統計量を増加させ、含まれない場合（miss）は減少させる。この際、ステップの増分は遺伝子と表現型の相関値の絶対値（べき乗 p で重み付け、本研究では p=1 を推奨）に比例させ、減分は非所属遺伝子数で均等に配分する。ES は、この random walk においてゼロから得られた最大絶対偏差として定義される。これは重み付き Kolmogorov-Smirnov 様統計量に相当する。
ステップ2：有意性の評価 得られた ES の統計的有意性（nominal P 値）を評価するため、表現型ラベルをランダムに並び替える permutation test（1000回実施）を行う。各 permutation において遺伝子を再順位付けし、ES を再計算することで経験的零分布を生成する。この際、遺伝子そのものを並び替えるのではなく、表現型ラベルを並び替えることで、遺伝子間の複雑な相互相関構造を完全に保存したまま、生物学的に妥当な零分布を得ることができる。
ステップ3：多重検定の補正 複数の遺伝子セットを同時に評価する場合、セットのサイズ（遺伝子数）による影響を補正するため、各遺伝子セットの ES を、permutation から得られた零分布の平均値で除して NES を算出する。さらに、全遺伝子セットにわたる多重仮説検定の補正として、Reiner et al に基づき FDR を計算し、偽陽性の割合を制御する。
さらに、ES が最大値に達するまでに寄与した遺伝子群を leading-edge subset として抽出し、生物学的なコアとなる遺伝子群を特定する。本研究では、1325 個の遺伝子セットを含む MSigDB v1.0 を構築した。これには、染色体位置に基づく C1（319 セット）、既知のパスウェイや化学的・遺伝学的摂動に基づく C2（522 セット）、転写因子結合モチーフに基づく C3（57 セット）、がん関連遺伝子の発現近傍に基づく C4（427 セット）が含まれる。解析は、R および Java で実装された GSEA-P ソフトウェアを用いて実行された。