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を実行し、その後フィルタリングを行ってから解析を行う必要がある。
GATK MarkDuplicatesの使い方 -GATK解説シリーズ-part 4
GATK MarkDuplicate
使用難易度★★☆☆☆
本記事は、GATK解説シリーズのPart 4です。
eupatho-bioinfomatics.hatenablog.com
今回は何をする？
- 前回の記事で得たBAM形式ファイルを使って、GATK MarkDuplicates/MarkDuplicateSparkにより重複したリードにタグを付けます。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
1. GATK MarkDupicatesの概要
このツールでは、DNAライブラリ中の単一DNA断片に由来する重複リードを検出してタグを付けることができます。
BAMファイルまたはSAMファイル内の重複リードを検出してタグ付けします。重複リードとは、単一のDNA断片に由来するリードと定義され、PCRを使用したライブラリ構築などのサンプル調製時に発生する可能性があります。また、PCRの重複とは別に光学的重複も検出します。光学的重複とは、1つの増幅クラスターがシーケンシング装置の光学センサーによって複数のクラスターとして誤検出されることによって生じる重複です。
PCRによる重複の検出は、EstimateLibraryComplexityに記載の方法に従って行われ、BAM形式ファイル内のペアリードの配列から重複を判定します。各リードのN塩基から5塩基（デフォルト）を使ってソートされ、ペア同士が互いにギャップなく一致し、全体のミスマッチ率がMAX_DIFF_RATE（デフォルトでは0.03）以下であれば重複しているとみなされます。重複リードを収集した後、base-quality scores でリードをランク付けするアルゴリズムを使用して、プライマルリードと重複リードを区別します。このツールの出力は、各リードの重複が識別された新しいBAMファイルです（重複リードを破棄するわけではない）。
以前は、重複かどうかがだけ判別されましたが、重複の種類を特定するものではありませんでした。このため、最近BAMファイルのオプション出力として、DT（duplicate type）という新しいタグが追加されました。TAGGING_POLICYオプションを使用すると、すべての重複をマークする（All）、光学的重複のみをマークする（OpticalOnly）、または重複をマークしない（DontTag）ようにプログラムに指示することができます。BAMファイルの出力は、ライブラリPCRで生成された重複（LB）、またはシーケンシングプラットフォームのアーティファクトの重複（SQ）のいずれかのDT値を持ちます。必要に応じてREMOVE_DUPLICATEまたはREMOVE_SEQUENCING_DUPLICATESを使うことで、重複リードを削除することができます。
2. MarkDuplicatesSparkとMarkDuplicatesの違い
細かな違いはいくつかあるようですが、簡単にいうと複数コアで実行可能になったものがSparkだと認識しています（間違っているかも？）。16GB以上のメモリを持つPCでの使用が推奨されているようです。時間があったら速度を比較してみて、ここに追記したいと思います。
公式HPのリンク
MarkDuplicates
https://gatk.broadinstitute.org/hc/en-us/articles/360057439771-MarkDuplicates-Picard-
MarkDuplicatesSpark
https://gatk.broadinstitute.org/hc/en-us/articles/360057438771-MarkDuplicatesSpark
Rのおすすめ書籍
バイオインフォマティクスを独学で勉強している中で、「参考になる本が知りたい！誰か教えて！」と思うことが頻繁にあります。専門的な分野なので書籍自体の数も多くない上に、Amazonのレビュー数も少ないです。そこで、私が読んだ範囲内で参考図書の紹介をします。
これら以外に、読者の方で初心者・中級者向けのRの参考書で良書を知っている方がいれば、コメントで教えて下さると大変嬉しいです。新しく読んだ本の情報を随時更新していきます
注意事項: みなさまの研究に必要としているものと、私が欲している情報が完全に一致することはないと思います。従って、おすすめ度は現時点の私が感じた評価です。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
全くの初心者向け
Rをはじめよう 生命科学のためのRStudio入門 （羊土社 2019年）
おすすめ度 ★★★★☆
本当は星5つけたいのですが、読んだ本の数が少なすぎるため保留しています。
この本は正にタイトルの通り、生物・医学のWet研究者で、これまでデータの解析にはExcelやコマーシャルなGUIの解析ソフトを使っていたけど、 なんとか独学でRを使えるようにこれからなりたい！という方にぶっ刺さります。つまり、二年前の私にぶっ刺さりました。 実際にこの本で勉強した内容だけでも、論文に投稿できる図を作成することができました。
Rの基礎から、演習データを使ったグラフの作成方法、２群間・多群間の統計解析を一通り習得することができます。
これからRを取り入れたい方、研究室に配属されたばかりの学生さんに大変おすすめします。
データの解釈についても、字数を割いて丁寧に説明されているのもGoodでした。見た目も綺麗で読みやすいレイアウトです。
ただし、難しい内容が記載された書籍ではありませんので、中級者以上の方には物足りないと思います。
バイオインフォマティクスのおすすめ書籍紹介
eupatho-bioinfomatics.hatenablog.com
今後紹介したい本
Rとグラフで実感する生命科学のための統計入門
演習で学ぶ生命科学
カエル教える生物統計コンサルテーション
実験でつかうとこどだけ生物統計２ キホンのホン
Rによるやさしい統計学
中級者以上 Rによる多変量解析入門 データ分析の実践と理論
今後のGATK解析で使用するWGSデータのマッピング -GATK解説シリーズ-part 3
FastQC, trimmomatic, bwa, samtools
導入難易度★☆☆☆☆
使用難易度★★★☆☆
今回は何をする？
前回の記事で収集したマラリア原虫のWGSデータを使ってマッピングを行い、BAM形式のファイルを生成します。内容は以前に投稿した下記の記事と重複しますが、これから解析で使うデータを扱いたいので、改めて解析を行います。
本記事では、トリミングからマッピングまでの工程をコンパクトに記載しているので、実際にご自身のデータを使って再現される場合には、この記事のスクリプトを置き替える形で使ってもらえると便利だと思います。
次回から、今回用意したデータを使ってGATK4を使用していきます。
ちなみに、以前まで使用していたWGSデータは、系統間の類似性が高すぎて使いづらかったため今後の解析には使用しない予定です。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
以前の記事と同様に、
① FastQCでデータのクオリティコントロールを行い、トリミングの閾値を設定する。
eupatho-bioinfomatics.hatenablog.com
② Trimmomaticで低クオリティーのリードを取り除くトリミングを実行する。
eupatho-bioinfomatics.hatenablog.com
③ BWA-memでリファレンスゲノムに対して、マッピングを行う。
eupatho-bioinfomatics.hatenablog.com
の３つを順番に行います。
