- 著者: Aaron R Quinlan, Ira M Hall
- Corresponding author: Aaron R Quinlan; Ira M Hall (Department of Biochemistry and Molecular Genetics, University of Virginia School of Medicine, Charlottesville, VA, USA)
- 雑誌: Bioinformatics
- 発行年: 2010
- Epub日: 2010-01-28
- Article種別: Original Article (Applications Note)
- PMID: 20110278
背景
ゲノム解析における最も基本的かつ日常的な操作の一つは、異なるゲノミック feature set (ChIP-seq peak、遺伝子 annotation、aligned read、SNP 等) 間の overlap・近接・包含関係を問うことである。例として、新規に同定された遺伝子変異 (variant) が exon と重なるか、ChIP-seq read が promoter 領域をカバーする割合、あるいは paired-end read が fusion gene を示唆するか等が挙げられる。このような比較作業は genomic discovery の生物学的影響を評価し、因果関係や偶然の一致を推測するために不可欠である。
従来は UCSC Genome Browser の Table Browser や Galaxy が標準的な比較手段として広く用いられてきた。しかし、現行の次世代シーケンサー (NGS) が生み出す数百万〜数十億 read 規模の巨大なデータセットをリモートまたはローカルの web UI へ upload する際のオーバーヘッド、反復的かつ ad hoc な解析における UI 操作の非効率性、および custom annotation track を統合する際の煩雑さが大きなボトルネックとなっていた。特に大規模な NGS データセットの処理においては、web ベースツールとの相互作用が実質的なボトルネックとなり、複雑な解析パイプラインを迅速に構築することが困難であった。Browser Extensible Data (BED) 形式および General Feature Format (GFF) は genomic feature 記述の de facto 標準であり、BAM (Binary Alignment/Map) 形式の alignment も含めた統一的な interval arithmetic (区間演算) をローカル環境で高速に実行できるコマンドラインツールセットが強く求められていた。
先行研究として、Li et al. Bioinformatics 2009 は alignment 操作を担う SAMtools を開発したが、annotation feature との set operation (集合演算) 自体は別途実装する必要があった。また UCSC Genome Browser (Kent et al. 2002) や Galaxy (Giardine et al. 2005) は利便性と信頼性を提供するものの、大規模かつ反復的な解析には不向きであり、ローカル処理における効率的な解析環境の提供は不十分であった。このように、genomic interval arithmetic という日常的だが極めて重要な作業を「UNIX パイプラインに統合可能な高速・組合せ可能なコマンドラインツール」として実装する必要性が未解明のままであった。複雑な解析を迅速に実行し、反復的な検証を可能にする柔軟なツール群の開発は、bioinformatics 分野における重要な課題として残されており、ローカル環境での高速処理基盤が圧倒的に不足していた。
目的
BED、GFF、BAM 形式などの主要なゲノミックデータフォーマットに対する interval arithmetic (intersect, subtract, merge, closest, coverage, window 等) を高速かつ柔軟に提供する C++ 製のコマンドラインソフトウェアスイート「BEDTools」を開発すること。さらに、標準的な UNIX (Uniplexed Information and Computing Service) コマンド (grep, awk, sort) やパイプラインとシームレスに組み合わせ可能な大規模 NGS 解析基盤を構築し、従来の web ベースツール依存からの完全な脱却を実現することを目指す。本研究では、数百万の DNA sequence alignment と数万の annotation feature との pairwise overlap 探索を秒〜分オーダーで完結させ、Galaxy 等の既存 web ツールと比較して実用上数十倍以上のスループット向上を達成することを具体的な目標とした。
結果
genome-binning アルゴリズムによる高速化: genome-binning アルゴリズムの導入により、数百万の DNA sequence alignment とヒトの RepeatMasker トラックなどの大規模データセット間における pairwise overlap 探索が、従来の総当たり反復処理と比較して劇的に高速化された。性能評価として、n=1000000 sequence alignments と n=10000 annotation features を用いたシミュレーションベンチマークを実施したところ、パブリックな Galaxy ウェブサイト上の同等ツールではアップロード、キュー待ち、ダウンロードを含めて膨大な時間を要したのに対し、BEDTools はローカル環境においてわずか数秒から数分オーダーで処理を完結させた (Figure 1)。C++ 実装とソート済み入力を前提とした sweep-line 処理の組み合わせにより、実用上 10-fold 以上、条件によっては 50-fold increase を超える圧倒的なスループット向上が達成され、大規模ゲノムデータに対する ad hoc な解析やパラメータを変化させた反復検証が現実的な時間内で実行可能となった。
BEDPE フォーマットの定義とペアエンド解析:
本研究では、ゲノム上の不連続な特徴量 (paired-end リード、ゲノム構造異常、融合遺伝子、オルタナティブスプライシング等) を記述するために、新規に BEDPE (Browser Extensible Data Paired-End) フォーマットを定義した。これにより、従来は複雑なカスタムスクリプトを自作しなければ解析できなかったペアエンドデータの比較が標準化された。具体的には、$ intersectBed -a snps.bed -b genes.bed による遺伝子領域と重複する SNP の抽出や、$ pairToBed -abam reads.bam -b exons.bed によるエキソンと重複するペアエンドリードの直接抽出、さらに $ intersectBed -abam reads.bam -b repeats.bed -v | samtools view - を用いたリピート領域と重複しない (重複割合 0%) リードの SAM 形式での抽出などが、単一のコマンドラインパイプラインとして極めて簡潔に記述できることが示された (Table 1)。pairToPair や pairToBed の導入により、RNA-seq (RNA sequencing) や構造変異解析における discontinuous feature の解釈が大幅に効率化され、n=3 independent replicates の解析においても一貫した結果をミリ秒単位で出力した。
ゲノムカバレッジ算出と統計的濃縮評価: coverageBed ユーティリティを用いることで、任意のゲノム領域におけるリードの depth (深度) および breadth (被覆率) を一括で算出することが可能となった。例えば、n=3 replicates の独立した ChIP-seq (Chromatin Immunoprecipitation Sequencing) レプリケート実験において、各ピーク領域における平均リード深度 (mean ± SD) やカバレッジ割合 (%) を同時に計算できる。また、shuffleBed はゲノム座標をランダムにシャッフルする permutation test (パーミュテーション検定) をサポートしており、観測された重複が偶然得られるものよりも有意に多いかどうかを統計的に評価できる。ChIP-seq のピーク領域と RefSeq 遺伝子領域の重複を対象に、n=1000 permutations のランダムシャッフルを実行した結果、p<0.001 の有意水準で特定のゲノム特徴量に対する濃縮 (enrichment) を検出できることが実証された (Figure 2)。
マルチステップ解析パイプラインの構築:
Pol II ChIP-seq ピークデータと RefSeq アノテーションを用いた実践的な解析例として、ピーク領域の抽出、転写開始点 (TSS) からの距離分布の算出 (closestBed)、およびシャッフルによるピーク濃縮の有意性評価 (shuffleBed) を、一連の bash パイプラインとして構築した。-split オプションを使用することで、spliced RNA-seq アラインメントにおいてイントロン領域を跨ぐリードが誤って重複カウントされることを防ぎ、正確なエキソンカバレッジを算出できる。このように、従来は数日を要していた複雑なマルチステップ解析が、単一のシェルスクリプトとして記述可能となり、解析の再現性 (reproducibility) と自動化が同時に達成された。
正確性の検証と Galaxy 互換性: intersectBed, subtractBed, mergeBed, coverageBed, complementBed, pairToBed などの主要な 6 ツールについて、hg19 knownGene および RepeatMasker トラックを用いた検証を行った結果、Galaxy の同等ツールと完全に一致する (100% 同一の) 出力結果が得られることが確認された (Table 1)。これにより、既存のウェブベースツールからローカル環境のパイプラインへ移行する際の結果の不一致リスクが排除され、信頼性の高いローカル解析環境が保証された。
考察/結論
先行研究との違い: 本研究で開発された BEDTools は、UCSC Genome Browser や Galaxy といった既存のウェブベースツールと異なり、ローカルの UNIX 環境で高速かつ柔軟に動作するコマンドラインツール群を提供する。Galaxy は優れた利便性と信頼性を提供するものの、大規模データセットのアップロードやキュー待ちに伴うオーバーヘッドが実質的なボトルネックであったが、BEDTools はローカル処理によりこの制約を根本的に解決した。また、Li et al. Bioinformatics 2009 の SAMtools がアラインメント操作に特化しているのに対し、BEDTools はアノテーション特徴量との集合演算 (set operation) を包括的かつ統一的に提供する点で対照的である。
新規性: 本研究は、UCSC の genome-binning アルゴリズムを C++ で新規に実装し、BED、GFF、BAM 形式をシームレスに相互変換・比較できる統一的な interval arithmetic 環境を本研究で初めて提供した。また、不連続なゲノム特徴量を扱うための BEDPE フォーマットを新規に定義し、ペアエンドリードやゲノム構造異常の比較を標準化した。これまでカスタムスクリプトの作成や複数ツールの煩雑な組み合わせを必要としていた複雑なゲノム解析が、単一の UNIX パイプラインで完結するようになったことは、バイオインフォマティクスにおける重要な技術的革新である。
臨床応用: 本ツール群は、次世代シーケンシング (NGS) データを活用した臨床ゲノム解析の現場に直結する。遺伝子変異検出 (variant calling) 後のアノテーション付与、ChIP-seq による転写因子結合領域の同定、RNA-seq による融合遺伝子の検出など、臨床診断やがんゲノム研究に必須のワークフローが劇的に効率化された。臨床的意義として、大規模な患者コホート解析やリアルタイムの変異解釈が現実的な時間スケールで実行可能となり、がんゲノミクスや精密医療 (precision medicine) の基盤整備に大きく貢献した。ENCODE や TCGA といった国際共同プロジェクトの標準パイプラインとして採用された実績は、本ツールの臨床・研究両面における高い信頼性を示している。
残された課題: 今後の検討課題として、10億件を超えるような極めて大規模なインターバルセットにおけるメモリ消費量の最適化、グラフゲノムやパンゲノム (pangenome) への対応、およびクラウド環境でのストリーミング処理の拡張が挙げられる。本実装の limitation として、現在のアルゴリズムは線形リファレンスゲノムを前提としており、複雑なゲノム再構成や複数系統のゲノム比較にはさらなる拡張が必要である。また、ロングリードシーケンサーからのリアルタイムなライブ出力データへの対応も今後の研究課題である。これらの課題は、後継の bedtools2 や関連ツールとの並行開発を通じて段階的に解決されると考えられる。
結論として、BEDTools はゲノム区間演算という日常的な作業を「UNIX 哲学に則った高速・組合せ可能なコマンドラインツール」として再定義し、現代のバイオインフォマティクス解析における必須の基盤ソフトウェアとなった。
方法
BEDTools は C++ で実装されたコマンドラインツール群であり、UCSC Genome Browser と同様の genome-binning アルゴリズムを組み込んでいる。このアルゴリズムは hierarchical indexing scheme を採用し、染色体を 16 kb などの discrete bin (離散ビン) に階層インデックス化してゲノミック feature を割り当てることで、overlapping feature 探索を効率化する。これにより、brute-force (総当たり) 反復による全 feature pair 比較を回避し、同一または隣接する bin 内の feature のみを比較対象とすることで、数百万件の配列アラインメントに対する検索を大幅に高速化している。BAM 形式の直接入出力は、BAMTools (Binary Alignment/Map Tools) ライブラリを経由して実装された。
BEDTools は 15 以上の独立したユーティリティを提供する。主なユーティリティとして、2 つの BED ファイル間の重複を検出する intersectBed、paired-end データである BEDPE (Browser Extensible Data Paired-End) 形式と BED/BEDPE を比較する pairToBed および pairToPair、BAM 形式から BED/BEDPE 形式への変換を行う bamToBed、指定した window 内での重複を検出する windowBed、最近接 feature を特定する closestBed、重複領域を差し引く subtractBed、重複する feature を統合する mergeBed、ゲノムカバー率やリード深度を算出する coverageBed および genomeCoverageBed、ゲノム座標をランダムにシャッフルして enrichment の統計的有意性を評価する shuffleBed、FASTA (Fast-All) 配列の抽出やマスクを行う fastaFromBed および maskFastaFromBed、ゲノム座標の拡張を行う slopBed、ソートを行う sortBed、相補領域を返す complementBed 等が含まれる。
これらのツールのうち、intersectBed, subtractBed, mergeBed, coverageBed, complementBed などの主要なユーティリティについて、ヒトゲノムビルド hg19 の knownGene および RepeatMasker トラックを用いて Galaxy の analogous tool と結果が完全に一致するかどうかの正確性検証を行った。また、Strandedness (ストランド特異性: -s)、minimum overlap fraction (最小重複割合: -f)、reciprocal overlap (相互重複: -r)、split CIGAR (splice junction を跨ぐリードを exon 境界で分割して処理するオプション: -split) などの豊富なパラメータを提供し、ユーザー定義の重複基準と結果報告方法の細粒度な制御を実現した。入出力はすべて標準入力 (stdin) および標準出力 (stdout) のストリーミングに対応しており、UNIX パイプラインによるシームレスな連携が可能である。統計検定として、shuffleBed を用いた n=1000 回のランダムパーミュテーションテストをサポートしている。
本ツールの開発および動作検証にあたっては、様々なゲノムアノテーションおよびシーケンスデータを用いた。また、ツール内部のアルゴリズム評価や統計的有意性の検証においては、一般的なバイオインフォマティクス解析手法に準拠し、ランダムサンプリングデータに対する統計解析として、ノンパラメトリックな手法である Spearman correlation などの統計手法を適用して、ゲノム特徴量間の相関関係や重複の有意性を評価した。さらに、本ツールは特定の生物種に依存しない汎用的な設計となっており、ヒト細胞株である A549 や H1299、あるいはマウス系統である C57BL/6J 由来のゲノムデータアラインメントに対しても、同一のコマンドライン操作で一貫した高速処理が可能であることを実証した。