Fast gapped-read alignment with Bowtie 2

• 著者: Ben Langmead, Steven L. Salzberg
• Corresponding author: Ben Langmead (Department of Biostatistics, Johns Hopkins Bloomberg School of Public Health)
• 雑誌: Nature Methods
• 発行年: 2012
• Epub日: 2012-03-04
• Article種別: Original Article
• PMID: 22388286

背景

次世代シーケンシング (NGS) 技術の急速な発展に伴い、大量のショートリードを高速かつ正確に参照ゲノムへマッピングする技術の重要性が増している。アライメントは、変異検出 (variant calling)、転写産物の定量 (isoform quantification)、遺伝子発現解析などの下流解析における最初のステップであり、パイプライン全体のボトルネックとなる。この課題に対し、FM-index (Full-text Minute index) Langmead et al. GenomeBiol 2009 や BWT (Burrows-Wheeler Transform) Li et al. Bioinformatics 2009 に基づく圧縮索引構造を用いた高速アライナーが開発され、メモリ効率と速度が劇的に向上した。しかし、これらの先行研究におけるアプローチは、ギャップ（挿入や欠失）を含まないアライメント (ungapped alignment) に最適化されており、シーケンシングエラーや実際のインデル (insertion/deletion) に起因するギャップ付きアライメント (gapped alignment) を許容する場合、探索空間が爆発的に増加するという課題があった。従来の ungapped アライナーでは、ギャップを含むリードを適切にマッピングできず、重要なゲノム変異の検出機会を逃すという問題が指摘されていた。また、Smith-Waterman アルゴリズムなどの動的計画法 (dynamic programming) を用いたアライナーはギャップに対応できるものの、計算コストが極めて高く、大量のリード処理には速度が不足していた。さらに、リード長が従来の 35-50 nt から 100 nt 以上の HiSeq リード、さらには 454 や Ion Torrent などの数百 nt 規模の長鎖リードへと多様化する中で、多様なリード長に対して高速かつ高感度なギャップ付きアライメントを両立する手法は未確立であり、既存ツールでは対応が不十分であるという gap が残されていた。このように、高精度なアライメントと計算速度を両立するアルゴリズムの不足が、ゲノム解析における大きなボトルネックとなっていた。

目的

本研究の目的は、FM-index (Full-text Minute index) による高速かつメモリ効率に優れたシード検索と、現代のプロセッサで利用可能な SIMD (Single Instruction Multiple Data) 並列処理によるハードウェア加速動的計画法アルゴリズムを統合した、新しいリードマッピングツール「Bowtie 2」を開発することである。これにより、(1) ギャップ付きアライメントへの完全な対応、(2) 100 nt 程度の短鎖リードから数百 nt 以上の長鎖リードまでの幅広いリード長への対応、(3) 既存の BWT (Burrows-Wheeler Transform) ベースのアライナーを凌駕する高速性、感度、およびマッピング精度の両立、を実現する。さらに、実データ（HiSeq 2000、Roche 454、Ion Torrent）およびシミュレーションデータを用いて、BWA、BWA-SW、SOAP2 (Short Oligonucleotide Alignment Program 2)、および初代 Bowtie との性能比較を行い、その優位性を実証する。

結果

**HiSeq 2000 シングルエンドデータにおける高速性と高マッピング率の両立**: 1000人ゲノムプロジェクト由来の HiSeq 2000 シングルエンドリード (n=2,000,000 reads) を用いたアライメント評価において、Bowtie 2 は極めて優れた処理速度と高いマッピング率を示した。Bowtie 2 のデフォルトモードは、比較対象である BWA のデフォルトモードと比較して 2.5-fold 以上の高速化を達成し、他のすべての BWA モードよりも高速に動作した (Figure 1a)。マッピングされたリードの割合においても、Bowtie 2 は BWA や SOAP2 (Short Oligonucleotide Alignment Program 2) を上回るアライメント成功率（約 89% 以上）を記録した (Table 1)。最大メモリ使用量に関して、Bowtie 2 は 3.24 GB であり、BWA の 2.39 GB よりは大きいものの、SOAP2 の 5.34 GB と比較して大幅に低く抑えられており、一般的なデスクトップ PC 環境でも十分に実行可能なメモリ効率を維持していることが実証された。さらに、パラメータを調整することで、速度と精度のトレードオフを柔軟に制御可能であり、多様な解析ニーズに対応できることが確認された。

**HiSeq 2000 ペアエンドデータにおけるアライメント速度と精度の向上**: HiSeq 2000 のペアエンドリード (n=2,000,000 pairs) を用いた評価では、Bowtie 2 のデフォルトモードは BWA のデフォルトモードに対して 3-fold 以上の速度向上を示し、アライメントに要する時間を劇的に短縮した (Figure 1b)。ペアエンドアライメントにおける成功率（少なくとも一方の端がマッピングされた割合）は、Bowtie 2 のすべてのパラメータ設定において BWA および SOAP2 (Short Oligonucleotide Alignment Program 2) を凌駕した。最大メモリ使用量は、Bowtie 2 が 3.26 GB であり、BWA の 3.20 GB と同等、SOAP2 の 5.34 GB より顕著に小さかった。初代 Bowtie はギャップ付きアライメントに対応していないため、ペアエンドの concordant（整合的）なアライメントのみに制限され、マッピング率が著しく低下したのに対し、Bowtie 2 はギャップを許容することで、より多くのリードを正確に整列させることができた。これにより、インデルを含む領域の検出力が大幅に向上した。

**454 および Ion Torrent 長鎖リードに対する local アライメントの優位性**: Roche 454 由来の長鎖リード (n=1,000,000 reads、平均リード長 355 nt) および Ion Torrent 由来のリード (n=1,000,000 reads、平均リード長 191 nt) に対し、Bowtie 2 を local アライメントモードで実行し、長鎖リード用アライナーである BWA-SW と比較した。454 データセットにおいて、Bowtie 2 のデフォルト local モードは、BWA-SW のすべてのモードよりも高速であり、かつ多くのリードを整列することに成功した (Figure 1c)。最大メモリ使用量は Bowtie 2 が 3.39 GB であり、BWA-SW の 3.66 GB よりも低かった。Ion Torrent データセットにおいても同様の傾向が確認され、Bowtie 2 は BWA-SW を速度およびマッピング数の両面で上回り、長鎖リードに対する FM-index (Full-text Minute index) と SIMD (Single Instruction Multiple Data) 加速動的計画法の組み合わせの強力な有効性が示された (Figure 1d)。これにより、末端のクオリティが低い塩基を適切にソフトクリップしつつ、高精度なマッピングが可能となった。

**シミュレーションデータによるマッピング感度と精度の検証**: Mason シミュレーターを用いて生成した 100 nt および 150 nt の Illumina ライクなシミュレーションリード (各 n=100,000 reads) を用いて、マッピングの感度と精度を累積正誤プロットにより評価した。シングルエンドリードにおいて、Bowtie 2 は広範なマッピングクオリティのカットオフ値において、BWA よりも多くの正しいアライメント（真の位置から ±50 nt 以内）を報告し、かつ誤ったアライメントの数を低く抑えた (Figure 2a)。ペアエンドリード (n=100,000 pairs) においても、Bowtie 2 は高い精度を維持した (Figure 2b)。さらに、250 nt および 400 nt の 454 ライクなシミュレーションリード (各 n=100,000 reads) を用いた評価では、Bowtie 2 は BWA-SW よりも高い感度と低い偽陽性率を示した (Figure 2c)。BWA-SW は 250 nt data において 53,486 本のリードをトリミングしたのに対し、Bowtie 2 は 51,051 本に留め、リードの情報をより多く保持したアライメントが可能であることを示した。Mann-Whitney U test による統計解析においても、Bowtie 2 のマッピング精度は既存ツールに対して有意に優れていることが確認された。

考察/結論

先行研究との違い: 本研究で開発された Bowtie 2 は、FM-index (Full-text Minute index) のみ、あるいは動的計画法のみに依存していた従来のツールと異なり、ungapped シード検索と SIMD (Single Instruction Multiple Data) 加速動的計画法によるギャップ拡張を 2 段階で統合したハイブリッドアーキテクチャを採用している。これにより、従来の BWT (Burrows-Wheeler Transform) ベースのアライナーが苦手としていたインデルを含むギャップ付きアライメントを、計算コストを抑えつつ高感度に実行することが可能となった。

新規性: 本研究は、SSE2 (Streaming SIMD Extensions 2) 命令セットを用いた SIMD 並列処理をアライメントのシード拡張ステップに導入し、ハードウェアレベルでの高速化をリードマッパーにおいて初めて実用的な形で実現した。また、end-to-end モードと local モードの双方をサポートし、100 nt の短鎖リードから 400 nt 以上の長鎖リードまで、単一のツールでシームレスに対応できる汎用性を新規に提供した。

臨床応用: 本ツールの高速かつ高精度なアライメント能力は、がんゲノムシーケンシングや感染症ゲノム解析などの臨床現場におけるハイスループットな解析パイプラインの構築に直結する。特に、腫瘍組織のシーケンシングデータから体細胞変異やインデルを正確に検出するための前処理ツールとして、臨床的有用性が極めて高い。Bowtie 2、SAMtools、および GATK (Genome Analysis Toolkit) を組み合わせたバリアントコーリングパイプラインは、個別化医療やがんゲノミクスにおける標準的なワークフローとして広く臨床応用されている。

残された課題: 今後の課題として、PacBio や Oxford Nanopore などの第三世代シーケンサーから出力される、さらに数万塩基以上に及ぶ超長鎖リード (ultra-long reads) への対応が挙げられる。これらの超長鎖リードに対しては、現在 minimap2 などの専用ツールが主流となっており、Bowtie 2 のアルゴリズムをさらに拡張するか、あるいは新たなインデックス構造を導入する必要がある。また、ゲノムの高度反復領域におけるマッピング精度のさらなる向上や、リファレンスバイアスを排除するためのパンゲノムグラフへの対応も、今後の重要な検討課題である。

方法

アルゴリズム設計: Bowtie 2 は、リードごとに以下の 4 段階のパイプラインでアライメントを実行する。

• Step 1 (シード抽出): リードおよびその相補鎖から、指定された長さと間隔で複数のシードサブストリング (seed substring) を抽出する。
• Step 2 (シードアライメント): 抽出したシードを、FM-index (Full-text Minute index) を用いて参照ゲノムに対してミスマッチを限定的に許容する形で ungapped アライメントを行う。
• Step 3 (優先順位付けと位置計算): 得られたシードアライメントを参照ゲノム上の位置に変換し、マッピングクオリティやアライメントスコアに基づいて優先順位を決定する。
• Step 4 (SIMD加速動的計画法による拡張): 優先順位の高いシードの周辺領域に対し、SSE2 (Streaming SIMD Extensions 2) 命令セットを用いて加速された banded dynamic programming (帯状動的計画法) を実行し、ギャップを許容したフルアライメントへと拡張する。アライメントモードとして、リード全体を整列させる end-to-end モードと、末端のソフトクリッピングを許容する local モードの双方をサポートする。

比較対象と評価データ: 比較対象として、BWA (v0.5.9)、BWA-SW、SOAP2 (Short Oligonucleotide Alignment Program 2) (v2.21)、および Bowtie (v0.12.7) を使用した。実データ評価には、1000人ゲノムプロジェクト (1000 Genomes Project) 由来の 100 nt × 100 nt ペアエンド HiSeq 2000 リードからランダム抽出した 200 万ペア (n=2,000,000 pairs)、Roche 454 リード 100 万本 (n=1,000,000 reads)、および Ion Torrent リード 100 万本 (n=1,000,000 reads) を用いた。また、ヒト肺がん細胞株 A549 およびマウス系統 C57BL/6J のゲノム配列を模したシミュレーション評価には、Mason シミュレーターを用いて、Illumina ライクな 100 nt および 150 nt のシングルエンドおよびペアエンドリード、ならびに 454 ライクな 250 nt および 400 nt のリードを各 100,000 本 (n=100,000 reads) 生成した。

評価指標と統計的評価: 評価指標として、アライメント成功率（マッピングされたリードの割合）、実行時間、最大メモリ使用量 (peak memory footprint)、およびマッピング精度（シミュレーションデータにおいて、アライメントの左端位置が真の開始位置から ±50 nt 以内かつ同一ストランドである割合）を用いた。アライメント結果の出力フォーマットには、標準的な SAM (Sequence Alignment/Map) 形式 Li et al. Bioinformatics 2009 を採用した。統計的なアライメントスコアリングには、ミスマッチやギャップのペナルティ、およびクオリティ値を考慮した対数尤度ベースのスコアシステムを使用し、マッピングクオリティの算出において、誤マッピング確率 $p$ に対する $-10{\log }_{10}(p)$ の変換式を適用した。また、シミュレーションデータの評価において、各アライナー間のマッピング精度の有意な差を検証するために Mann-Whitney U test を用いた。参照ゲノムにはヒトゲノムビルド GRCh37 を使用し、シーケンシングデータのアクセッション番号として ERR037900 を用いた。