バイオインフォマティクスでゲノムワイド関連解析（GWAS)

バイオインフォマティクスを頑張っている方が、本ブログの内容を真似することで、自分のデータで解析ができる情報を提供することが目標です! 今はGATKの解説をメインテーマにしています。

MENU
トップ > GATK > GATK BaseRecalibratorとApplyBQSRの使い方 -GATK解説シリーズ-part 9

GATK BaseRecalibratorとApplyBQSRの使い方 -GATK解説シリーズ-part 9

GATK その他

GATK BaseRecalibrator, ApplyBQSR

使用難易度★☆☆☆☆
本記事は、GATK解説シリーズのPart 9です。

GATK解説シリーズのリンクまとめは↓こちら eupatho-bioinfomatics.hatenablog.com

今回は何をする？

  • 前回のGATK VariantFilterationで出力された、merged_snps/indels_filtered.vcfを基にBQSRを行います。

  • これまであまり理解せずに使用していたBQSRですが、この機会に公式HPのGATK BaseRecalibratorのページを全部読んでみたので、内容をかいつまんで説明します。

  • さらに、 BaseRecalibratorとApplyBQSRの具体的な使用例を紹介します。


Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。

公式HPのリンク
Base Quality Score Recalibration (BQSR)　
BaseRecalibrator
ApplyBQSR

1. BaseRecalibratorの概要　

Base Quality Score Recalibration (BQSR)を一言で説明すると、
シーケンスの際に生じたエラーを補正することで、より精度の高いBAMファイルを得るためのステップです。

Base Quality Scoreは、シーケンシングマシンが生成する塩基ごとのエラーの推定値で、装置が毎回正しい塩基をCallしたという自信の度合いを表しています。例えば、Aという塩基を読み取った場合、Q20という品質スコアが割り当てられたとします。これは、Phredスケールでは、99％の確率で正しく塩基を特定したことを意味します。一見すると高精度と思われるかもしれませんが、100のうち1件は間違っていることが予想されることを意味しています。仮に、数十億回の塩基呼び出しがあったとすると（30倍のゲノムでは約900億回）、なんと9億回も塩基を間違えることになり、これは非常に多くのエラー塩基を含むこととなります。各ベースコールの品質スコアは、シーケンシングマシンのメーカーが秘匿しているアルゴリズムによって決定されます。

このことがなぜ重要かというと、、、

私たちのVariant Caliingのアルゴリズムは、個々に割り当てられた品質スコアに大きく依存しているためです。品質スコアが低いベースコールがあった場合、実際には別のものである可能性があることを意味します。そのため、品質の高い他のベースコールほどは信用できません。つまり、特定のサイトに存在するAlleleの有無を示す証拠を評価するために、このスコアを使用します。

では、BaseRecalibratorは何をしているかというと、、、

残念ながら、シーケンス装置が作り出すスコアには様々な系統的な（ランダムではない）技術的エラーがあり、データの塩基品質スコアが過大評価されたり過小評価されたりすることがあります。これらのエラーの中には、シーケンス反応がどのように機能するかという物理学や化学に起因するものもあれば、機器の製造上の欠陥に起因するものもあります。

BQSRは、シーケンサーによって様々な物理的・化学的・機器の構造上生じてしまう品質スコアに含まれるエラーを機械学習を使ってモデル化し、それに応じて品質スコアを調整するプロセスです。

例えば、あるプロセスにおいて、2つのA塩基をが連続してコールされた場合、次にコールした塩基のエラー率が1％高くなります。つまり、リードの中でAAの次に来る塩基のコールは、品質スコアを1％下げる必要があります。実際には、いくつかの異なる共変量に対して、相加的に行われます。これにより、全体としてより正確な塩基の品質スコアを得ることができ、その結果、Variant Callingの精度が向上します。つまり、低品質なAが実際にはTであるべきだったかどうかを判断することはできませんが、少なくともバリアントコールの担当者に、そのAをどこまで信頼できるかをより正確に伝えることはできます。この手順は、基本的な品質スコアを出力するあらゆるシーケンサーのデータを含むBAMファイルに適用することができます。


入力として

Markduplicatesで出力したmarkdup.bam
merged_snps_filtered.vcf
merged_indels_filtered.vcf

出力として

品質スコアによる再校正テーブル(*_recal_data.table)

なお、BQSR後のオプションとして、キャリブレーションの効果を可視化するためにBQSR前後のプロットを作成すると、品質管理に役立ちます。

2. ApplyBQSRの概要　

ApplyBQSRは、BaseRecalibratorによって作成された再校正テーブルに基づいて、入力リードの塩基品質を再校正し、再校正されたBAMを出力します。

第2段階では、第1段階で特定されたパターン（再較正テーブルの記録）に基づいて、個々のベースコールに数値補正を適用し、再較正されたデータを新しいBAMに書き出します。

入力として

Markduplicatesで出力した*.markdup.bam
BQSRで出力した*_recal_data.table

出力として

再校正されたbamファイル*.bqsr.bam　

3.　前提となるファイルについての説明。

Part-7-2で得たmerged_snps_filtered.vcfおよびmerged_indels_filtered.vcfファイルを使う。
また、以前の記事で出力した*.markdup.bamを全て使用します。 本ブログでいつも使用しているPlasmodium falciparum 3D7のリファレンスゲノムが必要となるので、GWASと名付けたフォルダの中にbam, bqsrと3D7_genomesというフォルダを用意して、*.markdup.bamファイル、 merged_snps/indels_filtered.vcfファイルとリファレンスゲノムをそれぞれ格納した状態で開始する。（リファレンスゲノムの入手方法は別の記事で紹介しています。）

また、解説シリーズ１の方法でGATK4を導入していることを前提とする。

2019年のScienceの論文で、アフリカのマラリア原虫の集団構造を2263株のWGSデータを使用して解析した大規模な研究のデータを一部拝借した。
science.sciencemag.org

4. BQSRを実行する。

cd /path to GWAS/bam
index=/path to GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta

for fpath in `ls *.markdup.bam`
do
fname=${fpath%.markdup.bam}

gatk BaseRecalibrator -R ${index} -I ${fname}.markdup.bam \
                        --known-sites ../bqsr/merged_snps_filtered.vcf \
                        --known-sites ../bqsr/merged_indels_filtered.vcf \
                        -O ${fname}_recal_data.table
  
gatk ApplyBQSR -R ${index} -I ${fname}.markdup.bam -bqsr ${fname}_recal_data.table -O ${fname}.bqsr.bam

done




お疲れ様でした。今回はこれで終わりです。

Part9までで新しい内容は終わりです。後はVariant Callからの内容をもう一度繰り返すだけです。
つまり、、、後一息であなたもGATK使いになれます！！！やったー！おめでとうございます👏👏👏
最後は、繰り返しになるので一気にコードだけを紹介しようと思います。

それが終わったら、ロードマップをリニューアルしてここまでの全体のスクリプトを記載したファイルを公開する予定です。
GATKに関する記事はこれからも投稿し続けます。


よければ他の記事のも見ていってください。

なお、本記事の執筆にあたり[GATK4 を使用した variant calling 時の前処理)(https://bi.biopapyrus.jp/gwas/gatk/preprocessing.html)と、
BWA および GATK4 を利用した SNPs/indels の検出方法のスクリプトを参考させていただきました。
とても参考になるサイトなので気になったらご自分で確認してみてください。

この他のおすすめ記事は↓こちら

eupatho-bioinfomatics.hatenablog.com

eupatho-bioinfomatics.hatenablog.com