recountmethylation_users_guide.R
In recountmethylation: Access and Analyze DNA Methylation Array Databases

## ----setup, echo = FALSE------------------------------------------------------
suppressMessages(library(knitr))
suppressMessages(library(GenomicRanges))
suppressMessages(library(limma))
suppressMessages(library(minfi))
suppressMessages(library(ExperimentHub))
knitr::opts_chunk$set(echo = TRUE, eval = TRUE, warning = FALSE, 
  message = FALSE)

## ----disclaimer, echo = FALSE, message = TRUE---------------------------------
library(recountmethylation)

## ---- echo = FALSE------------------------------------------------------------
url = "https://recount.bio/data/"
sslver <- FALSE
ftpuseopt <- FALSE
dirlistopt <- FALSE
dn <- RCurl::getURL(url, ftp.use.epsv = ftpuseopt, dirlistonly = dirlistopt,
                    .opts = list(ssl.verifypeer = sslver))
sm <- as.data.frame(servermatrix(dn, sslver), stringsAsFactors = FALSE)
sdf <- as.data.frame(sm)
tsv <- as.numeric(gsub("(.*_|\\.h5)", "", sdf[,1]))
sdff <- sdf[which(tsv == max(tsv, na.rm = TRUE)),]
rownames(sdff) <- NULL
knitr::kable(sdff, align = "c")

## -----------------------------------------------------------------------------
hub = ExperimentHub::ExperimentHub() # connect to the hubs
rmdat <- AnnotationHub::query(hub, "recountmethylation") # query the hubs
rmdat

## -----------------------------------------------------------------------------
eid <- "EH3778" # h5 test file id
fpath <- rmdat[[eid]] # download with default caching
rhdf5::h5ls(fpath)

## -----------------------------------------------------------------------------
dn <- "remethdb-h5se_gr-test_0-0-1_1590090412"
path <- system.file("extdata", dn, package = "recountmethylation")
h5se.test <- HDF5Array::loadHDF5SummarizedExperiment(path)

## -----------------------------------------------------------------------------
class(h5se.test)

## -----------------------------------------------------------------------------
dim(h5se.test)

## -----------------------------------------------------------------------------
summary(h5se.test)

## -----------------------------------------------------------------------------
h5se.md <- minfi::pData(h5se.test)
dim(h5se.md)

## -----------------------------------------------------------------------------
colnames(h5se.md)

## -----------------------------------------------------------------------------
h5se.bm <- minfi::getBeta(h5se.test)
dim(h5se.bm)

## -----------------------------------------------------------------------------
colnames(h5se.bm) <- h5se.test$gsm
knitr::kable(head(h5se.bm), align = "c")

## -----------------------------------------------------------------------------
an <- minfi::getAnnotation(h5se.test)
dim(an)

## -----------------------------------------------------------------------------
colnames(an)

## -----------------------------------------------------------------------------
ant <- as.matrix(t(an[c(1:4), c(1:3, 5:6, 9, 19, 24, 26)]))
knitr::kable(ant, align = "c")

## -----------------------------------------------------------------------------
dn <- "remethdb-h5_rg-test_0-0-1_1590090412.h5"
h5.test <- system.file("extdata", "h5test", dn, 
                    package = "recountmethylation")

## -----------------------------------------------------------------------------
h5.rg <- getrg(dbn = h5.test, all.gsm = TRUE)

## -----------------------------------------------------------------------------
h5.red <- minfi::getRed(h5.rg)
h5.green <- minfi::getGreen(h5.rg)
dim(h5.red)

## -----------------------------------------------------------------------------
knitr::kable(head(h5.red), align = "c")

## -----------------------------------------------------------------------------
knitr::kable(head(h5.green), align = "c")

## -----------------------------------------------------------------------------
identical(rownames(h5.red), rownames(h5.green))

## -----------------------------------------------------------------------------
dlpath <- tempdir()
gsmv <- c("GSM1038308", "GSM1038309")
geo.rg <- gds_idat2rg(gsmv, dfp = dlpath)
colnames(geo.rg) <- gsub("\\_.*", "", colnames(geo.rg))

## -----------------------------------------------------------------------------
geo.red <- minfi::getRed(geo.rg)
geo.green <- minfi::getGreen(geo.rg)

## -----------------------------------------------------------------------------
int.addr <- intersect(rownames(geo.red), rownames(h5.red))
geo.red <- geo.red[int.addr,]
geo.green <- geo.green[int.addr,]
geo.red <- geo.red[order(match(rownames(geo.red), rownames(h5.red))),]
geo.green <- geo.green[order(match(rownames(geo.green), rownames(h5.green))),]
identical(rownames(geo.red), rownames(h5.red))
identical(rownames(geo.green), rownames(h5.green))
class(h5.red) <- "integer"
class(h5.green) <- "integer"

## -----------------------------------------------------------------------------
identical(geo.red, h5.red)

## -----------------------------------------------------------------------------
identical(geo.green, h5.green)

## -----------------------------------------------------------------------------
geo.gr <- minfi::preprocessNoob(geo.rg)

## -----------------------------------------------------------------------------
geo.bm <- as.matrix(minfi::getBeta(geo.gr))

## -----------------------------------------------------------------------------
h5se.bm <- as.matrix(h5se.bm)
int.cg <- intersect(rownames(geo.bm), rownames(h5se.bm))
geo.bm <- geo.bm[int.cg,]
geo.bm <- geo.bm[order(match(rownames(geo.bm), rownames(h5se.bm))),]

## -----------------------------------------------------------------------------
identical(summary(geo.bm), summary(h5se.bm))

## -----------------------------------------------------------------------------
identical(rownames(geo.bm), rownames(h5se.bm))

## ----get_sessioninfo----------------------------------------------------------
sessionInfo()