2023年12月7日发(作者:)
生物基因数据文件——vcf格式详解
1.什么是vcf文件
VCF是用于描述SNP(单个碱基上的变异),INDEL(插入缺失标记)和SV(结构变异位点)结果的文本文件。在GATK软件中得到最好
的支持,当然SAMtools得到的结果也是VCF格式,和GATK的CVF格式有点差别。
的主体结构
1
##fileformat=VCFv4.2
2
##fileDate=20090805
3
##source=myImputationProgramV3.1
4
##reference=file:///seq/references/
5
##contig=
6
##phasing=partial
7
##INFO=
##INFO=
8
##INFO=
9
##INFO=
10
##INFO=
11
##INFO=
12
##FILTER=
13
##FILTER=
14
##FORMAT=
15
##FORMAT=
16
##FORMAT=
17
##FORMAT=
18
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
19
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
21
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
22
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
23
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
24
从范例上看,VCF文件分为两部分内容:以“#”开头的注释部分;没有“#”开头的主体部分。
值得注意的是,注释部分有很多对VCF的介绍信息。实际上不需要本文章,只是看看这个注释部分就完全明白了VCF各行各列代表的意义。
主体部分中每一行代表一个Variant的信息。
3.怎么解释Variation
CHROM:表示变异位点是在哪个contig 里call出来的,如果是人类全基因组的话那就是chr1…chr22,chrX,Y,M。
POS: 变异位点相对于参考基因组所在的位置,如果是indel,就是第一个碱基所在的位置。
ID: variant的ID。 如果call出来的SNP存在于dbSNP数据库里,就会显示相应的dbSNP里的rs编号;若没有,则用’.’表示其为一个novel
variant。
REF和ALT: 在这个变异位点处,参考基因组中所对应的碱基和研究对象基因组(Variant)中所对应的碱基。
QUAL: Phred格式(Phred_scaled)的质量值,可以理解为所call出来的变异位点的质量值。表 示在该位点存在variant的可能性;该值越
高,则variant的可能性越大;
计算方法:① Q=-10*lgP,Q表示质量值;P表示这个位点发生错误的概率。
②Phred值Q = -10 * lg (1-p) ,p为variant存在的概率;
通过计算公式可以看出值为10的表示错误概率为0.1,该位点为variant的概率为90%。
同理,当Q=20时,错误率就控制在了0.01。FILTER: 使用上一个QUAL值来进行过滤的话,是不够的。理想情况下,QUAL这个值应该是用所有的错误模型算出来的,这个值就可以代
表正确的变异位点了,但是事实是做不到的。因此,还需要对原始变异位点做进一步的过滤。无论你用什么方法对变异位点进行过滤,过滤
完了之后,在FILTER一栏都会留下过滤记录,如果是通过了过滤标准,那么这些通过标准的好的变异位点的FILTER一栏就会注释一个
PASS,如果没有通过过滤,就会在FILTER这一栏提示除了PASS的其他信息。如果这一栏是一个“.”的话,就说明没有进行过任何过滤。
INFO: 这一行是variant的详细信息,内容很多,以下再具体详述。
例子:
1
2
3
4
5
6
7
##fileformat=VCFv4.0
##FILTER= ##FORMAT= ##FORMAT= ##FORMAT= ##FORMAT= ##FORMAT= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO= ##INFO
#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
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
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
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
到现在,我们就可以解释上面的例子:
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”,在后续分析中可
以被过滤掉。
FORMAT 和 NA12878:这两行合起来提供了’NA12878′这个sample的基因型的信息。’NA12878′代表这该名称的样品,是由BAM文件中的
@RG下的 SM 标签决定的。
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)。这三种指
定的基因型(0/0,0/1,1/1)的概率总和为1。如果转换成支持该基因型概率(P)的话,由于L=-10lgP,那么P=10^(-L/10),因此,当L值为0
时,P=10^0=1。因此,这个值越小,支持概率就越大,也就是说是这个基因型的可能性越大。
GQ: 表示最可能的基因型的质量值。表示的意义同QUAL。Phred格式(Phred_scaled)的质量值,表示在该位点该基因型存在的可能性;该
值越高,则Genotype的可能性越 大;计算方法:Phred值 = -10 * log (1-p) p为基因型存在的概率。
举个例子说明一下:
1
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。
VCF第8列的信息
该列信息最多了,都是以 “TAG=Value”,并使用”;”分隔的形式。其中很多的注释信息在VCF文件的头部注释中给出。以下是这些TAG的解释
AC,AF 和 AN:AC(Allele Count) 表示该Allele的数目;AF(Allele Frequency) 表示Allele的频率; AN(Allele Number) 表示Allele的总数
目。对于1个diploid sample而言:则基因型 0/1 表示sample为杂合子,Allele数为1(双倍体的sample在该位点只有1个等位基因发生了突
变),Allele的频率为0.5(双倍体的 sample在该位点只有50%的等位基因发生了突变),总的Allele为2; 基因型 1/1 则表示sample为纯合
的,Allele数为2,Allele的频率为1,总的Allele为2。
DP: reads覆盖度。是一些reads被过滤掉后的覆盖度。
Dels: Fraction of Reads Containing Spanning Deletions。进行SNP和INDEL calling的结果中,有该TAG并且值为0表示该位点为SNP,
没有则为INDEL。
FS:使用Fisher’s精确检验来检测strand bias而得到的Fhred格式的p值。该值越小越好。一般进行filter的时候,可以设置 FS < 10~20。
HaplotypeScore: Consistency of the site with at most two segregating haplotypes
InbreedingCoeff: Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hard-
Weinberg expectation
MLEAC: 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
MLEAF: Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT alle in the
same order as listed
MQ: RMS Mapping Quality
MQ0: Total Mapping Quality Zero Reads
MQRankSum: Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities
QD: Variant Confidence/Quality by Depth
RPA: Number of times tandem repeat unit is repeated, for each allele (including reference)
RU: Tandem repeat unit (bases)
ReadPosRankSum: Z-score from Wilcoxon rank sum test of Alt vs. Ref read position biasSTR: Variant is a short tandem repeat
This is under continued development, please check the hts-specs page for the most recent specificationREF:
发布评论