The bumphunter package can be used for methylation analyses where you are interested in identifying differentially methylated regions. The vignette explains in greater detail the data set we are using in this example.
## Load bumphunter
library('bumphunter')## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colnames,
## do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, lengths, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff,
## sort, table, tapply, union, unique, unsplit##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:base':
##
## colMeans, colSums, expand.grid, rowMeans, rowSums## Loading required package: IRanges## Loading required package: GenomeInfoDb## Loading required package: GenomicRanges## Loading required package: foreach## foreach: simple, scalable parallel programming from Revolution Analytics
## Use Revolution R for scalability, fault tolerance and more.
## http://www.revolutionanalytics.com## Loading required package: iterators## Loading required package: locfit## locfit 1.5-9.1 2013-03-22## 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)## [bumphunterEngine] Using a single core (backend: doSEQ, version: 1.4.3).## [bumphunterEngine] Computing coefficients.## [bumphunterEngine] Finding regions.## [bumphunterEngine] Found 15 bumps.## Explore data
lapply(tab, head)## $table
## chr start end value area cluster indexStart indexEnd L
## 10 chr1 2316 2631 -1.5814747 15.8147473 2 39 48 10
## 7 chr2 451 551 1.5891293 4.7673878 3 68 70 3
## 2 chr1 456 526 1.0678828 3.2036485 1 14 16 3
## 5 chr1 2176 2211 0.7841794 1.5683589 2 35 36 2
## 6 chr1 2841 2841 1.2010184 1.2010184 2 54 54 1
## 4 chr1 771 771 0.7780902 0.7780902 1 23 23 1
## clusterL
## 10 29
## 7 20
## 2 29
## 5 29
## 6 29
## 4 29
##
## $coef
## [,1]
## [1,] 0.60960932
## [2,] -0.09052769
## [3,] -0.21482638
## [4,] 0.13053755
## [5,] -0.21723642
## [6,] 0.39934961
##
## $fitted
## [,1]
## [1,] 0.60960932
## [2,] -0.09052769
## [3,] -0.21482638
## [4,] 0.13053755
## [5,] -0.21723642
## [6,] 0.39934961
##
## $pvaluesMarginal
## [1] NAOnce we have the regions we can proceed to build the required GRanges object.
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
data(hg19Ideogram, package = 'biovizBase')
seqlengths(regions) <- seqlengths(hg19Ideogram)[names(seqlengths(regions))]
## Explore the regions
regions## GRanges object with 15 ranges and 5 metadata columns:
## seqnames ranges strand | value area
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr1 [2316, 2631] * | -1.58147472668612 15.8147472668612
## [2] chr2 [ 451, 551] * | 1.58912926544627 4.76738779633882
## [3] chr1 [ 456, 526] * | 1.06788283103948 3.20364849311844
## [4] chr1 [2176, 2211] * | 0.784179442360784 1.56835888472157
## [5] chr1 [2841, 2841] * | 1.20101842543356 1.20101842543356
## ... ... ... ... . ... ...
## [11] chr1 [ 631, 631] * | 0.618602646982714 0.618602646982714
## [12] chr1 [ 1, 1] * | 0.609609318450113 0.609609318450113
## [13] chr1 [2911, 2911] * | -0.576423040796044 0.576423040796044
## [14] chr2 [ 251, 251] * | -0.556159535339166 0.556159535339166
## [15] chr1 [2806, 2806] * | -0.521605772372843 0.521605772372843
## cluster L clusterL
## <numeric> <numeric> <integer>
## [1] 2 10 29
## [2] 3 3 20
## [3] 1 3 29
## [4] 2 2 29
## [5] 2 1 29
## ... ... ... ...
## [11] 1 1 29
## [12] 1 1 29
## [13] 2 1 29
## [14] 3 1 20
## [15] 2 1 29
## -------
## seqinfo: 2 sequences from an unspecified genomeNow 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 because we don’t have a p-value variable.
## Load regionReport
library('regionReport')## Make it so that the report will be available as a vignette
original <- readLines(system.file('regionExploration', 'regionExploration.Rmd',
package = 'regionReport'))
vignetteInfo <- c(
'vignette: >',
' %\\VignetteEngine{knitr::rmarkdown}',
' %\\VignetteIndexEntry{Basic genomic regions exploration}',
' %\\VignetteEncoding{UTF-8}'
)
new <- c(original[1:12], vignetteInfo, original[13:length(original)])
writeLines(new, '/Users/lcollado/Dropbox/JHSPH/Code/regionReportSupp/bumphunter-example/regionReportBumphunter.Rmd')
## Now create the report
report <- renderReport(regions, 'Example bumphunter', pvalueVars = NULL,
densityVars = c('Area' = 'area', 'Value' = 'value',
'Cluster Length' = 'clusterL'), significantVar = NULL,
output = 'index', outdir = 'bumphunter-example',
template = '/Users/lcollado/Dropbox/JHSPH/Code/regionReportSupp/bumphunter-example/regionReportBumphunter.Rmd', device = 'png')## Warning: replacing previous import 'ggplot2::Position' by
## 'BiocGenerics::Position' when loading 'ggbio'## Writing 11 Bibtex entries ...## OK
## Results written to file 'bumphunter-example/index.bib'##
##
## processing file: index.Rmd##
|
| | 0%
|
|.. | 2%
## inline R code fragments
##
##
|
|... | 5%
## label: docSetup (with options)
## List of 3
## $ bootstrap.show.code : logi FALSE
## $ dev : symbol device
## $ bootstrap.show.message: logi FALSE
##
##
|
|..... | 7%
## ordinary text without R code
##
##
|
|...... | 9%
## label: setup (with options)
## List of 1
## $ bootstrap.show.message: logi FALSE
##
##
|
|........ | 12%
## ordinary text without R code
##
##
|
|......... | 14%
## label: pvaluePlots (with options)
## List of 3
## $ echo : logi FALSE
## $ results: chr "asis"
## $ eval : symbol hasPvalueVars
##
##
|
|........... | 16%
## ordinary text without R code
##
##
|
|............ | 19%
## label: regLen (with options)
## List of 4
## $ fig.width : num 14
## $ fig.height: num 14
## $ eval : symbol hasSignificant
## $ echo : symbol hasSignificant
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|............... | 23%
## label: regLen2 (with options)
## List of 4
## $ fig.width : num 10
## $ fig.height: num 10
## $ eval : language !hasSignificant
## $ echo : language !hasSignificant##
|
|................. | 26%
## inline R code fragments
##
##
|
|.................. | 28%
## label: densityPlots (with options)
## List of 3
## $ echo : logi FALSE
## $ results: chr "asis"
## $ eval : symbol hasDensityVars##
|
|.................... | 30%
## inline R code fragments
##
##
|
|..................... | 33%
## label: genomeOverview1 (with options)
## List of 6
## $ message : logi FALSE
## $ fig.width : num 7
## $ fig.height: num 9
## $ dpi : num 300
## $ eval : symbol hasSignificant
## $ echo : symbol hasSignificant
##
##
|
|....................... | 35%
## inline R code fragments
##
##
|
|........................ | 37%
## label: manhattanPlots (with options)
## List of 3
## $ echo : logi FALSE
## $ results: chr "asis"
## $ eval : symbol hasPvalueVars
##
##
|
|.......................... | 40%
## ordinary text without R code
##
##
|
|........................... | 42%
## label: genomeOverview2 (with options)
## List of 4
## $ message : logi FALSE
## $ fig.width : num 7
## $ fig.height: num 9
## $ dpi : num 300## No annotationPackage supplied. Trying org.Hs.eg.db.## Loading required package: org.Hs.eg.db## Loading required package: DBI## ## Getting TSS and TSE.## Getting CSS and CSE.## Getting exons.## Annotating genes.## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.##
|
|............................. | 44%
## inline R code fragments
##
##
|
|.............................. | 47%
## label: annoReg (with options)
## List of 1
## $ results: chr "asis"
##
##
|
|................................ | 49%
## ordinary text without R code
##
##
|
|................................. | 51%
## label: genomeOverview3 (with options)
## List of 6
## $ message : logi FALSE
## $ fig.width : num 7
## $ fig.height: num 9
## $ dpi : num 300
## $ eval : symbol hasSignificant
## $ echo : symbol hasSignificant
##
##
|
|................................... | 53%
## inline R code fragments
##
##
|
|.................................... | 56%
## label: countTable (with options)
## List of 1
## $ results: chr "asis"
##
##
|
|...................................... | 58%
## ordinary text without R code
##
##
|
|....................................... | 60%
## label: vennDiagram (with options)
## List of 2
## $ fig.width : num 7
## $ fig.height: num 7##
|
|......................................... | 63%
## inline R code fragments
##
##
|
|.......................................... | 65%
## label: vennDiagramSignificant (with options)
## List of 4
## $ eval : symbol hasSignificant
## $ echo : symbol hasSignificant
## $ fig.width : num 7
## $ fig.height: num 7
##
##
|
|............................................ | 67%
## inline R code fragments
##
##
|
|............................................. | 70%
## label: bestRegionInfo (with options)
## List of 1
## $ results: chr "asis"
##
##
|
|............................................... | 72%
## ordinary text without R code
##
##
|
|................................................ | 74%
## label: unnamed-chunk-1 (with options)
## List of 2
## $ child: symbol customCode
## $ eval : symbol hasCustomCode
##
##
|
|.................................................. | 77%
## inline R code fragments
##
##
|
|................................................... | 79%
## label: thecall (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|..................................................... | 81%
## ordinary text without R code
##
##
|
|...................................................... | 84%
## label: reproducibility1 (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|........................................................ | 86%
## ordinary text without R code
##
##
|
|......................................................... | 88%
## label: reproducibility2 (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|........................................................... | 91%
## ordinary text without R code
##
##
|
|............................................................ | 93%
## label: reproducibility3 (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............................................................. | 95%
## inline R code fragments
##
##
|
|............................................................... | 98%
## label: bibliography (with options)
## List of 3
## $ results: chr "asis"
## $ echo : logi FALSE
## $ warning: logi FALSE
##
##
|
|.................................................................| 100%
## ordinary text without R code## output file: index.knit.md## /usr/local/bin/pandoc +RTS -K512m -RTS index.utf8.md --to html --from markdown+autolink_bare_uris+ascii_identifiers+tex_math_single_backslash-implicit_figures --output index.html --smart --email-obfuscation none --self-contained --standalone --section-divs --table-of-contents --toc-depth 3 --variable toc_float=1 --variable toc_selectors=h1,h2,h3 --variable toc_collapsed=1 --variable toc_smooth_scroll=1 --template /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rmarkdown/rmd/h/default.html --variable 'theme:spacelab' --include-in-header /var/folders/cx/n9s558kx6fb7jf5z_pgszgb80000gn/T//Rtmp41dyxK/rmarkdown-str3a9461833679.html --mathjax --variable 'mathjax-url:https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --no-highlight --variable highlightjs=/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rmarkdown/rmd/h/highlight --variable navigationjs=/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rmarkdown/rmd/h/navigation-1.0 --variable code_folding=hide##
## Output created: index.html## Clean up
file.remove('/Users/lcollado/Dropbox/JHSPH/Code/regionReportSupp/bumphunter-example/regionReportBumphunter.Rmd')## [1] TRUEYou can view the final report here.
## Date generated:
Sys.time()## [1] "2016-04-12 07:36:20 EDT"## Time spent making this page:
proc.time()## user system elapsed
## 172.555 7.336 200.023## R and packages info:
options(width = 120)
devtools::session_info()## Session info -----------------------------------------------------------------------------------------------------------## setting value
## version R version 3.3.0 alpha (2016-03-23 r70368)
## system x86_64, darwin13.4.0
## ui AQUA
## language (EN)
## collate en_US.UTF-8
## tz America/New_York
## date 2016-04-12## Packages ---------------------------------------------------------------------------------------------------------------## package * version date source
## acepack 1.3-3.3 2014-11-24 CRAN (R 3.3.0)
## annotate 1.49.1 2016-02-06 Bioconductor
## AnnotationDbi * 1.33.8 2016-04-10 Bioconductor
## AnnotationHub 2.3.16 2016-03-25 Bioconductor
## backports 1.0.2 2016-03-18 CRAN (R 3.3.0)
## bibtex 0.4.0 2014-12-31 CRAN (R 3.3.0)
## Biobase * 2.31.3 2016-01-14 Bioconductor
## BiocGenerics * 0.17.4 2016-04-07 Bioconductor
## BiocInstaller 1.21.4 2016-03-23 Bioconductor
## BiocParallel 1.5.21 2016-03-23 Bioconductor
## BiocStyle * 1.99.0 2016-04-05 Bioconductor
## biomaRt 2.27.2 2016-01-14 Bioconductor
## Biostrings 2.39.12 2016-02-21 Bioconductor
## biovizBase 1.19.6 2016-04-06 Bioconductor
## bitops 1.0-6 2013-08-17 CRAN (R 3.3.0)
## BSgenome 1.39.4 2016-02-21 Bioconductor
## bumphunter * 1.11.5 2016-03-29 Bioconductor
## checkmate 1.7.4 2016-04-08 CRAN (R 3.3.0)
## cluster 2.0.3 2015-07-21 CRAN (R 3.3.0)
## codetools 0.2-14 2015-07-15 CRAN (R 3.3.0)
## colorspace 1.2-6 2015-03-11 CRAN (R 3.3.0)
## DBI * 0.3.1 2014-09-24 CRAN (R 3.3.0)
## DEFormats 0.99.8 2016-03-31 Bioconductor
## derfinder * 1.5.30 2016-03-25 Bioconductor
## derfinderHelper 1.5.3 2016-03-23 Bioconductor
## derfinderPlot * 1.5.7 2016-03-23 Bioconductor
## DESeq2 1.11.42 2016-04-10 Bioconductor
## devtools * 1.10.0 2016-01-23 CRAN (R 3.3.0)
## dichromat 2.0-0 2013-01-24 CRAN (R 3.3.0)
## digest 0.6.9 2016-01-08 CRAN (R 3.3.0)
## doRNG 1.6 2014-03-07 CRAN (R 3.3.0)
## DT * 0.1 2015-06-09 CRAN (R 3.3.0)
## edgeR 3.13.8 2016-04-08 Bioconductor
## ensembldb 1.3.19 2016-04-03 Bioconductor
## evaluate 0.8.3 2016-03-05 CRAN (R 3.3.0)
## foreach * 1.4.3 2015-10-13 CRAN (R 3.3.0)
## foreign 0.8-66 2015-08-19 CRAN (R 3.3.0)
## formatR 1.3 2016-03-05 CRAN (R 3.3.0)
## Formula 1.2-1 2015-04-07 CRAN (R 3.3.0)
## genefilter 1.53.3 2016-03-23 Bioconductor
## geneplotter 1.49.0 2016-01-14 Bioconductor
## GenomeInfoDb * 1.7.6 2016-01-29 Bioconductor
## GenomicAlignments 1.7.20 2016-02-25 Bioconductor
## GenomicFeatures * 1.23.29 2016-04-05 Bioconductor
## GenomicFiles 1.7.9 2016-02-22 Bioconductor
## GenomicRanges * 1.23.25 2016-03-31 Bioconductor
## GGally 1.0.1 2016-01-14 CRAN (R 3.3.0)
## ggbio * 1.19.13 2016-04-03 Bioconductor
## ggplot2 * 2.1.0 2016-03-01 CRAN (R 3.3.0)
## graph 1.49.1 2016-01-14 Bioconductor
## gridExtra * 2.2.1 2016-02-29 CRAN (R 3.3.0)
## gtable 0.2.0 2016-02-26 CRAN (R 3.3.0)
## highr 0.5.1 2015-09-18 CRAN (R 3.3.0)
## Hmisc 3.17-3 2016-04-03 CRAN (R 3.3.0)
## htmltools 0.3.5 2016-03-21 CRAN (R 3.3.0)
## htmlwidgets 0.6 2016-02-25 CRAN (R 3.3.0)
## httpuv 1.3.3 2015-08-04 CRAN (R 3.3.0)
## httr 1.1.0 2016-01-28 CRAN (R 3.3.0)
## interactiveDisplayBase 1.9.0 2016-01-14 Bioconductor
## IRanges * 2.5.43 2016-04-10 Bioconductor
## iterators * 1.0.8 2015-10-13 CRAN (R 3.3.0)
## jsonlite 0.9.19 2015-11-28 CRAN (R 3.3.0)
## knitcitations 1.0.7 2015-10-28 CRAN (R 3.3.0)
## knitr * 1.12.3 2016-01-22 CRAN (R 3.3.0)
## knitrBootstrap 1.0.0 2016-03-24 Github (jimhester/knitrBootstrap@cdaa4a9)
## labeling 0.3 2014-08-23 CRAN (R 3.3.0)
## lattice 0.20-33 2015-07-14 CRAN (R 3.3.0)
## latticeExtra 0.6-28 2016-02-09 CRAN (R 3.3.0)
## limma 3.27.14 2016-03-23 Bioconductor
## locfit * 1.5-9.1 2013-04-20 CRAN (R 3.3.0)
## lubridate 1.5.6 2016-04-06 CRAN (R 3.3.0)
## magrittr 1.5 2014-11-22 CRAN (R 3.3.0)
## markdown 0.7.7 2015-04-22 CRAN (R 3.3.0)
## Matrix 1.2-4 2016-03-02 CRAN (R 3.3.0)
## matrixStats 0.50.1 2015-12-15 CRAN (R 3.3.0)
## memoise 1.0.0 2016-01-29 CRAN (R 3.3.0)
## mgcv * 1.8-12 2016-03-03 CRAN (R 3.3.0)
## mime 0.4 2015-09-03 CRAN (R 3.3.0)
## munsell 0.4.3 2016-02-13 CRAN (R 3.3.0)
## nlme * 3.1-126 2016-03-14 CRAN (R 3.3.0)
## nnet 7.3-12 2016-02-02 CRAN (R 3.3.0)
## org.Hs.eg.db * 3.3.0 2016-04-11 Bioconductor
## OrganismDbi 1.13.6 2016-04-05 Bioconductor
## pkgmaker 0.22 2014-05-14 CRAN (R 3.3.0)
## plyr 1.8.3 2015-06-12 CRAN (R 3.3.0)
## qvalue 2.3.2 2016-01-14 Bioconductor
## R6 2.1.2 2016-01-26 CRAN (R 3.3.0)
## RBGL 1.47.0 2016-01-14 Bioconductor
## RColorBrewer * 1.1-2 2014-12-07 CRAN (R 3.3.0)
## Rcpp 0.12.4 2016-03-26 CRAN (R 3.3.0)
## RCurl 1.95-4.8 2016-03-01 CRAN (R 3.3.0)
## RefManageR 0.10.13 2016-04-04 CRAN (R 3.3.0)
## regionReport * 1.5.47 2016-04-12 Bioconductor
## registry 0.3 2015-07-08 CRAN (R 3.3.0)
## reshape 0.8.5 2014-04-23 CRAN (R 3.3.0)
## reshape2 1.4.1 2014-12-06 CRAN (R 3.3.0)
## RJSONIO 1.3-0 2014-07-28 CRAN (R 3.3.0)
## rmarkdown 0.9.5 2016-02-22 CRAN (R 3.3.0)
## rngtools 1.2.4 2014-03-06 CRAN (R 3.3.0)
## rpart 4.1-10 2015-06-29 CRAN (R 3.3.0)
## Rsamtools 1.23.8 2016-04-10 Bioconductor
## RSQLite * 1.0.0 2014-10-25 CRAN (R 3.3.0)
## rstudioapi 0.5 2016-01-24 CRAN (R 3.3.0)
## rtracklayer 1.31.10 2016-04-07 Bioconductor
## S4Vectors * 0.9.46 2016-04-07 Bioconductor
## scales 0.4.0 2016-02-26 CRAN (R 3.3.0)
## shiny 0.13.2 2016-03-28 CRAN (R 3.3.0)
## stringi 1.0-1 2015-10-22 CRAN (R 3.3.0)
## stringr 1.0.0 2015-04-30 CRAN (R 3.3.0)
## SummarizedExperiment 1.1.23 2016-04-06 Bioconductor
## survival 2.38-3 2015-07-02 CRAN (R 3.3.0)
## TxDb.Hsapiens.UCSC.hg19.knownGene * 3.2.2 2016-03-24 Bioconductor
## VariantAnnotation 1.17.23 2016-04-07 Bioconductor
## whisker * 0.3-2 2013-04-28 CRAN (R 3.3.0)
## XML 3.98-1.4 2016-03-01 CRAN (R 3.3.0)
## xtable 1.8-2 2016-02-05 CRAN (R 3.3.0)
## XVector 0.11.8 2016-04-06 Bioconductor
## yaml 2.1.13 2014-06-12 CRAN (R 3.3.0)
## zlibbioc 1.17.1 2016-03-19 Bioconductor