4.1.1 Basic Syntax

1. Reading in a MAF file

# get the disk location for maf test file
testFileDir <- system.file("extdata", package = "GenVisR")
testFile <- Sys.glob(paste0(testFileDir, "/brca.maf"))

# define the objects for testing
mafObject <- MutationAnnotationFormat(testFile)

2. Reading in a VEP file

# get the disk location for test files
testFileDir <- system.file("extdata", package = "GenVisR")
testFile <- Sys.glob(paste0(testFileDir, "/*.vep"))

# define the object for testing
vepObject <- VEP(testFile)

3. Reading in a GMS file

# get the disk location for test files
testFileDir <- system.file("extdata", package = "GenVisR")
testFile <- Sys.glob(paste0(testFileDir, "/FL.gms"))

# define the objects for testing
gmsObject <- GMS(testFile)

## This function is part of the new S4 feature and is under active development

4.1.2 Viewing Data

In some cases you may want to view data after it has been read in, there are functions available to make viewing these data easier. Briefly these are getSample(), getPosition(), getMutation(), and getMeta() which are globally available for each object type (GMS, VEP, etc.). Let’s test one out by viewing the samples from the VEP file that was read in with VEP().

# view the samples from the VEP file
getSample(vepObject)

4.1.3 Additional Notes

All of these functions expect a path to a file of the appropriate type. Optionally a wildcard can be supplied with * as was done with VEP() causing the function to read in multiple files at once. With the exception of MAF files these functions expect the files being read in to have a column called “sample” (MAF files already have a sample designation within the file). If a sample column is not found one will be created based on the filenames. Each of these functions will attempt to infer a file specification version from the file header. If a version is not found one will need to be specified via the version parameter.

5.1.1 Using filtering parameters with Waterfall

Relevant Parameters: samples, recurrence, genes, geneMax, mutation

Often it is the case that the input data supplied to the Waterfall() function will contain thousands of genes and hundreds of samples. While Waterfall() can handle such scenarios the graphics device Waterfall() would neeed to output to would have to be enlarged to such a degree that the visualization may become unwieldy. To alleviate such issues and show the most relevant data Waterfall() provides a suite of filtering parameters. Let’s filter our plot using a few of these, suppose we only wanted to visualize those genes which occur in 50% of the cohort and we further didn’t care about sample “FLX0010Naive”. We could tell Waterfall() this by giving the samples we only wanted plotted via the samples parameter and setting the recurrence parameter to .50.

# show those genes which recur in 50% of the cohort for these 4 samples
drawPlot(Waterfall(vepObject, recurrence = 0.5, samples = c("FLX0040Naive", "FLX0070Naive", 
    "FLX0050Naive", "FLX0030Naive")))

## This function is part of the new S4 feature and is under active development, did you mean to use waterfall() with a lower case w?

5.1.2 Mutation sub-plot

Relevant Parameters: coverage, plotA, plotATally, plotALayers

You might have noticed that when filtering the waterfall plot the top sub-plot didn’t change i.e. for sample “FLX0040Naive” the Frequency of mutations remained around 50. You may also have noticed a warning when constructing a waterfall plot for the VEP data saying that duplicate genomic locations were detected. This is important to note and understand, regardless of any filtering that occurs the top sub-plot will always be based on the original input with one caveat. It will check for duplicate variants (i.e. same sample, allele, genomic position) and remove variants which are deemed to be duplicated. This occurs in the case of VEP as different transcripts can be annotated with different mutation types but still be the same variant which is why we see the warning. By default Waterfall() will simply output the frequencies observed in the cohort. However a mutation burden can be displayed as well by changing the value for parameter plotA from “frequency” to “burden”. In doing so you will also need to supply a value to coverage giving the aproximate space in base pairs for which a mutation could have been called. Doing so will plot a mutation burden instead of a mutation frequency calculated via the folowing formula:

Mutation Burden = (# of observed mutations per sample)/(# of bases for which a mutation could be called) * 1000000

For a more accurate mutation burden calculation we can also supply coverage values individual per sample via a named vector, let’s do that and tell Waterfall() to calculate a mutation burden by mutation type + sample rather than sample alone via the parameter plotATally. Lastly we will set the drop parameter to FALSE to avoid dropping mutations from the legend that are not in the main plot but are in the top sub-plot.

# define a coverage for each sample
sampCov <- c(FLX0040Naive = 14500000, FLX0070Naive = 13900000, FLX0050Naive = 12100000, 
    FLX0030Naive = 1.3e+07, FLX0010Naive = 1.1e+07)

drawPlot(Waterfall(vepObject, recurrence = 0.5, coverage = sampCov, plotA = "burden", 
    plotATally = "complex", drop = FALSE))

## This function is part of the new S4 feature and is under active development, did you mean to use waterfall() with a lower case w?

5.1.3 Altering the mutation colors and hierarchy

Relevant Parameters: mutationHierarchy

You may ask yourself what happens when there are two different mutations for the same gene/sample, which one is plotted? For each file type a pre-defined hierarchy of mutations has been established which follows the order of the “Mutation Type” legend from top to bottom. By default mutations are based on a hierarchy of being most deleterious however both the priority and color of mutations can be changed by supplying a data.table to the mutationHierarchy parameter. Let’s assume that we are particularly interested in splice variants, by default these might not show up because a “missense_variant” is deemed to be more deleterious but let’s change that by supplying a custom hierarchy along with colors to go along with it. To do that we will need to give the parameter a data.table object with column names “mutation” and “color”.

# find which mutations are in the input data
mutations <- getMutation(vepObject)$Conse

# define a new color and hierarchy for mutations
library(data.table)

## 
## Attaching package: 'data.table'

## The following object is masked from 'package:GenomicRanges':
## 
##     shift

## The following object is masked from 'package:IRanges':
## 
##     shift

## The following objects are masked from 'package:S4Vectors':
## 
##     first, second

## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt

newHierarchy <- data.table(mutation = c("splice_region_variant", "splice_acceptor_variant", 
    "splice_donor_variant", "missense_variant", "stop_gained"), color = c("tomato1", 
    "tomato2", "tomato3", "purple", "cyan"))

# draw the plot
drawPlot(Waterfall(vepObject, recurrence = 0.5, mutationHierarchy = newHierarchy, 
    plotATally = "complex", drop = FALSE))

## This function is part of the new S4 feature and is under active development, did you mean to use waterfall() with a lower case w?

Here we can see that for sample “FLX0040Naive” on gene BIRC6 there is actually a splice_region_variant as well as a missense_variant.

Note: If mutations are found in the input but are not listed in the mutationHierarchy parameter they will be addded and assigned random colors

Note: In a VEP formated file mutations can be comma delimited, if they are not specifically specified in the mutation hierarchy as comma delimited the mutations are split up and collapsed

5.2.1 Input

Input to the MutSpectra() function can be just one of the afore mentioned and supported file format object discussed in section 3.1.1. In cases where the file format does not provide a reference base, such as with vep, a BSgenome object will also need to be provided in order to annotate reference bases for a given position. These BSgenome objects are a core object class in bioconductor and can be installed through the instructions on the bioconductor site. Let’s go ahead and create a MutSpectra plot for the vep data. We first determine which assembly vep used for these annotations ussing the getHeader() function and can see that it is GRCh37. We next load in the BSgenome for GRCh37 (GRCh37 and hg19 are analogous assemblies with very minor differences) and supply the BSgenome object to the parameter BSgenome.

# determine the appropriate BSgenome object to use from the vep header
assembly <- getHeader(vepObject)
assembly <- assembly[grepl("assembly", assembly$Info), ]

# load in the correct BSgenome object
library(BSgenome.Hsapiens.UCSC.hg19)

# create a MutSpectra plot
drawPlot(MutSpectra(vepObject, BSgenome = BSgenome.Hsapiens.UCSC.hg19))

Note: In the example above you’ll see a warning related to chromosome mismatches, this is part of the difference between GRCh37 and hg19, chromosomes are designated differently between the two (i.e. chr1 and 1). MutSpectra() is smart enough to recognize this and append chr to all chromosomes.

Note: If you are unsure what assembly to use try running MutSpectra() without specifying one. It will attempt to match an assembly from bioconductor with the information in the header.

5.3.1 Input

As with the MutSpectra plots basic input consists of an object corresponding to one of the supported file types. Further a BSgenome object is required if the input format does not contain a reference column (see the MutSpectra section for more details). Let’s start by creating a Rainfall plot for our VEP object. You might notice in the plot that density estimates are missing for a few chromosomes, this is normal and is done intentionally to avoid artificially skewing the y-axis when there are too few points to construct a confident density plot.

# determine the appropriate BSgenome object to use from the vep header
assembly <- getHeader(vepObject)
assembly <- assembly[grepl("assembly", assembly$Info), ]

# load in the correct BSgenome object
library(BSgenome.Hsapiens.UCSC.hg19)

# create a Rainfall plot
drawPlot(Rainfall(vepObject, BSgenome = BSgenome.Hsapiens.UCSC.hg19))

5.3.2 Using filtering parameters with Rainfall

By default the Rainfall() function plots everything within the input, this means that it is not only plotting all chromosomes but also all samples. In some cases this may not be desireable but fortunately Rainfall() has two parameters to limit the data which is plotted, these are sample and chromosomes which will limit the samples and chromosomes ploted respectively. Let’s go ahead and try them out, we’ll limit our plot to only chromosomes 1-3 and only sample “FLX0010Naive”.

# create a Rainfall plot limiting the chromosomes and samples plotted
drawPlot(Rainfall(vepObject, BSgenome = BSgenome.Hsapiens.UCSC.hg19, sample = c("FLX0010Naive"), 
    chromosomes = c("chr1", "chr2", "chr3")))