- 著者: Sonja Hänzelmann, Robert Castelo, Justin Guinney
- Corresponding author: Robert Castelo (Hospital del Mar Medical Research Institute / Universitat Pompeu Fabra, Barcelona, Spain); Justin Guinney (Sage Bionetworks, Seattle, WA, USA)
- 雑誌: BMC Bioinformatics
- 発行年: 2013
- Epub日: 2013-01-16
- Article種別: Original Article (Methods/Software)
- PMID: 23323831
背景
遺伝子発現プロファイルの解釈において gene set enrichment (GSE) 法は標準的なアプローチとして定着している。GSEA (Gene Set Enrichment Analysis) をはじめとする多くの手法が開発されてきたが、これらは基本的に case-control 型の2群比較を前提とし、表現型ラベルを与えてランクリストを生成する設計である。competitive test (遺伝子セット内外の差を検定) と self-contained test (セット内遺伝子のみを検定) の区別があるが、いずれも事前のサンプル分類情報を必要とする点で共通している。
一方、TCGA (The Cancer Genome Atlas) に代表される大規模ゲノム研究は、多数のサンプルと複雑な多クラス・連続・打切りデータ構造を持つコホートを対象としている。このような heterogeneous なデータセットでは、表現型ラベルが事前に定まらないサブタイプ発見・生存時間解析・連続的バイオマーカーとの相関解析など、サンプル単位で経路活性を定量したいユースケースが増加していた。既存の single-sample 手法 (ssGSEA: Barbie et al. 2009, PLAGE: Tomfohr et al. 2005, combined z-score: Lee et al. 2008) は提案されていたが、多様な gene set サイズやデータ分布への頑健性が課題として残っていた。さらに RNA-Seq の count data に対応した GSE 手法の不足も問題であり、当時は有効な単一サンプル GSE 解析がほとんど存在しなかった。既存手法がどのサンプルにも適用可能な unsupervised な pathway-activity 定量化を未確立のまま残していたこと、および microarray と RNA-Seq の双方に対応した統一フレームワークが不足していたことが、何が足りなかったかという核心的なギャップであった。
目的
サンプル集団全体にわたる経路活性変動を非教師的 (unsupervised) かつサンプル単位で推定する新たな GSE 法 GSVA (Gene Set Variation Analysis) を提案し、以下を達成することを目的とした：(1) microarray と RNA-Seq の双方に対応する統一フレームワークの提供、(2) 既存 single-sample 手法 (ssGSEA・PLAGE・combined z-score) との比較によるロバストネスと検出パワーの実証、(3) 差次的経路活性解析・生存解析等の下流解析への応用例の提示、(4) Bioconductor R パッケージとしての公開による広い利用可能性の確保。
結果
所見1：シミュレーションにおける統計的パワーの優位性： n = 30 サンプルを仮定した第1シミュレーションで、subtle な経路活性差 (全体の50%の遺伝子のみが差次発現) の検出パワーを評価した結果、GSVA は ssGSEA・PLAGE・combined z-score のすべてを上回る検出感度を示した (significance threshold α = 0.05) (Fig. 2)。この優位性は4つのシミュレーション条件すべてで一貫しており、特に gene set size が小さい (10-50遺伝子) 場合に顕著であった。いずれの条件でも type-I error は同等に制御されており (α = 0.05 付近)、感度の向上がα膨張によるものでないことを示した。
所見2：AUC による比較 (第2シミュレーション)： p = 10,000 遺伝子・n = 60・500 DE gene sets を対象とした第2シミュレーションで、5% FDR での DE gene set 検出 AUC を100回独立シミュレーションで評価した結果、GSVA は12ペアワイズ比較のうち10 (P < 0.05) で他の手法より有意に高い平均 AUC 値を示した。1% FDR のより厳格な条件でも優位性は維持された。このことは、GSVA が FDR 制御のもとで real な pathway 変化を検出する能力が高いことを示す。
所見3：生存解析での concordance index (第3シミュレーション)： Cox 比例ハザードモデルを用いた生存予測 (training/test 独立データ) の評価では、GSVA は n = {25, 50, 75, 100} の全4サンプルサイズで他の手法より高い mean/median concordance index を達成した。n ≥ 50 では差の有意性が確認された (P < 0.05)。この結果は GSVA が生存解析などの下流予後モデリングにも適切であることを示す。
所見4：白血病実データ (ALL vs MLL) での分類性能： ALL 20例・MLL 17例の白血病データに対し、Broad C2 collection (833 canonical pathways + 2,392 chemical/gene perturbations) を適用し、limma による差次検定後の上位5 gene set の enrichment score で階層クラスタリングを実施した。細胞種で最小の fold change 三分位 (subtle な変化のみ) を用いた場合でも、GSVA の ARI (Adjusted Rand Index) は ssGSEA・PLAGE・combined z-score より有意に高い値を示した (t 検定 P < 2e-16)。最大 fold change 三分位 (strong な変化) では差が縮小したが、これは strong signal では全手法が同等に機能するためである。
所見5：TCGA卵巣がん生存解析： TCGA 卵巣漿液性腺がん (n=389 patients) に GSVA を適用し、Cox PHM で経路レベルの生存相関を評価した (Fig. 6, Table 1)。5分割交差検証で GSVA は training・test データセットの双方で他の手法より高い concordance index を達成した。上位5経路として同定されたのは SIMBULAN_UV_RESPONSE_NORMAL_DN (Cox P = 7.21×10^-6)、BIOCARTA_VIP_PATHWAY (P = 1.38×10^-5)、ZIRN_TRETINOIN_RESPONSE_WT1_UP (P = 3.38×10^-5)、DASU_IL6_SIGNALING_SCAR_UP (P = 3.46×10^-5)、WANG_HCP_PROSTATE_CANCER (P = 3.65×10^-5) で、DNA修復・免疫調節・tretinoin応答・EGF/RAS 関連経路が予後と関連した。FDR 評価では P < 10^-4 で FDR = 0.05 を達成し、P < 5×10^-3 では FDR = 0.2 と推定された。
所見6：RNA-Seq への適用 — microarrayとの高相関： HapMap Yoruba 個体由来リンフォブラストイド細胞株 (n=69 samples, Affymetrix microarray と RNA-Seq 両プラットフォーム) に GSVA を適用し、プラットフォーム間の相関を評価した。male-specific Y 染色体遺伝子セットでは Spearman ρ=0.82、X 不活性化逸脱遺伝子セットでは Spearman ρ=0.78 と高い相関が得られ (Fig. 7)、Poisson kernel を用いた count データ対応が有効であることを実証した。これにより microarray と RNA-Seq のプラットフォーム横断 meta-analysis が容易になることが示された。
考察/結論
本研究で初めて示されたことは、表現型ラベル不要でサンプル単位の pathway activity score を非教師的に推定できる統一フレームワーク GSVA を確立した点、および microarray と RNA-Seq の双方に対応する kernel density estimation により both platforms での高精度な gene set enrichment が実現可能であることである。GSVA は GSE を「バイオインフォマティクス解析の終点」から「pathway-centric な生物学モデリングの起点」へと位置づけ直す方法論的シフトを提示した。
先行研究 (Subramanian et al. ProcNatlAcadSciUSA 2005 の GSEA, ssGSEA, PLAGE, combined z-score) との比較では、PLAGE は SVD の不安定性により小サンプル (n < 20) で精度が低下し、combined z-score は large gene set で過剰なスコアを与える傾向があるのに対し、GSVA の KS-like statistic は gene set サイズに対して頑健でサイズ間で比較可能なスコアを生成する点で異なる。先行研究と異なり、本研究は supervised な 2群比較に限定されず、unsupervised な任意コホート構造に対応した点で根本的に新しい手法論を提供した。臨床応用として、n=389の TCGA卵巣がんコホートでの生存予測への適用や白血病 ALL/MLL 分類での有用性が示され、臨床的意義として pathway-level の予後バイオマーカー探索に直接使用可能なツールとしての価値が実証された。
Limitations として、(1) ランクベース手法のため絶対発現量情報が失われること、(2) batch effect への感受性が存在し CombBat などの前処理が推奨されること、(3) non-parametric kernel 推定には最低 n ≥ 10 程度のサンプルが必要であること (< 10 では精度低下) が挙げられる。また、遺伝子が独立でない (pathway 内相関がある) ことが多重検定補正を複雑にすることも指摘されており、著者は Benjamini-Hochberg 補正を保守的 FDR 推定として使用した。下流の差次発現解析には Robinson et al. Bioinformatics 2010 や limma との組み合わせが推奨され、gene set 解釈基盤として Liberzon et al. CellSyst 2015 との組み合わせが広く実用されている。
将来の展望として、GSVA を eQTL マッピングと類似した「pathway-QTL」解析に応用し、DNA多型が経路活性に与える影響を解析すること、さらにはゲノム → 遺伝子発現 → 表現型の因果推論モデルを pathway レベルに拡張することが提案されている。本ツールは公開後、Bioconductor download ランキング上位を継続して維持しており、がんトランスクリプトミクスにおける pathway-level 解析の事実上の標準基盤となっている。
方法
アルゴリズム設計：GSVA はサンプルごとに gene set enrichment score を生成する2段階手法。第1段階では各遺伝子の発現量をサンプル集団内で non-parametric にkernel density estimation (microarray: Gaussian kernel、RNA-Seq count: Poisson kernel) し、cumulative distribution function (CDF) を構築する。各 (遺伝子、サンプル) ペアの CDF 評価値からランクに変換し、ランクをさらに対称化 (ゼロ周辺の symmetric distribution) する。第2段階では gene set ごとに Kolmogorov-Smirnov (KS) like random walk statistic を計算し、gene set 内の遺伝子と外の遺伝子のランク分布の最大差から enrichment score を算出する。具体的には、ES の計算に2種類のアプローチを提供した：(1) 最大偏差スコア (ES_max、bimodal 分布) と (2) 正規化スコア (ES_diff = ES+ + ES-、unimodal 近似正規分布)。後者は concordant に活性化・抑制される遺伝子セットを検出するのに適する。
比較対象手法：ssGSEA (Barbie et al.)、PLAGE (SVD ベース)、combined z-score を同一フレームワークで実装して比較評価した。いずれも single-sample 手法であるが、パラメトリック仮定や pheno type 情報利用に差異がある。
シミュレーション：p = 1,000 遺伝子・n = {10, 20, 40, 60} サンプルの線形加法モデルによる microarray データをシミュレートし、DE gene set (30遺伝子) と non-DE gene set の2種類を設定。フラクション (50% vs 80%) とシグナル-ノイズ比 (weak vs strong) の4条件。別に p = 10,000 遺伝子・n = 60・1,000 gene sets (うち 500 DE) のシナリオで AUC を評価。さらに生存解析シミュレーションとして n = {25, 50, 75, 100} の4サイズで concordance index を評価 (100 反復)。
実データ解析：(1) 白血病データ (ALL 20例 vs MLL 17例、Broad C2 collection 適用、limma で差次検定、ARI で評価)、(2) TCGA卵巣がんコホート (n = 389、MSigDB C2 gene set、Cox PHM、5分割交差検証、concordance index で評価)、(3) HapMap リンフォブラストイド細胞株の microarray vs RNA-Seq paired data での Spearman 相関評価。