A comprehensive guide to using the minfi package for analyzing DNA methylation microarrays from Illumina.vignettes.
minfi 1.53.1
The 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.
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.
getSex(), quality control using getQC(), quantile normalization using preprocessQuantile(), bump hunting using bumphunter() and block finding using blockFinder().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.preprocessQuantile(), it would be considerate to also cite (Touleimat and Tost 2012), since this publication describes a method essentially identical to preprocessQuantile().bumphunter(), you should also cite the original bump hunter publication (Jaffe et al. 2012).preprocessSWAN() please cite (Maksimovic, Gordon, and Oshlack 2012).preprocessNoob(), please cite (Triche et al. 2013).preprocessFunnorm(), please cite (Fortin et al. 2014). The default in preprocessFunnorm() is to do noob background correction. If this is used, please also cite (Triche et al. 2013).compartments() and extractAB(), please cite (Fortin and Hansen 2015).If you’re using Bibtex, you can get the citations in this format by
toBibtex(citation("minfi"))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.
This document has the following dependencies
library(minfi)
library(minfiData)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 minfiData package. This is the the number of features and and classes as we move through the class hierarchy (code not run):
RGsetEx## class: RGChannelSet
## dim: 622399 6
## metadata(0):
## assays(2): Green Red
## rownames(622399): 10600313 10600322 ... 74810490 74810492
## rowData names(0):
## colnames(6): 5723646052_R02C02 5723646052_R04C01 ... 5723646053_R05C02
## 5723646053_R06C02
## colData names(13): Sample_Name Sample_Well ... Basename filenames
## Annotation
## array: IlluminaHumanMethylation450k
## annotation: ilmn12.hg19## RGsetEx: RGChannelSet, 622,399 features
MsetEx <- preprocessRaw(RGsetEx)
## MsetEx: MethylSet, 485,512 features
GMsetEx <- mapToGenome(MsetEx)
## GMsetEx: GenomicMethylSet, 485,512 featuresNote 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).
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)## [1] "5723646052" "5723646053" "SampleSheet.csv"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"))## [1] "5723646052_R02C02_Grn.idat" "5723646052_R02C02_Red.idat"
## [3] "5723646052_R04C01_Grn.idat" "5723646052_R04C01_Red.idat"
## [5] "5723646052_R05C02_Grn.idat" "5723646052_R05C02_Red.idat"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)## [read.metharray.sheet] Found the following CSV files:## [1] "/home/biocbuild/bbs-3.21-bioc/R/site-library/minfiData/extdata/SampleSheet.csv"targets## Sample_Name Sample_Well Sample_Plate Sample_Group Pool_ID person age sex
## 1 GroupA_3 H5 <NA> GroupA <NA> id3 83 M
## 2 GroupA_2 D5 <NA> GroupA <NA> id2 58 F
## 3 GroupB_3 C6 <NA> GroupB <NA> id3 83 M
## 4 GroupB_1 F7 <NA> GroupB <NA> id1 75 F
## 5 GroupA_1 G7 <NA> GroupA <NA> id1 75 F
## 6 GroupB_2 H7 <NA> GroupB <NA> id2 58 F
## status Array Slide
## 1 normal R02C02 5723646052
## 2 normal R04C01 5723646052
## 3 cancer R05C02 5723646052
## 4 cancer R04C02 5723646053
## 5 normal R05C02 5723646053
## 6 cancer R06C02 5723646053
## Basename
## 1 /home/biocbuild/bbs-3.21-bioc/R/site-library/minfiData/extdata/5723646052/5723646052_R02C02
## 2 /home/biocbuild/bbs-3.21-bioc/R/site-library/minfiData/extdata/5723646052/5723646052_R04C01
## 3 /home/biocbuild/bbs-3.21-bioc/R/site-library/minfiData/extdata/5723646052/5723646052_R05C02
## 4 /home/biocbuild/bbs-3.21-bioc/R/site-library/minfiData/extdata/5723646053/5723646053_R04C02
## 5 /home/biocbuild/bbs-3.21-bioc/R/site-library/minfiData/extdata/5723646053/5723646053_R05C02
## 6 /home/biocbuild/bbs-3.21-bioc/R/site-library/minfiData/extdata/5723646053/5723646053_R06C02First 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)## [1] "/5723646052/5723646052_R02C02" "/5723646052/5723646052_R04C01"
## [3] "/5723646052/5723646052_R05C02" "/5723646053/5723646053_R04C02"
## [5] "/5723646053/5723646053_R05C02" "/5723646053/5723646053_R06C02"(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## class: RGChannelSet
## dim: 622399 6
## metadata(0):
## assays(2): Green Red
## rownames(622399): 10600313 10600322 ... 74810490 74810492
## rowData names(0):
## colnames(6): 5723646052_R02C02 5723646052_R04C01 ... 5723646053_R05C02
## 5723646053_R06C02
## colData names(13): Sample_Name Sample_Well ... Basename filenames
## Annotation
## array: IlluminaHumanMethylation450k
## annotation: ilmn12.hg19pd <- pData(RGset)
pd[,1:4]## DataFrame with 6 rows and 4 columns
## Sample_Name Sample_Well Sample_Plate Sample_Group
## <character> <character> <character> <character>
## 5723646052_R02C02 GroupA_3 H5 NA GroupA
## 5723646052_R04C01 GroupA_2 D5 NA GroupA
## 5723646052_R05C02 GroupB_3 C6 NA GroupB
## 5723646053_R04C02 GroupB_1 F7 NA GroupB
## 5723646053_R05C02 GroupA_1 G7 NA GroupA
## 5723646053_R06C02 GroupB_2 H7 NA GroupBThe 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)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## Sample_Name Sample_Well Sample_Plate Sample_Group Pool_ID Sentrix_ID
## 1 GroupA_3 H5 NA GroupA NA 5723646052
## 2 GroupA_2 D5 NA GroupA NA 5723646052
## 3 GroupB_3 C6 NA GroupB NA 5723646052
## 4 GroupB_1 F7 NA GroupB NA 5723646053
## 5 GroupA_1 G7 NA GroupA NA 5723646053
## 6 GroupB_2 H7 NA GroupB NA 5723646053
## Sentrix_Position person age sex status
## 1 R02C02 id3 83 M normal
## 2 R04C01 id2 58 F normal
## 3 R05C02 id3 83 M cancer
## 4 R04C02 id1 75 F cancer
## 5 R05C02 id1 75 F normal
## 6 R06C02 id2 58 F cancerWe 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 Biobase is combine() that combines (“adds”) two sets of samples. This allows the user to manually build up an RGChannelSet.
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)## array annotation
## "IlluminaHumanMethylation450k" "ilmn12.hg19"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 affy package (an AffyBatch is preprocessed into an ExpressionSet using array design information stored in a CDF environment (package)).
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 (Maksimovic, Gordon, and Oshlack 2012).preprocessQuantile : Quantile normalization (adapted to DNA methylation arrays), described in Aryee et al. (2014)preprocessNoob : Noob preprocessing, described in (Triche et al. 2013).preprocessFunnorm : Functional normalization as described in (Fortin et al. 2014).FIXME: discuss literature