--- title: "The minfi User's Guide" shorttitle: "minfi guide" author: - Kasper D. Hansen - Jean-Phillipe Fortin package: minfi bibliography: minfi.bib abstract: > A comprehensive guide to using the minfi package for analyzing DNA methylation microarrays from Illumina.vignettes. vignette: > %\VignetteIndexEntry{minfi User's Guide} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} output: BiocStyle::html_document: toc_float: true --- # Introduction The `r Biocpkg("minfi")` package provides tools for analyzing Illumina's Methylation arrays, specifically the 450k and EPIC (also known as the 850k) arrays. We have partial support for the older 27k array. The tasks addressed in this package include preprocessing, QC assessments, identification of interesting methylation loci and plotting functionality. Analyzing these types of arrays is ongoing research in ours and other groups. The input data to this package are IDAT files, representing two different color channels prior to normalization. This is the most complete data type, because IDAT files includes measurements on control probes. It is possible to use Genome Studio files together with the data structures contained in this package, but only some functionality is available because Genome Studio output does not contain control probe information. In addition, usually Genome Studio output is normalized using the methods implemented in Genome Studio and these are generally considered inferior. ## Citing the minfi package The MINFI package contains methods which are described across multiple manuscripts, by different non-overlapping authors. This makes citing the package a bit difficult, so here is a guide. - If you are using MINFI in a publication, please cite [@minfi]. This publication includes details on sex estimation using `getSex()`, quality control using `getQC()`, quantile normalization using `preprocessQuantile()`, bump hunting using `bumphunter()` and block finding using `blockFinder()`. - If you are using MINFI to analyze EPIC or 27k arrays, please cite [@minfiEPIC]. The publication includes details on `convertArray()` and `combineArrays()`, extending NOOB to work in single-sample mode as well as using `estimateCellCounts()` with reference data from the 450k array to estimate cell type composition for EPIC data. - If you are using `preprocessQuantile()`, it would be considerate to also cite [@Touleimat:2012], since this publication describes a method essentially identical to `preprocessQuantile()`. - If you are using `bumphunter()`, you should also cite the original bump hunter publication [@bumphunter]. - If you are using SWAN normalization as implemented in `preprocessSWAN()` please cite [@SWAN]. - If you are using noob background correction as implemented in `preprocessNoob()`, please cite [@noob]. - If you are using functional normalization as implemented in `preprocessFunnorm()`, please cite [@funnorm]. The default in `preprocessFunnorm()` is to do noob background correction. If this is used, please also cite [@noob]. - If you are estimating A/B compartments as implemented in `compartments()` and `extractAB()`, please cite [@compartments]. If you're using Bibtex, you can get the citations in this format by ```{r citation,eval=FALSE} toBibtex(citation("minfi")) ``` ## Terminology The literature is often a bit unspecific wrt. the terminology for the DNA methylation microarrays. For the 450k microarray, each sample is measured on a single array, in two different color channels (red and green). Each array measures roughly 450,000 CpG positions. Each CpG is associated with two measurements: a methylated measurement and an "un"-methylated measurement. These two values can be measured in one of two ways: using a "Type I" design or a "Type II design". CpGs measured using a Type I design are measured using a single color, with two different probes in the same color channel providing the methylated and the unmethylated measurements. CpGs measured using a Type II design are measured using a single probe, and two different colors provide the methylated and the unmethylated measurements. Practically, this implies that on this array there is **not** a one-to-one correspondence between probes and CpG positions. We have therefore tried to be precise about this and we refer to a "locus" (or "CpG") when we refer to a single-base genomic locus, and we differentiate this from a "probe". The previous generation 27k methylation array uses only the Type I design, and the EPIC arrays uses both Type I and Type II. Differences in DNA methylation between samples can either be at a single CpG which is called a differentially methylated position (DMP), or at a regional level which is called a differentially methylated region (DMR). Physically, each sample is measured on a single "array". For the 450k design, there are 12 arrays on a single physical "slide" (organized in a 6 by 2 grid). Slides are organized into "plates" containing at most 8 slides (96 arrays). The EPIC array has 8 arrays per slide and 64 arrays per plate. ## Dependencies This document has the following dependencies ```{r dependencies, warning=FALSE, message=FALSE} library(minfi) library(minfiData) ``` # minfi classes The MINFI package is designed to be very flexible for methods developers. This flexibility comes at a cost for users; they need to understand a few different data classes: - `RGChannelSet` : raw data from the IDAT files; this data is organized at the probe (not CpG locus) level. This data has two channels: Red and Green. - `MethylSet` : data organized by the CpG locus level, but not mapped to a genome. This data has two channels: Meth (methylated) and Unmeth (unmethylated). - `RatioSet` : data organized by the CpG locus level, but not mapped to a genome. The data has at least one of two channels: Beta and/or M (logratio of Beta). It may optionally include a CN channel (copy number). - `GenomicMethylSet` : like a `MethylSet`, but mapped to a genome. - `GenomicRatioSet` : like a `RatioSet`, but mapped to the genome. A class hierarchy is as follows ![The class hierarchy of minfi](./minfiClasses.pdf) To make this more clear, let us look at the example data called `RGsetEx` from the `r Biocpkg("minfiData")` package. This is the the number of features and and classes as we move through the class hierarchy (code not run): ```{r RGsetEx} RGsetEx ## RGsetEx: RGChannelSet, 622,399 features MsetEx <- preprocessRaw(RGsetEx) ## MsetEx: MethylSet, 485,512 features GMsetEx <- mapToGenome(MsetEx) ## GMsetEx: GenomicMethylSet, 485,512 features ``` Note how the number of features changes. In the `RGChannelSet` a feature is a probe (which is different from a CpG locus, see the Terminology section). In the `MethylSet` each feature is now a methylation locus, and it has fewer features because some loci are measured using two probes. Finally, the `GenomicMethylSet` has the same size as the `MethylSet`, but it could in principle be smaller in case the annotation you use says that some probes do not map to the genome (in this case hg19). Finally we can convert to a `RatioSet` by `ratioConvert()`. The two functions `ratioConvert()` and `mapToGenome()` commute, as shown by the class hierarchy diagram above. Many preprocessing functions goes through several steps in the diagram, for example if the function needs to know the genomic location of the probes (several preprocessing functions handle probes measured on the sex chromosomes in a different way). # Reading data This package supports analysis of IDAT files, containing the summarized bead information. In our experience, most labs use a "Sample Sheet" CSV file to describe the layout of the experiment. This is based on a sample sheet file provided by Illumina. Our pipeline assumes the existence of such a file(s), but it is relatively easy to create such a file using for example Excel, if it is not available. We use an example dataset with 6 samples, spread across two slides. First we obtain the system path to the IDAT files; this requires a bit since the data comes from an installed package ```{r baseDir} baseDir <- system.file("extdata", package = "minfiData") list.files(baseDir) ``` This shows the typical layout of 450k data: each "slide" (containing 12 arrays, see Termiology) is stored in a separate directory, with a numeric name. The top level directory contains the sample sheet file. Inside the slide directories we find the IDAT files (and possible a number of JPG images or other files): ```{r baseDir2} list.files(file.path(baseDir, "5723646052")) ``` The files for each array has another numeric number and consists of a Red and a Grn (Green) IDAT file. Note that for this example data, each slide contains only 3 arrays and not 12. This was done because of file size limitations and because we only need 6 arrays to illustrate the package's functionality. First we read the sample sheet. We provide a convenience function for reading in this file `read.metharray.sheet()`. This function has a couple of attractive bells and whistles. Let us look at the output ```{r sheet} targets <- read.metharray.sheet(baseDir) targets ``` First the output: this is just a `data.frame`. It contains a column `Basename` that describes the location of the IDAT file corresponding to the sample, as well as two columns `Array` and `Slide`. In the sample sheet provided by Illumina, these two columns are named `Sentrix_Position` and `Sentrix_ID`, but we rename them. We provide more detail on the use of this function below. The `Basename` column tend to be too large for display, here it is simplified relative to `baseDir`: ```{r BasenameColumn>} sub(baseDir, "", targets$Basename) ``` (This is just for display purposes). With this `data.frame`, it is easy to read in the data ```{r readingTargets} RGset <- read.metharray.exp(targets = targets) ``` Let us look at the associated pheno data, which is really just the information contained in the targets object above. ```{r pData} RGset pd <- pData(RGset) pd[,1:4] ``` The `read.metharray.exp()` function also makes it possible to read in an entire directory or directory tree (with `recursive` set to `TRUE`) by using the function just with the argument `base` and `targets=NULL`, like ```{r read2} RGset2 <- read.metharray.exp(file.path(baseDir, "5723646052")) RGset3 <- read.metharray.exp(baseDir, recursive = TRUE) ``` ## Advanced notes on Reading Data The only important column in sheet `data.frame` used in the `targets` argument for the `read.metharray.exp()` function is a column names `Basename`. Typically, such an object would also have columns named `Array`, `Slide`, and (optionally) `Plate`. We used sheet data files build on top of the Sample Sheet data file provided by Illumina. This is a CSV file, with a header. In this case we assume that the phenotype data starts after a line beginning with `[Data]` (or that there is no header present). It is also easy to read a sample sheet manually, using the function `read.csv()`. Here, we know that we want to skip the first 7 lines of the file. ```{r sampleSheet2} targets2 <- read.csv(file.path(baseDir, "SampleSheet.csv"), stringsAsFactors = FALSE, skip = 7) targets2 ``` We now need to populate a `Basename` column. On possible approach is the following ```{r Basename} targets2$Basename <- file.path(baseDir, targets2$Sentrix_ID, paste0(targets2$Sentrix_ID, targets2$Sentrix_Position)) ``` Finally, MINFI contains a file-based parser: `read.metharray()`. The return object represents the red and the green channel measurements of the samples. A useful function that we get from the package `r Biocpkg("Biobase")` is `combine()` that combines ("adds") two sets of samples. This allows the user to manually build up an `RGChannelSet`. # Manifest / annotation ## What everyone needs to know For a methylation array, we have two types of annotation packages: "manifest" packages which contains the array design and "annotation" packages which contains information about where the methylation loci are located on the genome, which genomic features they map to and possible whether they overlap any known SNPs. You can see which packages are being used by ```{r annotation} annotation(RGsetEx) ``` ## Advanced discussion This discussion is intended for package developers or users who want to understand the internals of MINFI. A set of 450k data files will initially be read into an `RGChannelSet`, representing the raw intensities as two matrices: one being the green channel and one being the red channel. This is a class which is very similar to an `ExpressionSet` or an `NChannelSet`. The `RGChannelSet` is, together with a `IlluminaMethylationManifest` object, preprocessed into a `MethylSet`. The `IlluminaMethylationManifest` object contains the array design, and describes how probes and color channels are paired together to measure the methylation level at a specific CpG. The object also contains information about control probes (also known as QC probes). The `MethylSet` contains normalized data and essentially consists of two matrices containing the methylated and the unmethylated evidence for each CpG. Only the `RGChannelSet` contains information about the control probes. The process described in the previous paragraph is very similar to the paradigm for analyzing Affymetrix expression arrays using the `r Biocpkg("affy")` package (an `AffyBatch` is preprocessed into an `ExpressionSet` using array design information stored in a CDF environment (package)). # Quality control - use shinyMethyl - minfiQC - getSex checks - mds plots - plot conversion probes # Preprocessing Preprocessing in MINFI is done by a series of functions with names like `preprocessXX`. Different functions has different classes as output. Currently, we have - `preprocessRaw` : No processing. - `preprocessIllumina` : Illumina preprocessing, as performed by Genome Studio (reverse engineered by us). - `preprocessSWAN` : SWAN normalization, described in [@SWAN]. - `preprocessQuantile` : Quantile normalization (adapted to DNA methylation arrays), described in [@Touleimat:2012, @minfi] - `preprocessNoob` : Noob preprocessing, described in [@noob]. - `preprocessFunnorm` : Functional normalization as described in [@funnorm]. FIXME: discuss literature # SNPs and other issues - SNPs - cross reactive probes - remove bad probes - detection P - Gap Hunting # Identifying differentially methylated regions - DMP finding - Bump hunting - Block finding - plot DMRs with annotation. # Correcting for cell type heterogenity - estimateCellTypes - reference packages # Other stuff - Horvath age estimation - other packages # Sessioninfo ```{r sessionInfo, results='asis', echo=FALSE} toLatex(sessionInfo()) ``` # References