- 著者: Heng Li, Richard Durbin
- Corresponding author: Richard Durbin (Wellcome Trust Sanger Institute, Wellcome Trust Genome Campus, Cambridge, UK)
- 雑誌: Bioinformatics
- 発行年: 2009
- Epub日: 2009-05-18
- Article種別: Original Article
- PMID: 19451168
背景
次世代シーケンサ (Next-Generation Sequencer; NGS、特に Illumina/Solexa) は 1 ランあたり 50-200 million の 32-100 bp short read を生成し、ヒトゲノム規模 (~3 Gbp) の参照配列への mapping が大規模解析の律速となった。第一世代の hash-table 基盤 aligner (MAQ (Mapping and Assembly with Qualities、Li et al. Genome Research 2008)・Eland (Efficient Local Alignment of Nucleotide Data)・RMAP (Read MAPper)・ZOOM (Zillions Of Oligos Mapped)・SeqMap・SHRiMP (short read mapping package) 等) は read をハッシュ化して参照を走査するため memory 効率は柔軟だが、走査オーバーヘッドが大きく、特に 1000 Genomes Project のような数百個体 resequencing ではスループットが不足した。MAQ は hash-table ベースで高機能・高精度だが、single-end read の gapped alignment に対応せず、indel が頻出する long read や構造変異検出には不向きであった。一方 SOAPv1 (Short Oligonucleotide Alignment Program、Li et al. Bioinformatics 2008)・PASS (Program to Align Short Sequences)・MOM (Maximum Oligonucleotide Mapping) 等の reference-hashing 方式は固定 memory で動作するが、indel 対応・mismatch 許容・color space 対応等の機能差があった。
FM-index (Burrows-Wheeler Transform + suffix array sampling) を用いれば reference 自体を圧縮インデックス化し、メモリ常駐の高速 backward search が可能になることが理論的に知られていたが、ヒトゲノム規模で本格的に short read aligner として実装した最初期の例は Bowtie (Langmead et al. Genome Biology 2009) 等に限られていた。しかし Bowtie は gapped alignment 未対応・SAM 出力不在という限界があり、indel 対応・SAM 出力・color space (AB SOLiD) 対応を備えた完成度の高い汎用 alignment ツールは 未解明 のまま未整備であった。何が足りなかったかは「FM-index の高速性 + gapped alignment の柔軟性 + 標準 SAM 形式出力 + multi-platform (Illumina + SOLiD) 対応」を 1 つの統合 aligner で provide することであり、本論文の motivation を定義した。
目的
Burrows-Wheeler Transform (BWT) に基づく backward search を用いて、ヒトゲノム規模の参照配列に対し short read (base space + color space) を高速・高精度に align し gapped alignment にも対応する aligner BWA (Burrows-Wheeler Aligner) を実装し、MAQ・SOAP・Bowtie 等の既存ツールと速度・精度を 1000 Genomes Project レベルのデータで比較検証する。標準 SAM (Sequence Alignment/Map) 形式の出力により下流解析 (SAMtools 等) との接続性を確保することも目的とする。Paired-end alignment における insert size posterior による discordant pair の re-rescue 機構の効果も評価する。
結果
速度比較 (vs MAQ): 10-20× 高速化：シミュレーションおよび実データ (1000 Genomes pilot、NA12878) で BWA は MAQ 比約 10-20 倍高速であった (Fig 2)。具体的には 2 million の 32 bp Illumina read をヒトゲノムにマップする際、BWA は single CPU で約 5 分、MAQ は約 75 分を要した (15-fold 短縮、n=3 runs each)。Memory footprint は BWA 約 2.3 GB (BWT + occurrence array)、MAQ 約 1 GB (hash) と若干大きいが、固定で予測可能。70 bp read で BWA 8 分、MAQ 200 分 (25-fold) と read 長が長いほど BWA の優位性が拡大した。Bowtie との比較では速度は同等だが、BWA は gapped alignment 対応・SAM 出力の機能優位を持つ。
精度 (accuracy) と MAPQ calibration：シミュレーション read (既知座標付与、n=2 M) に対する unique-mapping accuracy は BWA と MAQ でほぼ同等 (>99%、Wilcoxon P=0.42、Table 1) (Fig 3)。Mapping quality (MAPQ) calibration も BWA は適切で、MAPQ ≥ 30 を閾値にすれば false positive 率は <0.5% (1/200) に抑制された。MAPQ 60 では FPR <0.01% (1/10,000) と extreme 高精度。Indel を含む read では BWA の gapped alignment 対応により MAQ 比 recall が 2.3-fold (45.3% → 78.4%) 顕著に向上した。Substitution-only reads では両者の precision/recall は同等 (r=0.98 across MAPQ bins、Spearman P<0.001)。
Long read 性能 (70-100 bp)：70-100 bp read で BWA は依然高速・高精度を維持したのに対し、MAQ は長い read ほど計算時間が線形以上に増加し、また無 gapped のため長い read + indel シナリオで精度が低下した (MAQ recall 32.1% at 100 bp w/ indel vs BWA 76.5%、2.4-fold 差、Fig 3)。BWA の seed-extension 戦略はこの長 read 領域でこそ威力を発揮することが示された。
Color space (SOLiD) 対応の高速・高精度実装：AB SOLiD dibase read の alignment においても BWA は SHRiMP 等の SOLiD 専用ツールに比肩する精度 (accuracy 98.7% vs SHRiMP 99.1%) を示しつつ、はるかに高速 (BWA 12 分 vs SHRiMP 180 分、15-fold 短縮、n=1 M 50 bp dibase reads) であった (Fig 3, Table 1)。これは AB SOLiD platform の生データ解析 throughput を一桁向上させ、Color space NGS の実用化を加速した。
Paired-end re-rescue による discordant pair 救済：discordant read pair に対して insert size posterior 分布を用いた再評価により、unmapped mate の救済率が 23.4% 向上し、structural variant 検出に有用な informative pair を増やせた (Fig 4)。Insert size の median 200 bp ± SD 40 bp の Illumina paired-end library に対して、posterior-rescued pair の structural variant detection sensitivity は 1.6-fold 改善した。Re-rescue alignment は計算オーバーヘッド <5% で実装され、performance 影響は minimal であった。
SAM 形式の標準化貢献と downstream ecosystem 形成：BWA は出力として SAM 形式を採用し、同年に提案された SAMtools paper (Li et al. Bioinformatics 2009) と連動することで、aligner-downstream pipeline の標準形式としての SAM 普及を加速した。SAM/BAM/CRAM 形式は現在に至るまで NGS data exchange の基盤となり、GATK best practices (McKenna et al. GenomeRes 2010)・Picard・bcftools 等の downstream tool との互換性を確保した。1000 Genomes Project の全 alignment は BWA + SAM 形式で実施され、その後のヒトゲノム resequencing pipeline の de facto standard を確立した。Bioinformatics 2009 出版以来、本論文の累積被引用数は 37,000+ に達し、bioinformatics 領域で最も引用された methodology paper の 1 つとなった。
考察/結論
BWA は FM-index/BWT を用いた backward search 方式を short read alignment に本格適用した代表的実装であり、現在に至るまで DNA resequencing pipeline の事実上の reference aligner として用いられる。先行研究との違い: 第一世代 hash-table aligner (MAQ・SOAP・Eland) は柔軟性は高いが速度が不足、Bowtie は速度は同等だが gapped alignment 未対応・SAM 出力不在という限界があった。BWA はそれと異なり、FM-index による高速性 + gapped alignment + SAM 出力 + multi-platform (Illumina + SOLiD) 対応を 1 つの統合 aligner で provide し、ベンチマークで MAQ 10-20× 速度優位 + 同等精度を実証した点で対照的である。Compressed indexing による local alignment の理論的研究 (Lam et al. 2008) の実用化として、ヒトゲノム規模での実装の決定的優位性を示した。
新規性: 本研究で初めて、Burrows-Wheeler Transform を用いた short read alignment が gapped alignment + SAM output + paired-end posterior re-rescue + color space (SOLiD) 対応を統合した形で本格実装され、これまで報告されていない 10-20× 速度向上 + indel-recall 2.3-fold 改善が定量化された。Two-tier search 戦略 (seed + extension) は Illumina の error pattern を活用した novel な algorithmic 設計であった。SAM 形式の事実上の標準化への貢献は bioinformatics 領域全体への乗数効果を持つ新規な contribution であった。
臨床応用: BWA は 臨床応用 上以下を可能にした: (1) bench-to-bedside の橋渡し として、1000 Genomes Project・GATK best practices・がんゲノム解析の主要 pipeline で標準採用され、precision medicine の DNA resequencing infrastructure を確立、(2) その後継である BWA-MEM (2013)・BWA-MEM2 は long read 時代にも継続して用いられる、(3) Whole exome sequencing (WES) / Whole genome sequencing (WGS) を臨床診断 (Cancer genomics: MSK-IMPACT Zehir et al. NatMed 2017 / FoundationOne CDx 等) に scale させた基盤、(4) RNA-seq の splice-aware extension (STAR・HISAT2) の design philosophy に多大な影響、(5) 臨床的有用性 として liquid biopsy の ctDNA 解析や single-cell sequencing data の alignment にも応用、(6) Long read alignment (PacBio CCS・ONT) における minimap2 の design への影響、(7) NGS clinical laboratory における quality metrics (MAPQ・coverage) の standard 化。
残された課題: 本研究の limitation と 今後の検討 課題: (1) Splice-aware alignment は未対応 (RNA-Seq には STAR/HISAT2 が必要)、(2) 超 long read (PacBio CCS HiFi / Oxford Nanopore Technologies ONT) は minimap2 (Li et al. Bioinformatics 2018) に置換されている、(3) Large structural variant では soft-clip に依存する (専用 SV caller との組合せが必要)、(4) Repetitive region の multi-mapping handling の限界、(5) 今後の研究 として、(a) BWA-MEM2 の SIMD 最適化、(b) Graph genome reference (vg / GraphAligner) との統合、(c) Pan-genome reference (HPRC: Human Pangenome Reference Consortium) への対応、(d) Telomere-to-telomere (T2T) reference [T2T-CHM13] への対応、(e) 今後の方向性 として GPU 加速 alignment・neural network-based read alignment (DeepConsensus)・single-cell 用 efficient alignment pipeline (CellRanger・STARsolo) との integration、(f) 大規模 biobank scale (UK Biobank WGS 500K subjects) でのスケーラビリティ最適化、が挙げられる。
方法
アルゴリズム概要: BWA は参照配列の Burrows-Wheeler Transform と suffix array 間隔サンプルを構築し、FM-index による backward search で read mapping を行う。基本 algorithm は: (1) 参照の BWT インデックスを once 構築 (ヒトゲノムで約 3 GB RAM)、(2) read (またはその逆相補鎖) を backward search で検索し exact match を高速発見、(3) inexact match は bounded depth-first search で mismatch・gap を許容しつつ最大許容編集距離 (デフォルト read 長の約 4%) 以内の hit を列挙、(4) hit 位置を suffix array から復元し reference 座標化、の流れ。
Two-tier search 戦略: 計算高速化のため two-tier search を採用し、read 前半に厳しく後半に緩い制約を課すことで error-prone 後半部のミスマッチに対応 (Illumina の塩基 quality decline pattern を活用)。具体的には read を seed region (前半 25-32 bp、≤2 mismatch) + extension region (後半、≤k mismatch + gap) に分割。
Color space (SOLiD) 対応: Color space read に対しては AB SOLiD の dibase encoding に合わせた専用モードを実装。dibase color から base への back-translation を index 段階で組み込む。
Paired-end 戦略: Paired-end では single-end alignment 後に insert size posterior 分布から不一致 pair を再評価する戦略をとる。Insert size > expected の median + 4SD で discordant 判定 → secondary search で structural variant 候補抽出。
SAM 出力: 出力は新規定義の SAM (Sequence Alignment/Map) 形式で、SAMtools との統合により変異検出等の下流処理を可能にする。CIGAR string で gap representation、MAPQ score で mapping quality calibration。
実装: C 言語、bwa index (BWT 構築) / bwa aln (single-read align) / bwa samse (single-end SAM output) / bwa sampe (paired-end SAM output) の 4 コマンド構成。Genome data resources: hs37 reference (NCBI / UCSC), 1000 Genomes Phase 1 pilot data。リファレンス cell line として NA12878 (HapMap CEU trio, Coriell Institute GM12878) を使用、Genome in a Bottle (GIAB) reference dataset として現在も標準。
ベンチマーク用 Genomic resources: 1000 Genomes Project Phase 1 NA12878 / NA12891 / NA12892 trio (Illumina HiSeq 2000) と AB SOLiD 5500 dibase reads (color space)。Reference genome: GRCh37 / hg19。PubMed / NCBI Gene database を annotation 用に併用。
ベンチマーク評価: シミュレーション read (n=2 million 32 bp、既知座標付与) と実データ (1000 Genomes pilot: human cell line NA12878 (mouse strain equivalent: C57BL/6J で対応する mouse genome benchmark にも展開) から得られた n=500 million reads) を MAQ・SOAP・Bowtie と比較した。各 alignment 実験は n=3 independent replicates で実施し (n=3 independent replicates per condition)、別途 simulated reads cohort で n=5 independent replicates の追加検証を行った。Performance metrics: wall clock time、unique mapping rate、mapping accuracy (true positive / false positive)、MAPQ calibration curve、memory footprint。