gVCFとはなんなのか

2020/02/14 一部修正

変異情報を格納するフォーマットといえばVCFファイルが一般的なわけですが、 似たようなフォーマットにgVCFというものがあります。 VCF自体自由度の高いフォーマットなこともあり、gVCFも使用するツールによって色々なフォーマットがありそうです。

今回はそんなgVCFについて調べてみました。

概要

gVCFgVCFってなんのこっちゃと思われるかもしれません。
Genomic Variant Call Formatの略です。基本的なフォーマットは通常のVCFと一緒です。

VCFについてはhts-specsから飛べる詳細PDFを参照してください。

VCFが単に変異を格納するフォーマットであったのに対し、 gVCFは、変異が検出されたされていないに関わらず、ゲノム上の全てのポジション(もしくは関心のある範囲)についてのレコードを含むことができるVCFを指します。

変異を格納するフォーマットなだけあって変異検出ソフトウェアの出力フォーマットや入力フォーマットとして規定されるわけですが、例えばGATK 3系と4系で出力されるgVCFは混合するのは避けるべきとのことです(下記で説明するGVCFBlockとアノテーション列が同様のフォーマットをしていれば、混合して使用しても恐らく問題ないはず)。

用途

通常のVCFは、変異が入った場所のみがレコードとして保持されます。一方で、臨床応用においては、変異が入っていない箇所(non-variant)の情報が重要になる場合もあるわけです。
この要因として、VCFファイルは1ポジションについて1レコードをとるため、情報量が大きくなってしまうことが挙げられます。

gVCFは変異の入っていない箇所に関しては情報量を圧縮しているため、 gVCFを使えば、典型的なヒト全ゲノムシーケンス結果は1ギガバイト未満に抑えることが可能となっています。

f:id:kimoppy126:20200213082044p:plain

また、コホートごとに変異を解析したい!といった用途が挙げられます。
個々のサンプルについてだけでなく統合的に解析を行いたい場合、gVCFとして出力しておき上記のフローのようにgVCFをJoint Callするステップを挟むことで、コホート内で有意な変異を抽出することが可能となります。

VCFとgVCFの違い

f:id:kimoppy126:20191030111946p:plain

上述した通り、重要な違いとしてgVCFは変異検出の有無に関わらず全てのポジションのレコードを含んでいる点が大きく異なります。
後続のステップでコホートの統合解析を行うために、ファイルに全サイトを表記することが必要となるわけです。
gVCFのレコードには、そのレコードがホモ接合であるかどうかの判断にどれだけ自信があるかの情報が含まれます。この推定は、HaplotypeCaller built-in のリファレンスモデルを用いて生成されます。

他の一部のツール(GATKのUnifiedGenotyperを含む)やBP_RESOLUTIONを用いて出力されるような個々のポジションに関して行を分割したVCFではこの情報が格納されないため、Joint Genotypingを行うことができません。

GVCF関連ツール

GATK HaplotypeCaller

GATK HaplotypeCallerでは、実行時に -ERC GVCF もしくは-ERC BP_RESOLUTIONのどちらかを指定することでgVCFを出力することができます。

-ERC BP_RESOLUTIONを用いた場合、変異検出の有無に関わらず全ての行が別々の行として表記されます。
一方-ERC GVCFを用いた場合、変異が検出されたポジションに関しては全て別々の行として表記されるが、変異が検出されなかったポジションに関しては、ポジションの領域としてGQごとにグループ化されて表記されます。

GQの範囲はgVCF中のヘッダー行##GVCFBlockで定義されます。 -ERC GVCF を使用することでファイルサイズを節約することができるため、こちらを使用することが推奨されています。

GATK UnifiedGenotyper

現在ではHaplotypeCallerに置き換わったことにより非推奨となっていますが、GATK UnifiedGenotyperでは--output_mode EMIT_ALL_SITESオプションを指定することでgVCFを出力することができます。
参照:Should I use UnifiedGenotyper or HaplotypeCaller to call variants on my data? — GATK-Forum

Illumina Strelka

Germline variant reporting uses the gVCF conventions to represent both variant and reference call confidence.

Illumina社のStrelkaではgermline変異検出に関して出力がgVCFとなっています。

Isaac Variant Caller

The output of the workflow is a single Genome VCF (gVCF) file describing all sites and indels in the genome. gVCF is a set of conventions for output in the VCF 4.1 format to represent all sites in the genome.

旧Strelkaです。この時から出力フォーマットはgVCFでした。

gvcftools

GATK Unified Genotyperの全サイトを含むGVCFからblockごとに分けられたGVCFを出力できたり(gatk_to_gvcf)、 複数のgVCFを1つにまとめることができたり(merge_variants)、gVCFをVCFに変換できたり(extract_variants)するツールです。

gVCFに特異的なヘッダー

FORMATフィールドのMIN_DP

gVCFの各ブロックにおける最低リード数を定義しています。

##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">

GVCFBlock

GVCFBlockは下記のようになっています。
GATKの最近のバージョンでは1-59までのGQについて、高い精度でブロックを分け、 60以上のGQについてはブロックサイズを大きくとることで情報量を圧縮しているようです。

[individual blocks from 1 to 55]
##GVCFBlock55-56=minGQ=55(inclusive),maxGQ=56(exclusive)
##GVCFBlock56-57=minGQ=56(inclusive),maxGQ=57(exclusive)
##GVCFBlock57-58=minGQ=57(inclusive),maxGQ=58(exclusive)
##GVCFBlock58-59=minGQ=58(inclusive),maxGQ=59(exclusive)
##GVCFBlock59-60=minGQ=59(inclusive),maxGQ=60(exclusive)
##GVCFBlock60-70=minGQ=60(inclusive),maxGQ=70(exclusive)
##GVCFBlock70-80=minGQ=70(inclusive),maxGQ=80(exclusive)
##GVCFBlock80-90=minGQ=80(inclusive),maxGQ=90(exclusive)
##GVCFBlock90-99=minGQ=90(inclusive),maxGQ=99(exclusive)
##GVCFBlock99-100=minGQ=99(inclusive),maxGQ=100(exclusive)

GVCFの例

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  NA12878
20  10001567    .   A   <NON_REF>   .   .   END=10001616    GT:DP:GQ:MIN_DP:PL  0/0:38:99:34:0,101,1114
20  10001617    .   C   A,<NON_REF> 493.77  .   BaseQRankSum=1.632;ClippingRankSum=0.000;DP=38;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.000;RAW_MQ=136800.00;ReadPosRankSum=0.170    GT:AD:DP:GQ:PL:SB   0/1:19,19,0:38:99:522,0,480,578,538,1116:11,8,13,6
20  10001618    .   T   <NON_REF>   .   .   END=10001627    GT:DP:GQ:MIN_DP:PL  0/0:39:99:37:0,105,1575
20  10001628    .   G   A,<NON_REF> 1223.77 .   DP=37;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=133200.00   GT:AD:DP:GQ:PL:SB   1/1:0,37,0:37:99:1252,111,0,1252,111,1252:0,0,21,16
20  10001629    .   G   <NON_REF>   .   .   END=10001660    GT:DP:GQ:MIN_DP:PL  0/0:43:99:38:0,102,1219
20  10001661    .   T   C,<NON_REF> 1779.77 .   DP=42;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=151200.00   GT:AD:DP:GQ:PGT:PID:PL:SB   1/1:0,42,0:42:99:0|1:10001661_T_C:1808,129,0,1808,129,1808:0,0,26,16
20  10001662    .   T   <NON_REF>   .   .   END=10001669    GT:DP:GQ:MIN_DP:PL  0/0:44:99:43:0,117,1755
20  10001670    .   T   G,<NON_REF> 1773.77 .   DP=42;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=151200.00   GT:AD:DP:GQ:PGT:PID:PL:SB   1/1:0,42,0:42:99:0|1:10001661_T_C:1802,129,0,1802,129,1802:0,0,25,17
20  10001671    .   G   <NON_REF>   .   .   END=10001673    GT:DP:GQ:MIN_DP:PL  0/0:43:99:42:0,120,1800
20  10001674    .   A   <NON_REF>   .   .   END=10001674    GT:DP:GQ:MIN_DP:PL  0/0:42:96:42:0,96,1197
20  10001675    .   A   <NON_REF>   .   .   END=10001695    GT:DP:GQ:MIN_DP:PL  0/0:41:99:39:0,105,1575
20  10001696    .   A   <NON_REF>   .   .   END=10001696    GT:DP:GQ:MIN_DP:PL  0/0:38:97:38:0,97,1220

ALT列に<NON_REF> という表現があります。これこそがGVCF特有の表記で、このポジションにおけるリファレンスがないアレルを表現しています。

また、変異の存在しない行のINFOフィールドにはENDタグが含まれています。これはどのポジションでそのブロックが終了しているかを示しており、例えば1行目の変異は染色体20番10001567番目から染色体20番10001616番目までのブロックとなっています。

最後の4行のように、変異を持たない領域で連続した行が存在しうるのは、各ブロックにおいて異なるGQ (Genotype Quality)を持っているためです。

gVCFに関するQ&A

Q. gVCFを入力にsnpEff等で直接アノテーションを行うことは可能なのか。
A. gVCFはあくまで中間ファイルです。GenotypeGVCFsを使ってVCFに変換してからアノテーションは行う必要があります。

Q. readが全くアライメントされなかった領域はどのように表記されるのか。gVCFに含まれるのか。
A. 下記の2つ目の変異のようにすべての値が0で表記されます。GenotypeGVCFsを使ってVCFに変換したときに 未検出を表すGT(./.)に変換されます。

1 668406 . G . . END=668423 GT:DP:GQ:MIN_DP:PL 0/0:11:21:8:0,21,281
1 668424 . T . . END=668553 GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0
1 668554 . A . . END=668563 GT:DP:GQ:MIN_DP:PL 0/0:7:21:7:0,21,217

参考