VCFファイルとはなにかを説明します-後編

VCF (Variant call format)ファイルの見方

今回は何をする？

本記事は、VCFファイルを解説する記事の後編です。
前半をご覧になっていない方はこちらからどうぞ。

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

eupatho-bioinfomatics.hatenablog.com

後編では各バリアントサイトのレコードに付随する項目やスコアについて説明します。

今回はこちらで用意したVCFファイルをもとに説明していきます。すべてを貼り付けると長大になってしまうため、一部を抜粋します。 全体を見たい方はこちらのリンクからダウンロードしてください。(400MBほど）

GATK公式HPのVCFについてのリンク

1. バリアントサイトの各項目について

各バリアントサイトのレコードは、次の10個の要素で構成されています。

#CHROM  POS  ID  REF  ALT  QUAL  FILTER  INFO  FORMAT  sample-nameが続く・・・ .

１行だけ抜粋したものを下に記載します。

#CHROM   POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  ERR015377
Pf3D7_01_v3 110 .   A   C   124.72  .   AC=1;AF=0.071;AN=14;DP=99;FS=0.000;MLEAC=2;MLEAF=0.143;MQ=36.03;QD=31.18;SOR=3.258  GT:AD:DP:GQ:PL  0:9,0:9:99:0,135

それでは順番に見ていきましょう。

CHROMとPOS

#CHROM   POS
Pf3D7_01_v3 110
Pf3D7_01_v3 260
Pf3D7_01_v3 264

バリアントがあるコンティグとゲノムの座標情報です。
欠失の場合は、実際にはその前の塩基が指定されることに注意。

ID

#CHROM   ID
Pf3D7_01_v3 .
Pf3D7_01_v3 .
Pf3D7_01_v3 .

バリアントの識別子。参照したデータベースに、このサイトのレコードが存在すると記載されます。
情報がない場合は、「・」で表記されます。 典型的な識別子はdbSNPのIDで、ヒトのデータではrs28548431のようになります。

REFとALT

#CHROM   REF ALT
Pf3D7_01_v3 A   C
Pf3D7_01_v3 C   CTCTTACTTACTTACT

参照遺伝子と対立遺伝子が記載されています。そのバリアントがSNPなのかインデルなのかを教えてくれます。上記の例の場合、一つ目がSNPで、二つ目がインデルです。

QUAL

#CHROM   QUAL
Pf3D7_01_v3 124.72
Pf3D7_01_v3 251.07
Pf3D7_01_v3 248.82

シーケンスデータから、この部位にREF/ALT多型が存在する可能性をPhred-scaledで表したものです。Phred-scaledについてはこちらを参照してください。ちなみに10=10分の1、20は100分の1のエラー確率を表します。ただし、この値は大量のサンプルを用いてバリアントコールを行うと、非常に大きくなる可能性があるため、QUALはバリアントコールの品質を評価するためにはあまり有用な特性ではありません。

VCFファイルとはなにかを説明します-前編

VCF (Variant call format)ファイルの見方

今回は何をする？
本記事では、これまでに全くVCFファイルに触れたことのない方に向けて、ファイルを構成する要素を詳しく解説しています。VCFは変異解析に不可欠な要素であり、必ず理解しなければなりません。初めて見た時の感想は、「うわ、ちんぷんかんぶん・・・」という印象を持つと思いますが、慣れてくると意外と単純な要素で構成されていることに気が付きます。VCFファイルの形式を理解しないと、GATKやvcf/bcftoolsを使ったフィルタリングで何をしているか理解できなくなってしまうので、この機会に一緒に勉強していきましょう。二つの記事を読み終わる頃には、きっとVCFファイルについて苦手意識は無くなっていると思います！
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。 フォローする @harrykun_blog
前編では、VCFを構成する要素について順番に解説します。後編では各バリアントサイトのレコードに付随する項目やスコアについて説明します。

今回はこちらで用意したVCFファイルをもとに説明していきます。すべてを貼り付けると長大になってしまうため、一部を抜粋します。 全体を見たい方はこちらのリンクからダウンロードしてください。(400MBほど）

GATK公式HPのVCFについてのリンク

1. VCFフォーマットとは

VCF（Variant Call Format）を一言で説明すると、SNP、インデル、構造的変異を記述するために標準化されたテキストファイル形式です。 フォーマットの詳細な仕様は、こちらのPDFで見ることができます。

VCFはGATKがバリアントコールに使用しているフォーマットで、主にバリアントの種類とシークエンス、および個々のバリアントに関する複数サンプルのジェノタイプについての情報が記述されます。とはいえ、GATK HaplotypeCallerのようなツールが生成するVCFファイルは少し複雑です。ここでは、GATKが出力するVCFファイルを理解するために知っておくべき点について説明します。

なお、VCFファイルはプレーンテキストファイルなので、Excelなどのテキストエディタで開いて編集することができますが、VCFファイルは大きな容量になりがちであるためファイルを読み込むのに時間がかかる場合があります。そのため、GATKのSelectVariantsのようなツールを使ってデータを編集するアプローチをとる必要がでてくるわけです。 また、Microsoft Wordなどのワープロでvcfを編集するとフォーマットが崩れてしまうので絶対にしないでください。

2. VCFファイルの構成

VCFファイルは、ヘッダーとバリアントコールレコードの2つの要素で構成されています。

GATK公式HPからの引用

ヘッダーには、データセットおよび関連するリファレンスソース（生物、ゲノムビルドバージョンなど）に関する情報に加えて、VCFファイルに含まれるバリアントコールの特性を修飾および定量化するために使用されるすべてのアノテーションの定義が含まれています。GATKツールで生成されたVCFのヘッダーには、生成に使用されたコマンドラインも含まれています。

実際のデータを順番に見ていきましょう。

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

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