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

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

MENU
トップ > GATK > GATK BQSR後のVariant Calling -GATK解説シリーズ-part 10

GATK BQSR後のVariant Calling -GATK解説シリーズ-part 10

GATK その他

GATK BQSR後のVariant Calling

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

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

今回は何をする？

  • 前回のGATK VariantFilterationで出力された、*.bqsr.bamを基に２回目のVariant Callingを行います。

  • 今回行う内容は、GATK解説シリーズ Part 5 からPart 8で行った仕事の繰り返しになります。新しく覚えることは特にありません。

  • 今回得たmerged_snps/indels_filtered.vcfを使って、多様なゲノムワイド解析を行うことができます。


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

1. 作業の概要

BQSR後のBAMファイルを使って、

 GATK  HaplotypeCaller
 GATK  GenomicDBImport
 GATK  GenotypeGVCFs
 GATK  SelectVariatns
 GATK  VariantFiltration


を一気に行います。

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

Part-9で得たBQSR後のBAM形式ファイルを使う。
本ブログでいつも使用しているPlasmodium falciparum 3D7のリファレンスゲノムが必要となるので、GWASと名付けたフォルダの中にbamと3D7_genomesというフォルダを用意して、bqsr.bamファイルとリファレンスゲノムをそれぞれ格納した状態で開始する。（リファレンスゲノムの入手方法は別の記事で紹介しています。）

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

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

3.　コマンドの実行

ステップ1 HaplotypeCaller（2回目）の実行。

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

cd ${workingdir}/bam

for fpath in `ls *.bqsr.bam`
do
    fname=${fpath%.bqsr.bam}
    gatk HaplotypeCaller -R ${index} --emit-ref-confidence GVCF \
                        -I ${fname}.bqsr.bam -O ${fname}.g.vcf --sample-ploidy 1
done

mv *.g.vcf ../vcf/
mv *.g.vcf.idx ../vcf/

ステップ 2 ジョイント・ジェノタイピング（２回目）の実行。

cd ${workingdir}/vcf

# list up all gVCF files
gvcf_files=""
for gvcf_file in `ls *.g.vcf`
do
  gvcf_files=${gvcf_files}"-V ${gvcf_file} "
done

##!!! I=にはどれでもいいので、サンプル名をいれること。
java -jar ${EBROOTPICARD}/picard.jar VcfToIntervalList I=ERS032649.g.vcf O=sample.interval_list
    gatk GenomicsDBImport -R ${index} ${gvcf_files} -L sample.interval_list \
                      --genomicsdb-workspace-path gvcfs_db

    gatk GenotypeGVCFs -R ${index} -V gendb://gvcfs_db -O merged.vcf

ステップ3 フィルタリング（２回目）の実行。

cd ${workingdir}/vcf

#snpsを抽出してフィルタリングを行う
    gatk SelectVariants -R ${index} -V merged.vcf --select-type-to-include SNP -O merged_snps.vcf

    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"

#indelsを抽出してフィルタリングを行う

cd ${workingdir}/vcf

    gatk SelectVariants -R ${index} -V merged.vcf --select-type-to-include INDEL -O merged_indels.vcf

    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"



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

Part１から２ヶ月に渡り、GATKの基本的な使い方について解説してきましたが、今回で一区切りつきました。
まだまだ不十分な部分も多くあると思いますが、これからも勉強し続けて少しずつアップデートしていきたいと考えています。

Part1からずっと読んでくれた方がいらっしゃいましたら、本当にありがとうございました。貴方の悩みの解決に貢献できたでしょうか。
みなさんの研究に、少しでもこのブログが役立てていれば幸いです。
もし、気に入っていただけたら、他の方にシェアしていただけるとブログ主がとても喜びます😄。

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

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

次回からは、最終的に取得したSNPファイルを公開して、系統解析をやっていきます。

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

eupatho-bioinfomatics.hatenablog.com

eupatho-bioinfomatics.hatenablog.com

eupatho-bioinfomatics.hatenablog.com