GATK HaplotypeCallerの使い方 後編 -GATK解説シリーズ-part 5
GATK HaplotypeCaller
使用難易度★★★☆☆
本記事は、GATK解説シリーズのPart 5の後編です。
eupatho-bioinfomatics.hatenablog.com
前編の記事は↓こちら eupatho-bioinfomatics.hatenablog.com
今回は何をする？
- GATK HaplotypeCallerを使って、Part 4の記事で得たBAM形式ファイルから、変異情報の記載されたVCFファイルを出力します。
- ものすごく時間のかかる（もし16GBのMacで実行すれば数日かかると思います）工程だったのですが、今回公式ページを確認したら、HaplotypeCallerにもマルチコアで解析するSparkモードがBeta版で搭載されていることに初めて気が付きました。開発中なのでまだ実際の解析に使用しないように注意喚起されていますが、今度試しに使ってみて感想をここに追記します。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
公式HPのHaplotypeCallerのリンク
https://gatk.broadinstitute.org/hc/en-us/articles/360056969012-HaplotypeCaller
2. 前提となるファイルについての説明。
Part-4で得たBAM形式ファイルを使う。
本ブログでいつも使用しているPlasmodium falciparum 3D7のリファレンスゲノムが必要となるので、GWASと名付けたフォルダの中にbamと3D7_genomesというフォルダを用意して、markdup.bamファイルとリファレンスゲノムをそれぞれ格納した状態で開始する。（リファレンスゲノムの入手方法は別の記事で紹介しています。）
また、解説シリーズ１の方法でGATK4を導入していることを前提とする。
2019年のScienceの論文で、アフリカのマラリア原虫の集団構造を2263株のWGSデータを使用して解析した大規模な研究のデータを一部拝借した。
science.sciencemag.org
3. Haplotype Callerを実行する。
workingdir=path/to/GWAS
cd ${workingdir}/bam
index=/path to GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta
for fpath in `ls *.markdup.bam`
do
fname=${fpath%.markdup.bam}
gatk HaplotypeCaller \
--reference ${workingdir}/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta \
--emit-ref-confidence GVCF \
--input ${fname}.markdup.bam \
--output ${fname}.g.vcf \
--sample-ploidy 1
done
オプションの説明
--input / -I : インプットファイル名
--output / -O: 出力ファイル名
--emit-ref-confidence / -ERC: 参照信頼度スコアを出力するモードを選択する。
--sample-ploidy / -ploidy :サンプルごとのPloidy（染色体の数）(デフォルト=2)。
#解析が終了するとg.vcfフォルダにファイルが生成される。
#作成したファイルをbqsrフォルダに移動する
mkdir bqsr
mv *.g.vcf ${workingdir}/bqsr/
mv *.g.vcf.idx ${workingdir}/bqsr/
お疲れ様でした。今回はこれで終わりです。
よければ他の記事のも見ていってください。
なお、本記事の執筆にあたりBWA および GATK4 を利用した SNPs/indels の検出方法のスクリプトを参考させていただきました。 とても参考になるサイトなので気になったらご自分で確認してみてください。
続きは↓こちら
eupatho-bioinfomatics.hatenablog.com
バイオインフォマティクス関連の書籍紹介記事 eupatho-bioinfomatics.hatenablog.com
ブログ主の自己紹介
