In kasperdanielhansen/minfi: Analyze Illumina Infinium DNA methylation arrays

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

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

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

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):

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

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):

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

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:

sub(baseDir, "", targets$Basename)

(This is just for display purposes).

With this data.frame, it is easy to read in the data

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.

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

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.

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

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

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

toLatex(sessionInfo())

References

kasperdanielhansen/minfi documentation built on May 4, 2024, 12:04 p.m.