当前，育种芯片数据主要是通过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: 0plot(chrom)
3.5 可视化数据
chromoqc()使用R函数layout()来制作数据的复合图。这些图包括有染色体坐标的散点图和barplot。
chromoqc(chrom,dp.alpha = 20)