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

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

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)