- 著者: Nicolas L Bray, Harold Pimentel, Páll Melsted, Lior Pachter
- Corresponding author: Lior Pachter (lpachter@math.berkeley.edu, UC Berkeley)
- 雑誌: Nature Biotechnology
- 発行年: 2016
- Epub日: 2016-04-04
- Article種別: Brief Communication (Methods)
- PMID: 27043002
背景
RNA-seqデータ処理の標準的な2ステップ (アライメント → 転写産物量定量) は計算量が膨大で、データ規模の拡大に伴い実用的なボトルネックとなっていた。代表的なパイプラインの TopHat2 (ref. 1) + Cufflinks では、20サンプル各3,000万リードの処理に20コアで28コア時間のアライメントとさらに14時間の定量が必要であり、次世代シーケンス技術の高スループット化とともに実行不能となりつつあった (Trapnell et al. NatBiotechnol 2014)。ストリーミングアルゴリズムやナイーブリードカウントによる高速化の試みは定量精度の低下を招いた。k-mer (k-塩基長のDNA断片) ハッシュテーブルを用いたSailfish (Patro et al. 2014) はアライメントステップを回避したが、リードをk-merへ分断することにより各k-merが元のリードより多くの転写体にアライメントされる情報損失が生じ、精度が著しく低下した (Supplementary Fig. 1)。核心的な問題は多重マッピングリードの取り扱いで、これには確率的統計モデルによる最尤推定が必要とされていたが、アライメントなしに同等の精度を達成するための数学的基盤が不足していた (Aibar et al. NatMethods 2017)。正確な定量に必要なのは「リードが転写体内のどこから来たか」でなく「どの転写体から来たか」という等価クラス情報で十分であることが理論的に示されており (ref. 7)、これを高速実装する計算幾何学的手法が未解明な課題として残されていた。
目的
T-DBG (transcriptome-derived Bruijn graph; 転写体deブライングラフ) とk-merハッシュを組み合わせた擬似アライメント (pseudoalignment; 個別塩基のアライメントを省略し等価クラスのみを同定) アルゴリズムを開発し、アライメントと同等の精度を維持しつつ2桁高速な RNA-seq 定量ツール「kallisto」を実装・検証すること。
結果
kallistoのアルゴリズム — T-DBGとk-mer擬似アライメントによる超高速等価クラス同定：
kallistoはまずトランスクリプトームからT-DBG (transcriptome-derived Bruijn graph) を構築する (Fig 1a, b)。各ノードはk-merに対応し、各転写体は有色パスとして表現される。パスカバリングによってk-compatibility class (k適合クラス) が各k-merに帰属し、同一適合クラスを持つk-merのコンティグが定義される。リードの擬似アライメントは、リードのk-merを順次ハッシュで参照してk-compatibility classを取得し、その積集合を取ることで等価クラス (equivalence class; そのリードが生起しうる転写体集合) を同定する (Fig 1c-e)。最重要な最適化は「skipping (スキップ)」で、同一コンティグ内のk-merは同一適合クラスを持つため、コンティグ境界のk-merだけを参照すれば全リードの大半でハッシュルックアップがわずか2回で完了する (Supplementary Fig. 11)。等価クラスが同定されれば、確率的RNAモデルに基づくEM (expectation-maximization; 期待値最大化) アルゴリズムで等価クラスカウントから転写体発現量を最尤推定し、TPM (transcripts per million; 百万転写産物当たりの転写産物数) として出力する。配列特異的バイアスはfragment 5’端周辺の6-merを100万リードから実測して補正する。
速度評価 — 従来法の100倍以上高速、標準ラップトップで実用的：
GEUVADIS (Genetic European Unified Variation-Analysis Dataset) コホートのサンプルNA12716_7をベースにRSEMで生成した30億リード (75 bp paired-end) のn=20 samples (RSEMシミュレーション) に対し、20コア・20スレッドで処理した (Fig 2b)。kallistoの合計実行時間は11.47分であり、TopHat2+Cufflinksの42コア時間に対して約220倍高速であった。同データセットの単純な文字数カウント (word count) が75秒という理論的下限値と比較してもkallistoの速度は「near-optimal (準最適)」水準にある (Fig 2b)。1.3-GHz (gigahertz) プロセッサの小型ラップトップでも3,000万リードを10分以内に処理でき (最大3.2 GB RAM)、大規模データセットでも汎用ハードウェア上での RNA-seq解析が現実的となった。
精度評価 — RSEMと同等、Cufflinks・Sailfishを上回る：
20シミュレーションのKallisto精度 (転写体推定カウントの中央相対差) はRSEMとほぼ同等で、Cufflinks・Sailfishを大幅に上回った (Fig 2a)。Cufflinksの精度劣化はゲノム跨ぎ多重マッピングリードへのEM適用が限定的なことに起因する。擬似アライメントの等価クラスの品質をBowtie2のアライメントと比較すると、70.7%のリードで完全一致し、乖離したリードではBowtie2が平均8.02転写体を報告したのに対しkallistoは4.96転写体 (より特異的) であった。感度は99.89%でBowtie2の99.99%に対して0.1%以下の差にとどまった。実データ (NA12716_7) でも同様の傾向で、66.22%で完全一致、アライメント率はBowtie2の86.5%に対しkallistoが90.8%と優れていた。対立遺伝子特異的発現定量ではSpearman r=0.833 (RSEM 0.848、eXpress 0.830、Sailfish 0.706) と競合ツールと同水準を達成した (Supplementary Fig. 2)。SEQC (Sequencing Experimental Quality Control) データセットの qPCR独立検証でもkallistoは他手法と同等のパフォーマンスを示した (Supplementary Table 1)。
ブートストラップによる発現量不確実性の定量 — 従来法では不可能だった計算を実用化：
高速な擬似アライメントにより、ブートストラップ (bootstrap; 復元抽出による反復再サンプリング解析) が実用時間で実施可能となった。等価クラスのカウントが十分な統計量であるため、各ブートストラップサンプルの生成は多項分布から等価クラスカウントをサンプリングし直すだけで完了し、1サンプルあたり約10秒で処理できる。216億リードの深いSEQC-MAQCIII (Microarray Analysis Quality Control III) ヒト RNA-seqデータ (Trapnell et al. NatBiotechnol 2014) からn=40 samples (各3,000万リードのbootstrapサブセット) を繰り返しサブサンプリングし、ブートストラップによる分散推定値がサブサンプル間の真の分散とR=0.933 (p<0.001) で高相関することを確認した (Supplementary Fig. 5, 6)。ポアソン分布を仮定したリードカウントの射撃雑音 (Poisson variance proxy) は真の転写体発現量の不確実性の貧弱な代理指標であり、kallistoブートストラップが下流の差次発現解析など統計解析に組み込まれるべき誤差構造を明示できることが示された (Supplementary Figs. 7, 8)。
考察/結論
① 先行研究との違い: TopHat2+Cufflinksなどの既存パイプラインと異なり、kallistoはリードを個別塩基レベルでアライメントすることなく、k-merハッシュのみで等価クラスを決定する。Sailfishがリードをk-merへ分断することで情報を損失していたのに対照的に、kallistoの擬似アライメントは完全なリードのk-merを使用することで等価クラスの精度を維持する。等価クラスが十分統計量であるという理論的根拠により、アライメントと同等の最尤推定が実現し、速度と精度のトレードオフが解消された。
② 新規性: 転写体deブラインGrafen由来のk-compatibility classと「skipping」最適化を組み合わせた擬似アライメントは本研究で初めて提案された手法であり、多くのリードで2回のハッシュルックアップのみで等価クラスが決定できるという新規な計算的洞察に基づいている。ブートストラップによる不確実性定量の実用化も、高速なEM計算と等価クラスの十分統計量性を組み合わせた新規な応用である。
③ 臨床応用: 汎用ラップトップで大規模RNA-seqを処理できることで、臨床現場での RNA-seqトランスクリプトーム解析のコスト・時間・情報技術インフラ要件が劇的に低下する。ブートストラップ不確実性の定量化は差次発現解析の統計的厳密性を向上させ、診断バイオマーカーの信頼区間推定や希少疾患の小コホート研究への臨床的意義をもたらす。
④ 残された課題: 本研究はRNA-seqに特化しているが、著者自身が配列センサスデータセット全般への適用可能性を示唆している。k-merの長さの最適化は転写体特異性と誤差耐性のバランスによって制約されるため、今後の研究によって更に詳細な検討が必要である。またRNAアイソフォームが極めて高い類似性をもつパラログ遺伝子 (同一ファミリー内) ではkallistoを含む全手法の精度が低下するという制約は、今後の課題として残される。単細胞RNA-seqへの適用拡張やスプライシング多様性の直接定量なども更なる検討が求められる。
方法
アルゴリズム: (1) 転写体k-merからT-DBG構築、コンティグ同定 (2) k-merハッシュテーブル → k-compatibility class付与 (3) リードk-merの積集合で等価クラス決定 (pseudoalignment) (4) EMアルゴリズムで等価クラスカウントから転写体発現量を最尤推定 (5) TPM単位出力。ブートストラップ: N リードの等価クラスカウントから多項分布で復元サンプリング → EM再計算、ブートストラップ標本数はオプション指定、出力は HDF5 (hierarchical data format 5; 階層型データフォーマット5) で圧縮。バイアス補正: 5’端6-merの実測頻度で有効長を調整 (Roberts et al. 方法準拠)。ベンチマーク: 20 RSEM simulation (GEUVADIS NA12716_7; 30 million 75-bp paired-end reads) + SEQC/MAQCIII consortium実データ + qPCR独立検証 (Supplementary Table 1)。比較手法: TopHat2+Cufflinks / Bowtie2+RSEM / Bowtie2+eXpress / Bowtie2+EMSAR (efficient multisplit alignment and RNA quantification) / HISAT2+Cufflinks / Sailfish。計算環境: 20コア/20スレッド (TopHat2/Cufflinks/RSEM/Bowtie2) または20並列プロセス (eXpress/kallisto)。解析再現: Snakefile at GitHub (pachterlab/kallisto_paper_analysis)。