- 著者: Michael I. Love, Wolfgang Huber, Simon Anders
- Corresponding author: Simon Anders (European Molecular Biology Laboratory, Heidelberg, Germany)
- 雑誌: Genome Biology
- 発行年: 2014
- Epub日: 2014-12-05
- Article種別: Original Article
- PMID: 25516281
背景
RNA-seq (RNAシーケンス) などのハイスループットシーケンス技術の普及に伴い、遺伝子ごとのリードカウントを群間で比較し、発現変動遺伝子を検出する統計的手法の確立が求められてきた。カウントデータは、(i) 非正規性、(ii) 平均値依存的な分散構造、(iii) わずか2〜3レプリケートという典型的な小標本サイズ、(iv) 大きなダイナミックレンジと外れ値の存在、といった特性を併せ持つため、適切な統計的枠組みが必要とされる。先行研究として、Robinson et al. Bioinformatics 2010が開発したedgeR（TMM正規化とtagwise dispersion shrinkage）や、Anders and Huber (2010) が開発したDESeq（median-of-ratios正規化とdispersion-mean trend）が標準的に使用されてきた。また、マイクロアレイデータの解析手法としてSmyth (2004) による線形モデルと経験ベイズ法も広く知られていた。しかし、これらの先行研究の手法にはいくつかの課題が残されていた。具体的には、(a) 分散推定の不安定さ、(b) 低カウント遺伝子においてLFC (logarithmic fold change; 対数フォールド変化) が極端な値となる問題、(c) サンプルサイズが大きい場合に統計的に有意であっても生物学的に取るに足らない小さな差を多数検出してしまう問題、(d) 外れ値の影響、などが挙げられる。これらの課題は、特にLFCによる遺伝子ランキングや可視化、ユーザー定義の生物学的閾値に対する仮説検定、品質評価のためのクラスタリングといったデータ解釈の場面で深刻な問題を引き起こした。例えば、LFCの不安定性は、独立した実験間での結果の再現性を低下させ、低発現遺伝子における偽陽性の原因となることが指摘されていた。また、大規模な研究では、統計的有意性のみに基づいて遺伝子リストを作成すると、生物学的に関連性の低い遺伝子が多数含まれることになり、その解釈が困難になるという課題も存在した。これらの問題に対処するため、より安定した推定と柔軟な解析機能を持つ新たなツールの開発が強く求められていた。特に、分散とLFCの両方に対する縮小推定 (shrinkage estimation) の導入は、小標本サイズでの推定の安定性を向上させる上で不可欠な要素であると考えられたが、この点において既存の手法では十分な対応が不足しており、低カウント領域におけるLFC推定の信頼性を担保する統計的アプローチが未確立であるという明確な知識ギャップが存在した。
目的
本研究の目的は、先行手法であるDESeqの後継として、RNA-seqカウントデータの発現変動解析における課題を克服する包括的なR/Bioconductorパッケージ「DESeq2」を開発することである。具体的には、以下の機能の統合と性能評価を目指した。(1) 分散（dispersion）とLFC (logarithmic fold change) の両方に経験ベイズ法に基づく縮小推定（shrinkage estimator）を統合的に適用し、推定の安定性と再現性を向上させること。(2) 特に低カウント遺伝子においてLFCが極端な値となる問題を抑制し、信頼性の高いLFC推定値を提供すること。(3) データ中の外れ値を自動的に検出し、頑健に処理するメカニズムを実装すること。(4) 分散安定化変換であるrlog (regularized log; 正則化対数) 変換を導入し、品質評価やクラスタリングなどの多変量解析に適したデータを提供すること。(5) ユーザー定義のLFC閾値（effect size threshold）を考慮した仮説検定を可能にし、生物学的に意味のある発現変動を検出すること。(6) 並列計算に対応し、大規模データセットの解析効率を向上させること。これらの新機能を統合したDESeq2の性能を、既存の主要な発現変動解析ツール（edgeR、DESeq、limma-voom、EBSeq、DSSなど）と比較し、感度、特異度、およびFDR (false discovery rate; 偽発見率) 制御の観点から優位性を示すことを目的とした。
結果
分散縮小による少サンプル推定の安定化: Bottomlyデータセット（n=4 replicates）を用いた解析では、遺伝子ごとのMLEによる分散推定値は広範囲に散らばり、特に低発現遺伝子で極めて不安定であった。これに対し、DESeq2の経験ベイズMAP推定は、分散-平均トレンド曲線に向けて適切に収束し、推定の安定性が大幅に向上した (Figure 1)。また、生物学的に真に高い分散を示す遺伝子（分散外れ値）は、トレンド曲線から大きく逸脱している場合、縮小されずにその遺伝子固有の分散値が維持される設計が機能した。これにより、少サンプル数でも信頼性の高い分散推定が可能となり、偽陽性のリスクが低減された。
LFC縮小による極端な値の抑制と再現性向上: 低カウント遺伝子において、MLEによるLFCが±10以上に膨らむという問題がRNA-seqデータで頻繁に観察された (Figure 2)。DESeq2の経験ベイズLFC縮小は、事前分布によりこれらの極端なLFC値をゼロ方向へ適切に縮小し、実用的な範囲に抑制した。Bottomlyデータセットにおいて、平均カウントが10未満の遺伝子のLFCは、縮小後には0付近に集中し、上位のLFCランキングは高カウントかつ統計的に頑健な遺伝子で占められるようになった。これにより、MAプロットの解釈性が劇的に改善し、LFCの再現性も向上した。例えば、Bottomlyデータセットを2つの独立したグループに分割してLFCを比較した結果、MLE LFCではRMSE (root-mean-square error; 平均二乗誤差平方根) が2.7であったのに対し、MAP LFCではRMSEが1.0に改善した (Figure 3)。さらに、上位100のLFC遺伝子の重複率は、MLE LFCで21%であったのに対し、MAP LFCでは81%に向上した。
シミュレーションによる感度とFDR制御の優位性: n=6 replicates vs n=6 replicates の実験デザインで、30%の遺伝子が発現変動し、効果量（fold change）が1.5x、2.0-fold、4.0-foldのいずれかであるという条件で1,000回反復シミュレーションを行った。DESeq2は、edgeRと同等以上の感度を維持しつつ、FDRを名目値（0.1）以下に適切に制御した (Figure 6)。特に、効果量が小さい（2.0-foldまたは3.0-fold）場合にDESeq2の感度が他のアルゴリズムよりも高かった。limma-voomは大規模サンプルで強力な性能を示したが、小規模サンプル（n=3 replicates）ではDESeq2の方が精度が高かった。EBSeqはDESeq2よりもLFC縮小が緩やかであり、低カウント遺伝子での偽陽性が多く見られた。LFC推定の精度に関して、DESeq2はGFOLDやedgeRのpredictive LFCと比較して、幅広いサンプルサイズとモデル条件で一貫して低いRMSEとMAE (mean absolute error; 平均絶対誤差) を示した。
独立フィルタリングによる検出力の改善: 正規化された平均カウントに基づく独立フィルタリングを適用することで、低発現遺伝子における多重検定のペナルティを回避し、真陽性検出力が約5%〜15%改善された。Benjamini-Hochberg法によるFDRは、名目レベル（5%）によく一致した。Pickrellデータセットを用いた模擬比較（n=5 replicates vs n=5 replicates、30回反復）では、DESeq2は名目FDR 0.01に対して実際のFPR (false positive rate; 偽陽性率) を0.01以下に適切に制御した (Figure 7)。
外れ値処理の頑健性と補完アルゴリズム: シミュレーションにおいて1%の遺伝子に外れ値を注入した条件下で、DESeq2のCook’s distanceに基づくフィルタリングは外れ値を確実に検出し、edgeRやvoomと比較して偽陽性率を大幅に低減した。レプリケート数がn=7 replicates以上の場合には、外れ値のカウントをトリム平均で補完し、モデルを再適合させることで、情報を完全に破棄することなく頑健な推定を可能とした。このアプローチにより、特定の遺伝子における異常なカウントが全体の統計的推論に与える影響を最小限に抑え、より信頼性の高い結果を得ることができた。
rlog変換による多変量解析とクラスタリングの最適化: GTExの組織横断データセットを用いてrlog変換後のサンプル間主成分分析（PCA）を行った結果、組織種別が明瞭に分離された。rlog変換は、従来の分散安定化変換よりも低発現遺伝子のノイズ抑制に優れており、n=3 replicatesのようなサンプル数の少ないグループでもクラスタリングの安定性が向上した (Figure 5)。これは、rlog変換が低カウント領域でのLFCの過剰な変動を抑制するためである。
生物学的閾値に基づく仮説検定の挙動: H0: |LFC| ≤ 0.5やH0: |LFC| ≤ 1.0のような生物学的に意味のある閾値に対する仮説検定が標準実装された (Figure 4)。これにより、大規模サンプル数の研究で「統計的有意性 ≠ 生物学的重要性」という問題に対処可能となり、ユーザーは生物学的に関連性の高い効果量を持つ遺伝子を直接検出できるようになった。
考察/結論
先行研究との違い: DESeq2は、先行するDESeqやRobinson et al. Bioinformatics 2010、limma-voomと対照的に、分散とLFCの両方に整合的な経験ベイズ縮小推定を適用した点、および生物学的閾値を直接組み込んだ仮説検定を可能にした点で大きく異なる。従来のツールでは低カウント遺伝子のLFC過大評価や小サンプルにおける分散の不安定性を十分に解決できなかったが、DESeq2はこれらを統計的に厳密な枠組みで克服した。
新規性: 本研究で初めて、負の二項一般化線形モデル（GLM）において、分散とLFCの双方に経験ベイズ縮小推定を統合的に導入した。また、Cookの距離を用いた外れ値の自動検出・置換アルゴリズムや、多変量解析のためのrlog変換を新規に開発・実装し、発現変動解析の安定性と再現性を劇的に向上させた。
臨床応用: 本知見は、がんゲノムアトラス（TCGA）などの大規模臨床RNA-seqデータ解析、患者由来異種移植（PDX）モデルやオルガノイドを用いた薬剤応答性評価、免疫チェックポイント阻害剤のバイオマーカー探索など、臨床現場におけるトランスクリプトームデータの解釈に広く臨床応用されている。特に、少数の臨床検体から得られるデータの信頼性を担保する上で、DESeq2は不可欠なデファクトスタンダードツールとなっている。
残された課題: 今後の検討課題として、シングルセルRNA-seqデータに特有のゼロインフレーション（zero-inflation）や細胞レベルのノイズに対する適応が挙げられる。また、10,000サンプルを超える極めて大規模なコホート解析においては、GLMの反復計算に伴う計算コストが課題となるため、計算効率のさらなる改善や、バッチ効果補正ツールとのシームレスな連携アルゴリズムの開発が望まれる。
方法
統計モデルと正規化: DESeq2では、遺伝子iのサンプルjにおけるリードカウントKijを、平均μijと分散αiを持つ負の二項分布NB(μij, αi)に従うとモデル化した。平均μijは、サンプルjのサイズファクターsijと遺伝子iの真の濃度qijの積（μij = sij × qij）として定義される。サイズファクターsijは、DESeqから継承されたmedian-of-ratios法を用いて推定された。qijは、logリンク関数を持つGLM (generalized linear model; 一般化線形モデル) により、log2(qij) = Σr xjr × βirとモデル化された。ここで、xjrはデザイン行列の要素、βirはモデル係数である。
分散推定の経験ベイズ縮小: 分散αiの推定には、経験ベイズ縮小法が用いられた。まず、各遺伝子についてデータのみに基づくMLE (maximum-likelihood estimate; 最大尤度推定) により遺伝子ごとの分散推定値（αigw）を算出した。次に、これらのαigwを平均発現量に対してプロットし、分散の平均依存性を示すパラメトリックなトレンド曲線（αtr = a₀ + a₁/μ̄）をGLM回帰により適合させた。最後に、個々の分散推定値αigwを、このトレンド曲線によって予測される値に向けて縮小し、最終的な分散推定値（αiMAP）を得た。トレンドから大きく外れる分散外れ値（log αigw > log αtr(μ̄i) + 2slr）は縮小されず、遺伝子ごとの推定値がそのまま採用された。
LFC (logarithmic fold change) 縮小推定: モデル係数βir（LFCを表す）に対しては、ゼロを中心とする正規事前分布を適用し、MAP (maximum a posteriori; 最大事後確率) 推定値として最終的なLFC推定値を取得した。縮小の強さは、平均カウント、分散、および自由度によって決定される情報量に依存する。
仮説検定: デフォルトでは、縮小されたMAP-LFCをその標準誤差で割ったWald統計量を用いて仮説検定を行った。p値はBenjamini-Hochberg法により多重検定調整され、FDRが制御された。
外れ値検出と処理: Cook’s distanceを用いて外れ値を検出した。Cook’s distanceがF分布の0.99分位点を超えるサンプルを持つ遺伝子は、外れ値としてフラグ付けされた。6レプリケート以下の条件では、外れ値を持つ遺伝子は解析から除外され、p値はNAとされた。7レプリケート以上の条件では、外れ値カウントはトリム平均で補完され、その後分散、LFC、p値が再推定された。
独立フィルタリング: 各遺伝子の正規化された平均カウントを発現変動検定のフィルタリング基準として使用した。これにより、検定力が低い低発現遺伝子を多重検定調整から除外し、FDRの制御を改善しつつ、真陽性検出力を向上させた。
rlog (regularized log) 変換: 分散安定化変換の改良版として、rlog変換を導入した。これは、高カウント遺伝子ではlog2変換と同様に振る舞う一方、低カウント遺伝子では異なるサンプル間の値を縮小し、ノイズの影響を抑制する。
性能評価とデータセット: DESeq2の性能は、Bottomly et al. (2011) のマウス系統データ（C57BL/6JおよびDBA/2J）、Pickrell et al. (2010) のヒトリンパ芽球様細胞株データ、GTExデータ、およびシミュレーションデータを用いて評価された。比較対象には、edgeR、voom (limma)、EBSeq (empirical Bayes RNA-seq)、DSS (dispersion shrinkage for RNA-seq)、SAMseq、Cuffdiff2などの既存ツールが含まれた。統計解析には R ソフトウェアおよび Bioconductor パッケージが使用され、グループ間の平均値比較や相関分析において Student t-test および Spearman correlation が適用された。