In Bioconductor/AnnotationHub: Client to access AnnotationHub resources

options(width=100)
suppressPackageStartupMessages({
    ## load here to avoid noise in the body of the vignette
    library(AnnotationHub)
    library(txdbmaker)
    library(Rsamtools)
    library(VariantAnnotation)
})
BiocStyle::markdown()

Package: r Biocpkg("AnnotationHub")
Authors: r packageDescription("AnnotationHub")[["Author"]]
Modified: Mon March 18 2024
Compiled: r date()

Accessing Genome-Scale Data

Non-model organism gene annotations

Bioconductor offers pre-built org.* annotation packages for model organisms, with their use described in the OrgDb section of the Annotation work flow. Here we discover available OrgDb objects for less-model organisms

library(AnnotationHub)
ah <- AnnotationHub()
query(ah, "OrgDb")
orgdb <- query(ah, c("OrgDb", "[e-mail]"))[[1]]

The object returned by AnnotationHub is directly usable with the select() interface, e.g., to discover the available keytypes for querying the object, the columns that these keytypes can map to, and finally selecting the SYMBOL and GENENAME corresponding to the first 6 ENTREZIDs

library(AnnotationDbi)
AnnotationDbi::keytypes(orgdb)
AnnotationDbi::columns(orgdb)
egid <- head(keys(orgdb, "ENTREZID"))
AnnotationDbi::select(orgdb, egid, c("SYMBOL", "GENENAME"), "ENTREZID")

Roadmap Epigenomics Project

All Roadmap Epigenomics files are hosted here. If one had to download these files on their own, one would navigate through the web interface to find useful files, then use something like the following R code.

url <- "http://egg2.wustl.edu/roadmap/data/byFileType/peaks/consolidated/broadPeak/E001-H3K4me1.broadPeak.gz"
filename <-  basename(url)
download.file(url, destfile=filename)
if (file.exists(filename))
   data <- import(filename, format="bed")

This would have to be repeated for all files, and the onus would lie on the user to identify, download, import, and manage the local disk location of these files.

r Biocpkg("AnnotationHub") reduces this task to just a few lines of R code

library(AnnotationHub)
ah = AnnotationHub()
epiFiles <- query(ah, "EpigenomeRoadMap")

A look at the value returned by epiFiles shows us that r length(epiFiles) roadmap resources are available via r Biocpkg("AnnotationHub"). Additional information about the files is also available, e.g., where the files came from (dataprovider), genome, species, sourceurl, sourcetypes.

epiFiles

A good sanity check to ensure that we have files only from the Roadmap Epigenomics project is to check that all the files in the returned smaller hub object come from Homo sapiens and the r unique(epiFiles$genome) genome

unique(epiFiles$species)
unique(epiFiles$genome)

Broadly, one can get an idea of the different files from this project looking at the sourcetype

table(epiFiles$sourcetype)

To get a more descriptive idea of these different files one can use:

sort(table(epiFiles$description), decreasing=TRUE)

The 'metadata' provided by the Roadmap Epigenomics Project is also available. Note that the information displayed about a hub with a single resource is quite different from the information displayed when the hub references more than one resource.

metadata.tab <- query(ah , c("EpigenomeRoadMap", "Metadata"))
metadata.tab

So far we have been exploring information about resources, without downloading the resource to a local cache and importing it into R. One can retrieve the resource using [[ as indicated at the end of the show method

metadata.tab <- ah[["AH41830"]]

metadata.tab <- ah[["AH41830"]]

The metadata.tab file is returned as a data.frame. The first 6 rows of the first 5 columns are shown here:

metadata.tab[1:6, 1:5]

One can keep constructing different queries using multiple arguments to trim down these r length(epiFiles) to get the files one wants. For example, to get the ChIP-Seq files for consolidated epigenomes, one could use

bpChipEpi <- query(ah , c("EpigenomeRoadMap", "broadPeak", "chip", "consolidated"))

To get all the bigWig signal files, one can query the hub using

allBigWigFiles <- query(ah, c("EpigenomeRoadMap", "BigWig"))

To access the 15 state chromatin segmentations, one can use

seg <- query(ah, c("EpigenomeRoadMap", "segmentations"))

If one is interested in getting all the files related to one sample

E126 <- query(ah , c("EpigenomeRoadMap", "E126", "H3K4ME2"))
E126

Hub resources can also be selected using $, subset(), and BiocHubsShiny(); see the main AnnotationHub vignette for additional detail.

Hub resources are imported as the appropriate Bioconductor object for use in further analysis. For example, peak files are returned as GRanges objects.

peaks <- E126[['AH29817']]

peaks <- E126[['AH29817']]
seqinfo(peaks)

BigWig files are returned as BigWigFile objects. A BigWigFile is a reference to a file on disk; the data in the file can be read in using rtracklayer::import(), perhaps querying these large files for particular genomic regions of interest as described on the help page ?import.bw.

Each record inside r Biocpkg("AnnotationHub") is associated with a unique identifier. Most GRanges objects returned by r Biocpkg("AnnotationHub") contain the unique AnnotationHub identifier of the resource from which the GRanges is derived. This can come handy when working with the GRanges object for a while, and additional information about the object (e.g., the name of the file in the cache, or the original sourceurl for the data underlying the resource) that is being worked with.

metadata(peaks)
ah[metadata(peaks)$AnnotationHubName]$sourceurl

Ensembl GTF and FASTA files for TxDb gene models and sequence queries

Bioconductor represents gene models using 'transcript' databases. These are available via packages such as r Biocannopkg("TxDb.Hsapiens.UCSC.hg38.knownGene") or can be constructed using functions such as r Biocpkg("txdbmaker")::makeTxDbFromBiomart().

AnnotationHub provides an easy way to work with gene models published by Ensembl. Let's see what Ensembl's Release-94 has in terms of data for pufferfish, Takifugu rubripes.

query(ah, c("Takifugu", "release-94"))

We see that there is a GTF file descrbing gene models, as well as various DNA sequences. Let's retrieve the GTF and top-level DNA sequence files. The GTF file is imported as a GRanges instance, the DNA sequence as a twobit file.

gtf <- ah[["AH64858"]]
dna <- ah[["AH66116"]]

head(gtf, 3)
dna
head(seqlevels(dna))

Let's identify the 25 longest DNA sequences, and keep just the annotations on these scaffolds.

keep <- names(tail(sort(seqlengths(dna)), 25))
gtf_subset <- gtf[seqnames(gtf) %in% keep]

It is trivial to make a TxDb instance of this subset (or of the entire gtf)

library(txdbmaker)         # for makeTxDbFromGRanges
txdb <- makeTxDbFromGRanges(gtf_subset)

and to use that in conjunction with the DNA sequences, e.g., to find exon sequences of all annotated genes.

library(Rsamtools)               # for getSeq,FaFile-method
exons <- exons(txdb)
length(exons)
getSeq(dna, exons)

There is a one-to-one mapping between the genomic ranges contained in exons and the DNA sequences returned by getSeq().

Some difficulties arise when working with this partly assembled genome that require more advanced GenomicRanges skills, see the r Biocpkg("GenomicRanges") vignettes, especially "GenomicRanges HOWTOs" and "An Introduction to GenomicRanges".

liftOver to map between genome builds

Suppose we wanted to lift features from one genome build to another, e.g., because annotations were generated for hg19 but our experimental analysis used hg18. We know that UCSC provides 'liftover' files for mapping between genome builds.

In this example, we will take our broad Peak GRanges from E126 which comes from the 'hg19' genome, and lift over these features to their 'hg38' coordinates.

chainfiles <- query(ah , c("hg38", "hg19", "chainfile"))
chainfiles

We are interested in the file that lifts over features from hg19 to hg38 so lets download that using

chain <- chainfiles[['AH14150']]

chain <- chainfiles[['AH14150']]
chain

Perform the liftOver operation using rtracklayer::liftOver():

library(rtracklayer)
gr38 <- liftOver(peaks, chain)

This returns a GRangeslist; update the genome of the result to get the final result

genome(gr38) <- "hg38"
gr38

Working with dbSNP Variants

One may also be interested in working with common germline variants with evidence of medical interest. This information is available at NCBI.

Query the dbDNP files in the hub:

query(ah, c("GRCh38", "dbSNP", "VCF" ))
vcf <- ah[['AH57960']]

This returns a VcfFile which can be read in using r Biocpkg("VariantAnnotation"); because VCF files can be large, readVcf() supports several strategies for importing only relevant parts of the file (e.g., particular genomic locations, particular features of the variants), see ?readVcf for additional information.

variants <- readVcf(vcf, genome="hg19")
variants

rowRanges() returns information from the CHROM, POS and ID fields of the VCF file, represented as a GRanges instance

rowRanges(variants)

Note that the broadPeaks files follow the UCSC chromosome naming convention, and the vcf data follows the NCBI style of chromosome naming convention. To bring these ranges in the same chromosome naming convention (ie UCSC), we would use

seqlevelsStyle(variants) <-seqlevelsStyle(peaks)

And then finally to find which variants overlap these broadPeaks we would use:

overlap <- findOverlaps(variants, peaks)
overlap

Some insight into how these results can be interpretted comes from looking a particular peak, e.g., the 3852nd peak

idx <- subjectHits(overlap) == 3852
overlap[idx]

There are three variants overlapping this peak; the coordinates of the peak and the overlapping variants are

peaks[3852]
rowRanges(variants)[queryHits(overlap[idx])]

sessionInfo

sessionInfo()

Bioconductor/AnnotationHub documentation built on April 19, 2024, 7:26 p.m.