6 min read

Introduction to VCF file

当前,育种芯片数据主要是通过VCF格式的文件提供,而且PLINK等质控软件,对于VCF文件的支持非常方便,因此了解VCF文件的格式和其中的各项参数,是非常必要的。

本文主要参考vcfR包的vignette:VCF data

1. 什么是VCF?

VCF (Variant Call Format)是储存变异信息的常用文件格式,通常为文本文件。VCF目前包括4.1-4.3三个正式版本,4.4为草稿版本。四个版本的具体格式,查看这里

2. VCF文件结构

vcfR包的作者把一个VCF文件分为三个区,分别命名为meta、fix和gt。参见这里的示意图。

2.1 meta region

位于文件顶端,主要是定义文件中用到的缩写信息。也存档了产生VCF文件的软件及其调用的参数信息。

看一下vcfR包中自带的vcf例子文件。可以看到,例子文件包括18个样本,一条染色体,共计有2533个SNP位点。

library(vcfR)
## Warning: 程辑包'vcfR'是用R版本4.1.1 来建造的
## 
##    *****       ***   vcfR   ***       *****
##    This is vcfR 1.12.0 
##      browseVignettes('vcfR') # Documentation
##      citation('vcfR') # Citation
##    *****       *****      *****       *****
data(vcfR_example)
vcf
## ***** Object of Class vcfR *****
## 18 samples
## 1 CHROMs
## 2,533 variants
## Object size: 3.2 Mb
## 8.497 percent missing data
## *****        *****         *****

接下来,查看一下meta区的前7行。strwrap()是自动换行函数。从输出行中,可以看出:

  • meta区中的每一行从两个##开始。
  • 第一行包括VCF文件格式版本。例子文件是4.1版本格式。该行是VCF文件必须要有的行
  • 第二行给出产生本VCF文件的软件,例子文件使用GATK产生。该行不是VCF文件必须要包括的行
  • VCF文件可能会针对每条染色体或supercontig或contig产生一行信息。
  • 剩余行定了在fix区中INFO列和在gt区中FORMAT列使用的各种缩写名称。
strwrap(vcf@meta[1:6])
## [1] "##fileformat=VCFv4.1"                                                   
## [2] "##source=\"GATK haplotype Caller, phased with beagle4\""                
## [3] "##FILTER=<ID=LowQual,Description=\"Low quality\">"                      
## [4] "##FORMAT=<ID=AD,Number=.,Type=Integer,Description=\"Allelic depths for" 
## [5] "the ref and alt alleles in the order listed\">"                         
## [6] "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Approximate read"   
## [7] "depth (reads with MQ=255 or with bad mates are filtered)\">"            
## [8] "##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype Quality\">"

meta区中的部分行可能会很长,不容易阅读。可以用vcfR包中的queryMETA()函数汇总meta部分的信息。但是并不是所有的信息,都在汇总结果中。譬如,“contig”信息就没有包括在其中。

queryMETA(vcf)
##  [1] "FILTER=ID=LowQual"                   
##  [2] "FORMAT=ID=AD"                        
##  [3] "FORMAT=ID=DP"                        
##  [4] "FORMAT=ID=GQ"                        
##  [5] "FORMAT=ID=GT"                        
##  [6] "FORMAT=ID=PL"                        
##  [7] "GATKCommandLine=ID=HaplotypeCaller"  
##  [8] "INFO=ID=AC"                          
##  [9] "INFO=ID=AF"                          
## [10] "INFO=ID=AN"                          
## [11] "INFO=ID=BaseQRankSum"                
## [12] "INFO=ID=ClippingRankSum"             
## [13] "INFO=ID=DP"                          
## [14] "INFO=ID=DS"                          
## [15] "INFO=ID=FS"                          
## [16] "INFO=ID=HaplotypeScore"              
## [17] "INFO=ID=InbreedingCoeff"             
## [18] "INFO=ID=MLEAC"                       
## [19] "INFO=ID=MLEAF"                       
## [20] "INFO=ID=MQ"                          
## [21] "INFO=ID=MQ0"                         
## [22] "INFO=ID=MQRankSum"                   
## [23] "INFO=ID=QD"                          
## [24] "INFO=ID=ReadPosRankSum"              
## [25] "INFO=ID=SOR"                         
## [26] "1 contig=<IDs omitted from queryMETA"

我们也可以通过设定参数element只查询某一个元素。譬如,我们查询DP (reads depth)-reads覆盖度。我们会在FORMAT和INFO中都查到DP的信息。

queryMETA(vcf, element = "DP")
## [[1]]
## [1] "FORMAT=ID=DP"                                                                         
## [2] "Number=1"                                                                             
## [3] "Type=Integer"                                                                         
## [4] "Description=Approximate read depth (reads with MQ=255 or with bad mates are filtered)"
## 
## [[2]]
## [1] "INFO=ID=DP"                                                           
## [2] "Number=1"                                                             
## [3] "Type=Integer"                                                         
## [4] "Description=Approximate read depth; some reads may have been filtered"

如果我们只想在FORMAT部分查询DP,可以进一步限定element参数的值为“FORMAT=<ID=DP”。

queryMETA(vcf, element = "FORMAT=<ID=DP")
## [[1]]
## [1] "FORMAT=ID=DP"                                                                         
## [2] "Number=1"                                                                             
## [3] "Type=Integer"                                                                         
## [4] "Description=Approximate read depth (reads with MQ=255 or with bad mates are filtered)"

queryMETA()中参数nice默认设置为true,返回整理后易于阅读的数据。譬如,将原始VCF文件中的尖括号删除,且分为四行。 如果想直接看VCF文件中的原始信息,可以将nice设置为FALSE。

queryMETA(vcf, element = "FORMAT=<ID=DP", nice = FALSE)
## [1] "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Approximate read depth (reads with MQ=255 or with bad mates are filtered)\">"

2.2 fix region

  • fix区包含了每个SNP位点的信息,有时会对所有样本进行汇总。
  • fix区的前八列,标题为CHROM, POS, ID, REF, ALT, QUAL, FILTER和INFO。这是每个变异位点的信息,是“固定的”,或者说对于所有的样本都是相同的。

通过getFIX()函数,查看实例数据vcf中fix区的信息,函数默认输出前7列,INFO列默认不输出,原因在下边讲。head()函数控制只输出前几行:

head(getFIX(vcf))
##      CHROM              POS   ID REF ALT QUAL     FILTER
## [1,] "Supercontig_1.50" "2"   NA "T" "A" "44.44"  NA    
## [2,] "Supercontig_1.50" "246" NA "C" "G" "144.21" NA    
## [3,] "Supercontig_1.50" "549" NA "A" "C" "68.49"  NA    
## [4,] "Supercontig_1.50" "668" NA "G" "C" "108.07" NA    
## [5,] "Supercontig_1.50" "765" NA "A" "C" "92.78"  NA    
## [6,] "Supercontig_1.50" "780" NA "G" "T" "58.38"  NA
  • CHROM:染色体名称;
  • POS:变异位点在染色体内的位置;
  • ID:变异位点名称,在例子中,该列信息缺失,用NA表示;
  • REF:参考基因组中对应的等位基因;
  • ALT:样本基因组中发生变异的等位基因,如果有多个变异等位基因,相互之间用英文逗号隔开;
  • QUAL:CALL出来的变异位点的质量。Phred格式(Phred_scaled)的质量值,表示在该位点存在变异的可能性;该值越高,则变异的可能性越大;计算方法:Phred值 = -10 * log10(1-p)。 1-p为call变异位点的错误概率; QUAl=10,表示错误概率为0.1,即正确率p=0.9;QUAL=20,表示错误概率为0.01(百分之一),即正确率p=0.99;QUAL=30,表示错误概率为0.001(千分之一),即正确率p=0.999,以此类推。
  • FILTER:该列信息缺失。但可以包含有关该位点是否通过某种形式的质量评估的信息。

第8列INFO包括由分号分割的信息。信息可能会相当长,且非常冗杂。INFO列中的信息所用到的缩写在meta区进行了定义。

我们看一下第一行中的INFO列的内容,包括的缩写参数非常多。

head(getFIX(vcf, getINFO=TRUE))[1,8]
##                                                                                                                                                                                              INFO 
## "AC=1;AF=0.042;AN=24;BaseQRankSum=1.332;ClippingRankSum=-0.097;DP=197;FS=2.536;InbreedingCoeff=-0.0879;MLEAC=1;MLEAF=0.042;MQ=21.70;MQ0=0;MQRankSum=0.148;QD=2.34;ReadPosRankSum=0.832;SOR=1.414"

如果我们对一些缩写,譬如AC、AF、AN含义不清楚,我们可以像前边写的那样,利用queryMETA()函数进行查询。其中,AC表示每一个变异等位基因在基因型中的数量,AF表示对应的频率,AN表示在call出来的基因型中的总数量。

queryMETA(vcf, element = "INFO=<ID=AC")
## [[1]]
## [1] "INFO=ID=AC"                           
## [2] "Number=A"                             
## [3] "Type=Integer"                         
## [4] "Description=Allele count in genotypes"
## [5] " for each ALT allele"                 
## [6] " in the same order as listed"
queryMETA(vcf, element = "INFO=<ID=AF")
## [[1]]
## [1] "INFO=ID=AF"                   "Number=A"                    
## [3] "Type=Float"                   "Description=Allele Frequency"
## [5] " for each ALT allele"         " in the same order as listed"
queryMETA(vcf, element = "INFO=<ID=AN")
## [[1]]
## [1] "INFO=ID=AN"                                             
## [2] "Number=1"                                               
## [3] "Type=Integer"                                           
## [4] "Description=Total number of alleles in called genotypes"

2.3 gt region

  • gt (genotype)区包括了每个样本每个变异位点的信息。每个变异位点和每个样本的值都用冒号分隔。每个基因型的多种类型的数据可以以这种方式存储。
  • 数据的格式由FORMAT列(即第9列)指定。这里我们看到了GT、AD、DP、GQ和PL的信息。如前所述,可以通过查询meta区来了解这些首字母缩略词的含义。
  • 每个变异位点不一定有相同的信息(例如,SNPs和indels可能有不同的处理方法),所以最好独立处理每一行。在gt区不同的变异位点可能包含不同的信息。

查看前6个变异位点,三个样本的分型信息。

vcf@gt[1:6,1:4]
##      FORMAT           BL2009P4_us23               DDR7602                  
## [1,] "GT:AD:DP:GQ:PL" "0|0:62,0:62:99:0,190,2835" "0|0:12,0:12:39:0,39,585"
## [2,] "GT:AD:DP:GQ:PL" "1|0:5,5:10:99:111,0,114"   NA                       
## [3,] "GT:AD:DP:GQ:PL" NA                          NA                       
## [4,] "GT:AD:DP:GQ:PL" "0|0:1,0:1:3:0,3,44"        NA                       
## [5,] "GT:AD:DP:GQ:PL" "0|0:2,0:2:6:0,6,49"        "0|0:1,0:1:3:0,3,34"     
## [6,] "GT:AD:DP:GQ:PL" "0|0:2,0:2:6:0,6,49"        "0|0:1,0:1:3:0,3,34"     
##      IN2009T1_us22              
## [1,] "0|0:37,0:37:99:0,114,1709"
## [2,] "0|1:2,1:3:16:16,0,48"     
## [3,] "0|0:2,0:2:6:0,6,51"       
## [4,] "1|1:0,1:1:3:25,3,0"       
## [5,] "0|0:1,0:1:3:0,3,31"       
## [6,] "0|0:3,0:3:9:0,9,85"

如果对GT、AD、DP、GQ和PL的含义不清楚,可以通过queryMETA()函数查询。

queryMETA(vcf,element = "FORMAT=<ID=GT")
## [[1]]
## [1] "FORMAT=ID=GT"         "Number=1"             "Type=String"         
## [4] "Description=Genotype"
queryMETA(vcf,element = "FORMAT=<ID=AD")
## [[1]]
## [1] "FORMAT=ID=AD"                                                              
## [2] "Number=."                                                                  
## [3] "Type=Integer"                                                              
## [4] "Description=Allelic depths for the ref and alt alleles in the order listed"
queryMETA(vcf,element = "FORMAT=<ID=DP")
## [[1]]
## [1] "FORMAT=ID=DP"                                                                         
## [2] "Number=1"                                                                             
## [3] "Type=Integer"                                                                         
## [4] "Description=Approximate read depth (reads with MQ=255 or with bad mates are filtered)"
queryMETA(vcf,element = "FORMAT=<ID=GQ")
## [[1]]
## [1] "FORMAT=ID=GQ"                 "Number=1"                    
## [3] "Type=Integer"                 "Description=Genotype Quality"
queryMETA(vcf,element = "FORMAT=<ID=PL")
## [[1]]
## [1] "FORMAT=ID=PL"                                                               
## [2] "Number=G"                                                                   
## [3] "Type=Integer"                                                               
## [4] "Description=Normalized"                                                     
## [5] " Phred-scaled likelihoods for genotypes as defined in the VCF specification"

试着解释一下BL2009P4_us23样品的第一个变异位点信息“0|0:62,0:62:99:0,190,2835”:

  • 0|0 对应GT,即基因型为纯合的ref即参考等位基因。两个数字中间用’/’或者’|’分开,后者表示基因型已经phased,已经明确等位基因属于父本还是母本。这两个数字表示双倍体样本的基因型。0表示参考等位基因; 1表示第一个变异等位基因; 2表示有第二个变异等位基因。因此: 0/0 表示该样本该位点基因型为纯合,均为参考等位基因; 0/1 表示基因型为杂合,有参考等位基因和第1个变异等位基因; 1/1 表示基因型为纯合,均为第1个变异等位基因。
  • 62,0 对应AD,即参考等位基因和变异等位基因的覆盖度。
  • 62 对应DP,即近似的测序覆盖度。
  • 99 对应GQ,即基因型的质量,类似于QUAL的计算公式,表示该基因型的错误率\(<10^{-9}\)
  • 0,190,2835 对应PL,即在VCF中定义的基因型的Phred尺度的似然率。具体而言,表示指定的三种基因型(0/0,0/1,1/1)的质量值。值越大,表明为该种基因型的可能性越小。Phred值= -10 * log (p),p为基因型存在的概率。

3. 利用vcfR包函数操作VCF文件

3.1 简单汇总

通过head()函数查看三个区的信息。

head(vcf)
## [1] "***** Object of class 'vcfR' *****"
## [1] "***** Meta section *****"
## [1] "##fileformat=VCFv4.1"
## [1] "##source=\"GATK haplotype Caller, phased with beagle4\""
## [1] "##FILTER=<ID=LowQual,Description=\"Low quality\">"
## [1] "##FORMAT=<ID=AD,Number=.,Type=Integer,Description=\"Allelic depths fo [Truncated]"
## [1] "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Approximate read  [Truncated]"
## [1] "##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype Quality\">"
## [1] "First 6 rows."
## [1] 
## [1] "***** Fixed section *****"
##      CHROM              POS   ID REF ALT QUAL     FILTER
## [1,] "Supercontig_1.50" "2"   NA "T" "A" "44.44"  NA    
## [2,] "Supercontig_1.50" "246" NA "C" "G" "144.21" NA    
## [3,] "Supercontig_1.50" "549" NA "A" "C" "68.49"  NA    
## [4,] "Supercontig_1.50" "668" NA "G" "C" "108.07" NA    
## [5,] "Supercontig_1.50" "765" NA "A" "C" "92.78"  NA    
## [6,] "Supercontig_1.50" "780" NA "G" "T" "58.38"  NA    
## [1] 
## [1] "***** Genotype section *****"
##      FORMAT           BL2009P4_us23               DDR7602                  
## [1,] "GT:AD:DP:GQ:PL" "0|0:62,0:62:99:0,190,2835" "0|0:12,0:12:39:0,39,585"
## [2,] "GT:AD:DP:GQ:PL" "1|0:5,5:10:99:111,0,114"   NA                       
## [3,] "GT:AD:DP:GQ:PL" NA                          NA                       
## [4,] "GT:AD:DP:GQ:PL" "0|0:1,0:1:3:0,3,44"        NA                       
## [5,] "GT:AD:DP:GQ:PL" "0|0:2,0:2:6:0,6,49"        "0|0:1,0:1:3:0,3,34"     
## [6,] "GT:AD:DP:GQ:PL" "0|0:2,0:2:6:0,6,49"        "0|0:1,0:1:3:0,3,34"     
##      IN2009T1_us22               LBUS5                     NL07434             
## [1,] "0|0:37,0:37:99:0,114,1709" "0|0:12,0:12:39:0,39,585" NA                  
## [2,] "0|1:2,1:3:16:16,0,48"      NA                        NA                  
## [3,] "0|0:2,0:2:6:0,6,51"        NA                        NA                  
## [4,] "1|1:0,1:1:3:25,3,0"        NA                        "0|0:1,0:1:3:0,3,28"
## [5,] "0|0:1,0:1:3:0,3,31"        "0|0:1,0:1:3:0,3,34"      "0|0:1,0:1:3:0,3,26"
## [6,] "0|0:3,0:3:9:0,9,85"        "0|0:1,0:1:3:0,3,34"      NA                  
## [1] "First 6 columns only."
## [1] 
## [1] "Unique GT formats:"
## [1] "GT:AD:DP:GQ:PL"
## [1]

3.2 快速全览VCF文件

查看示例文件中包括的染色体名称。可以看到,仅有一条染色体,名字为Supercontig_1.50。

unique(vcf@fix[,1])
## [1] "Supercontig_1.50"

创建chrom对象。

chrom <- create.chromR(name='Supercontig_1.50', vcf=vcf)
## Names in vcf:
##   Supercontig_1.50
## Initializing var.info slot.
## var.info slot initialized.

快速的看一下VCF文件中的信息。

plot(chrom)

  • reads覆盖度(DP)的分布非常突出。据推测,每个基因组的大部分区域都是在某个基础倍性水平进行测序的。在这里我们看到了一个峰值,这可能代表了那个基出倍性区域,但我们也看到了一个长尾,这可能代表了拷贝数变异。基因型通常期望一个恒定的倍性水平,因此在拷贝数变异区域中call 变异位点可信度不高。
  • 比对质量(MQ)都在60左右达到峰值。因此,如果我们想对这个参数进行过滤,控制在60左右较好。关于MQ的计算和解释,参见这里,讲的不错。
  • 对变异位点质量的解释(QUAL)似乎不那么直接,直方图显示,大量的位点堆积在0附近。因此,这可能不是一个理想的质控筛选参数。

3.3 筛选变异位点

  • 使用masker()函数来过滤掉我们不太相信的数据。该函数的作用是:使用QUAL质量、DP深度和MQ比对质量等参数来选择高质量的变异位点。
  • 使用masker()函数时,被认为是低质量的变异位点并不会从数据集中删除。相反,将创建一个逻辑向量来指示哪些变异位点已被过滤或未被过滤。这是为了在整个分析过程中维护数据矩阵的结构,并允许用户可以轻松地撤消任何更改。

代码中的参数,删除了位于拷贝数变异区域的位点,将比对质量控制在60左右。

chrom <- masker(chrom, min_QUAL = 1, min_DP = 300, max_DP = 700, min_MQ = 59.9,  max_MQ = 60.1)
plot(chrom)

3.4 处理染色体对象

使用proc.chromR()来处理chromR对象。这个函数调用几个辅助函数来处理变异位点、序列和注释数据以实现可视化。

注意右下角第四个图,表示每一个窗口的变异位点数量。

chrom <- proc.chromR(chrom, verbose=TRUE)
## Warning in proc.chromR(chrom, verbose = TRUE): seq slot is NULL.
## Warning in proc.chromR(chrom, verbose = TRUE): annotation slot has no rows.
## Warning in proc.chromR(chrom, verbose = TRUE): seq slot is NULL, chromosome
## representation not made (seq2rects).
## Warning in proc.chromR(chrom, verbose = TRUE): seq slot is NULL, chromosome
## representation not made (seq2rects, chars=n).
## Population summary complete.
##   elapsed time:  0.03
## window_init complete.
##   elapsed time:  0.02
## Warning in proc.chromR(chrom, verbose = TRUE): seq slot is NULL, windowize_fasta
## not run.
## Warning in proc.chromR(chrom, verbose = TRUE): ann slot has zero rows.
## windowize_variants complete.
##   elapsed time:  0
plot(chrom)

3.5 可视化数据

chromoqc()使用R函数layout()来制作数据的复合图。这些图包括有染色体坐标的散点图和barplot。

chromoqc(chrom,dp.alpha = 20)