1 Getting started

GenomicScores is an R package distributed as part of the Bioconductor project. To install the package, start R and enter:

source("http://bioconductor.org/biocLite.R")
biocLite("GenomicScores")

Once GenomicScores is installed, it can be loaded with the following command.

library(GenomicScores)

Often, however, GenomicScores will be automatically loaded when working with an annotation package that uses GenomicScores, such as phastCons100way.UCSC.hg19.

2 Genomewide position-specific scores

Genomewide scores assign each genomic position a numeric value denoting an estimated measure of constraint or impact on variation at that position. They are commonly used to filter single nucleotide variants or assess the degree of constraint or functionality of genomic features. Genomic scores are built on the basis of different sources of information such as sequence homology, functional domains, physical-chemical changes of amino acid residues, etc.

One particular example of genomic scores are phastCons scores. They provide a measure of conservation obtained from genomewide alignments using the program phast (Phylogenetic Analysis with Space/Time models) from Siepel et al. (2005). The GenomicScores package allows one to retrieve these scores through annotation packages (Section 4) or as AnnotationHub resources (Section 5).

Often, genomic scores such as phastCons are used within workflows running on top of R and Bioconductor. The purpose of the GenomicScores package is to enable an easy and interactive access to genomic scores within those workflows.

3 Lossy storage of genomic scores with compressed vectors

Storing and accessing genomic scores within R is challenging when their values cover large regions of the genome, resulting in gigabytes of double-precision numbers. This is the case, for instance, for phastCons (Siepel et al. 2005), CADD (Kircher et al. 2014) or M-CAP (Jagadeesh et al. 2016) scores.

We address this problem by using lossy compression, also called quantization, coupled with run-length encoding (Rle) vectors. Lossy compression attempts to trade off precision for compression without compromising the scientific integrity of the data (Zender 2016).

Sometimes, measurements and statistical estimates under certain models generate false precision. False precision is essentialy noise that wastes storage space and it is meaningless from the scientific point of view (Zender 2016). In those circumstances, lossy compression not only saves storage space, but also removes false precision.

The use of lossy compression leads to a subset of quantized values much smaller than the original set of genomic scores, resulting in long runs of identical values along the genome. These runs of identical values can be further compressed using the implementation of Rle vectors available in the S4Vectors Bioconductor package.

4 Retrieval of genomic scores through annotation packages

There are currently four different annotation packages that store genomic scores and can be accessed using the GenomicScores package (Table 1):


Table 1: Bioconductor annotation packages storing genomic scores
Annotation Package Description
phastCons100way.UCSC.hg19 phastCons scores derived from the alignment of the human genome (hg19) to other 99 vertebrate species.
phastCons100way.UCSC.hg38 phastCons scores derived from the alignment of the human genome (hg38) to other 99 vertebrate species.
phastCons7way.UCSC.hg38 phastCons scores derived from the alignment of the human genome (hg38) to other 6 mammal species.
fitCons.UCSC.hg19 fitCons scores: fitness consequences of functional annotation for the human genome (hg19).

This is an example of how genomic scores can be retrieved using the phastCons100way.UCSC.hg19 package. Here, a GScores object is created when the package is loaded.

library(phastCons100way.UCSC.hg19)
library(GenomicRanges)
gsco <- phastCons100way.UCSC.hg19

We should use the function scores() to retrive genomic scores for specific positions.

scores(gsco, GRanges(seqnames="chr7", IRanges(start=117232380, width=1)))
## GRanges object with 1 range and 1 metadata column:
##       seqnames                 ranges strand |    scores
##          <Rle>              <IRanges>  <Rle> | <numeric>
##   [1]     chr7 [117232380, 117232380]      * |       0.8
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

5 Retrieval of genomic scores through AnnotationHub resources

Another way to retrieve genomic scores is by using the AnnotationHub, which is a web resource that provides a central location where genomic files (e.g., VCF, bed, wig) and other resources from standard (e.g., UCSC, Ensembl) and distributed sites, can be found. A Bioconductor AnnotationHub web resource creates and manages a local cache of files retrieved by the user, helping with quick and reproducible access.

The first step to retrieve genomic scores is to check the ones available to download.

availableGScores()
## [1] "fitCons.UCSC.hg19"         "phastCons100way.UCSC.hg19"
## [3] "phastCons100way.UCSC.hg38" "phastCons7way.UCSC.hg38"

The selected resource can be downloaded with the function getGScores(). After the resource is downloaded the first time, the cached copy will enable quicker later retrieval.

gsco <- getGScores("phastCons100way.UCSC.hg19")

Finally, the phastCons score of a particular genomic position is retrieved as it has been seen before.

scores(gsco, GRanges(seqnames="chr7", IRanges(start=117232380, width=1)))
## GRanges object with 1 range and 1 metadata column:
##       seqnames                 ranges strand |    scores
##          <Rle>              <IRanges>  <Rle> | <numeric>
##   [1]     chr7 [117232380, 117232380]      * |       0.8
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

6 Annotating variants with genomic scores

A typical use case of the GenomicScores package is in the context of annotating variants with genomic scores, such as phastCons conservation scores. For this purpose, we load the VariantAnnotaiton and TxDb.Hsapiens.UCSC.hg19.knownGene packages. The former will allow us to read a VCF file and annotate it, and the latter contains the gene annotations from UCSC that will be used in this process.

library(VariantAnnotation)
## Loading required package: SummarizedExperiment
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: DelayedArray
## Loading required package: matrixStats
## 
## Attach