GATK JointGenotyping -GATK解説シリーズ-part 6

GATK GenomicsDBimport, GATK GenotypeGVCFs, Picard VcfToIntervalList

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

eupatho-bioinfomatics.hatenablog.com

今回は何をする？

• GATK GenomicsDBimport および GATK GenotypeGVCFs を使って、前回の記事で得たVCF形式ファイルから、変異情報を記述したローカルなデータベースを構築し、Joint Genotypingを実施して複数のvcfファイルをまとめたmerged.vcfファイルを出力します。

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

公式HPリンク
GATK GenomicsDBimport
GATK GenotypeGVCFs
Picard VcfToIntervalList

ジョイントジェノタイピングの概要

まず、対象とする染色体の名称および領域を指定する必要があるため、その情報を記載したインターバルリストを VcfToIntervalList (Picard) を使って用意します。

次に、ジョイントジェノタイピングを行う前に、GenomicsDBImportを使用してデータベースを構築します。 これにより、サンプルを中心としたバリアント情報をゲノム上の遺伝子座に移し替え、ツールがデータにアクセスしやすくします。

最後に、GenotypeGVCFs を使用して 、HaplotypeCallerで作成されたサンプルに対してジョイントジェノタイピングを行います。このツールは複数のサンプルを扱うことができます。入力サンプルは、HaplotypeCallerで -ERC GVCF または -ERC BP_RESOLUTION で生成した遺伝子型尤度を持っていることが前提となります。

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

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

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

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

3.　Joint Genotypingを実行する。

まずはインターバルリストを作成する。入力ファイルとして、適当なvcfファイルを選び入力する。 （ここで用いるVcfToIntervalListは、まだGATKには実装されていないようなので、必要に応じてPicardをインストールしてください。）

workingdir=path/to/GWAS
cd ${workingdir}/bqsr

java -jar ${EBROOTPICARD}/picard.jar VcfToIntervalList I=ERS032649.g.vcf O=sample.interval_list

注意! なぜか分からないが、interval_listはこの形にしないとエラーがでる。
可　interval_list
不可　interval.list

次に、 vcfファイルを連続で表記するための変数を定義してデータベースを構築し、ジョイントジェノタイピングを実施する。

gvcf_files=""
for gvcf_file in `ls *.g.vcf`
do
 gvcf_files=${gvcf_files}"-V ${gvcf_file} "
done

index=/path to GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta
gatk  GenomicsDBImport \
      --reference ${index} ${gvcf_files} \ 
      --intervals sample.interval_list \
      --genomicsdb-workspace-path gvcfs_db
               
gatk GenotypeGVCFs \
      --reference ${index} \ 
      --variant gendb://gvcfs_db \
      --output merged.vcf

オプションの説明
--variant / -V : インプットファイル名
--output / -O: 出力ファイル名
--reference / -R: リファレンスファイル
--intervals / -L :インターバルリスト名

解析が終了するとsample.interval_list、gvcfs_dbのフォルダ、merged.vcfが生成される。

お疲れ様でした。今回はこれで終わりです。
よければ他の記事のも見ていってください。

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

続きは↓こちら

eupatho-bioinfomatics.hatenablog.com

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