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

今回は何をする？

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

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

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

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つの要素で構成されています。

f:id:Harry-kun:20210728131822p:plain — GATK公式HPからの引用

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

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

VCFバージョンの情報

１行目は、ファイルが準拠しているVCFのバージョンを示します。バイオインフォマティクスは変化の激しい分野であり、ファイルのフォーマットも急速に進化しています。VCFファイルを解析しようとして予期せぬ問題に遭遇した場合、フォーマットの仕様変更がないか、バージョンを確認してください。

##fileformat=VCFv4.2

FILTERについて

FILTER行では、データにどのようなフィルターが適用されたかが記述されています。下の例では、FILTER=<ID=LowQual,Description="Low quality">というフィルターが適用されています。ここに記載されたフィルターに適合しなかったレコードのFILTERフィールドには、フィルタのID（ここではLowQual）が記載されます。羅列されたすべてのフィルター条件をすべて満たす場合にPASSと記載されます。

##FILTER=<ID=LowQual,Description="Low quality">
##FILTER=<ID=PASS,Description="All filters passed">

FORMAT および INFO について

これらの行は、VCFファイルのFORMATおよびINFO列に含まれるアノテーションが記述されています。
簡単に言うと、FORMATが個別サンプルに付随する情報で、INFOが全体に付随する情報になります。
各項目の詳細は、後半に記載します。

FORMAT

##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=RGQ,Number=1,Type=Integer,Description="Unconditional reference genotype confidence, encoded as a phred quality -10*log10 p(genotype call is wrong)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">

INFO

##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed">
##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed">
##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
##INFO=<ID=RAW_MQandDP,Number=2,Type=Integer,Description="Raw data (sum of squared MQ and total depth) for improved RMS Mapping Quality calculation. Incompatible with deprecated RAW_MQ formulation.">
##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">

GATKCommandLineについて

この行には、ファイルを生成したツールで使用されたすべてのパラメータが含まれています。ここでは、GATK HaplotypeCallerで使用したコマンドラインを指します。これらのパラメータには、ツールが受け付けたすべての引数と、適用された数値が含まれています。

##GATKCommandLine=<ID=GenomicsDBImport,CommandLine="GenomicsDBImport --genomicsdb-workspace-path gvcfs_db --variant ERR015377.g.vcf --variant ERR015425.g.vcf --variant ERR1035536.g.vcf --variant ERR1045266.g.vcf --variant ERR1106528.g.vcf --variant ERR1106529.g.vcf --variant ERR2000569.g.vcf --variant ERR211448.g.vcf --variant ERR343116.g.vcf --variant ERR405240.g.vcf --variant ERR405245.g.vcf --variant ERR450058.g.vcf --variant ERR450079.g.vcf --variant ERR484676.g.vcf --variant ERR562868.g.vcf --variant ERR580552.g.vcf --variant ERR636426.g.vcf --variant ERR636430.g.vcf --variant ERR701750.g.vcf --variant ERR701756.g.vcf --variant ERS032647.g.vcf --variant ERS032649.g.vcf --variant ERS347567.g.vcf --variant ERS347575.g.vcf --variant SRS399547.g.vcf --intervals sample.interval_list --reference /Hatena/GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta ・・・


##GATKCommandLine=<ID=GenotypeGVCFs,CommandLine="GenotypeGVCFs --output merged.vcf --variant gendb://gvcfs_db --reference /Hatena/GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta --include-non-variant-sites false --merge-input-intervals false --input-is-somatic false --tumor-lod-to-emit 3.5 --allele-fraction-error 0.001 --keep-combined-raw-annotations false --use-posteriors-to-calculate-qual false --dont-use-dragstr-priors false --use-new-qual-calculator true --annotate-with-num-discovered-alleles false ・・・


##GATKCommandLine=<ID=HaplotypeCaller,CommandLine="HaplotypeCaller --sample-ploidy 1 --emit-ref-confidence GVCF --output ERR015377.g.vcf --input ERR015377.markdup.bam --reference s /Hatena/GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta --use-posteriors-to-calculate-qual false --dont-use-dragstr-priors false --use-new-qual-calculator true --annotate-with-num-discovered-alleles false --heterozygosity 0.001 --indel-heterozygosity 1.25E-4 --heterozygosity-stdev 0.01 --standard-min-confidence-threshold-for-calling 30.0 --max-alternate-alleles 6 --max-genotype-count 1024 --num-reference-samples-if-no-call 0 --genotype-assignment-method USE_PLS_TO_ASSIGN ・・・

Contigについて

使用されたリファレンスアセンブリのコンティグ名、長さが記載されています。

##contig=<ID=Pf3D7_01_v3,length=640851>
##contig=<ID=Pf3D7_02_v3,length=947102>
##contig=<ID=Pf3D7_03_v3,length=1067971>
##contig=<ID=Pf3D7_04_v3,length=1200490>
##contig=<ID=Pf3D7_05_v3,length=1343557>
##contig=<ID=Pf3D7_06_v3,length=1418242>
##contig=<ID=Pf3D7_07_v3,length=1445207>
##contig=<ID=Pf3D7_08_v3,length=1472805>
##contig=<ID=Pf3D7_09_v3,length=1541735>
##contig=<ID=Pf3D7_10_v3,length=1687656>
##contig=<ID=Pf3D7_11_v3,length=2038340>
##contig=<ID=Pf3D7_12_v3,length=2271494>
##contig=<ID=Pf3D7_13_v3,length=2925236>
##contig=<ID=Pf3D7_14_v3,length=3291936>
##contig=<ID=Pf3D7_API_v3,length=34250>
##contig=<ID=Pf3D7_MIT_v3,length=5967>

3. バリアントレコードについて

各バリアントサイトのレコードは、次のように構成されています。

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

VCFレコードの最初の8列（INFOまで）は、バリアントサイトのレベルで測定されたスコアを表しています。
VCFファイルに複数のサンプルが含まれている場合、サイトレベルのアノテーションの中には、異なるサンプルからそのサイトについて得られた値の合計や平均値が含まれています。
遺伝子型や個々のサンプルレベルのアノテーション値などのサンプル固有の情報は、FORMAT列（9列目）およびsample-name列（10列目以降）に含まれています。ほとんどのプログラムでは、サンプル名のアルファベット順に並べていますが、必ずしもそうではありません。

#CHROM   POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  ERR015377   ERR015425   ERR1035536
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    0:2,0:2:49:0,49 .:0,0:0:.:0,0
Pf3D7_01_v3 260 .   A   G   251.07  .   AC=2;AF=0.087;AN=23;BaseQRankSum=0.699;DP=188;FS=4.150;MLEAC=2;MLEAF=0.087;MQ=31.11;MQRankSum=0.301;QD=11.41;ReadPosRankSum=0.886;SOR=1.981 GT:AD:DP:GQ:PL  0:7,0:7:99:0,118    0:3,0:3:99:0,131    0:3,0:3:90:0,90
Pf3D7_01_v3 264 .   C   CTCTTACTTACTTACT    248.82  .   AC=1;AF=0.048;AN=21;BaseQRankSum=-8.820e-01;DP=162;FS=6.990;MLEAC=1;MLEAF=0.048;MQ=36.03;MQRankSum=-1.400e-01;QD=24.88;ReadPosRankSum=-8.420e-01;SOR=4.407  GT:AD:DP:GQ:PL  0:7,0:7:99:0,118    0:3,0:3:99:0,131    0:3,0:3:90:0,90

お疲れ様でした。前半はここまでです。
後半は個別の項目について、詳しく見ていきましょう。

eupatho-bioinfomatics.hatenablog.com

よかったら他の記事もみていってください。なお本記事は、大部分をGATK公式HPの内容を元に、自分で用意したデータを使って執筆しています。

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

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

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