VCF
objects. A Shiny web-application, the Shiny Variant Explorer (tSVE), provides a convenient interface to demonstrate those functionalities integrated in a programming-free environment.
TVTB 1.4.0
The VCF Tool Box (TVTB) offers S4 classes and methods to filter, summarise and visualise genetic variation data stored in VCF files pre-processed by the Ensembl Variant Effect Predictor (VEP) (McLaren et al. 2010). An RStudio/Shiny web-application, the Shiny Variant Explorer (tSVE), provides a convenient interface to demonstrate those functionalities integrated in a programming-free environment.
Currently, major functionalities in the TVTB package include:
A class to store recurrent parameters of genetic analyses
Genotype counts and allele frequencies
ExpandedVCF
objects (i.e. bi-allelic records)Classes of VCF filter rules
fixed
slot of an VCF
objectinfo
slot of an VCF
objectVCF
objects using the above filter rulesThe VCF Tool Box can be installed using the following code:
source("http://bioconductor.org/biocLite.R")
biocLite("TVTB")
Once installed, the package can be loaded and attached as follows:
library(TVTB)
Most functionalities in TVTB require recurrent information such as:
<phenotype>_<level>_<suffix>
,<suffix>
.To reduce the burden of repetition during programming, and to facilitate analyses using consistent sets of parameters, TVTB implements the TVTBparam
class. The TVTBparam
class offer a container for parameters recurrently used across the package. A TVTBparam
object may be initialised as follows:
tparam <- TVTBparam(Genotypes(
ref = "0|0",
het = c("0|1", "1|0", "0|2", "2|0", "1|2", "2|1"),
alt = c("1|1", "2|2")),
ranges = GenomicRanges::GRangesList(
SLC24A5 = GenomicRanges::GRanges(
seqnames = "15",
IRanges::IRanges(
start = 48413170, end = 48434757)
)
)
)
TVTBparam
objects have a convenient summary view and accessor methods:
tparam
## class: TVTBparam
## @genos: class: Genotypes
## @ref (hom. ref.): "REF" {0|0}
## @het (heter.): "HET" {0|1, 1|0, 0|2, 2|0, 1|2, 2|1}
## @alt (hom. alt.): "ALT" {1|1, 2|2}
## @ranges: 1 GRanges on 1 sequence(s)
## @aaf (alt. allele freq.): "AAF"
## @maf (minor allele freq.): "MAF"
## @vep (Ensembl VEP key): "CSQ"
## @svp: <ScanVcfParam object>
## @bp: <SerialParam object>
In this example:
genos(x)
Genotypes
"0|0"
."REF"
."0|1"
, "1|0"
, "0|2"
, "2|0"
, "1|2"
, and "2|1"
."HET"
."1|1"
."ALT"
.ranges(x)
GRangesList
"15"
.aaf(x)
"AAF"
.maf(x)
"MAF"
.vep(x)
"CSQ"
.bp(x)
BiocParallelParam
svp(x)
ScanVcfParam
which
slot automatically populated with reduce(unlist(ranges(x)))
Default values are provided for all slots except genotypes, as these may vary more frequently from one data set to another (e.g. phased, unphased, imputed).
Functionalities in TVTB support CollapsedVCF
and ExpandedVCF
objects (both extending the virtual class VCF
) of the VariantAnnotation package.
Typically, CollapsedVCF
objects are produced by the VariantAnnotation readVcf
method after parsing a VCF file, and ExpandedVCF
objects result of the VariantAnnotation expand
method applied to a CollapsedVCF
object.
Any information that users deem relevant for the analysis may be imported from VCF files and stored in VCF
objects passed to TVTB methods. However, to enable the key functionalities of the package, the slots of a VCF
object should include at least the following information:
fixed(x)
"REF"
and "ALT"
.info(x)
<vep>
: where <vep>
stands for the INFO key where the Ensembl VEP predictions are stored in the VCF
object.geno(x)
GT
: genotypes.colData(x)
: phenotypes.In the near future, TVTB functionalities are expected to produce summary statistics and plots faceted by meta-features, each potentially composed of multiple genomic ranges.
For instance, burden tests may be performed on a set of transcripts, considering only variants in their respective sets of exons. The GenomicRanges GRangesList
class is an ideal container in example, as each GRanges
in the GRangesList
would represent a transcript, and each element in the GRanges
would represent an exon.
Furthermore, TVTBparam
objects may be supplied as the param
argument of the VariantAnnotation readVcf
. In this case, the TVTBparam
object is used to import only variants overlapping the relevant genomic regions. Moreover, the readVcf
method also ensured that the vep
slot of the TVTBparam
object is present in the header of the VCF file.
svp <- as(tparam, "ScanVcfParam")
svp
## class: ScanVcfParam
## vcfWhich: 1 elements
## vcfFixed: character() [All]
## vcfInfo:
## vcfGeno:
## vcfSamples:
Although VCF
objects may be constructed without attached phenotype data, phenotype information is critical to interpret and compare genetic variants between groups of samples (e.g. burden of damaging variants in different phenotype levels).
VCF
objects accept phenotype information (as S4Vectors DataFrame
) in the colData
slot. This practice has the key advantage of keeping phenotype and genetic information synchronised through operation such as subsetting and re-ordering, limiting workspace entropy and confusion.
An ExpandedVCF
object that contains the minimal data necessary for the rest of the vignette can be created as follows:
Step 1: Import phenotypes
phenoFile <- system.file(
"extdata", "integrated_samples.txt", package = "TVTB")
phenotypes <- S4Vectors::DataFrame(
read.table(file = phenoFile, header = TRUE, row.names = 1))
Step 2: Define the VCF file to parse
vcfFile <- system.file(
"extdata", "chr15.phase3_integrated.vcf.gz", package = "TVTB")
tabixVcf <- Rsamtools::TabixFile(file = vcfFile)
Step 3: Define VCF import parameters
VariantAnnotation::vcfInfo(svp(tparam)) <- vep(tparam)
VariantAnnotation::vcfGeno(svp(tparam)) <- "GT"
svp(tparam)
## class: ScanVcfParam
## vcfWhich: 1 elements
## vcfFixed: character() [All]
## vcfInfo: CSQ
## vcfGeno: GT
## vcfSamples:
Of particular interest in the above chunk of code:
TVTBparam
constructor previously populated the which
slot of svp
with “reduced” (i.e. non-overlapping) genomic ranges defined in the ranges
slot.vep
slot will be importedStep 4: Import and pre-process variants
# Import variants as a CollapsedVCF object
vcf <- VariantAnnotation::readVcf(
tabixVcf, param = tparam, colData = phenotypes)
# Expand into a ExpandedVCF object (bi-allelic records)
vcf <- VariantAnnotation::expand(x = vcf, row.names = TRUE)
Of particular interest in the above chunk of code, the readVcf
method is given:
TVTBparam
parameters, invoking the corresponding method signaturerownames
of those phenotypes defines the sample identifiers that are queried from the VCF file.colData
slot of the resulting VCF
object.The result is an ExpandedVCF
object that includes variants in the targeted genomic range(s) and samples:
## class: ExpandedVCF
## dim: 481 2504
## rowRanges(vcf):
## GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
## info(vcf):
## DataFrame with 1 column: CSQ
## info(header(vcf)):
## Number Type Description
## CSQ . String Consequence annotations from Ensembl VEP. Format: Al...
## geno(vcf):
## SimpleList of length 1: GT
## geno(header(vcf)):
## Number Type Description
## GT 1 String Genotype
Although interesting figures and summary tables may be obtained as soon as the first ExpandedVCF
object is created (see section Summarising Ensembl VEP predictions), those methods may benefit from information added to additional INFO keys after data import, either manually by the user, or through various methods implemented in the TVTB package.
For instance, the method addOverallFrequencies
uses the reference homozoygote (REF), heterozygote (HET), and homozygote alternate (ALT) genotypes defined in the TVTBparam
object stored in the VCF
metadata to obtain the count of each genotype in an ExpandedVCF
object. Immediately thereafter, the method uses those counts to calculate alternate allele frequency (AAF) and minor allele frequency (MAF). Finally, the method stores the five calculated values (REF, HET, ALT, AAF, and MAF) in INFO keys defined by suffixes also declared in the TVTBparam
object.
initialInfo <- colnames(info(vcf))
vcf <- addOverallFrequencies(vcf = vcf)
setdiff(colnames(info(vcf)), initialInfo)
## [1] "REF" "HET" "ALT" "AAF" "MAF"
Notably, the addOverallFrequencies
method is synonym to the addFrequencies
method missing the argument phenos
:
vcf <- addFrequencies(vcf = vcf, force = TRUE)
Similarly, the method addPhenoLevelFrequencies
obtains the count of each genotype in samples associated with given level(s) of given phenotype(s), and stores the calculated values in INFO keys defined as <pheno>_<level>_<suffix>
, with suffixes defined in the TVTBparam
object stored in the VCF
metadata.
initialInfo <- colnames(info(vcf))
vcf <- addPhenoLevelFrequencies(
vcf = vcf, pheno = "super_pop", level = "AFR")
setdiff(colnames(info(vcf)), initialInfo)
## [1] "super_pop_AFR_REF" "super_pop_AFR_HET" "super_pop_AFR_ALT"
## [4] "super_pop_AFR_AAF" "super_pop_AFR_MAF"
Notably, the addPhenoLevelFrequencies
method is synonym to the addFrequencies
method called with the argument phenos
given as a list where names
are phenotypes, and values are character
vectors of levels to process within each phenotype:
initialInfo <- colnames(info(vcf))
vcf <- addFrequencies(
vcf,
list(super_pop = c("EUR", "SAS", "EAS", "AMR"))
)
setdiff(colnames(info(vcf)), initialInfo)
## [1] "super_pop_EUR_REF" "super_pop_EUR_HET" "super_pop_EUR_ALT"
## [4] "super_pop_EUR_AAF" "super_pop_EUR_MAF" "super_pop_SAS_REF"
## [7] "super_pop_SAS_HET" "super_pop_SAS_ALT" "super_pop_SAS_AAF"
## [10] "super_pop_SAS_MAF" "super_pop_EAS_REF" "super_pop_EAS_HET"
## [13] "super_pop_EAS_ALT" "super_pop_EAS_AAF" "super_pop_EAS_MAF"
## [16] "super_pop_AMR_REF" "super_pop_AMR_HET" "super_pop_AMR_ALT"
## [19] "super_pop_AMR_AAF" "super_pop_AMR_MAF"
In addition, the addFrequencies
method can be given a character
vector of phenotypes as the phenos
argument, in which case frequencies are calculated for all levels of the given phenotypes:
vcf <- addFrequencies(vcf, "pop")
head(grep("^pop_[[:alpha:]]+_REF", colnames(info(vcf)), value = TRUE))
## [1] "pop_GBR_REF" "pop_FIN_REF" "pop_CHS_REF" "pop_PUR_REF" "pop_CDX_REF"
## [6] "pop_CLM_REF"
Although VCF
objects are straightforward to subset using either indices and row names (as they inherit from the SummarizedExperiment RangedSummarizedExperiment
class), users may wish to identify variants that pass combinations of criteria based on information in their fixed
slot, info
slot, and Ensembl VEP predictions, a non-trivial task due to those pieces of information being stored in different slots of the VCF
object, and the 1:N relationship between variants and EnsemblVEP predictions.
To facilitate the definition of VCF filter rules, and their application to VCF
objects, TVTB extends the S4Vectors FilterRules
class in four new classes of filter rules:
Class | Motivation |
---|---|
VcfFixedRules |
Filter rules applied to the fixed slot of a VCF object. |
VcfInfoRules |
Filter rules applied to the info slot of a VCF object. |
VcfVepRules |
Filter rules applied to the Ensembl VEP predictions stored in a given INFO key of a VCF object. |
VcfFilterRules |
Combination of VcfFixedRules , VcfInfoRules , and VcfVepRules applicable to a VCF object. |
Note that FilterRules
objects themselves are applicable to VCF
objects, with two important difference from the above specialised classes:
VCF
slotsVCF
slots, for instance:S4Vectors::FilterRules(list(
mixed = function(x){
VariantAnnotation::fixed(x)[,"FILTER"] == "PASS" &
VariantAnnotation::info(x)[,"MAF"] >= 0.05
}
))
## FilterRules of length 1
## names(1): mixed
Instances of those classes may be initialised as follows:
VcfFixedRules
fixedR <- VcfFixedRules(list(
pass = expression(FILTER == "PASS"),
qual = expression(QUAL > 20)
))
fixedR
## VcfFixedRules of length 2
## names(2): pass qual
VcfInfoRules
infoR <- VcfInfoRules(
exprs = list(
rare = expression(MAF < 0.01 & MAF > 0),
common = expression(MAF > 0.05),
mac_ge3 = expression(HET + 2*ALT >= 3)),
active = c(TRUE, TRUE, FALSE)
)
infoR
## VcfInfoRules of length 3
## names(3): rare common mac_ge3
The above code chunk illustrates useful features of FilterRules
:
FilterRules
are initialised in an active state by default (evaluating an inactive rule returns TRUE
for all items) The active
argument may be used to initialise specific filter rules in an inactive state.expression
(or function
) may refer to multiple columns of the relevant slot in the VCF
object.VCF
object.VcfVepRules
vepR <- VcfVepRules(exprs = list(
missense = expression(Consequence %in% c("missense_variant")),
CADD_gt15 = expression(CADD_PHRED > 15)
))
vepR
## VcfVepRules of length 2
## names(2): missense CADD_gt15
VcfFilterRules
VcfFilterRules
combine VCF filter rules of different types in a single object.
vcfRules <- VcfFilterRules(fixedR, infoR, vepR)
vcfRules
## VcfFilterRules of length 7
## names(7): pass qual rare common mac_ge3 missense CADD_gt15
This vignette offers only a brief peek into the utility and flexibility of VCF filter rules. More (complex) examples are given in a separate vignette, including filter rules using functions and pattern matching. The documentation of the S4Vectors package—where the parent class FilterRules
is defined—can also be a source of inspiration.
As the above classes of VCF filter rules inherit from the S4Vectors FilterRules
class, they also benefit from its accessors and methods. For instance, VCF filter rules can easily be toggled between active and inactive states:
active(vcfRules)["CADD_gt15"] <- FALSE
active(vcfRules)
## pass qual rare common mac_ge3 missense CADD_gt15
## TRUE TRUE TRUE TRUE FALSE TRUE FALSE
A separate vignette describes in greater detail the use of classes that contain VCF filter rules.
Once defined, the above filter rules can be applied to ExpandedVCF
objects, in the same way as FilterRules
are evaluated in a given environment (see the S4Vectors documentation):
summary(eval(expr = infoR, envir = vcf))
## Mode FALSE
## logical 481
summary(eval(expr = vcfRules, envir = vcf))
## Mode FALSE
## logical 481
summary(evalSeparately(expr = vcfRules, envir = vcf))
## pass qual rare common
## Mode:logical Mode:logical Mode :logical Mode :logical
## TRUE:481 TRUE:481 FALSE:45 FALSE:453
## TRUE :436 TRUE :28
## mac_ge3 missense CADD_gt15
## Mode:logical Mode :logical Mode:logical
## TRUE:481 FALSE:454 TRUE:481
## TRUE :27
Let us show the alternate allele frequency (AAF) of common variants, estimated in each super-population, in the context of the transcripts ovelapping the region of interest.
In the MAF
track:
plotInfo(
subsetByFilter(vcf, vcfRules["common"]), "AAF",
range(GenomicRanges::granges(vcf)),
EnsDb.Hsapiens.v75::EnsDb.Hsapiens.v75,
"super_pop",
zero.rm = FALSE
)
Alternatively, the minor allele frequency (MAF) of missense variants (as estimated from the entire data set) may be visualised in the same manner. However, due to the nature of those variants, the zero.rm
argument may be set to TRUE
to hide all data points showing a MAF of 0
; thereby variants actually detected in each super-population are emphasised even at low frequencies.
plotInfo(
subsetByFilter(vcf, vcfRules["missense"]), "MAF",
range(GenomicRanges::granges(vcf)),
EnsDb.Hsapiens.v75::EnsDb.Hsapiens.v75,
"super_pop",
zero.rm = TRUE
)
Using the GGally ggpairs
method, let us make a matrix of plots for common variants, showing:
pairsInfo(subsetByFilter(vcf, vcfRules["common"]), "AAF", "super_pop")
Note that the ellipsis ...
allows a high degree of customisation, as it passes additional arguments to the underlying ggpairs
method.
This section presents upcoming features.
As soon as genetic and phenotypic information are imported into an ExpandedVCF
object, or after the object was extended with additional information, the scientific value of the data may be revealed by a variety of summary statistics and graphical representations. This section will soon present several ideas being implemented in TVTB, for instance:
Dr. Stefan Gräf and Mr. Matthias Haimel for advice on the VCF file format and the Ensembl VEP script. Prof. Martin Wilkins for his trust and support. Dr. Michael Lawrence for his helpful code review and suggestions.
Last but not least, the amazing collaborative effort of the rep("many",n)
Bioconductor developers whose hard work appears through the dependencies of this package.
Here is the output of sessionInfo()
on the system on which this document was compiled:
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] TVTB_1.4.0 knitr_1.17 BiocStyle_2.6.0
##
## loaded via a namespace (and not attached):
## [1] ProtGenerics_1.10.0 bitops_1.0-6
## [3] matrixStats_0.52.2 EnsDb.Hsapiens.v75_2.99.0
## [5] bit64_0.9-7 RColorBrewer_1.1-2
## [7] progress_1.1.2 httr_1.3.1
## [9] rprojroot_1.2 GenomeInfoDb_1.14.0
## [11] tools_3.4.2 backports_1.1.1
## [13] R6_2.2.2 rpart_4.1-11
## [15] Hmisc_4.0-3 DBI_0.7
## [17] lazyeval_0.2.1 BiocGenerics_0.24.0
## [19] Gviz_1.22.0 colorspace_1.3-2
## [21] nnet_7.3-12 gridExtra_2.3
## [23] prettyunits_1.0.2 GGally_1.3.2
## [25] RMySQL_0.10.13 curl_3.0
## [27] bit_1.1-12 compiler_3.4.2
## [29] Biobase_2.38.0 htmlTable_1.9
## [31] DelayedArray_0.4.0 labeling_0.3
## [33] rtracklayer_1.38.0 bookdown_0.5
## [35] scales_0.5.0 checkmate_1.8.5
## [37] stringr_1.2.0 digest_0.6.12
## [39] Rsamtools_1.30.0 foreign_0.8-69
## [41] rmarkdown_1.6 XVector_0.18.0
## [43] pkgconfig_2.0.1 base64enc_0.1-3
## [45] dichromat_2.0-0 htmltools_0.3.6
## [47] ensembldb_2.2.0 limma_3.34.0
## [49] BSgenome_1.46.0 htmlwidgets_0.9
## [51] rlang_0.1.2 RSQLite_2.0
## [53] BiocInstaller_1.28.0 shiny_1.0.5
## [55] BiocParallel_1.12.0 acepack_1.4.1
## [57] VariantAnnotation_1.24.0 RCurl_1.95-4.8
## [59] magrittr_1.5 GenomeInfoDbData_0.99.1
## [61] Formula_1.2-2 Matrix_1.2-11
## [63] Rcpp_0.12.13 munsell_0.4.3
## [65] S4Vectors_0.16.0 ensemblVEP_1.20.0
## [67] stringi_1.1.5 yaml_2.1.14
## [69] SummarizedExperiment_1.8.0 zlibbioc_1.24.0
## [71] plyr_1.8.4 AnnotationHub_2.10.0
## [73] grid_3.4.2 blob_1.1.0
## [75] parallel_3.4.2 lattice_0.20-35
## [77] Biostrings_2.46.0 splines_3.4.2
## [79] GenomicFeatures_1.30.0 pander_0.6.1
## [81] GenomicRanges_1.30.0 reshape2_1.4.2
## [83] biomaRt_2.34.0 stats4_3.4.2
## [85] XML_3.98-1.9 evaluate_0.10.1
## [87] biovizBase_1.26.0 latticeExtra_0.6-28
## [89] data.table_1.10.4-3 httpuv_1.3.5
## [91] gtable_0.2.0 reshape_0.8.7
## [93] assertthat_0.2.0 ggplot2_2.2.1
## [95] mime_0.5 xtable_1.8-2
## [97] AnnotationFilter_1.2.0 survival_2.41-3
## [99] tibble_1.3.4 GenomicAlignments_1.14.0
## [101] AnnotationDbi_1.40.0 memoise_1.1.0
## [103] IRanges_2.12.0 cluster_2.0.6
## [105] interactiveDisplayBase_1.16.0
McLaren, W., B. Pritchard, D. Rios, Y. Chen, P. Flicek, and F. Cunningham. 2010. “Deriving the Consequences of Genomic Variants with the Ensembl API and SNP Effect Predictor.” Journal Article. Bioinformatics 26 (16): 2069–70. doi:10.1093/bioinformatics/btq330.