1 Getting started

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

install.packages("BiocManager")
BiocManager::install("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.hg38.

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) or CADD (Kircher et al. 2014).

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.

To enable a seamless access to genomic scores stored with quantized values in compressed vectors the GenomicScores defines the GScores class of objects. This class manages the location, loading and dequantization of genomic scores stored separately on each chromosome, while it also provides rich metadata on the provenance, citation and licensing of the original data.

4 Availability and retrieval of genomic scores

The access to genomic scores through GScores objects is available either through annotation packages or as AnnotationHub resources. To find out what kind of genomic scores are available, through which mechanism, and in which organism, we may use the function availableGScores().

avgs <- availableGScores()
avgs
                                                 Organism      Category
AlphaMissense.v2023.hg19                     Homo sapiens Pathogenicity
AlphaMissense.v2023.hg38                     Homo sapiens Pathogenicity
MafDb.1Kgenomes.phase1.GRCh38                Homo sapiens           MAF
MafDb.1Kgenomes.phase1.hs37d5                Homo sapiens           MAF
MafDb.1Kgenomes.phase3.GRCh38                Homo sapiens           MAF
MafDb.1Kgenomes.phase3.hs37d5                Homo sapiens           MAF
MafDb.ExAC.r1.0.GRCh38                       Homo sapiens           MAF
MafDb.ExAC.r1.0.hs37d5                       Homo sapiens           MAF
MafDb.ExAC.r1.0.nonTCGA.GRCh38               Homo sapiens           MAF
MafDb.ExAC.r1.0.nonTCGA.hs37d5               Homo sapiens           MAF
MafDb.TOPMed.freeze5.hg19                    Homo sapiens           MAF
MafDb.TOPMed.freeze5.hg38                    Homo sapiens           MAF
MafDb.gnomAD.r2.1.GRCh38                     Homo sapiens           MAF
MafDb.gnomAD.r2.1.hs37d5                     Homo sapiens           MAF
MafDb.gnomADex.r2.1.GRCh38                   Homo sapiens           MAF
MafDb.gnomADex.r2.1.hs37d5                   Homo sapiens           MAF
MafH5.gnomAD.v3.1.2.GRCh38                   Homo sapiens           MAF
cadd.v1.3.hg19                                       <NA>          <NA>
cadd.v1.6.hg19                               Homo sapiens Pathogenicity
cadd.v1.6.hg38                               Homo sapiens Pathogenicity
fitCons.UCSC.hg19                            Homo sapiens    Constraint
linsight.UCSC.hg19                           Homo sapiens    Constraint
mcap.v1.0.hg19                               Homo sapiens Pathogenicity
phastCons100way.UCSC.hg19                    Homo sapiens  Conservation
phastCons100way.UCSC.hg38                    Homo sapiens  Conservation
phastCons27way.UCSC.dm6           Drosophila melanogaster  Conservation
phastCons30way.UCSC.hg38                     Homo sapiens  Conservation
phastCons35way.UCSC.mm39                     Mus musculus  Conservation
phastCons46wayPlacental.UCSC.hg19            Homo sapiens  Conservation
phastCons46wayPrimates.UCSC.hg19             Homo sapiens  Conservation
phastCons60way.UCSC.mm10                     Mus musculus  Conservation
phastCons7way.UCSC.hg38                      Homo sapiens  Conservation
phyloP100way.UCSC.hg19                       Homo sapiens  Conservation
phyloP100way.UCSC.hg38                       Homo sapiens  Conservation
phyloP35way.UCSC.mm39                        Mus musculus  Conservation
phyloP60way.UCSC.mm10                        Mus musculus  Conservation
                                  Installed Cached BiocManagerInstall
AlphaMissense.v2023.hg19              FALSE   TRUE              FALSE
AlphaMissense.v2023.hg38              FALSE   TRUE              FALSE
MafDb.1Kgenomes.phase1.GRCh38         FALSE  FALSE               TRUE
MafDb.1Kgenomes.phase1.hs37d5          TRUE  FALSE               TRUE
MafDb.1Kgenomes.phase3.GRCh38          TRUE  FALSE               TRUE
MafDb.1Kgenomes.phase3.hs37d5          TRUE  FALSE               TRUE
MafDb.ExAC.r1.0.GRCh38                FALSE  FALSE               TRUE
MafDb.ExAC.r1.0.hs37d5                 TRUE  FALSE               TRUE
MafDb.ExAC.r1.0.nonTCGA.GRCh38        FALSE  FALSE               TRUE
MafDb.ExAC.r1.0.nonTCGA.hs37d5        FALSE  FALSE               TRUE
MafDb.TOPMed.freeze5.hg19             FALSE  FALSE               TRUE
MafDb.TOPMed.freeze5.hg38             FALSE  FALSE               TRUE
MafDb.gnomAD.r2.1.GRCh38              FALSE  FALSE               TRUE
MafDb.gnomAD.r2.1.hs37d5              FALSE  FALSE               TRUE
MafDb.gnomADex.r2.1.GRCh38            FALSE  FALSE               TRUE
MafDb.gnomADex.r2.1.hs37d5             TRUE  FALSE               TRUE
MafH5.gnomAD.v3.1.2.GRCh38             TRUE  FALSE               TRUE
cadd.v1.3.hg19                        FALSE  FALSE              FALSE
cadd.v1.6.hg19                        FALSE   TRUE              FALSE
cadd.v1.6.hg38                        FALSE   TRUE              FALSE
fitCons.UCSC.hg19                     FALSE  FALSE               TRUE
linsight.UCSC.hg19                    FALSE  FALSE              FALSE
mcap.v1.0.hg19                        FALSE  FALSE              FALSE
phastCons100way.UCSC.hg19              TRUE  FALSE               TRUE
phastCons100way.UCSC.hg38              TRUE  FALSE               TRUE
phastCons27way.UCSC.dm6               FALSE  FALSE              FALSE
phastCons30way.UCSC.hg38              FALSE  FALSE              FALSE
phastCons35way.UCSC.mm39              FALSE  FALSE              FALSE
phastCons46wayPlacental.UCSC.hg19     FALSE  FALSE              FALSE
phastCons46wayPrimates.UCSC.hg19      FALSE  FALSE              FALSE
phastCons60way.UCSC.mm10              FALSE  FALSE              FALSE
phastCons7way.UCSC.hg38               FALSE  FALSE               TRUE
phyloP100way.UCSC.hg19                FALSE  FALSE              FALSE
phyloP100way.UCSC.hg38                FALSE  FALSE              FALSE
phyloP35way.UCSC.mm39                 FALSE  FALSE              FALSE
phyloP60way.UCSC.mm10                 FALSE  FALSE              FALSE
                                  AnnotationHub
AlphaMissense.v2023.hg19                   TRUE
AlphaMissense.v2023.hg38                   TRUE
MafDb.1Kgenomes.phase1.GRCh38             FALSE
MafDb.1Kgenomes.phase1.hs37d5             FALSE
MafDb.1Kgenomes.phase3.GRCh38             FALSE
MafDb.1Kgenomes.phase3.hs37d5             FALSE
MafDb.ExAC.r1.0.GRCh38                    FALSE
MafDb.ExAC.r1.0.hs37d5                    FALSE
MafDb.ExAC.r1.0.nonTCGA.GRCh38            FALSE
MafDb.ExAC.r1.0.nonTCGA.hs37d5            FALSE
MafDb.TOPMed.freeze5.hg19                 FALSE
MafDb.TOPMed.freeze5.hg38                 FALSE
MafDb.gnomAD.r2.1.GRCh38                  FALSE
MafDb.gnomAD.r2.1.hs37d5                  FALSE
MafDb.gnomADex.r2.1.GRCh38                FALSE
MafDb.gnomADex.r2.1.hs37d5                FALSE
MafH5.gnomAD.v3.1.2.GRCh38                FALSE
cadd.v1.3.hg19                            FALSE
cadd.v1.6.hg19                             TRUE
cadd.v1.6.hg38                             TRUE
fitCons.UCSC.hg19                         FALSE
linsight.UCSC.hg19                        FALSE
mcap.v1.0.hg19                            FALSE
phastCons100way.UCSC.hg19                 FALSE
phastCons100way.UCSC.hg38                 FALSE
phastCons27way.UCSC.dm6                   FALSE
phastCons30way.UCSC.hg38                  FALSE
phastCons35way.UCSC.mm39                  FALSE
phastCons46wayPlacental.UCSC.hg19         FALSE
phastCons46wayPrimates.UCSC.hg19          FALSE
phastCons60way.UCSC.mm10                  FALSE
phastCons7way.UCSC.hg38                   FALSE
phyloP100way.UCSC.hg19                    FALSE
phyloP100way.UCSC.hg38                    FALSE
phyloP35way.UCSC.mm39                     FALSE
phyloP60way.UCSC.mm10                     FALSE

For example, if we want to use the phastCons conservation scores available through the annotation package phastCons100way.UCSC.hg38, we should first install it (we only need to do this once).

BiocManager::install("phastCons100way.UCSC.hg38")

Second, we should load the package, and a GScores object will be created and named after the package name, during the loading operation. It is often handy to shorten that name.

library(phastCons100way.UCSC.hg38)
phast <- phastCons100way.UCSC.hg38
class(phast)
[1] "GScores"
attr(,"package")
[1] "GenomicScores"

Typing the name of the GScores object shows a summary of its contents and some of its metadata.

phast
GScores object 
# organism: Homo sapiens (UCSC, hg38)
# provider: UCSC
# provider version: 11May2015
# download date: Apr 10, 2018
# loaded sequences: chr5_GL000208v1_random
# number of sites: 2943 millions
# maximum abs. error: 0.05
# use 'citation()' to cite these data in publications

The bibliographic reference to cite the genomic score data stored in a GScores object can be accessed using the citation() method either on the package name (in case of annotation packages), or on the GScores object.

citation(phast)
Adam Siepel, Gill Berejano, Jakob S. Pedersen, Angie S. Hinrichs,
Minmei Hou, Kate Rosenbloom, Hiram Clawson, John Spieth, LaDeana W.
Hillier, Stephen Richards, George M. Weinstock, Richard K. Wilson,
Richard A. Gibbs, W. James Kent, Webb Miller, David Haussler (2005).
"Evolutionarily conserved elements in vertebrate, insect, worm, and
yeast genomes." _Genome Research_, *15*, 1034-1050.
doi:10.1101/gr.3715005 <https://doi.org/10.1101/gr.3715005>.

Other methods tracing provenance and other metadata are provider(), providerVersion(), organism() and seqlevelsStyle(); please consult the help page of the GScores class for a comprehensive list of available methods.

provider(phast)
[1] "UCSC"
providerVersion(phast)
[1] "11May2015"
organism(phast)
[1] "Homo sapiens"
seqlevelsStyle(phast)
[1] "UCSC"

To retrieve genomic scores for specific consecutive positions we should use the method gscores(), as follows.

gscores(phast, GRanges(seqnames="chr22",
                       IRanges(start=50528591:50528596, width=1)))
GRanges object with 6 ranges and 1 metadata column:
      seqnames    ranges strand |   default
         <Rle> <IRanges>  <Rle> | <numeric>
  [1]    chr22  50528591      * |       1.0
  [2]    chr22  50528592      * |       1.0
  [3]    chr22  50528593      * |       0.8
  [4]    chr22  50528594      * |       1.0
  [5]    chr22  50528595      * |       1.0
  [6]    chr22  50528596      * |       0.0
  -------
  seqinfo: 455 sequences (1 circular) from Genome Reference Consortium GRCh38 genome

For a single position we may use this other GRanges() constructor.

gscores(phast, GRanges("chr22:50528593"))
GRanges object with 1 range and 1 metadata column:
      seqnames    ranges strand |   default
         <Rle> <IRanges>  <Rle> | <numeric>
  [1]    chr22  50528593      * |       0.8
  -------
  seqinfo: 455 sequences (1 circular) from Genome Reference Consortium GRCh38 genome

We may also retrieve the score values only with the method score().

score(phast, GRanges(seqnames="chr22",
                     IRanges(start=50528591:50528596, width=1)))
[1] 1.0 1.0 0.8 1.0 1.0 0.0
score(phast, GRanges("chr22:50528593"))
[1] 0.8

Let’s illustrate how to retrieve phastCons scores using data from the GWAS catalog available through the Bioconductor package gwascat. For the purpose of this vignette, we will filter the GWAS catalog data by (1) discarding entries with NA values in either chromosome name or position, or with multiple positions; (2) storing the data into a GRanges object, including the GWAS catalog columns STRONGEST SNP-RISK ALLELE and MAPPED_TRAIT, and the reference and alternate alleles, as metadata columns; (4) restricting variants to those located in chromosomes 20 to 22; and (3) excluding variants with multinucleotide alleles, or where reference and alternate alleles are identical.

library(BSgenome.Hsapiens.UCSC.hg38)
library(gwascat)

gwc <- get_cached_gwascat()
mask <- !is.na(gwc$CHR_ID) & !is.na(gwc$CHR_POS) &
        !is.na(as.integer(gwc$CHR_POS))
gwc <- gwc[mask, ]
grstr <- sprintf("%s:%s-%s", gwc$CHR_ID, gwc$CHR_POS, gwc$CHR_POS)
gwcgr <- GRanges(grstr, RISK_ALLELE=gwc[["STRONGEST SNP-RISK ALLELE"]],
                 MAPPED_TRAIT=gwc$MAPPED_TRAIT)
seqlevelsStyle(gwcgr) <- "UCSC"
mask <- seqnames(gwcgr) %in% c("chr20", "chr21", "chr22")
gwcgr <- gwcgr[mask]
ref <- as.character(getSeq(Hsapiens, gwcgr))
alt <- gsub("rs[0-9]+-", "", gwcgr$RISK_ALLELE)
mask <- (ref %in% c("A", "C", "G", "T")) & (alt %in% c("A", "C", "G", "T")) &
         nchar(alt) == 1 & ref != alt
gwcgr <- gwcgr[mask]
mcols(gwcgr)$REF <- ref[mask]
mcols(gwcgr)$ALT <- alt[mask]
gwcgr
GRanges object with 13540 ranges and 4 metadata columns:
          seqnames    ranges strand |   RISK_ALLELE           MAPPED_TRAIT
             <Rle> <IRanges>  <Rle> |   <character>            <character>
      [1]    chr20  35321981      * |   rs6088792-T            body height
      [2]    chr22  23250864      * |   rs5751614-A            body height
      [3]    chr20   6640246      * |    rs967417-C            body height
      [4]    chr20  33130847      * |   rs6059101-A     ulcerative colitis
      [5]    chr22  38148291      * |   rs2284063-G                  nevus
      ...      ...       ...    ... .           ...                    ...
  [13536]    chr20  12979237      * |   rs1321940-G sphingomyelin measur..
  [13537]    chr20  10148004      * |   rs2210584-C sphingomyelin measur..
  [13538]    chr20  11892294      * | rs397865364-C sphingomyelin measur..
  [13539]    chr20  12823511      * |   rs1413019-A sphingomyelin measur..
  [13540]    chr21  46284982      * |   rs9975588-A S100 calcium-binding..
                  REF         ALT
          <character> <character>
      [1]           C           T
      [2]           G           A
      [3]           G           C
      [4]           C           A
      [5]           A           G
      ...         ...         ...
  [13536]           A           G
  [13537]           T           C
  [13538]           T           C
  [13539]           C           A
  [13540]           G           A
  -------
  seqinfo: 24 sequences from an unspecified genome; no seqlengths

Finally, let’s obtain the phastCons scores for this GWAS catalog variant set, and examine their summary and cumulative distribution.

pcsco <- score(phast, gwcgr)
summary(pcsco)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0000  0.0000  0.0000  0.1217  0.0000  1.0000      38 
round(cumsum(table(na.omit(pcsco))) / sum(!is.na(pcsco)), digits=2)
   0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1 
0.81 0.85 0.87 0.88 0.88 0.89 0.89 0.90 0.90 0.91 1.00 

We can observe that only 10% of the variants in chromosomes 20 to 22 have a conservation phastCons score above 0.5. Let’s examine which traits have more fully conserved variants.

xtab <- table(gwcgr$MAPPED_TRAIT[pcsco == 1])
head(xtab[order(xtab, decreasing=TRUE)])

                                     body height 
                                              63 
                        neuroimaging measurement 
                                              34 
                         mean corpuscular volume 
                                              33 
high density lipoprotein cholesterol measurement 
                                              30 
                          hemoglobin measurement 
                                              26 
                    BMI-adjusted waist-hip ratio 
                                              23 

5 Genomic scores as AnnotationHub resources

The AnnotationHub (AH), is a Bioconductor 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. An AH web resource creates and manages a local cache of files retrieved by the user, helping with quick and reproducible access.

We can quickly check for the available AH resources by subsetting as follows the resources names from the previous table obtained with availableGScores().

rownames(avgs)[avgs$AnnotationHub]
[1] "AlphaMissense.v2023.hg19" "AlphaMissense.v2023.hg38"
[3] "cadd.v1.6.hg19"           "cadd.v1.6.hg38"          

The selected resource can be downloaded with the function getGScores(). After the resource is downloaded the first time, the cached copy will enable a quicker retrieval later. Let’s download other conservation scores, the phyloP scores (Pollard et al. 2010), for human genome version hg38.

phylop <- getGScores("phyloP100way.UCSC.hg38")
phylop
GScores object 
# organism: Homo sapiens (UCSC, hg38)
# provider: UCSC
# provider version: 11May2015
# download date: May 12, 2017
# loaded sequences: chr20
# maximum abs. error: 0.55
# use 'citation()' to cite these data in publications

Let’s retrieve the phyloP conservation scores for the previous set of GWAS catalog variants and compare them in Figure @(fig:phastvsphylop).

ppsco <- score(phylop, gwcgr)
plot(pcsco, ppsco, xlab="phastCons", ylab="phyloP",
     cex.axis=1.2, cex.lab=1.5, las=1)
Comparison between phastCons and phyloP conservation scores. On the y-axis, phyloP scores as function of phastCons scores on the x-axis, for a set of GWAS catalog variant in the human chromosome 22.

Figure 1: Comparison between phastCons and phyloP conservation scores
On the y-axis, phyloP scores as function of phastCons scores on the x-axis, for a set of GWAS catalog variant in the human chromosome 22.

We may observe that the values match in a rather discrete manner due to the quantization of the scores. In the case of the phastCons annotation package phastCons100way.UCSC.hg38, the GScore object gives access in fact to two score populations, the default one in which conservation scores are rounded to 1 decimal place, and an alternative one, named DP2, in which they are rounded to 2 decimal places. To figure out what are the available score populations in a GScores object, we should use the method populations().

populations(phast)
[1] "default" "DP2"    

Whenever one of these populations is called default, this is the one used by default. In other cases we can find out which is the default population as follows:

defaultPopulation(phast)
[1] "default"

To use one of the available score populations we should use the argument pop in the corresponding method, as follows.

pcsco2 <- score(phast, gwcgr, pop="DP2")
head(pcsco2)
[1] 0.00 0.00 0.00 0.14 0.00 0.00

Figure 2 below shows again the comparison of phastCons and phyloP conservation scores, this time at the higher resolution provided by the phastCons scores rounded at two decimal places.

plot(pcsco2, ppsco, xlab="phastCons", ylab="phyloP",
     cex.axis=1.2, cex.lab=1.5, las=1)