VCF data
Brian J. Knaus
2023-12-07
Most variant calling pipelines result in files containing variant
information. The variant call format (VCF) is a
popular format for this data. Variant callers typically attempt to
agressively call variants with the perspective that a downstream quality
control step will remove low quality variants. A first step in working
with this data is to understand their contents.
Three sections
A VCF file can be thought of as having three sections: a meta
region, a fix region and a gt
region. The meta region is at the top of the file. The
information in the meta region defines the abbreviations used elsewhere
in the file. It may also document software used to create the file as
well as parameters used by this software. Below the meta region, the
data are tabular. The first eight columns of this table contain
information about each variant. This data may be common over all
variants, such as its chromosomal position, or a summary over all
samples, such as quality metrics. These data are fixed, or the same,
over all samples. The fix region is required in a VCF file, subsequent
columns are optional but are common in my experience. Beginning at
column ten is a column for every sample. The values in these columns are
information for each sample and each variant. The organization of each
cell containing a genotype and associated information is specified in
column nine. The location of these three regions within a file can be
represented by the cartoon below.
The VCF file specification is flexible. This means that there are
slots for certain types of data, but any particular software which
creates a VCF file does not necessarily use them all. Similarly, authors
have the opportunity to include new forms of data, forms which may not
have been foreseen by the authors of the VCF specification. The result
is that all VCF files do not contain the same information.
For this example, we’ll use example data provided with vcfR.
library(vcfR)
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
## ***** ***** *****
The function library() loads libraries, in this case the
package vcfR. The function data() loads datasets that were
included with R and its packages. Our usage of data() loads
the objects ‘gff’, ‘dna’ and ‘vcf’ from the ‘vcfR_example’ dataset. Here
we’re only interested in the object ‘vcf’ which contains example VCF
data. When we call the object name with no function it invokes the
‘show’ method which prints some summary information to the console.
The meta region
The meta region contains information about the file, its creation, as
well as information to interpret abbreviations used elsewhere in the
file. Each line of the meta region begins with a double pound sign
(‘##’). The example which comes with vcfR is shown below. (Only the
first 10 lines are shown for brevity.)
## [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\">"
## [9] "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">"
The first line contains the version of the VCF format used in the
file. This line is required. The second line specifies the software
which created the VCF file. This is not required, so not all VCF files
include it. When they do, the file becomes self documenting. Note that
the alignment software is not included here because it was used upstream
of the VCF file’s creation (aligners typically create *.SAM or *.BAM
format files). Because the file can only include information about the
software that created it, the entire pipeline does not get documented.
Some VCF files may contain a line for every chromosome (or supercontig
or contig depending on your genome), so they may become rather long.
Here, the remaining lines contain INFO and FORMAT specifications which
define abbreviations used in the fix and gt portions of the file.
The meta region may include long lines that may not be easy to view.
In vcfR we’ve created a function to help prcess this data.
## [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"
When the function queryMETA() is called with only a vcfR
object as a parameter, it attempts to summarize the meta information.
Not all of the information is returned. For example, ‘contig’ elements
are not returned. This is an attempt to summarize information that may
be most useful for comprehension of the file’s contents.
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"
When an element parameter is included, only the information about
that element is returned. In this example the element ‘DP’ is returned.
We see that this acronym is defined as both a ‘FORMAT’ and ‘INFO’
acronym. We can narrow down our query by including more information in
the element parameter.
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)"
Here we’ve isolated the definition of ‘DP’ as a ‘FORMAT’ element.
Note that the function queryMETA() includes the parameter
nice which by default is TRUE and attempts to present the
data in a nicely formatted manner. However, our query is performed on
the actual information in the ‘meta’ region. It is therefore sometimes
appropriate to set nice = FALSE so that we can see the raw
data. In the above example the angled bracket (‘<’) is omitted from
the nice = TRUE representation but is essential to
distinguishing the ‘FORMAT’ element from the ‘INFO’ element.
The fix region
The fix region contains information for each variant which is
sometimes summarized over all samples. The first eight columns of the
fixed region and are titled CHROM, POS, ID, REF, ALT, QUAL, FILTER and
INFO. This is per variant information which is ‘fixed’, or the same,
over all samples. The first two columns indicate the location of the
variant by chromosome and position within that chromosome. Here, the ID
field has not been used, so it consists of missing data (NA). The REF
and ALT columns indicate the reference and alternate allelic states.
When multiple alternate allelic states are present they are delimited
with commas. The QUAL column attempts to summarize the quality of each
variant over all samples. The FILTER field is not used here but could
contain information on whether a variant has passed some form of quality
assessment.
## 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
The eigth column, titled INFO, is a semicolon delimited list of
information. It can be rather long and cumbersome. The function
getFIX() will suppress this column by default. Each
abbreviation in the INFO column should be defined in the meta section.
We can validate this by querying the meta portion, as we did in the
‘meta’ section above.
The gt region
The gt (genotype) region contains information about each variant for
each sample. The values for each variant and each sample are colon
delimited. Multiple types of data for each genotype may be stored in
this manner. The format of the data is specified by the FORMAT column
(column nine). Here we see that we have information for GT, AD, DP, GQ
and PL. The definition of these acronyms can be referenced by querying
the the meta region, as demonstrated previously. Every variant does not
necessarily have the same information (e.g., SNPs and indels may be
handled differently), so the rows are best treated independently.
Different variant callers may include different information in this
region.
## 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"
vcfR
Using the R package vcfR, we can read VCF format files into memory
using the function read.vcfR(). Once in memory we can use
the head() method to summarize the information in the three
VCF regions.
## [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]
We now have a summary of our VCF file which we can use to help
understand what forms of information are contained within it. This
information can be further explored with plotting functions and used to
filter the VCF file for high quality variants.
