全基因组(WGS) 思维导图
  1. 首页
  2. 生信组学教程
  3. 全基因组(WGS)
  4. VCF详细介绍

VCF详细介绍

最后发布时间 : 2023-08-11 17:44:15 浏览量 :

学习资料

VCF头部信息

##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele not already represented at this location by REF and ALT">
##FILTER=<ID=LowQual,Description="Low quality">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##GATKCommandLine=<ID=CombineGVCFs,CommandLine="CombineGVCFs --output output/gvcf/cohort.g.vcf.gz --variant output/gvcf/mother.sorted.markdup.g.vcf --referen>
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##contig=<ID=20,length=63025520>
  • 第一行是VCF的版本信息；
  • FILTER行是说过滤了什么内容；
  • FORMAT和INFO 相当于变异位点的注释信息；
  • CommandLine 是说使用的call variant工具信息；
  • contig和reference 是当重复别人数据时，恰好没告诉你数据来源，这时就可以参考这个。

VCF记录的信息

这一部分至少包括8列tab分割的常规信息（用来描述变异位点），第9列及以后表示各个样本的变异信息（可以包括成百上千个样本）。

生信小木屋

  • CHROM：变异发生的chromosome或者contig
  • POS：变异发生的基因组坐标（对缺失来讲，显示的是缺失开始的位置）
  • ID：一般是dbSNP的ID（可有可无）
  • REF：Forward strand（正链）上的参考等位基因
  • ALT：正链上对应的改变的等位基因（可以有多个）
  • QUAL：REF/ALT 变异位点存在的概率（和FASTQ的质量值、SAM的MAPQ一样，都是Phred值 -10 * log(1-p)）
  • FILTER：数据一般都要经过适当的过滤后才能继续使用variant callset（变异位点数据集），关于是否完成过滤会给出三种说明： 一是给出没有通过过滤的变异位点；二是PASS表示全部通过了过滤；三是. 表示这个位点没有任何过滤
  • INFO：以tag=value的形式给出位点注释信息；分号;分割
  • FORMAT：针对样本的注释；冒号:分割，并且对应后面各个sample中的信息

补充一下正链上的参考等位基因的知识：

生信小木屋

具体的 REF/ALT

开始想象场景：我们用软件发现，在60bp处发现了一个变异，是A变成了T，那么应该这么表示：

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	align.bam
AF086833	60	.	A	T	54	.	DP=43	GT:PL	0/1:51,0,48

如果同一个位置检测到有两种可能发生的变异呢？

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	align.bam
AF086833	60	.	A	T,C	43.2	.	DP=95	GT:PL	1/2:102,124,42,108,0,48

事情还没完，刚刚两个例子都只是一个变异位点，但实际上有可能多个位点发生缺失（而这才是VCF复杂的开始），例如：58、59、60碱基（GGA）发生了缺失

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	align.bam
AF086833	55	.	TGAGGA	TGA	182	.	INDEL;IDV=42	GT:PL	0/1:215,0,255
AF086833	60	.	A	C,T	39.8	.	DP=126	GT:PL	1/2:99,87,47,86,0,49

这里注意：虽然我们看到缺失是从58碱基开始发生的，但是记录的是从55碱基。你会不会好奇：为什么是55而不是准确的58？为什么要表示成TGAGGA->TGA，而不是GAGGA->GA或者GGA->空 呢？这个需要引入一个新词汇”variant normalizatoin“，也就是说，所有的结果展示都是有规定的。

VCF的标准化，为了方便交流，规定以下几点：

  • 用尽可能少的字母来表示变异位点；
  • 等位基因长度不为0；
  • 变异向左”贪婪比对“ （也就是说：一直向左比对，直到不匹配为止，然后以最左边的碱基位置表示变异的起始位置）

因此，这里写成TGAGGA->TGA就是表示缺失位点GGA向左最远可以匹配到第55号T碱基处

具体 FORMAT

GT:AD:DP:GQ:PL  0/1:17,15:32:99:399,0,439       0/1:11,12:23:99:291,0,292
GT:AD:DP:GQ:PL  0/1:13,10:23:99:243,0,341       0/1:21,14:35:99:355,0,526

从VCF的头部可以看到这些字段的解释如下：

  • GT：Genotype
  • AD: Allelic depths for the ref and alt alleles in the order listed
  • DP: Approximate read depth (reads with MQ=255 or with bad mates are filtered)
  • GQ: Genotype Quality
  • PL: List of Phred-scaled genotype likelihoods

关于PL与GQ的解释可以参考GATK4 HaplotypeCaller

生信小木屋

GT的解释

尽管我们可以根据REF和ALT知道了碱基发生了怎样的变化，但是我们想知道这个变化是只是发生在一个DNA拷贝中，还是两个拷贝都有？需要用一个指标来量化这种变异~基因型（Genotype）。GT就是用来表示样本中这个位点的基因型，其中0表示参考REF，1表示变异ALT的第一个entry，2表示ALT的第二个entry（以此类推）

对于二倍体生物，GT表示了一个样本中的等位基因：

  • 0/0 表示样本是纯合子，并且和参考的等位基因一样
  • 0/1表示样本是杂合子，有一个参考的等位基因，一个变异的等位基因
  • 1/2 样本是杂合子，两个都是变异的等位基因
  • 1/1 样本是纯合子，且两个都是变异的第一个等位基因
  • 2/2样本是纯合子，且两个都是变异的第二个等位基因（以此类推）

当然，如果不是二倍体，命名原理也是一样：单倍体（Haploid）只有一个GT值；多倍体有多个GT值

从INFO列也可以看出等位基因的相关信息（A开头的大多和等位基因Allele相关）

AC=2;AF=0.5;AN=4;
AC=3;AF=0.75;AN=4;
  • AC（Allele Count）：该Allele数目
  • AF(Allele Frequency)：该Allele频率
  • AN(Allel Number)：Allele总数目

对于二倍体样本（Diploid sample），基因型为0/1表示杂合子：AC=1表示Allele数为1（即该位点只有1个等位基因发生变异）；AF=0.5频率为0.5（该位点只有50%的等位基因发生变异）；AN=2表示总Allele为2

基因型为1/1表示纯合子，AC=2,AF=1,AN=2

正常人的二倍体基因组位点只有杂合和纯合两种情况，因此杂合AF一定是0.5，纯合AF一定是1，但实际的VCF数据是通过测序数据来的，随机性加上测序深度不够带来的系统误差，就使得结果变得不那么理想

PL的解释

List of Phred-scaled genotype likelihoods，直白来说就是一组等位基因的可能性，用来衡量不同基因型可能发生的概率，该值越小可能性越大，例如：

GT:PL  0/1:99:399,0,439

其中PL这一项有三个数值，分别对应三种可能的基因型（0/0，0/1，1/1）发生的概率：

  • 第一个数值291表示基因型为0/0的概率是Phred值291，表示-10*log10( 10^(-29) )=290，及基因型是0/0纯合子的概率大约是10^{-29}；
  • 第二个数值0表示基因型为0/1的概率是Phred值0，表示-10*log10( 1 )=0，及基因型是0/1杂合子的概率大约是1；
  • 第三个数值439表示基因型为1/1的概率是Phred值439，表示-10*log10( 10^(-44) )=440，及基因型是1/1纯合子的概率大约是10^{-44}。

AD与DP的解释

软件判断是那种基因型，到底是不是发生了变异，是需要一定的统计方法的，主体就是之前比对的结果BAM文件，其中包含了reads的比对信息，这里就是根据reads比对的数量进行判断

GT:AD:DP  0/1:17,15
  • AD（DepthPerAlleleBySample）：unfiltered allele depth 就是有多少reads出现了某个等位基因（其中也包含了没有经过variant caller过滤的reads），但是不包括没意义的reads（就是那些统计结果不过关的，没法说服软件相信这个等位基因）。在二倍体样本中表示为：逗号分隔的两个值，前者为对应的ref基因型，后者为对应alt的基因型
  • DP：（Approximate read depth ）：(Coverage，即reads覆盖度，指一些reads被过滤后的覆盖度)filtered depth, at the sample level 只有通过variant caller软件过滤后的reads才能计算入内，但是DP也纳入了那些经过过滤但没有意义的reads（uninformative reads），这一点又和AD不同

具体 Genotype Quality

GQ就是用Phred值来表示GT判断的准确性，它和PL相似，但是取值不同。PL最小值0表示最准确，GQ一般取PL的第二个小的值（除非第二小的PL大于99）。在GATK分析中，GQ最大就限制在99，因为超过99其实是没有什么意义的，并且还多占字符。因此，如果GATK中发现PL值中第二个小的值比99还要大，软件就将GQ标为99。用GQ值就可以得到，第一位和第二位之间到底差了多少，因此可以快速判断分析的准不准，选择第一个靠不靠谱

做一个小总结–VCF帮助判断基因型

1   899282  rs28548431  C   T   [CLIPPED] GT:AD:DP:GQ:PL    0/1:1,3:4:26:103,0,26

这个位点的GT=0/1，可以判断基因型是C/T；GQ=26表示排名第二基因型的可能性是0.0025，结果不是很好，因为虽然判断的第一位基因型的PL 为0很可靠，但是毕竟相差不多，很难推翻第二位（如果GQ值再大一些，我们就更有信心说明判断的C/T基因型是正确的）。当然，这里GQ的原因很有可能是取样太少，只有4条reads在这个位点作为参考（DP=4），这四条中有1条带参考的碱基信息，另外3条与参考不一致，存在变异（AD=1,3）

因此，重点的结论来啦！尽管我们相信，这个位点确实存在变异，但是假阳性依然存在，也就意味着基因型判断结果不一定是杂合子C/T，有一定的可能是变异纯合T/T(PL(1/1)=26)，但一定不可能是参考纯合C/C （PL(0/0)=103）

bcftools 工具介绍 变异注释