GATK HaplotypeCallerの使い方 前編 -GATK解説シリーズ-part 5
GATK HaplotypeCaller
使用難易度★★★☆☆
本記事は、GATK解説シリーズのPart 5の前編です。
eupatho-bioinfomatics.hatenablog.com
今回は何をする？
GATK HaplotypeCallerを使って、前回の記事で得たBAM形式ファイルから、変異情報の記載されたVCFファイルを出力します。
前編では、GATK HaplotypeCallerについて勉強しましょう。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
GATK HaplotypeCallerの概要
HaplotypeCallerは、WGSデータに含まれるSNPとインデルを、局所的なハプロタイプの再構築によって記述します。言いかえると、プログラムが変異を示す領域に遭遇すると既存のマッピング情報を破棄し、その領域を再構成します。さらに、HaplotypeCallerは、非二倍体生物のデータを扱うことができます。ただし、バリアントの可能性を計算するアルゴリズムは、（倍数性に比べて）極端な頻度の対立遺伝子には適していないので、体細胞性（癌）のバリアント検出には使用しないことが推奨されており、そのような場合には代わりにMutect2の使用が勧められています。
公式HPのHaplotypeCallerのリンク
https://gatk.broadinstitute.org/hc/en-us/articles/360056969012-HaplotypeCaller
HaplotypeCallerのワークフロー
ワークフローのイメージ図
1. アクティブな領域の定義
このステップでは、解析対象のサンプルがリファレンスと比較して変化しているゲノムの領域を選別する。これらの領域は「active regions」と定義される。バックグラウンドノイズのレベルを超す変動が見られない領域はスキップされる。この工程は、リファレンスと同一の領域で再構成を行わないことで、解析を高速化することを目的としている。
2. active regionsのアセンブリによるハプロタイプの決定
このステップでは、サンプルに存在する可能性のあるDNA配列を再構築する。active regionsにマッピングされたリードを使用して、ハプロタイプと呼ばれる配列を構築する。まず、可能性のあるハプロタイプのリストを作成するために、参照配列をテンプレートとしてactive regionsのアセンブリを構築する。
次に、各リードを順番に取得し、テンプレートとのマッチングを試みる。リードの一部がテンプレートと一致しない場合、ミスマッチを考慮してグラフに新しいノードを追加する。このプロセスを多くのリードで繰り返すと、多くの配列パターンを持つ複雑なアセンブリができあがる。しかし、GATKは各パターンを支持したリードの数を記録しているため、最も可能性の高い配列だけを選択することができる。
ハプロタイプが決定されると、それぞれのハプロタイプは元の参照配列に対して再調整され、潜在的なバリアントサイトが特定される。
3. リードデータからハプロタイプの尤度を計算する。
候補となるハプロタイプが揃ったところで、各ハプロタイプに対する各リードのペアワイズアライメントを行う。各リードを取り出し、PairHMMアルゴリズムを用いて、各ハプロタイプ（参照ハプロタイプを含む）に対して順番にアライメントを行う。PairHMMアルゴリズムでは、データの品質に関する情報（base quality scoreやindel quality scoreなど）を考慮し、各リードとハプロタイプのペアに対して、そのハプロタイプでそのリードが観測される可能性を表すスコアを出力する。
次に、これらのスコアを用いて、前のステップで特定された候補地において、個々の対立遺伝子の証拠がどれだけあるかを計算する。このプロセスはallelesに対するmarginalizationと呼ばれ、次のステップでサンプルに遺伝子型を割り当てるために最終的に使用される実際の数値を生成する。これにより、リードデータに対するハプロタイプの尤度が得られる。次に、これらの尤度データから各バリアントサイトの対立遺伝子の尤度を求める。
4. サンプルごとの遺伝子型の割り当て
前のステップでは、対象となる各候補のバリアントサイトについて、読み取りごとの対立遺伝子の尤度の表を作成した。最後に、これらの尤度を総合的に評価して、各サイトで最も可能性の高いサンプルの遺伝子型を決定する。これは、ベイズの定理を適用して、可能性のある各遺伝子型の尤度を計算し、最も可能性の高いものを選択することで行われる。これにより、ジェノタイプコールが生成されるとともに、バリアントコールが出力された場合には、出力VCFにアノテーションされる様々なメトリクスが計算される。
入力 :Variant call行うファイル (bam形式)
出力:フィルタリングされていない生のSNPおよびインデルを含むVCFまたはGVCFファイル。GVCFワークフローを利用する場合は、GVCFファイルが出力されるため、まずGenotypeGVCFsを実行し、その後フィルタリングを行ってから解析を行う必要がある。
