---
title: "Example report using bumphunter results"
author: 
  - name: Leonardo Collado-Torres
    affiliation:
    - &libd Lieber Institute for Brain Development, Johns Hopkins Medical Campus
    - &ccb Center for Computational Biology, Johns Hopkins University
    email: lcolladotor@gmail.com
output: 
  BiocStyle::html_document:
    self_contained: yes
    toc: true
    toc_float: true
    toc_depth: 2
    code_folding: show
date: "`r doc_date()`"
package: "`r pkg_ver('regionReport')`"
vignette: >
  %\VignetteIndexEntry{Example report using bumphunter results}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}  
---

`r Biocpkg('bumphunter')` example
====================

The `r Biocpkg('bumphunter')` package can be used for methylation analyses where you are interested in identifying differentially methylated regions. The [vignette](http://bioconductor.org/packages/release/bioc/vignettes/bumphunter/inst/doc/bumphunter.pdf) explains in greater detail the data set we are using in this example.

```{r 'findRegions'}
## Load bumphunter
library("bumphunter")

## Create data from the vignette
pos <- list(
    pos1 = seq(1, 1000, 35),
    pos2 = seq(2001, 3000, 35),
    pos3 = seq(1, 1000, 50)
)
chr <- rep(paste0("chr", c(1, 1, 2)), times = sapply(pos, length))
pos <- unlist(pos, use.names = FALSE)

## Find clusters
cl <- clusterMaker(chr, pos, maxGap = 300)

## Build simulated bumps
Indexes <- split(seq_along(cl), cl)
beta1 <- rep(0, length(pos))
for (i in seq(along = Indexes)) {
    ind <- Indexes[[i]]
    x <- pos[ind]
    z <- scale(x, median(x), max(x) / 12)
    beta1[ind] <- i * (-1)^(i + 1) * pmax(1 - abs(z)^3, 0)^3 ## multiply by i to vary size
}

## Build data
beta0 <- 3 * sin(2 * pi * pos / 720)
X <- cbind(rep(1, 20), rep(c(0, 1), each = 10))
set.seed(23852577)
error <- matrix(rnorm(20 * length(beta1), 0, 1), ncol = 20)
y <- t(X[, 1]) %x% beta0 + t(X[, 2]) %x% beta1 + error

## Perform bumphunting
tab <- bumphunter(y, X, chr, pos, cl, cutoff = .5)

## Explore data
lapply(tab, head)
```

Once we have the regions we can proceed to build the required `GRanges` object.

```{r 'buildGRanges'}
library("GenomicRanges")

## Build GRanges with sequence lengths
regions <- GRanges(
    seqnames = tab$table$chr,
    IRanges(start = tab$table$start, end = tab$table$end),
    strand = "*", value = tab$table$value, area = tab$table$area,
    cluster = tab$table$cluster, L = tab$table$L, clusterL = tab$table$clusterL
)

## Assign chr lengths
seqlengths(regions) <- seqlengths(
    getChromInfoFromUCSC("hg19", as.Seqinfo = TRUE)
)[
    names(seqlengths(regions))
]

## Explore the regions
regions
```

Now that we have identified a set of differentially methylated regions we can proceed to creating the HTML report. Note that this report has less information than the [DiffBind example](http://leekgroup.github.io/regionReportSupp/DiffBind.html) because we don't have a p-value variable.


```{r 'loadLib'}
## Load regionReport
library("regionReport")
```

```{r 'createReport', eval = FALSE}
## Now create the report
report <- renderReport(regions, "Example bumphunter",
    pvalueVars = NULL,
    densityVars = c(
        "Area" = "area", "Value" = "value",
        "Cluster Length" = "clusterL"
    ), significantVar = NULL,
    output = "bumphunter-example", outdir = "bumphunter-example",
    device = "png"
)
```

You can view the final report at `bumphunter-example/bumphunter-example.html` [or here](http://leekgroup.github.io/regionReport/reference/bumphunter-example/bumphunter-example.html).

In case the link does not work, a [pre-compiled version of this document](http://leekgroup.github.io/regionReportSupp/bumphunterExample.html) and its [corresponding report](http://leekgroup.github.io/regionReportSupp/bumphunter-example/index.html) are available at [leekgroup.github.io/regionReportSupp/](http://leekgroup.github.io/regionReportSupp/index.html).

# Reproducibility

```{r 'reproducibility'}
## Date generated:
Sys.time()

## Time spent making this page:
proc.time()

## R and packages info:
options(width = 120)
sessioninfo::session_info()
```