```{r front, child="front.Rmd", echo=FALSE}
```

## Dependencies

This document has the following dependencies:

```{r dependencies, warning=FALSE, message=FALSE}
library(GenomicRanges)
library(rtracklayer)
library(AnnotationHub)
```

Use the following commands to install these packages in R.

```{r biocLite, eval=FALSE}
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("GenomicRanges", "rtracklayer", "AnnotationHub"))
```

## Corrections

Improvements and corrections to this document can be submitted on its [GitHub](https://github.com/kasperdanielhansen/genbioconductor/blob/master/Rmd/Usecase_AnnotationHub_GRanges.Rmd) in its [repository](https://github.com/kasperdanielhansen/genbioconductor).

## Overview

We are going to use `r Biocpkg("AnnotationHub")` and `r Biocpkg("GenomicRanges")` to access ENCODE data on the H3K4me3 histone modification in a specific cell line.  This histone modification is believed to mark active promoters, and we are going to attempt to verify this statement.  This involves

1. Getting the ENCODE histone data using `r Biocpkg("AnnotationHub")`.
2. Getting promoters using `r Biocpkg("AnnotationHub")`.
3. Comparing the histone data and promoters using `findOverlaps` in `r Biocpkg("GenomicRanges")`.

## Getting the data

First we use `r Biocpkg("AnnotationHub")` to get data on homo sapiens.

```{r ahub_species}
ah <- AnnotationHub()
ah <- subset(ah, species == "Homo sapiens")
```

Next we search for two keywords: **H3K4me3** and **Gm12878** which is the name of the cell line we are interested in.

```{r ahub_histone}
qhs <- query(ah, "H3K4me3")
qhs <- query(qhs, "Gm12878")
```

(Note: I like to keep my full annotation hub around, so I can re-do my query with a different search term in case I end up with no hits.  This is why I start assigning output to the `qhs` object and not `ah`).

Lets have a look
```{r ahub_look}
qhs
```
Note how some of these hits don't contain **Gm12878** in their title.  This is a useful illustration of how `query` searches over multiple fields.

Lets have a closer look at this
```{r ahub_closerlook}
qhs$title
qhs$dataprovider
```
This result is a great illustration of the mess of public data.  It turns our that **E116** is a Roadmap Epigenetics code for the **Gm12878** cell line.  The first 5 hits are from ENCODE, hosted at UCSC and the last 6 hits are from Roadmap Epigenomics hosted at the Broad Institute.  The Roadmap data is different representation (and peaks) from the same underlying data.  For the ENCODE data, two different centers did the same experiment in the same cell line (`Broad`, hit 1) and (`Uw`, hit 2-5), where `Uw` exposed data on two replicates (whatever that means).  These two experiments seems to be analyzed using different algorithms.  It is even possible that the Roadmap data is from the same raw data but just analyzed using different algorithms.

Lets take a look at the `narrowPeak` data:
```{r ahub_twoGR}
gr1 <- subset(qhs, title == "wgEncodeUwHistoneGm12878H3k4me3StdPkRep1.narrowPeak.gz")[[1]]
gr1
gr2 <- subset(qhs, title == "E116-H3K4me3.narrowPeak.gz")[[1]]
gr2
```
(Note: I use this code, where I use `title` to refer to the different resources, to make this script more robust over time).

This gives us two `GRanges` objects.  Let us look at the distribution of peak widths:
```{r ahub_summary}
summary(width(gr1))
table(width(gr1))
summary(width(gr2))
```
Turns out that almost all peaks in `gr1` have a width of 150bp, whereas `gr2` is much more variable.  This is likely a product of the data processing algorithm; it can be very hard to figure out the details of this.

At this time one can spend a lot of time thinking about which datasets is best.  We will avoid this (important) discussion; since we referred to ENCODE data above (in the Overview section), we will stick with `gr1`.

Now we need to get some promoter coordinates.  There are multiple ways to do this in Bioconductor, but here I will do a quick lookup for RefSeq in my annotation hub.  RefSeq is a highly curated (aka conservative) collection of genes.  

Lets get started
```{r ahub_refseq}
qhs <- query(ah, "RefSeq")
qhs
```
So this gives us multiple datasets, all with very similar names.  We probably need the thing called `RefSeq Genes` and not `Other RefSeq` but why are there multiple resources with the same name?

Turns out the answer makes sense:
```{r ahub_refseq_genome}
qhs$genome
```
This looks like the same resources, but in different genome builds.  We have
```{r ahub_histone_genome}
genome(gr1)
```

so we know which one to get:
```{r ahub_get_refseq}
refseq <- qhs[qhs$genome == "hg19" & qhs$title == "RefSeq Genes"]
refseq
refseq <- refseq[[1]] ## Downloads
```

Lets have a look
```{r refseq}
refseq
```

Let us look at the number of isoforms per gene name:
```{r ahub_refseq_name}
table(table(refseq$name))
```

(the `table(table())` construction may seems weird at first, but its a great way to get a quick tabular summary of occurrences with the same name).  This shows that almost all genes have a single transcript, which reflects how conservative RefSeq is.

Notice that each transcript is a separate range.  Therefore, each transcript is represented as the "outer" coordinates of the gene; the ranges does not exclude introns.  Because we got this from UCSC I happen to know that the `$blocks` metadata really contains the coordinates of the different exons.  In a later example we will introduce a gene representation where we keep track of exons, introns and transcripts from the same gene, through something known as a `TxDb` (transcript database) object.

For now, we just need the promoters, so we don't really care about introns.  We need to keep track of which strand each transcript is on, to get the transcription start site.

Or we could just use the convenience function `promoters()`:
```{r promoters}
promoters <- promoters(refseq)
table(width(promoters))
args(promoters)
```

There are many definitions of promoters based on transcription start site.  The default in this function is to use 2kb upstream and 200bp downstream of the start site.  Let's keep this.

Now we compute which promoters have a H3K4me3 peak in them:
```{r findOverlaps}
ov <- findOverlaps(promoters, gr1)
ov
```

Let us compute how many percent of the peaks are in a promoter
```{r queryHits}
length(unique(queryHits(ov))) / length(gr1)
```

and how many percent of promoters have a peak in them
```{r subjectHits}
length(unique(subjectHits(ov))) / length(promoters)
```

My rule of thumb is that any cell type has at most 50% of genes expressed, which fits well with these numbers.  We also see that there are many H3K4me3 peaks which do not lie in a genic promoter.  This is actually expected.

Now, the overlap is non-empty, but is it interesting or "significant".  Answering this question thoroughly requires thinking deeply about background distributions etc.  We will avoid a careful statistical approach to this question, and will instead do a few back-of-the-envelope calculations to get a feel for it.

First we notice that both promoters and peaks are small compared to the size of the human genome:
```{r widthPercentage}
sum(width(reduce(gr1))) / 10^6
sum(width(reduce(promoters))) / 10^6
```

(result is in megabases; the human genome is 3000 megabases).  Just this fact alone should convince you that the overlap is highly unlikely happen purely by chance. Let us calculate the overlap in megabases:
```{r size}
sum(width(intersect(gr1, promoters))) / 10^6
```
That's weird; why is the overlap empty when we know (using `findOverlaps`) that there is plenty of overlap?  Turns out the answer is strand.  We need to ignore the strand when we do this calculation:
```{r size2}
sum(width(intersect(gr1, promoters, ignore.strand = TRUE))) / 10^6
```

This brings us back to the question of the size of the promoter set above; this is also affected by strand when there is a promoter on each strand which overlaps.  Contrast
```{r widthPercentage2}
sum(width(reduce(promoters))) / 10^6
sum(width(reduce(promoters, ignore.strand = TRUE))) / 10^6
```
Strand can be troublesome!

Let us compute a small 2x2 matrix for which bases are in promoters and/or peaks:
```{r promInOut}
prom <- reduce(promoters, ignore.strand = TRUE)
peaks <- reduce(gr1)
both <- intersect(prom, peaks)
only.prom <- setdiff(prom, both)
only.peaks <- setdiff(peaks, both)
overlapMat <- matrix(0,, ncol = 2, nrow = 2)
colnames(overlapMat) <- c("in.peaks", "out.peaks")
rownames(overlapMat) <- c("in.promoters", "out.promoter")
overlapMat[1,1] <- sum(width(both))
overlapMat[1,2] <- sum(width(only.prom))
overlapMat[2,1] <- sum(width(only.peaks))
overlapMat[2,2] <- 3*10^9 - sum(overlapMat)
round(overlapMat / 10^6, 2)
```
Here we have just used the genome size of 3 billion bases.  This is not correct.  There are many of these bases which have not been sequenced and in addition, there are many bases which cannot be mapped using short reads.  This will reduce the genome size.

Nevertheless, let us compute an odds-ratio for this table:
```{r oddsratio}
oddsRatio <- overlapMat[1,1] * overlapMat[2,2] / (overlapMat[2,1] * overlapMat[1,2])
oddsRatio
```

This odds-ratio shows an enrichment of peaks in promoters.  We can get a feel for how much the genome size (which we use incorrectly) affects our result by using a lower bound on the genome size.  Let us say 1.5 billion bases:
```{r oddsRatio2}
overlapMat[2,2] <- 1.5*10^9
oddsRatio <- overlapMat[1,1] * overlapMat[2,2] / (overlapMat[2,1] * overlapMat[1,2])
oddsRatio
```
The odds-ratio got smaller, but it is still bigger than 1.

Here we basically says that each base can be assigned to peaks or promoters independently, which is definitely false.  So there are many, many reasons why this calculation is not the "right" one.  Nevertheless, it appears that we have some enrichment, as expected.

## SessionInfo

\scriptsize

```{r sessionInfo, echo=FALSE}
sessionInfo()
```

