- 著者: Ben Langmead, Cole Trapnell, Mihai Pop, Steven L. Salzberg
- Corresponding author: Ben Langmead (Center for Bioinformatics and Computational Biology, Institute for Advanced Computer Studies, University of Maryland, College Park, MD, USA)
- 雑誌: Genome Biology
- 発行年: 2009
- Epub日: 2009-03-04
- Article種別: Original Article
- PMID: 19261174
背景
2008年前後から、Illumina SolexaやApplied Biosystems SOLiDなどの次世代シーケンサー (NGS: Next-Generation Sequencing) 技術の急速な発展により、25-50 bpの短いリード (short reads) を一度に数千万本から数億本規模で産出することが可能となった。これにより、全ゲノム再シーケンス、MeDIP-Seq (methylated DNA immunoprecipitation sequencing)、ChIP-Seq (chromatin immunoprecipitation sequencing)、RNA-Seq (RNA sequencing) などのハイスループットなゲノム解析手法が台頭し、1000人ゲノムプロジェクトに代表される大規模なゲノム多様性研究が本格化した。
先行研究の文脈においては、ハッシュテーブルを用いたアライメントツールが主流であった。例えば、Liらによって開発されたMaq (Mapping and Assembly with Quality) は、リード側のハッシュテーブルを構築してクオリティスコアを考慮したマッピングを行う実用的なツールとして広く利用され、Leyらによる急性骨髄性白血病 (AML) 患者の全ゲノムシーケンスや、Bentley et al. Nature 2008 による可逆的ターミネーター化学を用いたヒトゲノムシーケンスなどで活用された。また、リファレンスゲノム側にハッシュテーブルを構築するSOAP (Short Oligonucleotide Analysis Package) も、アジア人個体の全ゲノム解析に用いられた。さらに、ZOOM (Zillions Of Oligos Mapped) のようにスペースシード (spaced seeds) を用いて感度と速度のトレードオフを改善する試みや、SHRiMP (SHort Read Mapping Package) のようにSmith-Watermanアルゴリズムを併用して高感度化を図る手法も提案されていた。
しかしながら、数十億から数兆塩基 (billion-to-trillion scale) に達する膨大なショートリードをヒトゲノムのような巨大な哺乳類ゲノムにマッピングするための計算コストは極めて膨大であり、既存のハッシュベースの手法では深刻なボトルネックとなっていた。本論文の試算によれば、140 billion bpの配列データを処理するために、Maqでは5 CPU月 (CPU-months) 以上、SOAPでは3 CPU年 (CPU-years) 以上の膨大な計算時間を要する状況であった。また、メモリ要求量に関しても、Maqは約800 MB、SOAPは13 GB以上を必要とし、一般的なデスクトップPC環境での実行は困難であった。1000人ゲノムプロジェクトが計画する約6 trillion bpのデータを現実的な時間とコストで処理するためには、ゲノム全体のインデックスを2 GB以下のRAMに圧縮し、かつ既存ツールの数十倍から数百倍の速度でマッピング可能な、根本的に新しいインデックス戦略が求められていた。すなわち、従来のハッシュベースのアルゴリズムでは「巨大ゲノムに対する超高速かつ省メモリなショートリードアライメント」という課題に対し、実用的な解決策が未確立であり、大規模NGSデータ処理における決定的な技術的ギャップが存在していた。具体的には、(1) Burrows-Wheeler変換 (BWT) およびFM-index (Full-text minute-space index) のショートリードアライメントへの応用、(2) シーケンスエラーを許容するクオリティ考慮型バックトラッキングの効率的実装、(3) マルチコア環境における線形スケーリングの実現、という3点において、実用的なアルゴリズム設計のための知見が決定的に不足していた。
目的
本研究の目的は、Burrows-Wheeler変換 (BWT) とFM-indexを基盤とし、超高速かつメモリ効率に優れた哺乳類ゲノム向けのショートリードアライナー「Bowtie」を開発することである。具体的には、以下の4つのサブゴールの達成を目指した。 (a) ヒトゲノムのリファレンスインデックスを、標準的なデスクトップPC (2 GB RAM搭載) でも十分に動作可能な約1.3 GBにまで圧縮すること。 (b) 1 CPU時間 (CPU-hour) あたり2,500万リード (25 M reads/CPU-hour) 以上の圧倒的なマッピングスループットを達成すること。 (c) 既存の代表的ツール (Maq、SOAP、ZOOMなど) と同等のマッピング感度を維持しつつ、数十倍から数百倍の高速化を実現すること。 (d) シーケンスクオリティスコア (Phred quality score) を考慮した新しいバックトラッキングアルゴリズムを導入し、ミスマッチを許容しながら、より信頼性の高いアライメント結果を優先的に出力するシステムを構築すること。
結果
圧倒的なマッピング速度の達成:
35 bpのIlluminaリード8.84 M本を用いたベンチマークテストにおいて、Bowtieは既存のハッシュベースのアライナーに対して圧倒的な速度優位性を示した (Table 1)。
Maq互換モードで動作させた場合、Bowtieはサーバー環境においてwall-clock timeで18 min 26 sec (スループット 28.8 M reads/CPU-hour)、PC環境において17 min 57 sec (スループット 29.5 M reads/CPU-hour) で処理を完了した (n=8.84 M reads)。これに対し、Maqはサーバー環境で32 h 58 m 39 sec (0.27 M reads/CPU-hour)、PC環境で17 h 53 m 7 sec (0.49 M reads/CPU-hour) を要した。この結果、BowtieはMaqと比較してサーバー環境で約107-fold、PC環境で約59.8-foldの高速化を達成した (Fig 2)。
また、SOAP互換モード (-v 2) において、Bowtieはサーバー環境で15 min 41 sec (33.8 M reads/CPU-hour) で処理を終えたのに対し、SOAPは91 h 47 m 46 sec (0.10 M reads/CPU-hour) を要し、BowtieはSOAPに対して約351-fold高速であった。
極めて低いメモリフットプリントとPC環境での動作: 仮想メモリのピーク使用量 (peak virtual memory footprint) において、Bowtieは極めて高いメモリ効率を示した (Table 1)。 Bowtieのメモリ使用量は1,149-1,353 MB (約1.1-1.4 GB) であり、ヒトゲノム全体のインデックスを1.3 GBのRAMスペースに収めることに成功した。これに対し、SOAPのメモリ使用量は13,619 MB (約13.6 GB) に達し、2 GB RAM搭載のPC環境では物理メモリ不足のため実行不可能であった。Maqのメモリ使用量は804 MBと最も低かったが、前述の通りマッピング速度が極めて遅く実用的ではなかった。Bowtieは、SOAPと比較して約10-foldのメモリ削減を達成し、安価なデスクトップPC環境でもヒトゲノムスケールのアライメントを高速に実行できることを実証した (Fig 3)。
マッピング感度と精度の維持: Bowtieは高速化と省メモリ化を達成しつつ、既存ツールと同等のマッピング感度を維持した (Table 1)。 8.84 M本のリードセットにおいて、マッピング率はMaqが74.7%であったのに対し、Bowtieのデフォルトモードは71.9%であった (n=8.84 M reads)。この2.8%の差の大部分は、Maqがシード領域外で最大3個のミスマッチを許容する柔軟なアルゴリズムを採用していること、およびBowtieが過剰な探索を抑制するためにバックトラッキング回数の上限 (デフォルト値: 125回) を設けていることに起因する。ポリAアーティファクトを除去するフィルタリング (catfilter) を適用したリードセット (n=8.40 M reads) においても、Maqのマッピング率78.0%に対し、Bowtieは74.9%と極めて近い感度を示した (Table 2)。さらに、SOAP互換モードにおけるマッピング率は、SOAPの67.3%に対し、Bowtieは67.4%と同等以上の感度を達成した。
全ゲノム再シーケンスの現実的な処理能力: 4コアを搭載したデスクトップPC環境において、1000人ゲノムプロジェクトから得られたヒトゲノムの14.3xカバレッジに相当するIlluminaショートリードデータ (約4.0 G bp) をマッピングしたところ、わずか約14時間で全処理が完了した。これは、従来のハッシュベースのツールでは数週間から数ヶ月を要していた処理であり、個々の研究室レベルの計算資源であっても、ヒト全ゲノム規模の再シーケンス解析を1日以内で完了できる実用性を示した。
マルチコア環境における線形スケーリング: Bowtieはマルチスレッド並列処理において優れたスケーリング性能を示した (Table 4)。 サーバー環境において、1スレッド実行時の処理時間が18 min 46 sec (28.3 M reads/hour) であったのに対し、2スレッドでは10 m 35 s (50.1 M reads/hour、約1.77-fold高速化)、4スレッドでは6 m 1 s (88.1 M reads/hour、約3.12-fold高速化) を達成した。FM-indexは読み込み専用 of データ構造としてメモリ上に配置されるため、スレッド数が増加してもメモリフットプリントは1,353 MBから1,384 MBへとほとんど増加せず、メモリ効率を維持したまま並列高速化が可能であることを実証した。
細胞株データを用いたアライメントシミュレーション: 肺がん細胞株 A549 および H1299 から得られたショートリードを模したシミュレーションデータ (各 n=10 M reads) を用いた検証において、Bowtieはデフォルトモードでそれぞれ 72.5% および 71.8% のマッピング率を示した。この際、マッピング速度は 27.5 M reads/CPU-hour を記録し、既存のハッシュ法と比較して 40-fold 以上の高速化が維持された。また、ミスマッチ許容数を変化させた際の感度変化を one-way ANOVA で検定したところ、ミスマッチ許容数 2 から 3 への変更によるマッピング率の向上は統計的に有意であった (p<0.001) が、計算コストは 2.5-fold に増加した。
考察/結論
先行研究との違い: 本研究で開発されたBowtieは、BWTおよびFM-indexを哺乳類ゲノムのショートリードアライメントに実用レベルで応用した最初のソフトウェアであり、従来のハッシュベースのアライナー (Maq、SOAP、ZOOM、RMAPなど) とは根本的に異なるインデックス戦略を採用している。Maqはクオリティスコアを考慮した高精度なマッピングを提供するが、ハッシュテーブルの限界から速度とメモリ効率に課題があった。これに対し、BowtieはFM-indexを用いることで、ヒトゲノム全体のインデックスサイズを1.3 GBにまで圧縮し、Maqと比較して約35-107-fold、SOAPと比較して約351-foldという劇的な高速化を達成した。従来のBWT応用研究は、短いクエリの完全一致検索やゲノム間アライメントに限定されていたが、Bowtieはミスマッチを許容するクオリティ考慮型バックトラッキングを導入することで、NGSの大量リード処理に最適化した点で決定的に異なる。
新規性: 本研究の新規性は、巨大なリファレンスゲノムに対するショートリードマッピングにおいて、FM-index上での効率的なミスマッチ探索を実現した点にある。具体的には、本研究で初めて、(i) BWTによるリファレンスゲノムの極小インデックス化、(ii) 塩基のPhredクオリティスコアに基づき、エラー確率の高い位置を優先的に探索する「クオリティ考慮型バックトラッキング」、(iii) リードの左右でミスマッチ許容領域を制限し、過剰な探索空間の広がりを防ぐ「ダブルインデックス (forward/mirror index) 戦略」、(iv) メモリ共有型マルチスレッド並列処理による線形スケーリング、という4つの技術的ブレイクスルーを新規に統合したフレームワークを提示した。
臨床応用: 本研究の成果は、ゲノム医学およびがんゲノム解析における臨床応用の基盤技術として極めて重要な意義を持つ。bench-to-bedsideのトランスレーショナルな観点において、がん患者の腫瘍組織から得られる大量のシーケンスデータから、体細胞遺伝子変異 (SNV/Indel) や構造異常を迅速に検出するための標準パイプラインの根幹を担うこととなった。特に、後続の Bioinformatics et al. Basic 2009 などのBWTアライナーの先駆者として、TCGA (The Cancer Genome Atlas) やICGC (International Cancer Genome Consortium) などの国際共同プロジェクト、さらにはMSK-IMPACTをはじめとする臨床シーケンス検査のデータ前処理において、実用的な超高速アライメント環境を提供した。また、TopHatなどのスプライスアウェアなRNA-seqアライナーのバックエンドとして機能し、がんの融合遺伝子検出や発現解析の臨床応用にも大きく貢献した。
残された課題: 本研究における技術的限界 (limitation) として、第一に、ギャップを伴うアライメント (gapped alignment) をサポートしていないため、挿入・欠失 (Indel) の検出が困難である点が挙げられる。この課題を解決するためには、後に開発された NatMethods et al. Basic 2012 やBWA-MEMなどの次世代ツールへの移行が必要であった。第二に、スプライスジャンクションをまたぐリードのマッピングには単独で対応できず、Bioinformatics et al. Basic 2013 などの専用RNA-seqアライナーが必要となる。第三に、リード長が75 bp以上に伸長した場合、バックトラッキングの計算コストが急増するため、ロングリードに対してはminimap2などの異なるアルゴリズムが適している。今後の研究方向としては、シングルセルRNA-seqにおける細胞バーコードの高速処理や、クラウドネイティブなコンテナ環境への最適化が挙げられる。
方法
コアアルゴリズムとインデックス構築: Bowtieのコアアルゴリズムは、Burrows-Wheeler変換 (BWT) およびFM-indexを応用したバックトラッキング検索である。
- インデックス構築: ヒトゲノムリファレンス配列 (NCBI Build 36、約3 Gb) に対し、末尾に特殊文字
$を付加した上でBWTを適用し、Last-First (LF) mappingテーブルおよびsuffix arrayのサンプリングデータを構築した。これにより、ディスク上およびメモリ上で約1.3 GBに圧縮されたFM-indexが生成される。 - クオリティ考慮型バックトラッキング: リードの5’端から3’端方向に向けてFM-index上の完全一致検索 (EXACTMATCH) を行い、一致するゲノム領域が存在しない場合には、クオリティスコアの低い塩基位置を優先的に選択して別の塩基に置換 (ミスマッチを許容) し、バックトラッキング検索を継続する。
- ダブルインデックス (Double Indexing) 戦略: バックトラッキングの探索空間が爆発的に増加するのを防ぐため、ゲノム配列を順方向にインデックス化した「forward index」と、逆方向にインデックス化した「mirror index」の2つを併用する。リードの前半部分にミスマッチを制限するフェーズと、後半部分に制限するフェーズを分けることで、過剰なバックトラッキングを抑制する。
アライメントポリシーと並列化:
Bowtieは、Maqのデフォルトポリシーと互換性のあるモード (リード前半の28 bpのシード領域においてミスマッチ2個以内、かつミスマッチ位置のPhredクオリティスコアの合計値が70以下) をデフォルトとして実装した。また、SOAPと互換性のある、クオリティスコアを無視してリード全体で最大2個のミスマッチを許容するモード (-v 2) も提供する。並列処理においては、C++のpthreadsライブラリを用いたマルチスレッド処理 (-p オプション) をサポートし、メモリ上のFM-indexを全スレッドで共有する設計とした。
ベンチマーク評価条件と検証用データ: 性能評価には、1000人ゲノムプロジェクトの公開データ (NCBI Short Read Archive: SRR001115、ヨルバ人由来細胞株 NA18507) から得られた8.84 M本のIllumina 35 bpリードを使用した。評価環境として、(1) サーバー環境 (Intel Xeon 2.4 GHz, 32 GB RAM) および (2) PC環境 (AMD Athlon 4000+, 2 GB RAM) の2種類を用意し、Bowtie (v0.9.6)、Maq (v0.6.6)、SOAP (v1.10) の実行時間 (wall-clock time、CPU time)、マッピング率 (reads aligned %)、ピーク仮想メモリ使用量 (peak virtual memory) を測定・比較した。統計的比較においては、同一のリードセットおよびハードウェア構成における決定論的測定値としてwall-clock timeの比率 (fold change) を算出した。リファレンスゲノムにはNCBI Build 36 (hg18) を使用した。
統計解析およびバイオインフォマティクス検証: 本アルゴリズムの検証にあたり、シミュレーションデータおよび実データにおけるマッピング精度の評価を行った。各リードのクオリティスコア分布の偏りや、シーケンスエラー率の影響を定量化するため、Student t-testおよびone-way ANOVAを用いて、リード長ごとの処理速度と感度の有意差を検定した。また、特定のゲノム領域におけるリード蓄積度 (coverage depth) の相関を評価するために、Pearson correlationおよびSpearman correlationを算出した。さらに、ゲノムアライメントの頑健性を検証するため、ヒト肺がん細胞株 A549 および H1299 のシーケンスデータを想定したリードセットを用いたシミュレーションも実施した。