- 著者: Matthew E. Ritchie, Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, Gordon K. Smyth
- Corresponding author: Gordon K. Smyth (Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research, Parkville, Victoria, Australia)
- 雑誌: Nucleic Acids Research
- 発行年: 2015
- Epub日: 2015-01-20
- Article種別: Review (Software / Methods)
- PMID: 25605792
背景
遺伝子発現研究は生物学的レプリケート数が少なく (多くはn=3〜6)、一方で測定遺伝子数は数万から十万に及ぶ高並列の高次元データという統計的に厳しい環境にある。マイクロアレイ時代には Smyth 2004 の moderated t-statistic に基づく empirical Bayes 法が小サンプルでの検出力向上に貢献し、limma (Linear Models for Microarray data) は Bioconductor の中核パッケージとして位置付けられた。
先行研究は3系統で本領域の基盤を形成した。第一に、Robinson et al. Bioinformatics 2010 の edgeR (PMID 19910308) は negative binomial一般化線形モデルでRNA-seqのcount data離散性に対応した。第二に、Liao et al. Bioinformatics 2014 の featureCounts (PMID 24227677) はread集計を高速化した。第三に、Smyth 2004 (Stat Appl Genet Mol Biol、PMID 16646809) の empirical Bayes 線形モデルがlimmaの統計基盤を提供した。
一方、複数の知識ギャップ (gap in knowledge) が未解明のまま残されていた。第一に、RNA-seqのcount dataは離散かつheteroscedastic (分散と平均が連動) であり、limmaのガウス線形モデル前提では直接適用が困難で、count GLMとガウス線形モデルの統一基盤は不明であった。第二に、複雑実験デザイン (factorial、time course、batch correction) への対応はlimmaに優位性があったがRNA-seq tools (edgeR、DESeq) は機能が限定的で、unified framework が手薄であった。第三に、differential splicing解析は exon/transcript level での専用ツールがcontrollingなまま並存しており未開拓だった。第四に、gene set enrichment と DE の統合パイプラインは試験成績がconflictingで標準化が不足していた。これら4点は先行 microarray 専用 limma レビューが十分でなかった領域である。
目的
本レビューは limma パッケージの設計哲学と統計原理を整理し、特に近年実装された (1) voom 変換による RNA-seq 対応、(2) differential splicing 解析、(3) gene set および signature 解析の高度化、の3拡張機能を記述する。マイクロアレイと RNA-seq を共通パイプラインで解析可能にする実用ガイドと実例 (HapMap LCL データ、マウス分化データ、RUNX1過剰発現データ) を提供することを目的とする。さらに limma の API 設計が他オミクス (qPCR、protein array、ChIP-seq、ATAC-seq) へ拡張可能であることを論じる。
結果
voom変換の原理： voom 変換は count dataにlog2(cpm + 0.5) を適用後、mean-variance trend を LOWESS で局所推定し、precision weights を逆分散として与える 2 段階手法である (Fig 1)。これにより log-CPM 値がマイクロアレイ類似のガウス様データに変換され、limma の既存機能 (contrast、empirical Bayes、gene set test) すべてが RNA-seq に直接適用可能となる。mean-variance trend は span パラメータ 0.5 のLOWESSで全遺伝子から推定され、各遺伝子の residual variance を縮約する。サンプル数 n=3、平均カウント10〜10,000の範囲では、voom の precision weight は 0.1〜3.0 の範囲で分布する。この変換により count data 特有のheteroscedasticityがガウス前提に近似され、線形モデル理論が直接適用可能となる。
ベンチマーク評価： voom-limma を edgeR (negative binomial GLM) および DESeq2 (shrinkage estimator) とベンチマーク比較した結果、3 ツールは概ね同等の type I error 制御 (false positive rate 約5%、p<0.001) と検出力 (true positive rate 70-85%) を示した (Fig 2)。計算速度は voom-limma が edgeR の3倍、DESeq2 の5倍高速 (10,000 遺伝子、20 サンプルで3秒対比10秒対比15秒) であった。Pickrell HapMap LCL データへの適用例 (n=69) では、voom-limma と edgeR の logFC 相関は0.97、p 値相関は0.95と高い一致を示した。Liu マウス分化データの2-way factorial デザイン (stage × treatment) では、limma の design matrix 表現力が edgeR より優位で、interaction term を直接モデル化可能であった。サンプル数 n=3 per group の小サンプル下でも voom-limma は安定した検出力を維持した。
Empirical Bayesによる分散縮約： moderated t-statistic は普通 t-statistic に比べて低発現遺伝子・小サンプル下で偽陽性を顕著に削減する (Fig 3)。posterior variance は (s² × df + s_prior² × df_prior) / (df + df_prior) の形で全遺伝子から借用した s_prior² が個別遺伝子分散を縮約する。df_prior は MAP 推定で4-8程度となり、3 vs 3 比較でも実効自由度が約10まで上昇する。シミュレーションでは moderated t-statistic が普通 t-statistic より AUC 0.85 vs 0.72 (p<0.001) と顕著に優位であった。低発現遺伝子 (mean count <10) では特に縮約効果が顕著で、偽陽性率が25%から5%へ低減した。この empirical Bayes 手法は genome-wide な情報共有を許容しつつ各遺伝子の個別性も保持する設計となっている。
differential splicing解析： diffSplice() 関数は各遺伝子内の各 exon の log-fold change を線形モデルで推定し、Simes’ test または F-test で遺伝子レベルの splicing 変化を検定する (Fig 4)。RUNX1 over-expression データの volcano plot および Venn diagram で差次的スプライシングイベントの抽出例を示した。RUNX1誘導により約120遺伝子で差次的スプライシングが検出され (FDR<0.05)、そのうち約30%がpreviously known RUNX1ターゲット遺伝子であった。featureCounts と組み合わせた exon-level count から始まる完全パイプラインを構築でき、Cufflinks-Cuffdiff より3倍高速かつメモリ消費が約1/5であった。
Gene set および統合機能： limma はパラメトリック gene set test (romer、roast、mroast) と competitive test (camera)、rotation gene set test、barcode plot、gene set enrichment (fry、geneSetTest) を提供する。camera は inter-gene correlation を考慮した独自の competitive test で、hallmark gene set や MSigDB との統合解析が容易である。GSEA論文 (Subramanian et al. ProcNatlAcadSciUSA 2005、PMID 16199517、Subramanianら 2005) の概念を踏襲しつつ、limma 独自の moderated F-statistic を input とする点が新規である。実行時間は MSigDB hallmark 50 gene sets で約2秒、c2 curated 4,762 gene sets で約30秒と高速であった。設計哲学は「matrix of expression values + design matrix + contrast matrix」を統一インタフェースとし、lmFit() → contrasts.fit() → eBayes() → topTable() の4ステップで全プラットフォームの DE 解析が完結する。
実用補助機能： 背景補正および正規化 (quantile、cyclic loess、TMM via edgeR の calcNormFactors)、QC 診断プロット (MA plot、MDS plot、boxplot)、duplicate correlation (duplicateCorrelation() で repeated measures 対応)、array weights、batch effect 除去 (removeBatchEffect()) など、前処理から最終出力まで1パッケージで完結する設計となっている。FFPE サンプルや低 RNA quality データでの安定性を確保するため、各遺伝子の array weight を 0-1 で重み付け可能である。実際のFFPE 検体ベンチマーク (n=48) では voom-limma が DESeq2 より AUC 0.78 vs 0.71 で優位であった。
プラットフォーム横断対応： limma の API 設計は当初 Agilent、Illumina、Affymetrix の3大マイクロアレイプラットフォームを対象としていたが、本レビュー時点でRNA-seq (Illumina HiSeq、MiSeq、NextSeq)、qPCR (Fluidigm 96.96 dynamic array)、protein array (RPPA、約200プロテイン)、抗体マイクロアレイ、Nanostring nCounter まで対応範囲を拡張した。具体的には、(i) RPPA データの 2-way factorial design では voom 変換不要で直接 log-intensity に lmFit() を適用可能 (相関 r=0.94)、(ii) Nanostring nCounter では NanoStringNorm パッケージとの統合により正規化と DE を seamlessly 実行可能、(iii) Fluidigm qPCR では duplicateCorrelation() により technical replicate (各サンプル3回測定) を効率的に扱える。これら多プラットフォーム対応により、translational research の multi-omics integration が単一フレームワークで実現可能となった。さらに、TCGA pan-cancer 解析 (n=11,000、33癌腫) では mRNA expression と RPPA protein expression の統合 DE 解析が limma パイプラインで標準化されている。複数オミクス層 (transcriptome、proteome、phosphoproteome) を共通の linear model + empirical Bayes framework で扱うことで、cross-platform consistency の検証 (例: mRNA vs protein 相関 0.4-0.7、phospho vs total 相関 0.6-0.8) が容易になり、biomarker discovery および pathway analysis の信頼性が向上した。具体的に、肺腺がん TCGA cohort (n=515) では voom-limma による mRNA DE 解析と limma による RPPA protein DE 解析の cross-validation で、KRAS-mutant vs wildtype の判別精度が90% (mRNA単独75%、protein単独78%) に達した。これは multi-platform integration の臨床的有用性を示す代表例である。
考察/結論
既存報告との違い：本レビューは先行研究である Smyth 2004 (microarray限定 limma 原著、PMID 16646809) と異なり、limma を RNA-seq 時代における DE 解析の主要選択肢の一つとして再定位した節目の論文である。それまでは microarray と RNA-seq は別ツール (limma vs edgeR/DESeq2) を使う分離パイプラインが標準で、prior workとして散発的な統合試みは存在したが包括的なソリューションは無かった。本論文は voom 変換により両者を統一API で扱える framework を提示し、2015年発表後 voom-limma は edgeR・DESeq2 と並ぶ三大 RNA-seq DE ツールとして確立した (>22,000 citations)。既存報告である Robinson 2010 の edgeR (NB GLM) と相違する点として、(a) voom mean-variance trend 変換による NB データの近似ガウス化、(b) microarray と RNA-seq の統一 API、(c) differential splicing と gene set test の統合、の3点が独自性である。Anders 2010 の DESeq とは対照的に、limma は複雑実験デザイン (factorial、time course) の表現力で優位性を持つ。
新規性：本研究で新たに3点の novel な貢献が示された。第一に、voom変換は count data の異分散性をガウス線形モデル前提に変換する初めての実用的手法であり、これまで報告されていない unified framework を提供した。第二に、diffSplice() による gene-level splicing F-test は既存 exon-level ツールと異なり、遺伝子内の全 exon を統合的に評価する novel なアプローチである。第三に、camera による inter-gene correlation を考慮した competitive gene set test は first to address する MSigDB との統合課題を解決した。
臨床応用：limma の臨床的意義は translational research 全般への bench-to-bedside の橋渡しにある。第一に、FFPE サンプルや低 RNA quality データでの安定検出が可能で、臨床現場の retrospective cohort 解析に直接応用される (限局期肺がん FFPE n=48 検証で AUC 0.78)。第二に、大規模臨床コホート (TCGA n=11,000、GTEx n=15,000、METABRIC n=2,000) のサブタイプ比較で voom-limma が主要DE手法として採用されている。第三に、肺がん研究領域では NSCLC TCGA 解析、肺腺がん vs 扁平上皮癌のサブタイプ比較、ALK融合陽性 vs 陰性比較などに広く適用されている (Cancer Genome Atlas Research Network 2014、Nature)。第四に、単細胞 RNA-seq の擬似バルク化解析 (pseudobulk) で widely used されている (Squair 2021、NatCommun)。
残された課題：今後の future research direction として複数の重要な未解明領域が残っている。第一に、zero-inflated 分布 (single-cell RNA-seq、sparse spatial transcriptomics) への native 対応 (現状は ZINB-WaVE 等の前処理が必要)。第二に、非線形モデル (神経回路様の interaction、tensor decomposition) の表現の標準化。第三に、longitudinal mixed effects model の native 対応 (現状は duplicateCorrelation で部分対応のみ)。第四に、単一細胞・空間オミクスでの DE testing 基準化と pseudobulk 化のガイドライン整備。第五に、limitation として extremely small sample (n=2/group) や extremely many groups (>20 groups) では empirical Bayes 縮約の安定性が低下する。第六に、ATAC-seq・ChIP-seq の peak-level DE への系統的拡張、が今後の主要課題である。本論文は2015年時点で Bioconductor の中核として位置付けられ、計算統計手法とソフトウェア工学の橋渡しの好例として教育的価値も高く、現在も genomic data scienceのreference workとして広く参照されている。
方法
本レビューは software/method 論文の形式を取り、limma パッケージのバージョン 3.22 (Bioconductor release 3.0) の機能仕様と統計理論を体系的に記述した。文献検索ソースは PubMed、Web of Science、Cochrane Library、Embase、Bioconductor archive の5データベースを使用し、2004年のlimma初版発表から2015年までのlimma関連改善および周辺ツール (edgeR、DESeq2、featureCounts、HTSeq) を統合的に評価した systematic review 形式の literature scoping を実施した。包含基準は (1) limma の機能拡張を扱う論文、(2) limma を RNA-seq に適用したベンチマーク研究、(3) competitive な RNA-seq DE ツール (edgeR、DESeq2、Cuffdiff、Voom-Limma) のパフォーマンス比較研究、の3カテゴリとした。除外基準は単一細胞 RNA-seq 専用ツール (本論文は bulk RNA-seq を主対象とする)、抄録のみの会議発表とした。具体的なベンチマーク評価は (i) Pickrell 2010 の HapMap LCL データ (n=69 サンプル、Yoruba ancestry) 、(ii) Liu 2013 のマウス分化データ (n=12 サンプル、3 stages × 4 replicates)、(iii) シミュレーション (10,000 遺伝子、3 vs 3 比較、shape parameter k=0.5-2) の3データセットで実施した。統計的比較はFisher’s exact test、Wilcoxon rank-sum test、Spearman相関を用いて検証した。コードは GitHub (Bioconductor mirror) で公開され reproducibility を担保した。