05_chow_hfd_full.Rmd

---
author: "Ryan Dale"
output:
  html_document:
    keep_md: true
---

# R component of "genomics tools"


## Comparing H3K27ac data with differential expression


Load the H3K27ac data:

```{r}
setwd('tmp/data')
h3k27ac <- read.table('tsses_with_gained_h3k27ac.bed')
head(h3k27ac)
```

This copies over some stuff from the ggplot part to add up/down genes, in case
we didn't get to it by this point. Otherwise, we could load the table created
from that analysis.

```{r}
setwd('tmp/data')
df <- read.table('GSE77625.txt')
head(df)
valid <- !is.na(df$padj) & !is.na(df$log2FoldChange)
sig <- df$padj < 0.1
up <- df$log2FoldChange > 0 & valid & sig
dn <- df$log2FoldChange < 0 & valid & sig
df[up, 'direction'] <- 'up'
df[dn, 'direction'] <- 'dn'
df[!(up | dn), 'direction'] <- 'un'
head(df)
```

Similar to how we added a "direction" variable to the dataframe, let's add
a column indicating if a gene has H3K27ac at any of its TSSes. This uses the
``%in%`` operator. `x %in% y` gives a vector of TRUE or FALSE for every value
in `x` indicating if it is in `y` or not. We want to get a vector for every
gene in `df` which is TRUE where that gene's name is in the list of genes with
H3K27ac. Which column was that again?

```{r}
head(h3k27ac)
```
Ah, `V4`. So:

```{r}
has_h3k27ac <- rownames(df) %in% h3k27ac$V4
```

Spot checks . . . how many TRUE?

```{r}
table(has_h3k27ac)
```

Wait, how long was `h3k27ac`?

```{r}
nrow(h3k27ac)
```

Recall that we generated that H3K27ac data using trancript TSSes, so a gene can
show up multiple times. Since the DESeq2 results are at the gene level, we are
"collapsing" transcripts to genes by saying TRUE when *any* TSS of
a transcript has H3K27ac.

So we would hope that the number of *unique* genes in `h3k27ac$V4` is close to
the number of TRUE in `has_h3k27ac`:

```{r}
length(unique(h3k27ac$V4))
```

Great! Let's add the column to the dataframe so we can use it for plotting.

```{r}
df$h3k27ac <- has_h3k27ac
head(df)
```

### Plotting

Here's the MA plot from the ggplot example:

```{r}
ggplot(df) +
    geom_point(mapping=aes(x=log10(baseMean), y=log2FoldChange, color=direction)) +
    scale_color_manual(
        values=c('un'='#888888', 'up'='firebrick', 'dn'='dodgerblue4'),
        limits=c('up', 'dn', 'un')
    )
```

Now let's color by H3k27ac instead of direction.

```{r}
ggplot(df) +
    geom_point(mapping=aes(x=log10(baseMean), y=log2FoldChange, color=h3k27ac))
```

Hmm, it looks like the TRUE ones are being hidden. Here is how we can plot
*layers* in ggplot, to make sure the TRUE gets plotted *after* the FALSE. The
key is in the `data=subset(...)` part for the geom:

```{r}
ggplot(df) +
    geom_point(data=subset(df, !df$h3k27ac), mapping=aes(x=log10(baseMean), y=log2FoldChange)) +
    geom_point(data=subset(df, df$h3k27ac),  mapping=aes(x=log10(baseMean), y=log2FoldChange), color='red')
```

It might be nice to combine the previous plots together:

```{r}
ggplot(df) +
    geom_point(mapping=aes(x=log10(baseMean), y=log2FoldChange, color=direction)) +
    scale_color_manual(
        values=c('un'='#888888', 'up'='firebrick', 'dn'='dodgerblue4'),
        limits=c('up', 'dn', 'un')
    ) +
    geom_point(
      data=subset(df, df$h3k27ac),
      mapping=aes(x=log10(baseMean), y=log2FoldChange),
      color='orange', size=.5)
```

### Fisher's exact tests

Let's assign some statistics to this. Are upregulated genes enriched for
H3K27ac in at least one of their TSSes?

R comes with a built-in Fisher's exact test. However, it is very particular
about how it wants the data to be provided. There are two ways of doing it;
here is the more verbose way (which is more explicit). We build the table that
follows the pattern:

```
         with_condition   without_condition
selected      .                   .
not           .                   .

```

Here's how to do that. Note that we're doing boolean operations, and then
taking the sum. That tells us how many TRUE. For example, `x` is TRUE wherever
both `df$direction == 'up'` is TRUE *and* `df$h3k27ac` is TRUE. So it's TRUE
wherever a gene went up AND has H3K27ac.

```{r}
x <- df$direction == 'up' & df$h3k27ac
```

And then we can get the *number* of upregulated genes that have H2K27ac:

```{r}
sum(x)
```

Here is how to prepare the data for a Fisher's exact test:

```{r}
m <-matrix(
           c(
             sum(df$direction == 'up' & df$h3k27ac),
             sum(df$direction == 'up' & !df$h3k27ac),
             sum(df$direction != 'up' & df$h3k27ac),
             sum(df$direction != 'up' & !df$h3k27ac)
             ),
           nrow=2,
           dimnames=list(
                         c('up', 'not'),
                         c('h3k27ac', 'not')
                         )
           )
m
```

Now that we have the data prepared, all we have to do is:

```{r}
fisher.test(m)
```

The odds ratio is greater than 1, which means the top-left corner of the table
is higher than you'd expect, given the other cells of the table. This is very
highly significant.

**Question:** How would we check if *downregulated* genes are enriched?

So we can conclude the following:

    "Genes upregulated under high-fat diet were enriched for gained H3K27ac at
    their promoters (p<2.2e-16, odds ratio 2.0, Fisher's exact test).

We also have all the numbers to calculate the various percentages:

```{r}
m
rowSums(m)
colSums(m)
```

It turns out that the *downregulated* genes are *also* enriched for H3K27ac at
their promoters (though not quite as strongly). So it would be good to report
this as well.

**Importantly**, in the methods we would add the details about how we handled
genes and transcripts and how we defined "gain" of H3K27ac.

## Inspecting enhancer distances

This is an optional section in case we still have time.

### Load the closest transcript data

```{r}
setwd('tmp/data')
closest_transcripts <- read.table('closest_transcripts_to_enhancer_chow.bed')
```
Let's inspect what we brought in:

```{r}
head(closest_transcripts)
```

No header. Notice how R adds ``V*`` column names if there's no header. One
option would be to manually type in a header row. However, for now we just want
the distance (``V14``), and the gene symbol (``V9``).

```{r}
closest_transcripts <- closest_transcripts[, c('V9', 'V14')]
head(closest_transcripts)

colnames(closest_transcripts) <- c('gene', 'dist')
head(closest_transcripts)
```

### Basic plotting
- How far away are these enhancers from genes, on average?
 - What kind of plot is appropriate for this?

```{r}
library(ggplot2)
ggplot(closest_transcripts) +
    aes(x=dist) +
    geom_histogram()
```

Looks like there's something pretty far away. What gene is it?

```{r}
# order small to large. Recall we have negative distances, but the one we're
# trying to find is positive. So we want the last thing.
o <- order(closest_transcripts$dist)
tail(o)
closest_transcripts[tail(o),]
```

Hmm. Looks like ``Ttll7`` is the farthest. But why are there repeats of
``Bcor`` and ``Rnf144a``?

Remember the `transcripts.bed` file? It had all transcripts of all genes. If
multiple transcripts have the same TSS, *and* an enhancer is upstream of that
transcript (say 10kb away), then we will get multiple hits for that gene. So we
need to interpret these results as "number of **transcripts**.

Just to demonstrate, here's the same thing we just did, but now including the
Ensembl transcript ID so we can verify that those genes are repeated because
they are different transcripts:

```{r}
setwd('tmp/data')
demo <- read.table('closest_transcripts_to_enhancer_chow.bed')
demo <- demo[, c('V9', 'V12', 'V14')]
colnames(demo) <- c('gene', 'transcript', 'dist')
o <- order(demo$dist)
tail(o)
demo[tail(o),]
```