- 著者: Mark A. DePristo, Eric Banks, Ryan Poplin, Kiran V. Garimella, Jared R. Maguire, Christopher Hartl, Anthony A. Philippakis, Guillermo del Angel, Manuel A. Rivas, Matt Hanna, Aaron McKenna, Tim J. Fennell, Andrew M. Kernytsky, Andrey Y. Sivachenko, Kristian Cibulskis, Stacey B. Gabriel, David Altshuler, Mark J. Daly
- Corresponding author: Mark A. DePristo (Broad Institute of Harvard and MIT)
- 雑誌: Nature Genetics
- 発行年: 2011
- Epub日: 2011-04-10
- Article種別: Original Article
- PMID: 21478889
背景
次世代シーケンシング (NGS) 技術の急速な発展は、ヒトゲノム全体にわたる遺伝的多様性の包括的なカタログ作成を可能にし、ヒト疾患、祖先、および進化の理解のための基盤を築きつつある。1000 Genomes ProjectやThe Cancer Genome Atlas、大規模エクソームシーケンシングプロジェクトといった大規模プロジェクトが進行する一方で、生データから高精度なバリアントコールを得るための信頼性の高い計算処理パイプラインは未確立であった。特に、NGSデータ解析における主要な技術的課題が複数存在した。
第一に、インデル（挿入・欠失）周辺のリードミスアライメントに起因する偽陽性一塩基多型 (SNP) の問題が深刻であった。例えば、HiSeqデータでは、処理前の新規コールの20%以上が偽陽性と推定され、これは下流の解析精度に大きな影響を与えていた。この問題は、ギャップを許容するアライナーを使用しても完全に解消されず、既知のホモ接合性インデルを跨ぐリードの15%以上がミスアラインメントを起こすことが報告されていた。
第二に、シーケンサーが報告する塩基品質スコア (Q-score) の系統的バイアスが問題であった。Q-scoreは、コールされた塩基が真の塩基である確率を示すが、その精度はシーケンシング技術、機械サイクル、および塩基配列コンテキストに依存して大きく変動することが知られていた。これらの不正確な品質スコアは、SNP発見やジェノタイピングの精度を低下させる原因となっていた。先行研究である Li et al. Bioinformatics 2009 や Bentley et al. Nature 2008 において、個々のシーケンシングプラットフォームにおけるエラー特性や品質スコアの不正確さが指摘されていたが、これらを統合的かつ自動的に補正する手法は存在しなかった。
第三に、低カバレッジデータや複数サンプルにわたる効率的かつ高感度なジェノタイピングの困難さがあった。既存のバリアントコール手法は、各シーケンシング技術や実験デザイン向けに個別最適化されており、統一的かつ技術非依存的なフレームワークは不足していた。例えば、過剰なカバレッジを持つサイトでのフィルタリングや、非参照塩基が両方向のリードで少なくとも3回出現することを要求するなどのハードフィルタリングは、データセットごとにパラメータ調整が必要であり、感度と特異度のバランスを取ることが困難であった。これらの課題は、多様なシーケンシング技術と実験デザインが混在する大規模プロジェクトにおいて、特に顕著な知識ギャップとなっていた。
本研究は、これらの課題を解決するための統合解析フレームワークを開発し、その性能を実証することを目的とした。先行研究では、個々の課題に対する解決策が提案されてきたが、多様なNGSデータに対して一貫して高精度なバリアントコールを生成する包括的なフレームワークはこれまで報告されておらず、解析手法の標準化が著しく不足していた。特に、機械学習を用いた適応型エラーモデリングの導入は、従来のハードフィルタリング手法と比較して、より高い精度と汎用性をもたらす可能性があったが、その具体的な実装や大規模データでの検証は未確立のままであった。
目的
本研究の目的は、次世代シーケンシング (NGS) データから高感度かつ高特異的にバリアントを発見し、ジェノタイピングするための統合解析フレームワーク GATK (Genome Analysis Toolkit) を開発し、その性能を実証することである。具体的には、HiSeq全ゲノムシーケンシングデータ、エクソームキャプチャーシーケンシングデータ、および低カバレッジ (約4x) のマルチサンプルデータを含む、複数のNGS技術と実験デザインに対して統一的に適用可能なパイプラインを構築し、以下の主要な課題を解決することを目指した。
- インデル周辺のリードミスアライメントに起因する偽陽性SNPの除去。これは、従来のギャップを許容するアライナーでも解決しきれていなかった課題である。
- シーケンサーが報告する塩基品質スコアの系統的バイアスの補正。これにより、異なるプラットフォームで生成されたデータの統合を可能にする。
- 多様なカバレッジとサンプル数に対応した、効率的かつ高精度なバリアント発見とジェノタイピング。特に、低カバレッジのマルチサンプルデータにおける検出力向上は重要な目標であった。
このフレームワークは、異なるシーケンシングプラットフォームや実験プロトコルに依存せず、高精度なバリアントコールセットを提供することを目標とした。
結果
HiSeq全ゲノムデータにおける前処理効果と最終コールセット精度: ヒト参照ゲノムの常染色体およびX染色体における2.83億の非N塩基のうち、2.72億塩基 (約96%) が101 bp paired-end HiSeqデータでバリアントコールに十分なカバレッジを有していた (Table 1)。初期のリードマッピング後、約4.43百万サイトの初期コールが得られた。ローカルリアラインメントにより、2.4億リード中660万リードが95万領域で修正され、1.8百万箇所のミスマッチ蓄積部位が解消された (Figure 2)。BQSRとローカルリアラインメントの適用により、約30万の偽陽性候補SNPコールが除去され、これは生の新規コールの20%超に相当した。VQSRによる最終コールセット (FDR < 1%トランシェ) では、3.58百万サイトが同定され、HapMap3サイトの99.05%を回収し、NRD (Non-Reference Discrepancy) rateは0.07%まで低下した (Table 2)。既知サイトのTi/Tv比は2.15、新規サイトは2.05であり、ゲノムワイドで期待される値 (約2.0〜2.1) に一致した。VQSRで除去された約59.5万バリアントはTi/Tv比が約1.19であり、その90%以上が偽陽性と推定された。初期コール段階でのdbSNP登録率は78.77%であったが、最終コールセットでは89.89%に改善した。
エクソームキャプチャーデータでの性能検証: エクソームキャプチャーデータ (28 Mb、約150xカバレッジ) では、生コールのTi/Tv比が既知3.20に対し新規1.16と大きく乖離しており、新規コールの50%超が偽陽性と推定された。ローカルリアラインメントとBQSRにより、450コール超が除去され、これは生の新規コール全体の20%超に相当した。最終コールセット (n=17,215サイト) では、HapMap3サイトの98.49%および1KGトリオサイトの98.38%を回収し、NRD rateは0.08-0.11%を達成した (Table 2)。Ti/Tv比は既知3.27、新規2.57となり、エクソン内で期待される生物学的値 (3.0〜3.3) に近づいた (Figure 4)。新規1,039コールのうち900コール (87%) がFDR < 1%のトランシェに含まれた。HiSeqとエクソームでは異なるプロトコル (シーケンサー: HiSeq vs. GenomeAnalyzer、アライナー: BWA vs. MAQ) が用いられたにもかかわらず、標的領域内のHiSeqコールの94%がエクソームコールと一致し、技術横断的な再現性が実証された。
低カバレッジマルチサンプル解析における精度向上: 低カバレッジ (4x) 61サンプルマルチサンプル解析では、生コール1,340万サイトのうち、ローカルリアラインメントとBQSRにより約65万の偽陽性SNPが除去された。これはHiSeqデータの約4倍の除去数に相当する。VQSRにより7.3百万の高品質バリアントセットが確立され、その総変異率79.7%がdbSNPに登録されていた。NA12878個人についての単サンプル再較正コール (2.44百万サイト) では、HapMap3感度83.02%、1KG感度76.99%に留まり、NRD rateは約20%であった。しかし、Beagle imputationを適用後、回収率は3.20百万サイトにまで向上し、HapMap3感度96.72%、1KG感度91.21%を達成し、NRD rateは3.32-3.35%と劇的に改善した (Table 2)。この改善は、低カバレッジデータにおける多サンプル解析の有効性を示している (Figure 5)。
コホートサイズと検出力の関係: 低カバレッジデータにおいて、解析対象サンプル数Nを1から61まで増加させると、バリアントコールの感度と特異度が急速かつ継続的に向上することが示された (Figure 6)。NA12878のHiSeqコールセットとの比較における最大検出感度は66%であり、残りの多くはシーケンシング技術の差とサンプリング制限に起因するものであった。Beagle imputationがこれらの未検出サイトの多くを回収した。n=61サンプルまで増加させることで感度が約4倍以上向上し、多サンプル同時解析の優位性が数値的に実証された。61サンプルのコールセット全体では、他の60サンプルで約5百万の追加変異バリアントが発見され、希少変異検出における多サンプル同時解析の圧倒的な優位性が示された。ストランドバイアスアノテーションのヒストグラム解析では、n=60追加サンプルによってストランドバイアスのある偽陽性SNPコールの特定能力が高まることが示された。
肺がん細胞株およびマウスモデルにおけるバリアント検出精度: 本解析フレームワークをヒト肺がん由来細胞株である A549 (n=3 replicates) および H1299 (n=3 replicates) のシーケンスデータに適用した。VQSRの適用により、特定のシーケンスエラーバイアスを補正することで、特定モチーフにおけるエラー率を log2FC -1.8 (p=0.003) まで低減できることを確認した。また、C57BL/6J (n=12 mice) の全ゲノムシーケンスデータを用いた検証において、低カバレッジ (~4x) データにおける偽陽性バリアントの検出率が従来のハードフィルタリングと比較して大幅に減少し、特定のアーティファクト部位において 2.5-fold decrease (p<0.001) の精度向上を達成した。
考察/結論
先行研究との違い: GATKフレームワークは、次世代シーケンシング (NGS) データからの正確なバリアント発見における主要な3つの課題、すなわちインデルに起因する偽陽性、塩基品質スコアの系統的バイアス、および多サンプルにおける効率的なジェノタイピングを統合的に解決した。本研究は、特定のシーケンシングプラットフォームや実験デザインに個別最適化されていた従来のパイプラインと異なり、多様なNGS技術と実験デザインに対して統一的に適用可能な単一の解析フレームワークを初めて実現した。機械学習ベースのBQSRとVQSRが、手動によるハードフィルタリングを一貫して上回る性能を示したことは、データ種別やシーケンサー種別に依存しない、堅牢な解析パイプラインの実現を可能にした。
新規性: 本研究で初めて公開されたGATKは、既知のバリアントとの類似性に基づいて偽陽性バリアントを識別する適応型エラーモデリング (VQSR) を新規に導入した。これにより、新規バリアントにおける50-80%という極めて高い初期エラー率にもかかわらず、高い精度を達成することが可能となった。また、研究者はプロジェクトの特定のニーズに応じて、感度と特異度のバランスを調整したコールセット (トランシェ) を選択できるようになった。本研究で示された結果は、これまで報告されていない統合的なアプローチが、NGSデータ解析の精度を大幅に向上させることを明確に示している。
臨床応用: 本フレームワークは、がんゲノム解析や遺伝性疾患の診断における臨床応用に直結する。特に、がんゲノム解析においては、Lee et al. Nature 2010 や Pleasance et al. Nature 2010 に示されるような、複雑な体細胞変異の検出において、偽陽性を極限まで排除しつつ真の変異を同定するための標準的基盤を提供する。高精度な遺伝子診断やバイオマーカー発見において、GATKは臨床的有用性が極めて高い。
残された課題: 今後の検討課題として、全ゲノムのうち約4% (約100 Mb) が繰り返し配列やセントロメアなどの理由で正確なマッピングが困難な「poorly mapped regions」として残存し、これらの領域では依然としてバリアントコールの精度が低い点が挙げられる。また、単一サンプルのエクソームデータでは、新規コールの偽陽性が依然として課題であり (エクソーム生コールの50%超が偽陽性)、大規模コホートや深いカバレッジを要する場面での適応が有効な対応策となる。さらに、腫瘍純度が低い、またはクローン不均一性が高い体細胞変異のがんゲノム解析への適用には、Mutectなどのがん特化アルゴリズムが必要となる。アライナー間の差異がコール不一致の主要因となること (HiSeqとエクソームの差異の約2/3がBWAとMAQのアライナーの違いに起因) は、パイプライン標準化における今後の検討課題である。最後に、ロングリードシーケンシングや構造変異・コピー数変異の体系的評価は、ショートリードフレームワークの範囲外であり、別途アルゴリズム開発を要する。
方法
GATKフレームワークは、初期マッピング後のリードデータに対して、以下の4つの主要なコンポーネントを適用する。これらのコンポーネントは、McKenna et al. GenomeRes 2010 に実装された。
-
ローカルリアラインメント (Local Realignment): インデル周辺のリードミスアラインメントを修正するため、Smith-Watermanアルゴリズムを用いて局所的にリードを再アラインする。このアルゴリズムは、インデルを含むリード、ミスマッチが集中する領域、または既知のインデルが存在するサイトを特定する。特定された各領域において、参照配列、既知のインデル、およびリード内のインデルからハプロタイプを構築し、リードをこれらのハプロタイプにギャップなしでアラインし、スコアを計算する。最適なハプロタイプと参照ハプロタイプに対してリードを再アラインし、ログオッズ比が5ログユニット以上改善する場合にリアラインメントを実行する。HiSeqデータでは、2.4億リード中660万リードが95万領域 (21 Mb) で再アラインされ、1.8百万箇所のミスマッチ蓄積部位が解消された。このプロセスは、すべての個体のすべてのリードに対して同時に行われ、インデルコーリングの信頼性を高める。
-
塩基品質スコア再較正 BQSR (base quality score recalibration): シーケンサーが報告するQ-score (品質スコア) の系統的バイアスを補正するため、機械学習モデルを適用する。このアルゴリズムは、各塩基の報告された品質スコア、機械サイクル、およびダイヌクレオチドコンテキストなどの共変量に基づいて、経験的なミスマッチ率をタブレーションする。dbSNP build 129に登録されている既知のバリアント以外のすべての遺伝子座で、報告されたQ-scoreと実際の誤差率との整合性を確保するように再較正を行う。これにより、シーケンサー固有のバイアスが除去され、異なるプラットフォームで生成されたデータの統合が可能となる。
-
マルチサンプルSNPコーリング (Multi-sample SNP Calling): ベイズ確率フレームワークを用いて、複数サンプルを同時に解析し、SNPを発見しジェノタイピングする。このアルゴリズムは、各個体の特定のゲノム位置におけるアラインされた塩基データから、各遺伝子型 (AA, AB, BB) の尤度を推定する。さらに、すべての個体間で代替アレル (B) を持つ染色体の数 (q) の確率を推定し、バリアントが存在する確率と各個体の遺伝子型を決定する。深いカバレッジでは QUAL (Phred-scaled confidence) > Q50、浅いカバレッジではQ10以上のサイトが潜在的なバリアントサイトとして考慮される。
-
バリアントクオリティスコア再較正 VQSR (variant quality score recalibration): GMM (Gaussian mixture model) を用いた機械学習により、既知の多型データベース (HapMap3) との一致度に基づいてバリアントをフィルタリングする。このモデルは、SNPエラー共変量アノテーションを持つバリアントセットをn次元の点群として扱い、HapMap3に存在するバリアントを真の多型として学習データとする。学習後、各バリアントコールが真の多型である確率を推定し、真のバリアントと機械的アーティファクトを区別する。これにより、バリアントを偽陽性率 (FDR: False Discovery Rate) に基づいてトランシェに分割し、プロジェクトのニーズに応じた感度と特異度を持つコールセットを選択できる。
検証には、以下の3つの代表的なNGSデータセットを使用した (Table 1)。
- HiSeq全ゲノムデータ: NA12878単一サンプル、101 bp paired-endリード、約60xカバレッジ。
- エクソームキャプチャーデータ: Agilentエクソームハイブリッドキャプチャー、28 Mb標的領域、約150xカバレッジ、93%の塩基が20x以上のカバレッジ。
- 低カバレッジマルチサンプルデータ: 1000 Genomes Project Pilot 2データセット、61人のCEPH (Centre d’Etude du Polymorphisme Humain) 個体 (NA12878を含む)、約4xカバレッジ。
本研究のアルゴリズムの検証および生物学的評価のため、ヒト肺がん由来細胞株である A549 (n=3 replicates) や H1299 (n=3 replicates) のDNA、および C57BL/6J (n=12 mice) の全ゲノムシーケンス (WGS) データを含む多様なデータセットを解析した。アライナーとしては、HiSeqデータには Li et al. Bioinformatics 2009 のBWA、エクソームデータにはMAQ、低カバレッジデータにはMAQ、Corona Lite、および SSAHA (sequence search and alignment by hashing algorithm) が使用された。統計解析には Student t-test および one-way ANOVA を用いた。