---
title: "Profile Plots & Bootstrapping"
package: BRGenomics
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
  BiocStyle::pdf_document:
    toc: true
vignette: |
  %\VignetteIndexEntry{Profile Plots and Bootstrapping}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

# Conventional Profiles

As an example, we'll plot PRO-seq signal in promoter-proximal regions by 
calculating the mean signal intensity at each base in the first 100 bases of a 
transcript.

First, calculate signal at each base for all promoter-proximal regions:

```{r, message = FALSE}
library(BRGenomics)
data("PROseq")
data("txs_dm6_chr4")
txs_pr <- promoters(txs_dm6_chr4, 0, 100)
```

```{r, collapse = TRUE}
countmatrix_pr <- getCountsByPositions(PROseq, txs_pr, ncores = 1)
dim(countmatrix_pr)
dim(countmatrix_pr) == c(length(txs_pr), unique(width(txs_pr)))
```

For each position (each column of the matrix), calculate the mean, and plot:

```{r}
plot(x = 1:ncol(countmatrix_pr),
     y = colMeans(countmatrix_pr),
     type = "l", 
     xlab = "Distance to TSS (bp)",
     ylab = "Mean PRO-seq Reads")
```

\ 

## Drawbacks

One drawback of using the arithmetic mean is that means are not robust to 
outliers. In other words, the mean signal at any position is liable to be 
determined by a small number of highly influential points. This is especially 
problematic for high dynamic range data like PRO-seq.

It's common to see this issue addressed through the use medians/quantiles in 
place of arithmetic means. However, observe what happens as we plot the median 
signal across our gene list using several different gene-filtering thresholds:

```{r, fig.height=7, fig.width=6}
plot_meds <- function(sig_thresh) {
  idx <- which(rowSums(countmatrix_pr) > sig_thresh)
  plot(x = 1:ncol(countmatrix_pr),
       y = apply(countmatrix_pr[idx, ], 2, median),
       type = "l", 
       main = sprintf("Regions with >%s reads", sig_thresh),
       xlab = "Distance to TSS (bp)",
       ylab = "Median PRO-seq Reads")
}

par(mfrow = c(3, 2))

for (i in c(0, 30*2^(0:4))) {
  plot_meds(i)
}
```

The paucity of transcription on Drosophila chromosome 4 makes this a somewhat 
extreme example, and you might find that the above situation is a non-issue for 
your data. However, it's important to keep in mind how unexpressed genes can 
affect the median. A common rule of thumb is that a given cell is only apt to 
express around half of its genes; if that holds true, the median of an 
unfiltered genelist will be close to zero, and could be insensitive to changes 
occurring at expressed genes. 

# The Rationale for Bootstrapping 

A robust alternative to plotting mean or median signal profiles is to plot 
__bootstrapped__ mean signal profiles, known as metaprofile plots or metaplots. 

To bootstrap the mean signal by position, a small number of genes are randomly 
sampled from a genelist, and the mean signal at each position is calculated for 
that group of genes. This process of randomly sampling genes and calculating 
mean signals is repeated over many iterations, and the final bootstrapped mean 
for each position is the median of the sampled means.

One feature of bootstrapping is that it is robust to outliers. But more than 
that, the a bootstrapped mean provides an _expectation_ for what the mean signal 
would be for any arbitrary group of genes, and associated with that expectation 
is a measure of uncertainty. By taking quantiles of the sub-sampled means _other 
than_ the median (the expected value), we can estimate the extent to which the 
mean varies across arbitrary groups of genes. 

For example, the 75th percentile of the sub-sampled means is a number for which, 
25% of the time, we calculated a higher mean. Similarly, a 90% confidence 
interval about the bootstrapped mean encompasses all values between the 5th and 
95th percentiles of the sub-sampled means.

# Generating and Plotting Metaprofiles

We can use the `metaSubsample()` function to bootstrap mean values by position 
for a genelist. The function accepts the same arguments as 
`getCountsByPositions()`, in addition to other arguments related to the 
bootstrapping. 

By default, 10% of the genelist is randomly sampled 1000 times, and confidence 
bands are returned for the 12.5th and 87.5th percentiles (i.e. a 75% confidence 
interval). 

Given the small size our dataset, we'll reduce this to a 30% confidence 
interval, and we'll additionally use 5 bp bins:

```{r}
bootmeans.df <- metaSubsample(PROseq, txs_pr, binsize = 5, 
                              lower = 0.35, upper = 0.65, ncores = 1)
head(bootmeans.df)
```

A dataframe is returned for plotting, and notice how the x-values have been 
automatically adjusted to be the center of the bins.

Below, we show how to plot the confidence bands using base R plotting, as well 
as `ggplot2`.

```{r}
plot(mean ~ x, data = bootmeans.df, type = "l", 
     main = "PRO-seq Signal", ylim = c(0, 1.4),
     xlab = "Distance from TSS",
     ylab = "Mean Signal + 30% CI")

# draw a polygon to add confidence bands,
# and use adjustcolor() to add transparency
polygon(c(bootmeans.df$x, rev(bootmeans.df$x)), 
        c(bootmeans.df$lower, rev(bootmeans.df$upper)),
        col = adjustcolor("black", 0.1), border = FALSE)
```

```{r}
require(ggplot2)
ggplot(bootmeans.df, aes(x, mean)) + 
  geom_line() + 
  geom_ribbon(aes(x, ymin = lower, ymax = upper),
              alpha = 0.1) + 
  labs(title = "PRO-seq Signal", 
       x = "Distance from TSS", 
       y = "Mean Signal + 30% CI") + 
  theme_bw()
```

## Example: Comparative Metaplots

Like other functions in BRGenomics, we can pass a list of GRanges to 
`metaSubsample()`, and the output is conveniently combined for plotting.

```{r}
# make 3 datasets
ps_list <- list(ps1 = PROseq[seq(1, length(PROseq), 3)], 
                ps2 = PROseq[seq(2, length(PROseq), 3)], 
                ps3 = PROseq[seq(3, length(PROseq), 3)])
```

```{r}
bm_list.df <- metaSubsample(ps_list, txs_pr, binsize = 5, 
                            lower = 0.35, upper = 0.65, ncores = 1)
head(bm_list.df)
```

```{r}
require(ggplot2)
ggplot(bm_list.df, aes(x, mean, color = sample.name)) + 
  geom_line() + 
  geom_ribbon(aes(x, ymin = lower, ymax = upper,
                  color = NULL, fill = sample.name),
              alpha = 0.2) + 
  labs(title = "PRO-seq Signal", 
       x = "Distance from TSS", 
       y = "Mean Signal + 30% CI") + 
  theme_bw()
```

_Bear in mind that the above confidence intervals are ignoring 70% of the 
sub-sampling experiments, which is excessive. More reasonable parameters would 
reveal how lacking this data is, i.e. how un-confident we are that there is a 
robust difference in the mean signal at various positions._ 
