Contents

1 Introduction

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

Classes of VCF filter rules

2 Installation

The 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)

3 Recurrent settings: TVTBparam

Most functionalities in TVTB require recurrent information such as:

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:

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).

4 Data import

4.1 Genetic variants

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)
    • fields "REF" and "ALT".
  • info(x)
    • field <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.

4.2 List of genomic ranges

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:

4.3 Phenotypes

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.

4.4 Example

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:

  • The TVTBparam constructor previously populated the which slot of svp with “reduced” (i.e. non-overlapping) genomic ranges defined in the ranges slot.
  • Only the INFO key defined in the vep slot will be imported
  • Only the matrix of called genotypes will be imported

Step 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 signature
  • phenotypes
    • The rownames of those phenotypes defines the sample identifiers that are queried from the VCF file.
    • Those phenotypes are stored in the 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

5 Adding allele frequencies

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.

5.1 Adding overall frequencies

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)

5.2 Adding frequencies within phenotype level(s)

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"

6 Filtering variants

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.

6.1 Definition of VCF filter rules

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:

Motivation for each of the new classes extending FilterRules, to define VCF 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:

  • Expressions must explicitely refer to the VCF slots
  • As a consequence, a single expression can refer to fields from different VCF 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.
  • A single rule expression (or function) may refer to multiple columns of the relevant slot in the VCF object.
  • Rules may include calculations, allowing filtering on values not necessarily stored as such in any slot of the 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.

6.2 Control of VCF filter rules

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.

6.3 Evaluation of 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

7 Visualising data

7.1 Visualise INFO data by phenotype level on a genomic axis

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:

  • the points represent the MAF in each super-population on a common Y axis.
  • the heatmap represents the MAF as the color intensity, given a row for each super-population.
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
    )

8 Pairwise comparison of INFO data between phenotype levels

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.

9 A taste of future…

This section presents upcoming features.

9.1 Summarising Ensembl VEP predictions

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:

  • Count of discrete Ensembl VEP predictions
    • by phenotype level
    • by genomic feature affected (i.e. transcript)
  • Distribution of continuous Ensembl VEP predictions
    • by phenotype level
    • by genomic feature affected (i.e. transcript)

10 Acknowledgement

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.

11 Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

## R version 3.4.0 (2017-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.5-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.5-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.2.0      knitr_1.15.1    BiocStyle_2.4.0
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.2.1                    Biobase_2.36.0               
##  [3] AnnotationHub_2.8.0           splines_3.4.0                
##  [5] shiny_1.0.2                   Formula_1.2-1                
##  [7] interactiveDisplayBase_1.14.0 stats4_3.4.0                 
##  [9] latticeExtra_0.6-28           pander_0.6.0                 
## [11] BSgenome_1.44.0               GenomeInfoDbData_0.99.0      
## [13] Rsamtools_1.28.0              yaml_2.1.14                  
## [15] RSQLite_1.1-2                 backports_1.0.5              
## [17] lattice_0.20-35               biovizBase_1.24.0            
## [19] limma_3.32.0                  digest_0.6.12                
## [21] GenomicRanges_1.28.0          RColorBrewer_1.1-2           
## [23] XVector_0.16.0                checkmate_1.8.2              
## [25] colorspace_1.3-2              httpuv_1.3.3                 
## [27] htmltools_0.3.5               Matrix_1.2-9                 
## [29] plyr_1.8.4                    XML_3.98-1.6                 
## [31] ensemblVEP_1.16.0             biomaRt_2.32.0               
## [33] bookdown_0.3                  zlibbioc_1.22.0              
## [35] xtable_1.8-2                  scales_0.4.1                 
## [37] BiocParallel_1.10.0           EnsDb.Hsapiens.v75_2.1.0     
## [39] tibble_1.3.0                  htmlTable_1.9                
## [41] AnnotationFilter_1.0.0        IRanges_2.10.0               
## [43] ggplot2_2.2.1                 SummarizedExperiment_1.6.0   
## [45] GenomicFeatures_1.28.0        nnet_7.3-12                  
## [47] Gviz_1.20.0                   BiocGenerics_0.22.0          
## [49] lazyeval_0.2.0                mime_0.5                     
## [51] survival_2.41-3               magrittr_1.5                 
## [53] memoise_1.1.0                 evaluate_0.10                
## [55] GGally_1.3.0                  foreign_0.8-67               
## [57] BiocInstaller_1.26.0          data.table_1.10.4            
## [59] tools_3.4.0                   matrixStats_0.52.2           
## [61] stringr_1.2.0                 S4Vectors_0.14.0             
## [63] munsell_0.4.3                 cluster_2.0.6                
## [65] DelayedArray_0.2.0            ensembldb_2.0.0              
## [67] AnnotationDbi_1.38.0          Biostrings_2.44.0            
## [69] compiler_3.4.0                GenomeInfoDb_1.12.0          
## [71] grid_3.4.0                    RCurl_1.95-4.8               
## [73] dichromat_2.0-0               VariantAnnotation_1.22.0     
## [75] htmlwidgets_0.8               labeling_0.3                 
## [77] bitops_1.0-6                  base64enc_0.1-3              
## [79] rmarkdown_1.4                 gtable_0.2.0                 
## [81] DBI_0.6-1                     reshape_0.8.6                
## [83] reshape2_1.4.2                R6_2.2.0                     
## [85] GenomicAlignments_1.12.0      gridExtra_2.2.1              
## [87] rtracklayer_1.36.0            Hmisc_4.0-2                  
## [89] rprojroot_1.2                 ProtGenerics_1.8.0           
## [91] stringi_1.1.5                 parallel_3.4.0               
## [93] Rcpp_0.12.10                  rpart_4.1-11                 
## [95] acepack_1.4.1

References

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.