1 Introduction

Illumina Infinium HumanMethylation 450K BeadChip assay has become a standard tool to analyse methylation in human samples. Developed in 2011, it has already been used in projects such as The Cancer Genome Atlas (TCGA). Their 450.000 probes provide a good overall image of the methylation state of the genome, being one of the reasons of its success.

Given its complex design1 More information can be found at this minfi tutorial, many Bioconductor packages have been developed to assess normalization and pre-processing issues (e.g. minfi (Aryee et al. 2014) or lumi (Du, Kibbe, and Lin 2008)). In addition, these packages can detect differentially methylated probes (DMPs) and differentially methylated regions (DMRs). However, the interfaces are not very intuitive and several scripting steps are usually required.

MEAL aims to facilitate the analysis of Illumina Methylation 450K chips. We have included two methods to analyze DMPs (Differentially Methylated Probes), that test differences in means (limma) or differences in variance (DiffVar). We have included three DMRs (Differentially Methylated Regions) detection algorithms (bumphunter, blockFinder and DMRcate) and a new method to test differences in methylation in a target region (RDA). Finally, we have prepared plots for all these analyses as well as a wrapper to run all the analyses in the same dataset.

2 Input data

MEAL is meant to analyze methylation data already preprocessed. All our functions accept a GenomicRatioSet as input, which is a class from minfi package designed to manage preprocessed methylation data. Users willing to preprocess their own data are encouraged to take a look to minfi’s vignette

In this vignette, we will use methylation data from minfiData package.

library(MEAL)
library(MultiDataSet)
library(minfiData)
library(minfi)
library(ggplot2)

data("MsetEx")

MsetEx is a MethylationRatioSet that contains measurements for 485512 CpGs and 6 samples, as well as some phenotypic variables such as age or sex. The first step will be to convert it to a GenomicRatioSet. Then, we will add some extra features annotation. Finally, we will remove any CpG with NAs:

meth <- mapToGenome(ratioConvert(MsetEx))
rowData(meth) <- getAnnotation(meth)[, -c(1:3)]
meth <- meth[!apply(getBeta(meth), 1, function(x) any(is.na(x))), ]

3 Analyzing Methylation data

3.1 Pipeline

The function runPipeline run all methods included in MEAL to the same dataset. We only need to pass to this function a GenomicRatioSet and the name of our variable of interest. In our case, we will analyze the effect of cancer on methylation:

res <- runPipeline(set = meth, variable_names = "status")

runPipeline includes several parameters to customize the analyses. The most important parameters are covariable_names, betas and sva. covariable_names is used to include covariates in our models. betas allows the user choosing between running the analyis with beta (TRUE) or M-values (FALSE). If sva is TRUE, Surrogate Variable Analysis is run and surrogate variables are included in the models. Finally, some parameters modify the behaviour of the methods included in the wrapper and they will be covered later on. More information about the parameters can be found in the documentation (by typing ?runPipeline).

We will run a new analysis including age as covariate:

resAdj <- runPipeline(set = meth, variable_names = "status", 
                      covariable_names = "age")
resAdj
## Object of class 'ResultSet'
##  . created with: runPipeline 
##  . sva:  no 
##  . #results: 5 ( error: 0 )
##  . featureData: 485500 probes x 35 variables

3.2 Managing the results

runPipeline generates a ResultSet object. ResultSet is a class designed to encapsulate different results from the same dataset. It contains the results of the different methods, the feature data and other data required to get tables or plots. We can examine the analyses included in a ResultSet with the function names:

names(resAdj)
## [1] "DiffMean"    "DiffVar"     "bumphunter"  "blockFinder" "dmrcate"

Both objects contains five analyses. DiffMean is an analysis of difference of means performed with limma while the others are named with the method name (DiffVar, bumphunter, blockFinder and dmrcate).

We can use the function getAssociation to get a data.frame with the results, independent of the original method. This function has two main arguments: object and rid. object is the ResultSet with our data and rid is the name or the index of the analysis we want to extract.

head(getAssociation(resAdj, "DiffMean"))
##                 logFC   AveExpr         t      P.Value  adj.P.Val        B
## cg09383816 -0.5938196 0.4885486 -42.60661 7.389032e-07 0.01865675 5.544293
## cg27651090  0.5433331 0.5411453  40.28169 9.421312e-07 0.01865675 5.459338
## cg21938148 -0.6659934 0.4712352 -39.20530 1.059345e-06 0.01865675 5.415866
## cg25104555 -0.5254995 0.3906114 -39.05219 1.077443e-06 0.01865675 5.409451
## cg25937714 -0.5906359 0.4350042 -38.97928 1.086194e-06 0.01865675 5.406375
## cg15732851 -0.5760397 0.3869865 -38.03189 1.208270e-06 0.01865675 5.365130
head(getAssociation(resAdj, "DiffVar"))
##                logFC  AveExpr         t      P.Value adj.P.Val         B
## cg02939019 -2.166225 1.203729 -6.808503 0.0003492050 0.1828971 0.1498925
## cg11847929 -2.327910 1.465679 -6.780464 0.0003576462 0.1828971 0.1352083
## cg22976979 -1.910704 1.173032 -6.780138 0.0003577456 0.1828971 0.1350370
## cg25385529 -2.223012 1.347720 -6.760792 0.0003637068 0.1828971 0.1248493
## cg22676401 -1.914225 1.233126 -6.718301 0.0003772001 0.1828971 0.1023137
## cg27402591  2.710765 1.706821  6.715763 0.0003780238 0.1828971 0.1009608
head(getAssociation(resAdj, "bumphunter"))
##         chr     start       end      value     area cluster indexStart
## 89924  chr6 133561647 133562776 -0.4128635 15.68881  169280     180989
## 76143 chr10 118030848 118034357 -0.4244696 12.73409   28532     267188
## 80548 chr16  51183988  51187807 -0.3357711 11.75199   76311     380944
## 89156  chr6  29520698  29521803 -0.3009453 11.73687  163137     158788
## 78728 chr13  78492568  78493590 -0.4345396 11.73257   55851     333318
## 85384 chr20  61050560  61051915 -0.4558366 11.39591  122286     459589
##       indexEnd  L clusterL
## 89924   181026 38       43
## 76143   267217 30       30
## 80548   380978 35       42
## 89156   158826 39       40
## 78728   333344 27       50
## 85384   459613 25       26
head(getAssociation(resAdj, "blockFinder"))
##        chr     start       end     value     area cluster indexStart
## 431   chr2 217468708 219096237 0.1692934 29.96494      93      36345
## 1255  chr7  40098119  42446212 0.1706023 20.47228     331      91718
## 118   chr1 152056869 153234423 0.1674880 19.93108      32      13744
## 519   chr3  54169983  56448270 0.1621731 19.29860     105      43843
## 1834 chr11  55066025  56956400 0.1703365 17.71500     507     133354
## 1900 chr11 131278641 132757156 0.1475531 17.70638     520     141047
##      indexEnd   L clusterL
## 431     36561 177      527
## 1255    91861 120     1845
## 118     13891 119     1646
## 519     43985 119     1761
## 1834   133468 104     1803
## 1900   141176 120     1754
head(getAssociation(resAdj, "dmrcate"))
##                          coord no.cpgs        minfdr     Stouffer  maxbetafc
## 3503    chr6:33130696-33148812     145  1.584571e-89 1.862417e-48  0.4429870
## 3657  chr6:133561224-133564578      53  0.000000e+00 2.043200e-30 -0.6600767
## 1565   chr16:51183363-51190201      48 6.573279e-189 3.493944e-30 -0.5757717
## 457  chr10:118030292-118034357      31  0.000000e+00 3.028362e-29 -0.6980408
## 1570   chr16:54964677-54973966      42 9.691969e-112 4.130315e-29 -0.6023573
## 3411    chr6:29520527-29521803      40 2.138728e-204 1.422662e-27 -0.5406367
##      meanbetafc
## 3503  0.2063092
## 3657 -0.3541019
## 1565 -0.3008899
## 457  -0.4162689
## 1570 -0.2518663
## 3411 -0.2956325

DiffMean and DiffVar are internally stored as a MArrayLM, the class from limma results. This class allows testing different constrasts or evaluating different variables simultaneously. The function getProbeResults helps the user performing these operations. It also has the arguments object and rid from getAssociation. coef is a numeric with the index of the coefficient from which we want the results. If we did not pass a custom model to runPipeline, the first coefficient (coef = 1) is the intercept and the second coefficient (coef = 2) is the first variable that we included in variable_names. We can evaluate different coefficients simultaneously by passing a vector to coef. contrast is a matrix with the contrasts that we want to evaluate. This option is useful when our variable of interest is a factor with several levels and we want to do all the different comparisons. Finally, the argument fNames is used to select the variables from features annotation that will be added to the tables.

To exemplify the use of this function, we will evaluate our whole adjusted model, including age coefficient. We will also add some annotation of the CpGs:

head(getProbeResults(resAdj, rid = 1, coef = 2:3, 
                     fNames = c("chromosome", "start")))
##            statusnormal           age   AveExpr        F      P.Value
## cg09383816   -0.5938196 -0.0026657333 0.4885486 915.6123 2.002722e-06
## cg27651090    0.5433331 -0.0009235097 0.5411453 812.3259 2.594906e-06
## cg21938148   -0.6659934 -0.0028869823 0.4712352 774.8049 2.874494e-06
## cg25104555   -0.5254995 -0.0031493548 0.3906114 774.4414 2.877414e-06
## cg25937714   -0.5906359  0.0009181122 0.4350042 760.4902 2.992841e-06
## cg15732851   -0.5760397 -0.0050629750 0.3869865 747.4969 3.106531e-06
##             adj.P.Val chromosome     start
## cg09383816 0.04930423       chr8  67344556
## cg27651090 0.04930423      chr13 109270071
## cg21938148 0.04930423      chr13 110958977
## cg25104555 0.04930423      chr10  16562998
## cg25937714 0.04930423       chr3  44036492
## cg15732851 0.04930423       chr1 203598761

When more than one coefficient is evaluated, a estimate for each coefficient is returned and the t-statistic is substituted by a F-statistic. More information about linear models, including a detailed section of how to create a constrast matrix can be found in limma users’ guide.

Finally, we can obtain the results of CpGs mapped to some genes with the function getGeneVals. This function accepts the same arguments than getProbeResults but includes the arguments gene and genecol to pass the names of the genes to be selected and the column name of feature data containing gene names.

We will retrieve the difference in variance results for all CpGs mapped to ARMS2. We can see in the rowData of meth that gene names are in the column ‘UCSC_RefGene_Name’:

getGeneVals(resAdj, "ARMS2", genecol = "UCSC_RefGene_Name", fNames = c("chromosome", "start"))
##                  logFC   AveExpr          t     P.Value  adj.P.Val         B
## cg24296920  0.32587475 0.5876988  5.2591402 0.004978747 0.06840951 -1.899236
## cg24884230  0.19254764 0.6064401  4.1256910 0.012289860 0.11028866 -2.938109
## cg00676728  0.11068528 0.8206485  3.0338330 0.034758240 0.19634111 -4.124759
## cg03623097  0.12550116 0.5342317  2.5137794 0.060976354 0.26724923 -4.753218
## cg18222240  0.18747014 0.5432686  2.4828934 0.063131502 0.27251533 -4.791562
## cg13265583  0.19110580 0.6897652  2.0494239 0.104280251 0.35390133 -5.336792
## cg25542438 -0.01966479 0.8436013 -0.7907323 0.470081928 0.75622351 -6.770906
##            chromosome     start UCSC_RefGene_Name
## cg24296920      chr10 124214120             ARMS2
## cg24884230      chr10 124216658             ARMS2
## cg00676728      chr10 124213760             ARMS2
## cg03623097      chr10 124213466             ARMS2
## cg18222240      chr10 124213527             ARMS2
## cg13265583      chr10 124214151             ARMS2
## cg25542438      chr10 124214250             ARMS2

3.3 Plotting the results

We can easily get Manhattan plots, Volcano plots and QQ-plots for the probes results (DiffMean and DiffVar) using plot method. Our extension of plot method to ResultSet includes the arguments rid or coef that were already present in getProbeResult. In addition, the argument type allows choosing between a Manhattan plot (“manhattan”), a Volcano plot (“volcano”) or a qq-plot (“qq”).

We will get these three plots for the unadjusted analysis of differential means:

plot(resAdj, rid = "DiffMean", type = "manhattan")

plot(resAdj, rid = "DiffMean", type = "volcano")

plot(resAdj, rid = "DiffMean", type = "qq")

3.3.1 Manhattan plot

We can customize different aspects of a Manhattan plot. We can highlight the CpGs of a target region by passing a GenomicRanges to the argument highlight. Similarly, we can get a Manhattan plot with only the CpGs of our target region passing a GenomicRanges to the argument subset. It should be noticed that the GenomicRange should have the chromosome as a number (1-24).

We will show these capabilities by highlighting and subsetting a region of ten Mb in chromosome X:

targetRange <- GRanges("23:13000000-23000000")
plot(resAdj, rid = "DiffMean", type = "manhattan", highlight = targetRange)

plot(resAdj, rid = "DiffMean", type = "manhattan", subset = targetRange)

We can also change the height of lines marking different levels of significance. Height of blue line can be set with suggestiveline parameter and red line with genomewideline parameter. It should be noticed that these values are expressed as -log10 of p-value:

plot(resAdj, rid = "DiffMean", type = "manhattan", suggestiveline = 3, 
     genomewideline = 6)

Finally, as our Manhattan plot is done with base framework, we can customize the plot using base plotting functions such as points, lines or text or arguments of plot function like main:

plot(resAdj, rid = "DiffMean", type = "manhattan", 
     main = "My custom Manhattan")
abline(h = 13, col = "yellow")

3.3.2 Volcano plot

In our Volcano plot, we can also customize the thresholds for statistical significance and magnitude of the effect using the arguments tPV and tFC. As in the previous case, tPV is expressed as -log10 of p-value. On the other hand, tFC units will change depending if we used beta or M-values. Finally, show.labels can turn on and turn off the labelling of significant features:

plot(resAdj, rid = "DiffMean", type = "volcano", tPV = 14, tFC = 0.4, 
     show.labels = FALSE)

Volcano plot is based on ggplot2 so we can further customize the plot adding new layers:

plot(resAdj, rid = "DiffMean", type = "volcano", tPV = 30, tFC = 0.1, 
     show.labels = FALSE) + ggtitle("My custom Volcano")

3.3.3 QQplot

Our QQplot include the computation of the lambda, a measure of the inflation of the p-values. We can remove this value with the parameter show.lambda.

Our qqplot is also based on ggplot2 so we will add a title to customize it:

plot(resAdj, rid = "DiffMean", type = "qq") + ggtitle("My custom QQplot")

3.3.4 Features

MEAL incorporates the function plotFeature to plot the beta values distribution of a CpG. plotFeature has three main arguments. set is the GenomicRatioSet with the methylation data. feat is the index or name of our target CpG. variables is a character vector with the names of the variables used in the plot. We can include two variables in our plot.

In the next line, we will plot a CpG with high difference in means between male and female (cg17547524) and a CpG with high difference in variance (cg05111645) vs sex. As plotFeature is based on ggplot2, we can customize it:

plotFeature(set = meth, feat = "cg17547524", variables = "sex") + 
  ggtitle("Diff Means")

plotFeature(set = meth, feat = "cg05111645", variables = "sex") + 
  ggtitle("Diff Vars")

3.3.5 Regional plotting

We can simultaneously plot the different results in a target region along with gene and CpG annotation with the function plotRegion. This function has two main arguments. rset is the ResultSet and range is a GenomicRanges with our target region.

We will plot a region of 1 Mb in chromosome X:

targetRange <- GRanges("chrX:13000000-14000000")
plotRegion(resAdj, targetRange)

Our plot has three main parts. The top contains the annotation of the regional genes and the CpGs included in the analysis. The middle part contains the results of the DMR detection methods (Bumphunter, blockFinder and DMRcate). The bottom part contains the results of the single probe analyses (differential mean and differential variance). Each analysis has two parts: the coefficients and the p-values. The line in the p-values plot marks the significance threshold.

By default,plotRegion includes all analyses run in the plot. However, we can plot only few analyses with the parameter results. We can also modify the height of the p-value line with the parameter tPV (units are -log10 of p-value):

plotRegion(resAdj, targetRange, results = c("DiffMean", "bumphunter"), 
           tPV = 10)

3.4 Methods wrappers

MEAL includes wrappers to run the different methods of the pipeline individually. All these functions accept a GenomicRatioSet as input and can return the results in a ResultSet. Consequently, functionalities described in the above section for the results of the pipeline also apply for the results of a single method.

3.4.1 Differences of mean analysis

We can test if a phenotype causes changes in methylation means using the runDiffMeanAnalysis. This function is a wrapper of lmFit function from limma and requires two arguments: set and model. set contains the methylation data, either in a GenomicRatioSet or a matrix. model can be a matrix with the linear model or a formula indicating the model. In the former case, set must be a GenomicRatioSet and the variables included in the model must be present in the colData of our set.

We exemplify the use of this function by running the same linear model than in our pipeline:

resDM <- runDiffMeanAnalysis(set = meth, model = ~ status)

runDiffMeanAnalysis also has other parameters to customize the analysis. If set is a GenomicRatioSet, the parameter betas allows us choosing between betas (TRUE) and M-values (FALSE). We can also run a robust linear model changing the parameter method to “robust”. Finally, resultSet indicates if the function will return a ResultSet (TRUE) or a MArrayLM (FALSE).

All these parameters can be set in the runPipeline function with the argument DiffMean_params.

3.4.2 Differences of Variance analysis

We can test if a phenotype causes changes in methylation variance using the runDiffVarAnalysis. This function is a wrapper of varFit function from missMethyl and requires three arguments: set, model and coefficient. set contains the methylation data in a GenomicRatioSet. model can be a matrix with the linear model or a formula indicating the model. In the former case, the variables included in the model must be present in the colData of our set. coefficient indicates the variables of the linear model for which the difference of variance will be computed. By default, all discrete variables will be included.

We exemplify the use of this function by running the same model than in our pipeline:

resDV <- runDiffVarAnalysis(set = meth, model = ~ status, coefficient = 2)

runDiffVarAnalysis also has the parameter resultSet that allows returning a MArrayLM object instead of a ResultSet. Finally, we can change other parameters of varFit function using the ... argument. These parameters can also be set in the runPipeline function passing them to the argument DiffVar_params.

3.4.3 Bumphunter

We can detect DMRs using Bumphunter from minfi with the function runBumphunter. This function requires three arguments: set, model and coefficient. set contains the methylation data in a GenomicRatioSet. model can be a matrix with the linear model or a formula indicating the model. In the former case, the variables included in the model must be present in the colData of our set. coefficient indicates the variable used to detect the DMRs.

We exemplify the use of this function by running bumphunter as in the pipeline:

resBH <- runBumphunter(set = meth, model = ~ status, coefficient = 2)

runBumphunter also has other parameters to customize the analysis. The parameter betas allows us choosing between betas (TRUE) and M-values (FALSE). bumphunter_cutoff specifies the minimum beta change to include a probe in a bump. num_permutations indicates the number of permutations run to compute bumps p-values (by default is 0 so no permutations are run and no p-values are returned). resultSet allows returning a data.frame object instead of a ResultSet. Finally, we can change other parameters of bumphunter function using the ... argument. These parameters can also be set in the runPipeline function passing them to the argument bumphunter_params.

3.4.4 Blockfinder

blockFinder is an adaptation of Bumphunter to detect DMRs from open sea probes. The function runBlockFinder has essentially the same arguments than runBumphunter.

We exemplify the use of this function by running blockFinder as in the pipeline:

resBF <- runBlockFinder(set = meth, model = ~ status, coefficient = 2)

To change the parameters in the runPipeline function, we can pass them to the argument blockFinder_params.

3.4.5 DMRcate

We can detect DMRs using DMRcate with the function runDMRcate. This function only has four parameters. set is the GenomicRatioSet, model is the linear model or a formula, coefficient is the variable used to detect the DMRs and resultSet to change the class of the output.

We exemplify the use of this function by running DMRcate as in the pipeline:

resDC <- runDMRcate(set = meth, model = ~ status, coefficient = 2)

We can change other parameters of DMRcate functions (cpg.annotate and dmrcate) passing them to the ... argument. These parameters can also be set in the runPipeline function passing them to the argument dmrcate_params.

3.4.6 RDA

We can determine if a genomic region is differentially methylated with RDA (Redundancy Analysis). This analysis can be run with the function runRDA that requires three arguments: set, model and range. As in the previous functions, set is a GenomicRatioSet with the methylation data and model contains the linear model either in a matrix or in a formula. range is a GenomicRanges with the coordinates of our target region.

We will exemplify the use of this function by running RDA in a region of chromosome X:

targetRange <- GRanges("chrX:13000000-23000000")
resRDA <- runRDA(set = meth, model = ~ status, range = targetRange)

runRDA also has other parameters to customize the analysis. The parameter betas allows us choosing between betas (TRUE) and M-values (FALSE). num_vars selects the number of columns in model matrix considered as variables. The remaining columns will be considered as covariates. num_permutations indicates the number of permutations run to compute p-values. resultSet allows returning a rda object from vegan package instead of a ResultSet.

We can run RDA in our pipeline when we are a priori interested in a target genomic range. In this case, we will pass our target region to the argument range of runPipeline. We can pass other parameters of runRDA using the argument rda_params.

3.4.6.1 Managing and plotting RDA results

We can retrieve RDA results using the function getAssociation:

getAssociation(resRDA, rid = "RDA")
## Call: rda(X = t(mat), Y = varsmodel, Z = covarsmodel)
## 
##                 Inertia Proportion Rank
## Total         1.403e+03  1.000e+00     
## Constrained   1.285e+02  9.157e-02    1
## Unconstrained 1.274e+03  9.084e-01    4
## Inertia is variance 
## Some constraints were aliased because they were collinear (redundant)
## 
## Eigenvalues for constrained axes:
##   RDA1 
## 128.45 
## 
## Eigenvalues for unconstrained axes:
##    PC1    PC2    PC3    PC4 
## 1090.0   96.6   62.8   24.9

RDA results are encapsulated in a rda object from vegan package. We can get a summary of RDA results with the function getRDAresults:

getRDAresults(resRDA)
##          R2        pval   global.R2 global.pval 
##  0.09157019  0.50000000  0.41556242  0.99010099

This function returns four values: R2, pval, global.R2 and global.pval. R2 is the ammount of variance that the model explains in our target region. pval is the probability of finding this ammount of variance of higher by change. global.R2 is the ammount of variance that our model explains in the whole genome. global.pval is the probability of finding a region with the same number of probes explaining the same or more variance than our target region. With these values, we can determine if our target region is differentially methylated and if this phenomena is local or global.

The function topRDAhits returns a data.frame with features associated to first two RDA components. This functions computes a Pearson correlation test between the methylation values and the RDA components. Only CpGs with a p-value lower than tPV parameter (by default 0.05) with any of the components are included in the data.frame:

topRDAhits(resRDA, tPV = 1e-17)
## [1] feat      RDA       cor       P.Value   adj.P.Val
## <0 rows> (or 0-length row.names)

Finally, we can plot the first two dimensions of our RDA with the function plotRDA. This function makes a biplot of samples and features. We can color the samples using categorical variables by passing in a data.frame to argument pheno.

We will plot RDA using status variable of our sets colData:

plotRDA(object = resRDA, pheno = colData(meth)[, "status", drop = FALSE])

The RDA plot prints a label at the center of each group and the summary of RDA results (R2 and p-value) in the legend. plotRDA has two additional arguments. main is a character vector with the plot’s title. n_feat is a numeric with the number of feats that will have a label in the text. Only the n_feat features most associated to each of the components will be displayed.

plotRDA relies on base paradigm, so we can add layers using functions from this infrastructure (e.g. lines, points…):

plotRDA(object = resRDA, pheno = colData(meth)[, "status", drop = FALSE])
abline(h = -1)

4 Session Info

sessionInfo()
## R version 3.5.0 (2018-04-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.4 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.7-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] doRNG_1.6.6                                       
##  [2] rngtools_1.3.1                                    
##  [3] pkgmaker_0.22                                     
##  [4] registry_0.5                                      
##  [5] ggplot2_2.2.1                                     
##  [6] minfiData_0.26.0                                  
##  [7] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0
##  [8] IlluminaHumanMethylation450kmanifest_0.4.0        
##  [9] minfi_1.26.0                                      
## [10] bumphunter_1.22.0                                 
## [11] locfit_1.5-9.1                                    
## [12] iterators_1.0.9                                   
## [13] foreach_1.4.4                                     
## [14] Biostrings_2.48.0                                 
## [15] XVector_0.20.0                                    
## [16] SummarizedExperiment_1.10.1                       
## [17] DelayedArray_0.6.0                                
## [18] BiocParallel_1.14.1                               
## [19] matrixStats_0.53.1                                
## [20] GenomicRanges_1.32.3                              
## [21] GenomeInfoDb_1.16.0                               
## [22] IRanges_2.14.10                                   
## [23] S4Vectors_0.18.2                                  
## [24] MEAL_1.10.1                                       
## [25] MultiDataSet_1.8.0                                
## [26] Biobase_2.40.0                                    
## [27] BiocGenerics_0.26.0                               
## [28] BiocStyle_2.8.1                                   
## 
## loaded via a namespace (and not attached):
##   [1] rms_5.1-2                                          
##   [2] R.utils_2.6.0                                      
##   [3] tidyselect_0.2.4                                   
##   [4] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
##   [5] RSQLite_2.1.1                                      
##   [6] AnnotationDbi_1.42.1                               
##   [7] htmlwidgets_1.2                                    
##   [8] grid_3.5.0                                         
##   [9] munsell_0.4.3                                      
##  [10] codetools_0.2-15                                   
##  [11] preprocessCore_1.42.0                              
##  [12] statmod_1.4.30                                     
##  [13] colorspace_1.3-2                                   
##  [14] knitr_1.20                                         
##  [15] rstudioapi_0.7                                     
##  [16] labeling_0.3                                       
##  [17] GenomeInfoDbData_1.1.0                             
##  [18] bit64_0.9-7                                        
##  [19] rhdf5_2.24.0                                       
##  [20] rprojroot_1.3-2                                    
##  [21] TH.data_1.0-8                                      
##  [22] xfun_0.1                                           
##  [23] qqman_0.1.4                                        
##  [24] biovizBase_1.28.0                                  
##  [25] R6_2.2.2                                           
##  [26] illuminaio_0.22.0                                  
##  [27] AnnotationFilter_1.4.0                             
##  [28] bitops_1.0-6                                       
##  [29] reshape_0.8.7                                      
##  [30] assertthat_0.2.0                                   
##  [31] scales_0.5.0                                       
##  [32] bsseq_1.16.0                                       
##  [33] multcomp_1.4-8                                     
##  [34] nnet_7.3-12                                        
##  [35] gtable_0.2.0                                       
##  [36] methylumi_2.26.0                                   
##  [37] sandwich_2.4-0                                     
##  [38] ensembldb_2.4.1                                    
##  [39] rlang_0.2.0                                        
##  [40] MatrixModels_0.4-1                                 
##  [41] genefilter_1.62.0                                  
##  [42] calibrate_1.7.2                                    
##  [43] splines_3.5.0                                      
##  [44] rtracklayer_1.40.2                                 
##  [45] lazyeval_0.2.1                                     
##  [46] acepack_1.4.1                                      
##  [47] DSS_2.28.0                                         
##  [48] GEOquery_2.48.0                                    
##  [49] dichromat_2.0-0                                    
##  [50] checkmate_1.8.5                                    
##  [51] yaml_2.1.19                                        
##  [52] snpStats_1.30.0                                    
##  [53] GenomicFeatures_1.32.0                             
##  [54] backports_1.1.2                                    
##  [55] Hmisc_4.1-1                                        
##  [56] tools_3.5.0                                        
##  [57] bookdown_0.7                                       
##  [58] nor1mix_1.2-3                                      
##  [59] RColorBrewer_1.1-2                                 
##  [60] siggenes_1.54.0                                    
##  [61] SNPassoc_1.9-2                                     
##  [62] Rcpp_0.12.17                                       
##  [63] plyr_1.8.4                                         
##  [64] base64enc_0.1-3                                    
##  [65] progress_1.1.2                                     
##  [66] zlibbioc_1.26.0                                    
##  [67] purrr_0.2.4                                        
##  [68] RCurl_1.95-4.10                                    
##  [69] BiasedUrn_1.07                                     
##  [70] prettyunits_1.0.2                                  
##  [71] rpart_4.1-13                                       
##  [72] openssl_1.0.1                                      
##  [73] zoo_1.8-1                                          
##  [74] ggrepel_0.8.0                                      
##  [75] cluster_2.0.7-1                                    
##  [76] magrittr_1.5                                       
##  [77] data.table_1.11.2                                  
##  [78] SparseM_1.77                                       
##  [79] mvtnorm_1.0-7                                      
##  [80] ProtGenerics_1.12.0                                
##  [81] missMethyl_1.14.0                                  
##  [82] hms_0.4.2                                          
##  [83] evaluate_0.10.1                                    
##  [84] xtable_1.8-2                                       
##  [85] XML_3.98-1.11                                      
##  [86] mclust_5.4                                         
##  [87] gridExtra_2.3                                      
##  [88] compiler_3.5.0                                     
##  [89] biomaRt_2.36.1                                     
##  [90] tibble_1.4.2                                       
##  [91] R.oo_1.22.0                                        
##  [92] htmltools_0.3.6                                    
##  [93] mgcv_1.8-23                                        
##  [94] Formula_1.2-3                                      
##  [95] tidyr_0.8.1                                        
##  [96] DBI_1.0.0                                          
##  [97] MASS_7.3-50                                        
##  [98] DMRcate_1.16.0                                     
##  [99] Matrix_1.2-14                                      
## [100] readr_1.1.1                                        
## [101] permute_0.9-4                                      
## [102] quadprog_1.5-5                                     
## [103] R.methodsS3_1.7.1                                  
## [104] Gviz_1.24.0                                        
## [105] bindr_0.1.1                                        
## [106] pkgconfig_2.0.1                                    
## [107] GenomicAlignments_1.16.0                           
## [108] haplo.stats_1.7.9                                  
## [109] foreign_0.8-70                                     
## [110] xml2_1.2.0                                         
## [111] annotate_1.58.0                                    
## [112] multtest_2.36.0                                    
## [113] beanplot_1.2                                       
## [114] ruv_0.9.7                                          
## [115] stringr_1.3.1                                      
## [116] VariantAnnotation_1.26.0                           
## [117] digest_0.6.15                                      
## [118] vegan_2.5-2                                        
## [119] rmarkdown_1.9                                      
## [120] base64_2.0                                         
## [121] htmlTable_1.11.2                                   
## [122] DelayedMatrixStats_1.2.0                           
## [123] curl_3.2                                           
## [124] Rsamtools_1.32.0                                   
## [125] gtools_3.5.0                                       
## [126] quantreg_5.35                                      
## [127] nlme_3.1-137                                       
## [128] bindrcpp_0.2.2                                     
## [129] Rhdf5lib_1.2.1                                     
## [130] limma_3.36.1                                       
## [131] BSgenome_1.48.0                                    
## [132] pillar_1.2.2                                       
## [133] lattice_0.20-35                                    
## [134] httr_1.3.1                                         
## [135] survival_2.42-3                                    
## [136] GO.db_3.6.0                                        
## [137] glue_1.2.0                                         
## [138] DMRcatedata_1.16.0                                 
## [139] bit_1.1-13                                         
## [140] stringi_1.2.2                                      
## [141] HDF5Array_1.8.0                                    
## [142] blob_1.1.1                                         
## [143] polspline_1.1.12                                   
## [144] org.Hs.eg.db_3.6.0                                 
## [145] latticeExtra_0.6-28                                
## [146] memoise_1.1.0                                      
## [147] dplyr_0.7.5

References

Aryee, Martin J, Andrew E Jaffe, Hector Corrada-Bravo, Christine Ladd-Acosta, Andrew P Feinberg, Kasper D Hansen, and Rafael A Irizarry. 2014. “Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays.” Bioinformatics (Oxford, England) 30 (10):1363–9. https://doi.org/10.1093/bioinformatics/btu049.

Du, Pan, Warren A Kibbe, and Simon M Lin. 2008. “lumi: a pipeline for processing Illumina microarray.” Bioinformatics (Oxford, England) 24 (13):1547–8. https://doi.org/10.1093/bioinformatics/btn224.