---
title: "recountmethylation User's Guide"
author:
- Sean K. Maden
- Reid F. Thompson
- Kasper D. Hansen
- Abhinav Nellore
date: "`r format(Sys.time(), '%d %B, %Y')`"
bibliography: bibliography.bib
package: recountmethylation
vignette: > 
  %\VignetteIndexEntry{recountmethylation User's Guide}
  %\VignetteDepends{RCurl}
  %\usepackage[UTF-8]{inputenc} 
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
output:
  BiocStyle::html_document:
    code_folding: show
    toc: yes
    tocfloat: yes
  BiocStyle::pdf_document: 
    toc: yes
    toc_depth: 2
---

```{r 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)
```

# Introduction and overview

The `recountmethylation` package provides access to databases of DNA 
methylation (DNAm) data from over 35,000 sample records with IDATs in the Gene 
Expression Omnibus (GEO, available through March 31, 2019) run using the 
Illumina HM450K BeadArray platform, including metadata, raw/unnormalized red 
and green signals, raw/unnormalized methylated and unmethylated signals, and 
normalized DNAm fractions or Beta-values (@maden_human_2020). Normalization 
was performed using the out-of-band signal correction (a.k.a. "noob") method, 
a type of within-sample normalization (@triche_low-level_2013).

This User's Guide shows how to use the `recountmethylation` package to obtain, 
load, and query the DNAm databases with 2 small example files. Background about 
DNAm arrays, DNAm measurement, `SummarizedExperiment` objects, database file 
types, and samples metadata is also provided. Further analysis examples are 
contained in the `data_analyses` vignette.

## Disclaimer

```{r disclaimer, echo = FALSE, message = TRUE}
library(recountmethylation)
```

## Database files and access

Database access, including downloads and file loads, are managed by the 
`get_db` functions. These download and access the latest available database 
files (see `?get_db` and below examples for details). Note you will need between 
90-120 Gb of disk space to store a single database file. Files pair sample 
metadata with assays including red and green channel signals, methylated and 
unmethylated level signals, and DNAm fractions in 3 `HDF5-SummarizedExperiment` 
entities, and red and green signals in an `HDF5` `h5` database. 

The database files are contained at the file server, located at the URL: [https://recount.bio/data/](https://recount.bio/data/). Details about 
the latest files are as follows.

```{r, 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")
```

## ExperimentHub integration

The DNAm array database files are indexed on `ExperimentHub`, and are viewable as follows.

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

In addition to using the `getdb` functions, the `HDF5` (`.h5` extension) files may alternatively be downloaded from the hubs, as follows.

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

Note that whether downloads use the hubs or `getdb` functions, caching is implemented to check for previously downloaded database files.

# Background

This section includes essential background about DNAm array platforms, assays 
and file types, and sample metadata.

## DNAm arrays

Databases include human samples run on the Illumina Infinium HM450K BeadArray 
platform. HM450K is a popular 2-channel platform that probes over 480,000 CpG 
loci genome-wide, with enriched coverage at CG islands, genes, and enhancers 
(@sandoval_validation_2011). Array processing generates 2 intensity files 
(IDATs) per sample, one each for the red and green color channels. These raw 
files also contain control signals useful for quality evaluations 
@noauthor_illumina_2010. 

HM450K probes use either of 2 bead technologies, known as Type I and Type II, 
where the majority (72%) of probes use the latter. For Type II probes, a single 
bead assay informs a single probe, while Type I probes use 2 beads each. 
Practically, this means the bead-specific matrices found in `RGChannelSet` 
objects are larger than the probe-specific matrices found in derived object 
types (e.g. 622,399 assays for red/green signal matrices versus 485,512 assays 
for methylated/unmethylated signal, DNAm fractions matrices, see below).

## SummarizedExperiment object classes

DNAm array sample IDATs can be read into an R session as an object of class 
`RGChannelSet`, a type of `SummarizedExperiment`. These objects support 
analyses of high-throughput genomics datasets, and they include slots for 
assay matrices, sample metadata, and experiment metadata. During a typical 
workflow, normalization and preprocessing convert `RGChannelSet` objects into 
new types like `MethylSet` and `RatioSet`. While not all IDAT information is 
accessible from every object type (e.g. only `RGChannelSet`s can contain 
control assays), derived objects like `MethylSet`s and `RatioSet`s may be 
smaller and/or faster to access.

Three `SummarizedExperiment` databases are provided as 
`HDF5-SummarizedExperiment` files, including an unnormalized `RGChannelSet` 
(red/green signals), an unnormalized `MethylSet` (methylated/unmethylated 
signals) and a normalized `GenomicRatioSet` (DNAm fractions). For the latter, 
DNAm fractions (logit2 Beta-values, or M-values) were normalized using the 
out-of-band signal or "noob" method, an effective within-sample normalization 
that removes signal artifacts (@triche_low-level_2013).

## Database file types

Database files are stored as either `HDF5` or `HDF5-SummarizedExperiment`. For 
most R users, the latter files will be most convenient to work with. `HDF5`, or 
hierarchical data format 5, combines compression and chunking for convenient 
handling of large datasets. `HDF5-SummarizedExperiment` files combine the 
benefits of `HDF5` and `SummarizedExperiment` entities using a 
DelayedArray-powered backend. Once an `HDF5-SummarizedExperiment` file is 
loaded, it can be treated similarly to a `SummarizedExperiment` object in 
active memory. That is, summary and subset operations execute rapidly, and 
realization of large data chunks in active memory is delayed until called for 
by the script (see examples).

## Sample metadata

Sample metadata are included with DNAm assays in the database files. Currently, 
metadata variables include GEO record IDs for samples (GSM) and studies (GSE), 
sample record titles, learned labels for tissue and disease, sample type 
predictions from the MetaSRA-pipeline, and DNAm model-based predictions for 
age, sex, and blood cell types. Access sample metadata from 
`SummarizedExperiment` objects using the `pData` minfi function (see examples). 
Examples in the `data_analyses` vignette illustrate some ways to utilize the 
provided sample metadata.

Provided metadata derives from the GSE-specific SOFT files, which contain 
experiment, sample, and platform metadata. Considerable efforts were made to 
learn, harmonize, and predict metadata labels. Certain types of info lacking 
in the `recountmethylation` metadata may be available in the SOFT files, 
especially if it is sample non-specific (e.g. methods text, PubMed ID, etc.) 
or redundant with DNAm-derived metrics (e.g. DNAm summaries, predicted sex, 
etc.).

It is good practice to validate the harmonized metadata with original metadata 
records, especially where labels are ambiguous or there is insufficient 
information for a given query. GEO GSM and GSE records can be viewed from a 
browser, or SOFT files may be downloaded directly. Packages like GEOmetadb and 
GEOquery are also useful to query and summarize GEO metadata.

# HDF5-SummarizedExperiment example

This example shows basic handling for `HDF5-SummarizedExperiment` (a.k.a. 
"h5se") files. For these files, the `getdb` function returns the loaded file. 
Thanks to a `DelayedArray` backend, even full-sized `h5se` databases can be 
treated as if they were fully loaded into active memory.

## Obtain the test database

The test h5se dataset includes sample metadata and noob-normalized 
DNAm fractions (Beta-values) for chromosome 22 probes for 2 samples. 
Datasets can be downloaded using the `getdb` series of functions 
(see `?getdb` for details), where the `dfp` argument specifies the 
download destination. The test h5se file is included in the package 
"inst" directory, and can be loaded as follows.

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

## Inspect and summarize the database

Common characterization functions can be used on the dataset after it has been 
loaded. These include functions for `SummarizedExperiment`-like objects, such 
as the `getBeta`, `pData`, and `getAnnotation` minfi functions. First, inspect 
the dataset using standard functions like `class`, `dim`, and `summary` as 
follows.

```{r}
class(h5se.test)
```

```{r}
dim(h5se.test)
```

```{r}
summary(h5se.test)
```

Access the sample metadata for the 2 available samples using `pData`. 

```{r}
h5se.md <- minfi::pData(h5se.test)
dim(h5se.md)
```
```{r}
colnames(h5se.md)
```

Next get CpG probe-specific DNAm fractions, or "Beta-values", with `getBeta` 
(rows are probes, columns are samples).

```{r}
h5se.bm <- minfi::getBeta(h5se.test)
dim(h5se.bm)
```
```{r}
colnames(h5se.bm) <- h5se.test$gsm
knitr::kable(head(h5se.bm), align = "c")
```

Access manifest information for probes with `getAnnotation`. This includes the 
bead addresses, probe type, and genome coordinates and regions. For full details 
about the probe annotations, consult the minfi and Illumina platform documentation.

```{r}
an <- minfi::getAnnotation(h5se.test)
dim(an)
```

```{r}
colnames(an)
```

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

# HDF5 database and example

To provide more workflow options, bead-specific red and green signal data have 
been provided with sample metadata in an `HDF5`/h5 file. This example shows 
how to handle objects of this type with `recountmethylation`.

## Obtain the test database

The test h5 file includes metadata and bead-specific signals from 
chromosome 22 for the same 2 samples as in the h5se test file. 
Note `getdb` functions for h5 files simply return the database path.
Since the test h5 file has also been included in the package "inst" folder,
get the path to load the file as follows.

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

## Inspect and summarize the database

Use the file path to read data into an `RGChannelSet` with the `getrg` 
function. Setting `all.gsm = TRUE` obtains data for all samples in the
database files, while passing a vector of GSM IDs to `gsmv` argument 
will query a subset of available samples. Signals from all available probes are retrieved by default, and probe subsets can be obtained by passing a vector of 
valid bead addresses to the `cgv` argument.

```{r}
h5.rg <- getrg(dbn = h5.test, all.gsm = TRUE)
```

To avoid exhausting active memory with the full-sized `h5` dataset, provide 
either `gsmv` or `cgv` to `getrg`, and set either `all.cg` or `all.gsm` to 
FALSE (see `?getrg` for details).

As in the previous example, use `pData` and `getAnnotation` to get sample 
metadata and array manifest information, respectively. Access the green and 
red signal matrices in the `RGChannelSet` with the `getRed` and `getGreen` 
minfi functions.

```{r}
h5.red <- minfi::getRed(h5.rg)
h5.green <- minfi::getGreen(h5.rg)
dim(h5.red)
```
```{r}
knitr::kable(head(h5.red), align = "c")
```
```{r}
knitr::kable(head(h5.green), align = "c")
```
```{r}
identical(rownames(h5.red), rownames(h5.green))
```

Rows in these signal matrices map to bead addresses rather than probe IDs. 
These matrices have more rows than the h5se test Beta-value matrix because 
any type I probes use data from 2 beads each.

# Validate DNAm datasets

This section demonstrates validation using the test databases. As the 
disclaimer notes, it is good practice to validate data against the latest 
available GEO files. This step may be most useful for newer samples published 
close to the end compilation date (end of March 2019 for current version), 
which may be more prone to revisions at initial publication.

## Download and read IDATs from the GEO database server

Use the `gds_idat2rg` function to download IDATs for the 2 test samples 
and load these into a new `RGChannelSet` object. Do this by passing a vector
of GSM IDs to `gsmv` and the download
destination to `dfp`.

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

## Compare DNAm signals

Extract the red and green signal matrices from `geo.rg`.

```{r}
geo.red <- minfi::getRed(geo.rg)
geo.green <- minfi::getGreen(geo.rg)
```

Match indices and labels between the GEO and `h5` test signal matrices.

```{r}
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"
```

Finally, compare the signal matrix data.

```{r}
identical(geo.red, h5.red)
```
```{r}
identical(geo.green, h5.green)
```

## Compare DNAm Beta-values

Before comparing the GEO-downloaded data to data from the `h5se.test` database, normalize the data using the same out-of-band or "noob" normalization technique 
that was used to generate data in the h5se database.

```{r}
geo.gr <- minfi::preprocessNoob(geo.rg)
```

Next, extract the Beta-values.

```{r}
geo.bm <- as.matrix(minfi::getBeta(geo.gr))
```

Now match row and column labels and indices.

```{r}
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))),]
```

Finally, compare the two datasets.

```{r}
identical(summary(geo.bm), summary(h5se.bm))
```
```{r}
identical(rownames(geo.bm), rownames(h5se.bm))
```

# Get more help

Consult the Data Analyses [vignette](link.url) and main [manuscript](link.url) 
for analysis examples and details about data compilations.

# Session info

```{r get_sessioninfo}
sessionInfo()
```

# Works Cited