VCF文件全称为Variant Call Format,通常为GATK和Samtools软件处理所得到。在没接触过的人乍一看确实很吓人,但是只要理解了其中的一些关键tag再进行过滤就可以有的放矢了,下面以如下例子做讲解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##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=.,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 GTCT G,GTACT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3

VCF文件大致可分为三个部分:

  • 带有##的注释部分
  • 表头信息
  • INFO(变异详细信息)

注释部分

这一部分的信息做大致了解即可,通常包括VCF版本、参考基因组信息、运行命令以及介绍各个参数所代表的意义。此部分的格式分为四种:

1.key=value格式

1
2
3
4
5
##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial

介绍了VCF版本,日期,参考基因组等信息
另一种为如下形式

2.INFO

1
##INFO=<ID=ID,Number=number,Type=type,Description=”description”>

type通常为整数、浮点数、标识符、字符串

3.FILTERs

1
##FILTER=<ID=ID,Description=”description”>

4.FORMAT

1
##FORMAT=<ID=ID,Number=number,Type=type,Description=”description”>

表头信息

表头信息分为以下:

  • CHROM

    1
    2
    3

    - POS
    ```SNP或者INDEL所在染色体的位置,通常INDEL为第一个碱基所在染色体的位置
  • ID

    1
    2
    3

    - REF
    ```此位置碱基在参考基因组的序列,通常为A、T、C、G、N里面的一种
  • ALT

    1
    2
    3

    - QUAL
    ```变异碱基的可信度,该值越高,变异的可信度越高;Q=-10lgP,Q表示质量值;P表示这个位点发生错误的概率。如果想把错误率从控制在90%以上,P的阈值就是0.1,那lg(0.1)=-1,Q=(-10)*(-1)=10。同理,当Q=20时,错误率就控制在了0.01。
  • FILTER
    过滤情况,PASS说明质量较好,通过过滤;q10表示过滤质量低于10;s50表示低于总样本数的50%;'.'代表没有经过任何过滤

例子

20 14370 rs6054257 G A 29 PASS
表示这个变异位于20号染色体,14370这个碱基位置,在SNP数据库中对应ID为rs6054257,参考基因组碱基为G,变异碱基为A,质量值为29,通过过滤

20 17330 . T A 3 q10
表示这个变异位于20号染色体,17330这个碱基位置,在SNP数据库中无对应ID,参考基因组碱基为T,变异碱基为A,质量值为3,未通过过滤

20 1110696 rs6040355 A G,T 67 PASS
表示这个变异位于20号染色体,1110696这个碱基位置,在SNP数据库中对应ID为1110696,参考基因组碱基为A,变异碱基为A或T,质量值为67,通过过滤

变异详细信息

  • GT 表示这个样本的基因型,对于一个二倍体生物,GT值表示的是这个样本在这个位点所携带的两个等位基因。0表示跟REF一样;1表示表示跟ALT一样;2表示第二个ALT。
    20 14370 rs6054257 G A 29 PASS
    0/0–>G/G
    0/1–>G/A
    1/1–>A/A
  • AD 这两个值分别表示覆盖到REF和ALT碱基的reads数,相当于支持REF和支持ALT的测序深度。
  • DP 覆盖到这个位点的总的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。

  • AC 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。

  • 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。
  • MQ 表示覆盖序列质量的均方值RMS Mapping Quality
  • FQ phred值关于所有样本相似的可能性
  • AF1 AF(Allele Frequency) 表示Allele的频率,AF1为第一个ALT allele 发生频率的可能性评估
  • AC1 AC表示Allele(等位基因)的数目,AC1为对第一个ALT allele count的最大可能性评估
  • IS 插入缺失或部分插入缺失的reads允许的最大数量
  • G3 ML 评估基因型出现的频率
  • HWE chi^2基于HWE的测试p值和G3
  • CLR 在受到或者不受限制的情况下基因型出现可能性log值
  • UGT 最可能不受限制的三种基因型结构
  • CGT 最可能受限制三种基因型的结构
  • PV4 四种P值得误差,分别是(strand、baseQ、mapQ、tail distance bias)
  • INDEL 表示该位置的变异是插入缺失
  • PC2 非参考等位基因的phred(变异的可能性)值在两个分组中大小不同
  • PCHI2 后加权chi^2,根据p值来测试两组样本之间的联系
  • QCHI2 Phred scaled PCHI2.
  • PR 置换产生的一个较小的PCHI2
  • QBD Quality by Depth,测序深度对质量的影响
  • RPB 序列的误差位置(Read Position Bias)
  • MDV样本中高质量非参考序列的最大数目
  • VDBVariant Distance Bias,RNA序列中过滤人工拼接序列的变异误差范围

参考资料
http://www.internationalgenome.org/wiki/Analysis/vcf4.0/
http://samtools.github.io/hts-specs/VCFv4.2.pdf
http://samtools.github.io/bcftools/bcftools.html
http://www.cnblogs.com/emanlee/p/4562064.html
http://www.bio-info-trainee.com/863.html