2020/02/14 一部修正
2020/04/27 NON_REFについて加筆
変異情報を格納するフォーマットといえば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は変異の入っていない箇所に関しては情報量を圧縮しているため、
gVCFを使えば、典型的なヒト全ゲノムシーケンス結果は1ギガバイト未満に抑えることが可能となっています。
また、コホートごとに変異を解析したい！といった用途にも有用です。
個々のサンプルについてだけでなく統合的に解析を行いたい場合、gVCFとして出力しておき上記のフローのようにgVCFをJoint Callするステップを挟むことで、コホート内で有意な変異を抽出することが可能となります。
VCFとgVCFの違い
上述した通り、重要な違いとして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> という表現があります。
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
ヘッダーの説明文からもわかるように、<NON_REF> は、このPOSにおけるあらゆる可能なアレルをとりうることを示しています。
因みに<NON_REF>変異はSelectVariantsを使用して除去することが可能です。
また、変異の存在しない行のINFOフィールドにはENDタグが含まれています。これはどのポジションでそのブロックが終了しているかを示しており、例えば1行目の変異は染色体20番10001567番目から染色体20番10001616番目までのブロックとなっています。
最後の4行のように、変異を持たない領域で連続した行が存在しうるのは、それぞれが異なるGVCFBlockに含まれるためです。
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
