文章来源:http://www.cnblogs.com/emanlee/p/4562064.htmlhtml
VCF文件示例(VCFv4.2)git
##fileformat=VCFv4.2 ##fileDate=20090805 ##source=myImputationProgramV3.1 ##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x> ##phasing=partial ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data"> ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth"> ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency"> ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele"> ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129"> ##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership"> ##FILTER=<ID=q10,Description="Quality below 10"> ##FILTER=<ID=s50,Description="Less than 50% of samples have data"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality"> ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth"> ##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003 20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,. 20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3 0/0:41:3 20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2 2/2:35:4 20 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2 20 1234567 microsat1 GTC G,GTCT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3
CHROM: 表示变异位点是在哪一个contig 里call出来的,若是是人类全基因组的话那就是chr1…chr22,chrX,Y,M。github
POS: 变异位点相对于参考基因组所在的位置,若是是indel,就是第一个碱基所在的位置。数据库
ID: 若是call出来的SNP存在于dbSNP数据库里,就会显示相应的dbSNP里的rs编号。api
REF和REF: 在这个变异位点处,参考基因组中所对应的碱基和研究对象基因组中所对应的碱基。ui
QUAL: 能够理解为所call出来的变异位点的质量值。Q=-10lgP,Q表示质量值;P表示这个位点发生错误的几率。所以,若是想把错误率从控制在90%以上,P的阈值就是1/10,那lg(1/10)=-1,Q=(-10)*(-1)=10。同理,当Q=20时,错误率就控制在了0.01。spa
FILTER: 理想状况下,QUAL这个值应该是用全部的错误模型算出来的,这个值就能够表明正确的变异位点了,可是事实是作不到的。所以,还须要对原始变异位点作进一步的过滤。不管你用什么方法对变异位点进行过滤,过滤完了以后,在FILTER一栏都会留下过滤记录,若是是经过了过滤标准,那么这些经过标准的好的变异位点的FILTER一栏就会注释一个PASS,若是没有经过过滤,就会在FILTER这一栏提示除了PASS的其余信息。若是这一栏是一个“.”的话,就说明没有进行过任何过滤。.net
例子:3d
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878 chr1 873762 . T G 5231.78 PASS AC=1;AF=0.50;AN=2;DP=315;Dels=0.00;HRun=2;HaplotypeScore=15.11;MQ=91.05;MQ0=15;QD=16.61;SB=-1533.02;VQSLOD=-1.5473 GT:AD:DP:GQ:PL 0/1:173,141:282:99:255,0,255 chr1 877664 rs3828047 A G 3931.66 PASS AC=2;AF=1.00;AN=2;DB;DP=105;Dels=0.00;HRun=1;HaplotypeScore=1.59;MQ=92.52;MQ0=4;QD=37.44;SB=-1152.13;VQSLOD= 0.1185 GT:AD:DP:GQ:PL 1/1:0,105:94:99:255,255,0 chr1 899282 rs28548431 C T 71.77 PASS AC=1;AF=0.50;AN=2;DB;DP=4;Dels=0.00;HRun=0;HaplotypeScore=0.00;MQ=99.00;MQ0=0;QD=17.94;SB=-46.55;VQSLOD=-1.9148 GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26 chr1 974165 rs9442391 T C 29.84 LowQual AC=1;AF=0.50;AN=2;DB;DP=18;Dels=0.00;HRun=1;HaplotypeScore=0.16;MQ=95.26;MQ0=0;QD=1.66;SB=-0.98 GT:AD:DP:GQ:PL 0/1:14,4:14:60.91:61,0,255
到如今,咱们就能够解释上面的例子:code
chr1:873762 是一个新发现的T/G变异,而且有很高的可信度(qual=5231.78)。
chr1:877664 是一个已知的变异为A/G 的SNP位点,名字rs3828047,而且具备很高的可信度(qual=3931.66)。
chr1:899282 是一个已知的变异为C/T的SNP位点,名字rs28548431,但可信度较低(qual=71.77)。
chr1:974165 是一个已知的变异为T/C的SNP位点,名字rs9442391,可是这个位点的质量值很低,被标 成了“LowQual”,在后续分析中能够被过滤掉。
Vcf文件看起来很复杂,挺吓人的样子,可是里面大部分都是一些tags,而这些tags基本上都是在VASR中过滤用的,可以理解每一个tags的意思最好,若是实在不理解也就不用管了。其实最关键的信息也就是那么几列:
chr1 873762 . T G [CLIPPED] GT:AD:DP:GQ:PL 0/1:173,141:282:99:255,0,255
chr1 877664 rs3828047 A G [CLIPPED] GT:AD:DP:GQ:PL 1/1:0,105:94:99:255,255,0
chr1 899282 rs28548431 C T [CLIPPED] GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26
其中最后面两列是相对应的,每个tag对应一个或者一组值,如:
chr1:873762,GT对应0/1;AD对应173,141;DP对应282;GQ对应99;PL对应255,0,255。
GT: 表示这个样本的基因型,对于一个二倍体生物,GT值表示的是这个样本在这个位点所携带的两个等位基因。0表示跟REF同样;1表示表示跟ALT同样;2表示第二个ALT。当只有一个ALT 等位基因的时候,0/0表示纯和且跟REF一致;0/1表示杂合,两个allele一个是ALT一个是REF;1/1表示纯和且都为ALT; The most common format subfield is GT (genotype) data. If the GT subfield is present, it must be the first subfield. In the sample data, genotype alleles are numeric: the REF allele is 0, the first ALT allele is 1, and so on. The allele separator is '/' for unphased genotypes and '|' for phased genotypes.
0 - reference call
1 - alternative call 1
2 - alternative call 2
AD: 对应两个以逗号隔开的值,这两个值分别表示覆盖到REF和ALT碱基的reads数,至关于支持REF和支持ALT的测序深度。
DP: 覆盖到这个位点的总的reads数量,至关于这个位点的深度(并非多有的reads数量,而是大概必定质量值要求的reads数)。
PL: 对应3个以逗号隔开的值,这三个值分别表示该位点基因型是0/0,0/1,1/1的没通过先验的标准化Phred-scaled似然值(L)。若是转换成支持该基因型几率(P)的话,因为L=-10lgP,那么P=10^(-L/10),所以,当L值为0时,P=10^0=1。所以,这个值越小,支持几率就越大,也就是说是这个基因型的可能性越大。
GQ: 表示最可能的基因型的质量值。表示的意义同QUAL。
举个例子说明一下:
chr1 899282 rs28548431 C T [CLIPPED] GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26
在这个位点,GT=0/1,也就是说这个位点的基因型是C/T;GQ=25.92,质量值并不算过高,多是由于cover到这个位点的reads数太少,DP=4,也就是说只有4条reads支持这个地方的变异;AD=1,3,也就是说支持REF的read有一条,支持ALT的有3条;在PL里,这个位点基因型的不肯定性就表现的更突出了,0/1的PL值为0,虽然支持0/1的几率很高;可是1/1的PL值只有26,也就是说还有10^(-2.6)=0.25%的可能性是1/1;但几乎不多是0/0,由于支持0/0的几率只有10^(-10.3)=5*10-11。
The VCF specification is no longer maintained by the 1000 Genomes Project. The group leading the management and expansion of the format is the Global Alliance for Genomics and Health Data Working group file format team, http://ga4gh.org/#/fileformats-team
The main version of the specification can be found on https://github.com/samtools/hts-specs
This is under continued development, please check the hts-specs page for the most recent specification
A PDF of the v4.1 spec is http://samtools.github.io/hts-specs/VCFv4.1.pdf
A PDF of the v4.2 spec is http://samtools.github.io/hts-specs/VCFv4.2.pdf
VCFTools host a discussion list about the specification called vcf-spec http://sourceforge.net/p/vcftools/mailman/
REF:
http://blog.sina.com.cn/s/blog_12d5e3d3c0101qv1u.html