查看原文
其他

【直播】我的基因组28-必须要理解vcf格式记录的变异位点信息

2016-12-27 生信菜鸟团 生信技能树

根据前面的讲解我们也对NGS数据分析基本的数据格式和数据流有了一个基本的理解了,不管是fasta、fastq、sam、bam、vcf、gtf 以及 biom 等格式都是具有相对标准的数据组织形式的,含有特定的信息。其中vcf就是一个尤为重要的基本数据格式,里面描述的是变异位点的信息,是一个表格形式,不同的表头下面有着不同的信息。

VCF格式本来由千人基因组计划提出来,方便描述他们找到的海量(当时是海量)变异位点。本质上也是个文本文件而已,普通编辑器打开即可。但是它对每一行每一列有具体的定义,包括文件最前面一些#开头的注释信息(这个非常重要,后面每一个位点的描述的tag都在这个注释信息里面可以找到)

详细说明可查看



1.完整地看一下vcf的头文件部分

vcf的文件表头是很长的,我们先完整地来看一遍vcf文件的head


##fileformat 首先是VCF文件的版本,该文件是GATK跑出的结果。所以此处展示的是使用GATK跑出来的数据,显示的是4.1版本。

##FILTER下面是过滤的信息,这里显示这个文件已经进行过了过滤,低质量的SNV/indel被标记上了"lowqual"。

##FORMAT 下面是,几个标记的解释,分别是:AD、DP、GQ、GT、PL的含义。

##GATKCommandLine 就是你在跑软件的时候使用的那些命令,此处略过,那些命令还是很长的。


##INFO是一些ID的信息 AC、AF、AN等,此处的#INFO的ID可以很多,后面都有解释。

##contig 是参考基因组的contig信息和参考基因组的信息,此处使用的是hg19。



vcf格式的文件内容还挺多的 ,文件头是对文件情况的一个描述和正文中一些ID的解释,他们是如何对应的?下面就让我们看一下正文的内容。

2.vcf文件的正文部分

vcf的正文部分,必须要有的是前面8列,一般来说可以有10列,分别是:

  1. #CHROM

  2. POS

  3. ID

  4. REF

  5. ALT

  6. QUAL

  7. FILTER [来自于##FILTER]

  8. INFO

  9. FORMAT

  10. 可能会有样本的名称


CHROM 和 POS:参考序列名和variant的位置;如果是INDEL的话,位置是INDEL的第一个碱基位置。

ID:variant的ID。比如在dbSNP中有该SNP的id,则会在此行给出;若没有,则用’."表示其为一个novel variant。

REF 和 ALT:参考序列的碱基 和 Variant的碱基。

QUAL:Phred格式(Phred_scaled)的质量值,表 示在该位点存在variant的可能性;该值越高,则variant的可能性越大;计算方法:Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1,该位点为variant的概率为90%。

FILTER:使用上一个QUAL值来进行过滤的话,是不够的。GATK能使用其它的方法来进行过滤,过滤结果中通过则该值为”PASS”;若variant不可靠,则该项不为”PASS”或”.”。

INFO:这一行是variant的详细信息,内容很多,以下再具体详述。

FORMAT 和 TTG11B:这两行合起来提供了’TTG11B′这个sample的基因型的信息。’TTG11B′代表这该名称的样品,是由BAM文件中的@RG下的 SM 标签决定的。


前面7列都很简单,顾名思义,分别就是该变异位点位于参考基因组的哪条染色体,哪个位置,是否被一下数据库给标记了ID(通常说的是dbSNP),该位置的参考基因组是什么碱基,这个变异位点变异成了什么碱基。找到这个变异的软件给它的质量值是多少,是否合格。下面这个表格里面我们可以看到第十列就是'realign',可以看到比对时候@RG留下来的sam的样本名称,就可以知道这个vcf是经过realign的那个bam里面call出来的突变。


vcf只学七列是远远不够的,我们有必要下功夫把较为复杂的第8列和第9列的内容好好学习一下!



第8列 INFO 就非常复杂了,该列信息最多了,看起来是一列,但是里面可以无限包容,可以根据字段拆分成多列,都是以 “TAG=Value”,并使用”;”分隔的形式。其中很多的TAG含义在VCF文件的头部注释息##INFO中已给出。

通常我们熟悉的tag有:

ACAF 和,AN[A开头的多和等位基因有关]:

    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被过滤掉后的覆盖度。[注意,第八列和第九列都有DP,都表示该位点覆盖深度的信息,但是详细意义可能是不同的大家可以探究一下,在head里面就可以找到相应信息]

Dels:Fraction of Reads Containing Spanning Deletions。进行SNP和INDEL calling的结果中,有该TAG并且值为0表示该位点为SNV,没有则为INDEL。[这个值很重要,可以根据这个tag分离indel和snv]

如果你觉得call变异的软件默认给出的tag不符合你的要求,你可以继续用其它软件在该列里面不停的增加tag,我见过给该列直接添加到180个tag的,我们后面主要讲如何来添加tag。

有了这8列,已经是标准的vcf文件了,但是大家肯定会奇怪,还没有关于这个位点的基因型,测序深度的描述的信息。

这就是属于后面的第9列FORMAT规定的了,如果有多个样本,就会按照第九列的格式不停的增加下去。

第九列可以是GT,DP,FT,GL,PL,GP等等,都可以在该vcf文件的表头里面找到关于它们的解释。前面所讲的 ##FORMAT 表头部分  便是对第九列的解释

第九列相对于第八列来说没有那么复杂的信息,数据格式是比较固定的,其中包含的信息也很重要,主要是某一个特定位点基因型,测序深度的描述,因此有必要弄清楚。


第9列数据,包含两列内容,两列内容是对应的,前者为格式,后者为格式对应的数据。

GT:样品的基因型(genotype)。两个数字中间用’/"分 开,这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele; 1 表示样品中variant的allele; 2表示有第二个variant的allele。因此: 0/0 表示sample中该位点为纯合的,和ref一致; 0/1 表示sample中该位点为杂合的,有ref和variant两个基因型; 1/1 表示sample中该位点为纯合的,和variant一致。


AD 和 DP:AD(Allele Depth)为sample中每一种allele的reads覆盖度,在diploid中则是用逗号分割的两个值,前者对应ref基因型,后者对应variant基因型; DP(Depth)为sample中该位点的覆盖度。

GQ:基因型的质量值(Genotype Quality)。Phred格式(Phred_scaled)的质量值,表示在该位点该基因型存在的可能性;该值越高,则Genotype的可能性越 大;计算方法:Phred值 = -10 * log (1-p) p为基因型存在的概率。

PL:指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1),这三种基因型的概率总和为1。和之前不一致,该值越大,表明为该种基因型的可能 性越小。 Phred值 = -10 * log (p) p为基因型存在的概率。


最需要理解的就是DP4和GT了:

第十列的话就是样本的信息 可以在比对的时候使用@RG来做一个标记


vcf格式记录变异位点的信息大该就讲到了这里,有什么不清楚的地方留言讨论哦!

文:Jimmy、阿尔的太阳

图文编辑:吃瓜群众




您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存