在用于存储reads比对结果的SAM/BAM格式提出之后（2009年），用于存储变异检测结果的格式VCF( variant call format )也被提出（2010年），与此同时作者还提供了配套的处理工具——vcftools。
1. vcf文件
vcf主要包括两部分，以##和#开头的header信息，以及主体variation信息。
1.1 header部分
以#开头的一行（header部分的最后一行）是主体部分的列名，以##开头的行是一些描述解释信息，比如主体中的"FILTER", "FORMAT", "INFO"都是什么意思，另外还能看到一些历史命令，通过这些命令可以知道这个vcf文件是如何得到的，比如，通过上面的图片可知，这个vcf是GATK和bcftools分别call variation之后取交集得到。
要完全理解vcf的每一列，每一个英文大写的缩写名词其实挺困难的，仅INFO就有很多名词，有的可能还需要高深的统计学知识，所以目前我只掌握的一些简单的。
下面我以我的vcf文件为例，尽可能多的理解一下vcf中字段的含义。
1.2 variation部分
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample1 sample2 sample3 sample4 sample5 sample6 sample7 sample8
chrA01 895 . T A 101.53 PASS AC=2;AF=0.125;AN=16;BaseQRankSum=-2.261;ClippingRankSum=0;DP=239;ExcessHet=3.3099;FS=4.905;MLEAC=2;MLEAF=0.125;MQ=80.42;MQRankSum=0;QD=2.21;ReadPosRankSum=-0.033;SOR=1.609 GT:AD:DP:GQ:PGT:PID:PL 0/0:22,0:22:66:0|1:895_T_A:0,66,1689 0/0:36,2:38:65:.:.:0,65,1576 0/1:10,3:13:67:.:.:67,0,396 0/0:19,0:19:51:.:.:0,51,765 0/0:31,0:31:84:.:.:0,84,1260 0/1:28,5:33:74:.:.:74,0,1087 0/0:29,0:29:81:.:.:0,81,1215 0/0:32,0:32:90:.:.:0,90,1350
chrA01 912 . G A 181.85 PASS AC=3;AF=0.188;AN=16;BaseQRankSum=-1.042;ClippingRankSum=0;DP=234;ExcessHet=3.9794;FS=2.43;MLEAC=3;MLEAF=0.188;MQ=74.43;MQRankSum=0;QD=2.64;ReadPosRankSum=-0.339;SOR=0.857 GT:AD:DP:GQ:PGT:PID:PL 0/0:20,0:20:60:0|1:895_T_A:0,60,1619 0/0:36,1:37:66:.:.:0,66,1496 0/1:10,3:13:41:.:.:41,0,385 0/0:16,0:16:42:.:.:0,42,630 0/0:33,0:33:90:.:.:0,90,1350 0/1:27,5:32:82:.:.:82,0,1046 0/1:19,5:24:99:.:.:102,0,631 0/0:35,0:35:99:.:.:0,99,1485
chrA01 921 . T C 2635.31 PASS AC=6;AF=0.375;AN=16;BaseQRankSum=1.9;ClippingRankSum=0;DP=221;ExcessHet=2.98;FS=0;MLEAC=6;MLEAF=0.375;MQ=59.33;MQRankSum=0;QD=24.86;ReadPosRankSum=1.15;SOR=0.73 GT:AD:DP:GQ:PGT:PID:PL 0/1:10,7:17:99:1|0:895_T_A:264,0,724 0/0:36,0:36:99:.:.:0,99,1485 0/1:3,9:12:73:.:.:314,0,73 0/1:10,7:17:99:0|1:921_T_C:286,0,361 1/1:0,39:39:99:1|1:921_T_C:1751,117,0 0/0:28,0:28:75:.:.:0,75,1125 0/1:18,3:21:72:.:.:72,0,1015 0/0:35,0:35:99:.:.:0,99,1485
chrA01 931 . A G 3003.26 PASS AC=8;AF=0.5;AN=16;BaseQRankSum=0.645;ClippingRankSum=0;DP=219;ExcessHet=2.5225;FS=2.735;MLEAC=8;MLEAF=0.5;MQ=66.13;MQRankSum=0;QD=22.75;ReadPosRankSum=-0.772;SOR=0.929 GT:AD:DP:GQ:PGT:PID:PL 0/1:10,7:17:99:0|1:895_T_A:261,0,767 0/0:34,1:35:60:.:.:0,60,1399 1/1:0,11:11:33:.:.:464,33,0 0/1:9,7:16:99:0|1:921_T_C:267,0,357 1/1:0,41:41:99:1|1:921_T_C:1867,126,0 0/1:24,3:27:54:.:.:54,0,974 0/1:15,5:20:99:.:.:144,0,4730/0:35,0:35:90:.:.:0,90,1350
chrA01 1093 . A G 1117.16 PASS AC=3;AF=0.188;AN=16;BaseQRankSum=-1.18;ClippingRankSum=0;DP=163;ExcessHet=0.4576;FS=11.522;MLEAC=3;MLEAF=0.188;MQ=49.72;MQRankSum=-1.559;QD=29.67;ReadPosRankSum=0.015;SOR=1.008 GT:AD:DP:GQ:PGT:PID:PL 0/1:3,6:9:99:0|1:1093_A_G:243,0,108 0/0:32,0:32:90:.:.:0,90,13500/0:9,0:9:27:.:.:0,27,405 0/0:13,0:13:36:.:.:0,36,540 1/1:0,21:21:63:1|1:1093_A_G:923,63,0 0/0:34,0:34:90:.:.:0,90,1350 0/0:17,0:17:9:.:.:0,9,527 0/0:25,0:25:72:.:.:0,72,1080
chrA01 1103 . C T 1122.24 PASS AC=3;AF=0.188;AN=16;BaseQRankSum=0.903;ClippingRankSum=0;DP=164;ExcessHet=0.4576;FS=10.474;MLEAC=3;MLEAF=0.188;MQ=51.31;MQRankSum=-1.383;QD=34.04;ReadPosRankSum=-0.319;SOR=0.727 GT:AD:DP:GQ:PGT:PID:PL 0/1:3,7:10:99:0|1:1093_A_G:285,0,105 0/0:29,0:29:81:.:.:0,81,12150/0:11,0:11:33:.:.:0,33,456 0/0:12,0:12:33:.:.:0,33,495 1/1:0,20:20:60:1|1:1093_A_G:886,60,0 0/0:35,0:35:99:.:.:0,99,1485 0/0:18,0:18:8:.:.:0,8,557 0/0:26,0:26:66:.:.:0,66,990
| 字段 | 含义 |
|---|---|
| CHROM | 染色体名称 |
| POS | 变异位点在参考基因组上的位置，染色体开始的位置编号为1 |
| ID | 变异检测过程中如果引入外部数据库，可能会有数据库中的位点ID |
| REF | 参考碱基（序列），在Genotype中用0表示 |
| ALT | 变异碱基（序列），不止一个时，逗号分隔，表示为1,2,3等 |
| QUAL | 检测质量值，越大越好，等于-10*log10(该变异位点检测错误的概率) |
| FILTER | 根据自定义阈值设置字符串标记，再根据标记挑选（比如使用GATK的VariantFiltration和SelectVariants模块） |
| INFO | 注释信息 |
| FORMAT | 该列与其之后的样本Genotype信息一一对应，比如GT:AD:DP:GQ:PGT:PID:PL 0/0:22,0:22:66:0|1:895_T_A:0,66,1689
|
1.2.1 INFO
AC=3;AF=0.188;AN=16;BaseQRankSum=-0.601;ClippingRankSum=0;DP=209;ExcessHet=0.4576;FS=0;MLEAC=3;MLEAF=0.188;MQ=56.28;MQRankSum=0.319;QD=34.04;ReadPosRankSum=1.67;SOR=0.569
| 键 | 值 | 解释 |
|---|---|---|
| AC | 3 | 二倍体，8个样本，该位点总共出现16次，该碱基出现3次 |
| AF | 0.188 | 该碱基的频率，3/16 |
| AN | 16 | 该位点总数 |
| DP | 209 | 8个样本的该位点深度之和 |
| MQ | 56.28 | RMS Mapping Quality |
1.2.3 FORMAT
#CHROM POS ID REF ALT QUAL FILTER FORMAT sample1
chrA01 895 . T A 101.53 PASS GT:AD:DP:GQ:PGT:PID:PL 0/0:22,0:22:66:0|1:895_T_A:0,66,1689
chrA01 912 . G A 181.85 PASS GT:AD:DP:GQ:PGT:PID:PL 0/0:20,0:20:60:0|1:895_T_A:0,60,1619
chrA01 921 . T C 2635.31 PASS GT:AD:DP:GQ:PGT:PID:PL 0/1:10,7:17:99:1|0:895_T_A:264,0,724
chrA01 931 . A G 3003.26 PASS GT:AD:DP:GQ:PGT:PID:PL 0/1:10,7:17:99:0|1:895_T_A:261,0,767
chrA01 1093 . A G 1117.16 PASS GT:AD:DP:GQ:PGT:PID:PL 0/1:3,6:9:99:0|1:1093_A_G:243,0,108
chrA01 1103 . C T 1122.24 PASS GT:AD:DP:GQ:PGT:PID:PL 0/1:3,7:10:99:0|1:1093_A_G:285,0,105
GT: Genotype, /分隔表示未定相, |表示已定相
PGT: phased Genotype, 邻近的位点是如何定相排列的，以最后两列为例：
chrA01 1093 A G 0|1
chrA01 1103 C T 0|1
竖线左边的碱基来源于一条染色体，竖线右边的碱基来源于另一条染色体，即A,C来自相同的染色体；G,T来源于另一条染色体。定相这一步需要结合多个样本的信息！
PID: 一个phasing group会有一个共用的group ID，并且这个ID是以第一个phased位点来定义的，比如上面两个phased位点组成的group的ID是1093_A_G。
AD: 每一个allele的深度
DP: 所有allele的总深度
GQ: Genotype Quality, -10*log10(基因型错误的概率)
PL: 每种基因型的似然值，越小越可信，比如最后一列285,0,105，分别对应基因型0/0,0/1,1/1，说明0/1为可能的基因型。
reference
The variant call format and VCFtools: https://academic.oup.com/bioinformatics/article/27/15/2156/402296
VCFtools--A set of tools written in Perl and C++ for working with VCF files:
https://vcftools.github.io/man_latest.html#EXAMPLES
VCFv4.2详解（英文版）
https://samtools.github.io/hts-specs/VCFv4.2.pdf
