---
title: "AnnotationHub How-To's"
output:
  BiocStyle::html_document:
    toc: true
vignette: >
  % \VignetteIndexEntry{AnnotationHub: AnnotationHub HOW TO's}
  % \VignetteDepends{AnnotationHub, GenomicFeatures, Rsamtools}
  % \VignetteEngine{knitr::rmarkdown}
---

```{r style, echo = FALSE, results = 'asis', warning=FALSE}
options(width=100)
suppressPackageStartupMessages({
    ## load here to avoid noise in the body of the vignette
    library(AnnotationHub)
    library(GenomicFeatures)
    library(Rsamtools)
    library(VariantAnnotation)
})
BiocStyle::markdown()
```

**Package**: `r Biocpkg("AnnotationHub")`<br />
**Authors**: `r packageDescription("AnnotationHub")[["Author"]] `<br />
**Modified**: Sun Jun 28 10:41:23 2015<br />
**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](http://bioconductor.org/help/workflows/annotation/Annotation_Resources/#OrgDb)
section of the Annotation work flow. Here we discover available `OrgDb`
objects (more than 1000!) for less-model organisms

```{r less-model-org}
library(AnnotationHub)
ah <- AnnotationHub()
query(ah, "OrgDb")
query(ah, c("OrgDb", "Takifugu"))
takifugu <- ah[["AH48147"]]
```

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

```{r less-model-org-select}
keytypes(takifugu)
columns(takifugu)
egid <- head(keys(takifugu, "ENTREZID"))
select(takifugu, egid, c("SYMBOL", "GENENAME"), "ENTREZID")
```

## Roadmap Epigenomics Project 

All Roadmap Epigenomics files are hosted
[here](http://egg2.wustl.edu/roadmap/data/byFileType/). 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.

```{r, eval=FALSE}
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 
```{r results='hide'}
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.

```{r}
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 
```{r}
unique(epiFiles$species)
unique(epiFiles$genome)
```
Broadly, one can get an idea of the different files from this project 
looking at the sourcetype
```{r}
table(epiFiles$sourcetype)
```
To get a more descriptive idea of these different files one can use:
```{r}
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.
```{r}
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

```{r echo=FALSE, results='hide'}
metadata.tab <- ah[["AH41830"]]
```
```{r}
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:

```{r}
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
```{r}
bpChipEpi <- query(ah , c("EpigenomeRoadMap", "broadPeak", "chip", "consolidated"))
```
To get all the bigWig signal files, one can query the hub using 
```{r}
allBigWigFiles <- query(ah, c("EpigenomeRoadMap", "BigWig"))
```
To access the 15 state chromatin segmentations, one can use
```{r}
seg <- query(ah, c("EpigenomeRoadMap", "segmentations"))
```
If one is interested in getting all the files related to one sample
```{r}
E126 <- query(ah , c("EpigenomeRoadMap", "E126", "H3K4ME2"))
E126
```
Hub resources can also be selected using `$`, `subset()`, and
`display()`; see the main
[_AnnotationHub_ vignette](AnnotationHub.html) 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.

```{r echo=FALSE, results='hide'}
peaks <- E126[['AH29817']]
```
```{r}
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.

```{r}
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("GenomicFeatures")`::`makeTxDbFromBiomart()`.

_AnnotationHub_ provides an easy way to work with gene models
published by Ensembl. Let's see what Ensembl's Release-80 has in terms
of data for pufferfish, _Takifugu rubripes_.

```{r takifugu-gene-models}
query(ah, c("Takifugu", "release-80"))
```

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 compressed, indexed Fasta file


```{r takifugu-data}
gtf <- ah[["AH47101"]]
dna <- ah[["AH47477"]]

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

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

```{r takifugu-seqlengths}
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)

```{r takifugu-txdb}
library(GenomicFeatures)         # 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.

```{r takifugu-exons}
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.

```{r}
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

```{r echo=FALSE, results='hide'}
chain <- chainfiles[['AH14150']]
```
```{r}
chain <- chainfiles[['AH14150']]
chain
```
Perform the liftOver operation using `rtracklayer::liftOver()`:

```{r}
library(rtracklayer)
gr38 <- liftOver(peaks, chain)
```
This returns a _GRangeslist_; update the genome of the result to get
the final result

```{r}
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](http://www.ncbi.nlm.nih.gov/variation/docs/human_variation_vcf/).

Query the dbDNP files in the hub:

```{r echo=FALSE, results='hide', message=FALSE}
query(ah, c("GRCh37", "dbSNP", "VCF" ))
vcf <- ah[['AH50424']]
```
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.

```{r message=FALSE}
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

```{r}
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

```{r}
seqlevelsStyle(variants) <-seqlevelsStyle(peaks)
```

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

```{r}
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

```{r}
idx <- subjectHits(overlap) == 3852
overlap[idx]
```

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

```{r}
peaks[3852]
rowRanges(variants)[queryHits(overlap[idx])]
```

# sessionInfo

```{r}
sessionInfo()
```