---
title: Differential peak calling in forebrain/hindbrain of mice
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: show
    theme: readable
    df_print: paged
    fig_width: 8.5
    fig_height: 5
---

```{r set-publication-theme, echo=FALSE, cache=FALSE}
ggplot2::theme_set(ggpubr::theme_pubr(base_size=10))
```

```{r setup, echo = FALSE, cache = FALSE}
options(max.print="75")
library(knitr)
library(rmdformats)
knitr::opts_chunk$set(echo=FALSE,
                      tidy=TRUE,
                      prompt=FALSE,
                      cache=FALSE,
                      comment=NA,
                      message=FALSE,
                      warning=FALSE,
                      dev = c('png', 'cairo_pdf'),
                      fig.align = 'center', 
                      pdf.options(encoding = "ISOLatin9.enc"),
                      fig.path='figures/', warning=FALSE, message=FALSE)
opts_knit$set(width=75)
set.seed(12345)
```

# Overview

## Goal
This is an example peak affinity analysis, looking for differences in peak affinity in the hindbrain
and the forebrain of P0 mice, using two datasets from ENCODE 
[ENCSR312LQX](https://www.encodeproject.org/experiments/ENCSR312LQX) from the hindbrain and
[ENCSR310MLB](https://www.encodeproject.org/experiments/ENCSR310MLB/) from the forebrain. Each of
these datasets has a single replicate.

```{r load-libraries}
library(tidyverse)
library(DESeq2)
library(knitr)
library(kableExtra)
library(plyranges)
library(tidyverse)
library(stringr)
library(biobroom)
library(pheatmap)
library(ggrepel)
```


```{r load-metadata}
metadata = readr::read_csv("metadata/metadata.csv") %>%
  rename(X1="sample") %>%
  dplyr::select(-batch, -phenotype) %>%
  mutate(region=relevel(factor(region), ref="hindbrain")) %>%
  as.data.frame()
metadata = as.data.frame(metadata)
rownames(metadata) = metadata$sample
counts = readr::read_tsv("data/consensus-counts.tsv") %>%
  tibble::column_to_rownames("id")
```

```{r metadata-table}
metadata %>%
  kable() %>%
  kable_styling()
```

```{r get-symbols}
library(biomaRt)
ensembl = useMart("ensembl", dataset="mmusculus_gene_ensembl")
symbols = getBM(attributes=c("entrezgene_id", "mgi_symbol"), mart=ensembl) %>%
  rename(entrezgene_id="geneId") %>%
  mutate(geneId=as.character(geneId))
```

# Cleanup
We need to do a little it of cleanup before contining. We need to remove peaks
that appear in regions known to be false positive machines, called the blacklist
regions. We also are interested in regions that could be affecting expression
via chromatin accessability, so we'll only consider peaks within a window around
the transcription start sites of genes. 

## Blacklist region removal

```{r remove-blacklist}
expand_region_string = function(df, column) {
  df = as.data.frame(df)
  tokens = str_split_fixed(df[,column], ":", 2)
  chrom = tokens[,1]
  start = as.numeric(str_split_fixed(tokens[,2], "-", 2)[, 1])
  end = as.numeric(str_split_fixed(tokens[,2], "-", 2)[, 2])
  regions = data.frame(seqnames=chrom, start=start, end=end)
  df %>%
    bind_cols(regions)
}
blacklist = readr::read_tsv("metadata/mm10.blacklist.bed.gz", col_names=c("seqnames", "start", "end")) %>%
  as_granges()
blacklist_peaks = counts %>%
  tibble::rownames_to_column("peak") %>%
  expand_region_string("peak") %>%
  as_granges() %>%
  join_overlap_inner(blacklist) %>%
  tidy() %>%
  pull(peak)
counts = counts[!rownames(counts) %in% blacklist_peaks, ]
```

A few peaks overlap the ENCODE blacklist regions, so we removed them from the analysis. There 
were `r length(blacklist_peaks)` that overlapped with the blacklist regions.

## Open region annotation
Here we annotate the peaks that we called with contextual genomic information. We can also see that
some peaks are annotated in multiple regions, this can't be helped as genes overlap so there
are areas where multiple features of genes overlap.

```{r open-region-annotation}
library(ChIPseeker)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
annotatedobj = ChIPseeker::annotatePeak(counts %>%
                                        tibble::rownames_to_column("peak") %>%
                                        expand_region_string("peak") %>%
                                        as_granges(), TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene)
plotAnnoPie(annotatedobj)
upsetplot(annotatedobj)
```

## Remove peaks not near TSS
Generally we only care about peaks that are close to gene, here we are pretty lenient and consider
peaks that are within 1000 bases of a transcription start site of a gene.

```{r tss-peaks}
annotated = annotatedobj %>%
  as.GRanges() %>%
  tidy() 
ggplot(annotated, aes(distanceToTSS)) +
  stat_density(geom="line") +
  xlim(c(-10000,10000)) 
```

```{r keep-tss-peaks}
KEEP_DISTANCE = 1000
close_to_tss = annotated %>%
  filter(abs(distanceToTSS) < KEEP_DISTANCE) %>%
  pull(peak)
counts = counts[rownames(counts) %in% close_to_tss,]
```

Most peaks fall close to the TSS, but there are some that are very far away from the TSS. We'll 
remove peaks that aren't anywhere near the TSS for a gene from the analysis. We'll keep only
peaks within `r KEEP_DISTANCE` bases of a TSS. This leaves us with `r length(close_to_tss)` peaks to consider.

# PCA
We can see there is quite a bit of variation between these samples that is not associated with the **region**
the sample is from. **hindbrain_rep2** is clustered far away from the other samples, and has almost double
the input reads as the other samples. **forebrain_rep1** is also off on its own, and that sample was an
outlier for the percent of alignments in peaks, around 40% vs 20-30% for the other samples. 

```{r sample-clustering}
pc = prcomp(t(counts), scale=TRUE)
pcdev = data.frame(variation=100 * pc$sdev/sum(pc$sdev), component=paste0("PC", 1:length(pc$sdev)))
ggplot(pcdev, aes(component, variation)) +
  geom_bar(stat='identity') +
  xlab("") +
  ylab("% variation explained")
pcdat = pc$x %>%
  as.data.frame() %>%
  tibble::rownames_to_column("sample") %>%
  left_join(metadata, by="sample") 
ggplot(pcdat, aes(PC1, PC2, color=region, label=sample)) +
  geom_point() +
  geom_text_repel()
```
## mean-variance relationship
Here we can observe a mean-variance relationship in the data, similar to other count-based NGS
methods. At low counts we have a higher dispersion, meaning it is harder to estimate the true
counts at low count levels.

Since we have a small number of replicates, we borrow information from peaks close to the same
counts as other peaks to better estimate the variation of individual peaks. This plot shows 
the variation of individual peaks as the black dots, the red line fitting the mean-variance relationship
and the final result of weighting the individual peak variation with the overall mean-variance line,
which are the blue dots. 

```{r model-fit}
dds = DESeqDataSetFromMatrix(counts, metadata, design=~region)
dds = DESeq(dds)
plotDispEsts(dds)
CUTOFF = 0.1
```

# Differential affinity analysis

Here we look at the **region_forebrain_vs_hindbrain** coefficient, which will show us the peak 
affinity differences between the mouse forebrain and the mouse hindbrain samples.

```{r forebrain-hindbrain}
markup_deseq2 = function(res) {
  resanno = res %>%
    expand_region_string("peak") %>%
    as_granges()
  ChIPseeker::annotatePeak(resanno, TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene) %>%
    as.GRanges() %>%
    tidy() %>%
    left_join(symbols, by="geneId")
}
res = results(dds, name="region_forebrain_vs_hindbrain") %>%
  tidy() %>%
  rename(gene="peak")
dir.create("results")
readr::write_csv(res %>% markup_deseq2(), "results/forebrain-hindbrain.csv")
```

There are `r res %>% filter(p.adjusted < CUTOFF) %>% nrow()` peaks using an adjusted p-value cutoff of 
`r CUTOFF`. Here positive log2 fold changes are peaks with more counts in the **forebrain** than the
**hindbrain** samples.

#### MA-plot
```{r forebrain-hindbrain-plots-MA}
ggplot(res, aes(baseMean, estimate, color=p.adjusted < CUTOFF)) +
  geom_point(alpha=0.8, size=0.5) +
  scale_x_log10() +
  scale_color_manual(values=c("gray", "red")) +
  theme(legend.position='none') +
  xlab("mean reads per peak") +
  ylab(expression(log[2]~fold~change)) +
  ggtitle("forebrain vs hindbrain")
```

#### Volcano plot
```{r forebrain-hindbrain-plots-volcano}
ggplot(res, aes(estimate, -10*log10(p.value), color=p.adjusted < CUTOFF)) +
  geom_point(alpha=0.8, size=0.5) +
  scale_color_manual(values=c("gray", "red")) +
  theme(legend.position='none') +
  ylab(expression("-10"~log[10]~pvalue)) +
  xlab(expression(log[2]~fold~change)) +
  ggtitle("forebrain vs hindbrain")
```

#### Heatmap
```{r forebrain-hindbrain-plots-heatmap}
res_top25 = res %>%
  filter(p.adjusted < CUTOFF) %>%
  arrange(p.value) %>%
  head(25) %>%
  pull(peak)
ncounts = counts(dds, normalized=TRUE)
res_sig_counts = ncounts[rownames(ncounts) %in% res_top25,]
p = pheatmap(res_sig_counts, scale="row", annotation_col=metadata, main="forebrain vs hindbrain, top 25")
p
```

# Description of output files

- **start**: start of peak
- **end**: end of peak
- **width**: width of peak
- **strand**: strand of peak
- **seqname**: chromosome of peak
- **peak**: ID of peak
- **baseMean**: mean count of reads at peak for all samples
- **estimate**: log2 fold change of comparison
- **stderror**: standard error of the log2 fold change
- **statistic**: value of test statistic (Wald test)
- **p.value**: unadjusted p-value
- **p.adjusted**: adjusted p-value by BH correction
- **annotation**: proximity to nearest gene feature
- **geneChr**: chromosome of gene
- **geneStart**: start coordinate of gene
- **geneEnd**: end coordinate of gene
- **geneId**: Entrez gene ID of gene
- **transcriptId**: transcript ID
- **distanceToTSS**: distance to nearest transcription start site
- **mgi_symbol**: gene symbol from MGI
