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

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

MENU

GATK VariantFiltrationの使い方 -GATK解説シリーズ-part 8

GATK VariantFiltration

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

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

今回は何をする?

  • 前回のGATK SelectVariantsの使い方の続きになります。

  • GATK VariantFiltrationを使って、vcfファイルから低クオリティのバリアントを除外します。

  • FilteringのステップはSelectVarinatsとVariantFiltrationをセットで実行することになります
    GATK SelectVariatnsについては前回の記事を参照してください。

  • 今回のVariantFiltrationの条件は解析結果に大きな影響を及ぼしますが、実用レベルでもよく使用するオプションの種類が多く、最適な設定の検討が難しいことなどから奥が深くて難解なステップだと感じています。しかし、ここを乗り切れば難所はおしまいです!!

  • FIlteringの条件として、GATK Legacy Forumの値を参考にしました。リンクはこちら
    少し情報が古いので、もっといい参考資料を発見したらアップデートします。

  • VCFファイルについて理解していない方には難しい内容だと思います。
    VCFファイルについて知りたい方は、先に↓こちらの記事を読むことをお勧めします

eupatho-bioinfomatics.hatenablog.com

eupatho-bioinfomatics.hatenablog.com


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

公式HPの VariantFiltrationのリンク

https://gatk.broadinstitute.org/hc/en-us/articles/4404604763547-VariantFiltration#--set-filtered-genotype-to-no-call

1. VariantFiltrationの概要

VariantFiltrationでは、INFOおよび/またはFORMATアノテーションに設定した任意の基準に基づいて、バリアントのフィルタリングを行います。
フィルタリングされたレコードは、コマンドラインで指定しない限り出力ファイルにも保存されます。

入力として

フィルタリングを行うVCFを入力に持ちます。 今回は、Part7-1で作成したmerged_snps.vcfおよびmerged=indels.vcfを使用します。

出力として

基準を通過した variant は PASS と表記され、失敗した variant はそのフィルタ名がFILTERに記載された新しいVCFを得る。

2. よく使うオプションの説明

  • --filter-expression: INFO フィールドのフィルタリングに適用される条件
    注意!! ひとつの--filter-expressionに複数の条件を記載している例を見かけるかもしれませんが、公式HPでは1つの--filter-expressionに一つにフィルターを記載することを推奨しています。

  • --filter-name :フィルターのリストに使用する名前
    失敗した variant のFILTERに記載されるフィルタ名。

  • --genotype-filter-expression: FORMAT フィールドのフィルタリングに適用される条件
    INFOの代わりにFORMAT(遺伝子型)フィールドに適用されます。

  • ---genotype-filter-name : フィルターのリストに使用する名前
    失敗した variant のFILTERに記載されるフィルタ名。

  • --set-filtered-genotype-to-no-call : フィルタリングされたジェノタイプをノーコールに設定する
    この引数が与えられた場合、フィルタリングされたジェノタイプをノーコール./.に設定します。デフォルトはFALSE

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

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

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

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

3. GATK VariantFiltrationを実行する。


SNPと INDEL を異なる条件でフィルタリングしていることに注意!

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

gatk VariantFiltration  \
                        -R ${index}  \
                        -V merged_snps.vcf  \
                       -O merged_snps_filtered.vcf \
                       -filter "MQ < 20.0" --filter-name "MQ20"     \
                       -filter "QD < 2.0" --filter-name "QD2"       \
                       -filter "FS > 60.0" --filter-name "FS60"     \
                       -filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
                       -filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8"



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

gatk VariantFiltration  \
       -R ${index}  \
                        -V merged_indels.vcf  \
                       -O merged_indels_filtered.vcf \
                       -filter "QD < 2.0" --filter-name "QD2"       \
                       -filter "FS > 200.0" --filter-name "FS200"   \
                       -filter "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20"



出力されたファイルの詳しい解説は別で行います。


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


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

なお、本記事の執筆にあたりBWA および GATK4 を利用した SNPs/indels の検出方法のスクリプトを参考させていただきました。 とても参考になるサイトなので気になったらご自分で確認してみてください。


バイオインフォマティクス関連の書籍紹介は↓こちら eupatho-bioinfomatics.hatenablog.com