---
title: "Interval statistics"
author: "Jay Hesselberth"
date: '`r format(Sys.Date(), "%B %d %Y")`'
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
    vignette: >
      %\VignetteIndexEntry{Interval Statistics}
      %\VignetteEngine{knitr::rmarkdown}
      %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "img/interval-stats-",
  fig.height = 4,
  fig.align = "center",
  fig.width = 4
)
```

## Overview

`valr` can be used to explore relationships between sets of intervals. Here we explore the relationship between transcription start sites and repetitive elements in the human genome. 

```{r load-data, message = FALSE, warning = FALSE}
library(valr)
library(dplyr)
library(ggplot2)
library(tidyr)

# load repeats and genes. Data in the valr package is restricted to chr22; the entire
# files can be downloaded from UCSC.
rpts <- read_bed(valr_example('hg19.rmsk.chr22.bed.gz'), n_fields = 6) 
genes <- read_bed(valr_example('hg19.refGene.chr22.bed.gz'), n_fields = 12)

# load chrom sizes
genome <- read_genome(valr_example('hg19.chrom.sizes.gz'))

# create 1bp intervals representing transcription start sites
tss <-
  mutate(genes,
         .start = ifelse(strand == '+', start, end),
         .end = .start + 1) %>%
  select(chrom, start = .start, end = .end, name)

tss
```

## Distance metrics

First we define a function that takes `x` and `y` intervals and computes distance statistics (using `bed_reldist()` and `bed_absdist()`) for specified groups. The value of each statistic is assigned to a `.value` column.

```{r stats}
distance_stats <- function(x, y, genome, grp, type = NA) {
  group_by_(x, .dots = grp) %>%
    do(
       reldist = bed_reldist(., y, detail = TRUE) %>%
         select(.value = .reldist),
       absdist = bed_absdist(., y, genome) %>%
         select(.value = .absdist)
       ) %>%
    tidyr::gather_('stat', 'value', setdiff(names(.), list(grp))) %>%
    mutate(type = type)
} 
```

We use the `distance_stats()` function to apply the `bed_absdist()` function to each group of data.

```{r compute_obs}
obs_stats <- distance_stats(rpts, tss, genome, 'name', 'obs')
obs_stats
```

And the same is done for a set of  shuffled group of data. `bed_shuffle()` is used to shuffle coordinates of the repeats within each chromosome (i.e., the coordinates change, but the chromosome stays the same.)

```{r compute_shf}
shfs <- bed_shuffle(rpts, genome, within = TRUE)
shf_stats <- distance_stats(shfs, tss, genome, 'name', 'shuf')
```

Now we can bind the observed and shuffled data together, and do some tidying to put the data into a format appropriate for a statistical test. This involves:

1. `unnest()`ing the data frames
1. creating groups for each repeat (`name`), stat (`reldist` or `absdist`) and type (`obs` or `shf`)
1. adding unique surrogate row numbers for each group
1. using `tidyr::spread()` to create two new `obs` and `shuf` columns
1. removing rows with `NA` values.

```{r bind_res}
res <- bind_rows(obs_stats, shf_stats) %>%
  tidyr::unnest(value) %>% 
  group_by(name, stat, type) %>%
  mutate(.id = row_number()) %>%
  tidyr::spread(type, .value) %>%
  na.omit()

res
```

Now that the data are formatted, we can use `ks.test()` to determine whether there are significant differences between the observed and shuffled data for each group. `broom::tidy()` is used to reformat the results of each test into a `data_frame`, and the results of each test are `gather`ed to into a `type` column for each test type.

```{r pvalues}
library(broom)

# silence the tests with ties ...
ks.test.quiet <- function(...) {suppressWarnings(ks.test(...))}

pvals <- res %>% do(twosided = broom::tidy(ks.test.quiet(.$obs, .$shuf)),
                    less = broom::tidy(ks.test.quiet(.$obs, .$shuf, alternative = 'less')),
                    greater = broom::tidy(ks.test.quiet(.$obs, .$shuf, alternative = 'greater'))) %>%
  tidyr::gather(alt, type, -name, -stat) %>%
  unnest(type) %>%
  select(name:p.value) %>%
  arrange(p.value)
```

Histgrams of the different stats help visulaize the distribution of p.values.

```{r pvalue_viz}
ggplot(pvals, aes(p.value)) +
  geom_histogram(binwidth = 0.05) +
  facet_grid(stat ~ alt) + theme_bw()
```

We can also assess false discovery rates (q.values) using `p.adjust()`.

```{r qvalues}
pvals <-
  group_by(pvals, stat, alt) %>%
  mutate(q.value = p.adjust(p.value)) %>%
  ungroup() %>%
  arrange(q.value)
```

Finally we can visualize these results using `stat_ecdf()`.

```{r ecfs}
res_gather <- tidyr::gather(res, type, value, -name, -stat, -.id)

signif <- head(pvals, 5)

res_signif <-
  signif %>%
  left_join(res_gather, by = c('name', 'stat'))

ggplot(res_signif, aes(x = value, color = type)) +
  stat_ecdf() + 
  facet_grid(stat ~ name) + theme_classic() + scale_x_log10()
```

## Projection test

`bed_projection()` is a statistical approach to assess the relationship between two intervals based on the  binomial distribution. Here, we examine the distribution of repetitive elements within the promoters of coding or non-coding genes.

First we'll extract 5 kb regions upstream of the transcription start sites to represent the promoter regions for coding and non-coding genes. 

```{r get_promoters}
# create intervals 5kb upstream of tss representing promoters
promoters <-
  bed_flank(genes, genome, left = 5000, strand = TRUE) %>%
  mutate(name = ifelse(grepl('NR_', name), 'non-coding', 'coding')) %>%
  select(chrom:strand)
  
# select coding and non-coding promoters
promoters_coding <- filter(promoters, name == 'coding')
promoters_ncoding <- filter(promoters, name == 'non-coding')

promoters_coding

promoters_ncoding
```

Next we'll apply the `bed_projection()` test for each repeat class for both coding and non-coding regions. 

```{r get_projections}

# function to apply bed_projection to groups
projection_stats <- function(x, y, genome, grp, type = NA) {
  group_by_(x, .dots = grp) %>%
    do(n_repeats = nrow(.),
       projection = bed_projection(., y, genome)) %>%
    mutate(type = type)
} 

pvals_coding <- projection_stats(rpts, promoters_coding, genome, 'name', 'coding')
pvals_ncoding <- projection_stats(rpts, promoters_ncoding, genome, 'name', 'non_coding')

pvals <-
  bind_rows(pvals_ncoding, pvals_coding) %>% 
  tidyr::unnest() %>%
  select(-chrom)

# filter for repeat classes with at least 10 intervals
pvals <- filter(pvals, n_repeats > 10,
                obs_exp_ratio != 0)

# adjust pvalues 
pvals <- mutate(pvals, q.value = p.adjust(p.value))

pvals
```

The projection test is a two-tailed statistical test. A significant p-value indicates either enrichment or depletion of query intervals compared to the reference interval sets. A value of `lower_tail = TRUE` column indicates that the query intervals are depleted, whereas `lower_tail = FALSE` indicates that the query intervals are enriched. 

```{r table}
library(knitr)

# find and show top 5 most significant repeats
signif_tests <-
  pvals %>% 
  arrange(q.value) %>%
  group_by(type) %>%
  top_n(-5, q.value) %>%
  arrange(type)

knitr::kable(signif_tests)
```
