1 Introduction

The missMethyl package contains functions to analyse methylation data from Illumina’s HumanMethylation450 and MethylationEPIC beadchip. These arrays are a cost-effective alternative to whole genome bisulphite sequencing, and as such are widely used to profile DNA methylation. Specifically, missMethyl contains functions to perform SWAN normalisation (Maksimovic, Gordon, and Oshlack 2012), perform differential methylation analysis using RUVm (Maksimovic et al. 2015), differential variability analysis (Phipson and Oshlack 2014) and gene set analysis (Phipson, Maksimovic, and Oshlack 2016). As our lab’s research into specialised analyses of these arrays continues we anticipate that the package will be continuously updated with new functions.

Raw data files are in IDAT format, which can be read into R using the minfi package (Aryee et al. 2014). Statistical analyses are usually performed on M-values, and \(\beta\) values are used for visualisation, both of which can be extracted from objects, which is a class of object created by minfi. For detecting differentially variable CpGs we recommend that the analysis is performed on M-values. All analyses described here are performed at the CpG site level.

2 Reading data into R

We will use the data in the minfiData package to demonstrate the functions in missMethyl. The example dataset has 6 samples across two slides. The sample information is in the targets file. An essential column in the targets file is the Basename column which tells where the idat files to be read in are located. The R commands to read in the data are taken from the minfi User’s Guide. For additional details on how to read the IDAT files into R, as well as information regarding quality control please refer to the minfi User’s Guide.

library(missMethyl)
library(limma)
library(minfi)

library(minfiData)
baseDir <- system.file("extdata", package = "minfiData")
targets <- read.metharray.sheet(baseDir)

## [1] "/home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/SampleSheet.csv"

targets[,1:9]

##   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
## 1 normal
## 2 normal
## 3 cancer
## 4 cancer
## 5 normal
## 6 cancer

targets[,10:12]

##    Array      Slide
## 1 R02C02 5723646052
## 2 R04C01 5723646052
## 3 R05C02 5723646052
## 4 R04C02 5723646053
## 5 R05C02 5723646053
## 6 R06C02 5723646053
##                                                                                Basename
## 1 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646052/5723646052_R02C02
## 2 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646052/5723646052_R04C01
## 3 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646052/5723646052_R05C02
## 4 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646053/5723646053_R04C02
## 5 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646053/5723646053_R05C02
## 6 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646053/5723646053_R06C02

rgSet <- read.metharray.exp(targets = targets)

The data is now an RGChannelSet object and needs to be normalised and converted to a MethylSet object.

3 Subset-quantile within array normalization (SWAN)

SWAN (subset-quantile within array normalization) is a within-array normalization method for Illumina 450k & EPIC BeadChips. Technical differencs have been demonstrated to exist between the Infinium I and Infinium II assays on a single Illumina HumanMethylation array (Bibikova et al. 2011, Dedeurwaerder, Defrance, and Calonne (2011)). Using the SWAN method substantially reduces the technical variability between the assay designs whilst maintaining important biological differences. The SWAN method makes the assumption that the number of CpGs within the 50bp probe sequence reflects the underlying biology of the region being interrogated. Hence, the overall distribution of intensities of probes with the same number of CpGs in the probe body should be the same regardless of assay type. The method then uses a subset quantile normalization approach to adjust the intensities of each array (Maksimovic, Gordon, and Oshlack 2012).

SWAN can take a MethylSet, RGChannelSet or MethyLumiSet as input. It should be noted that, in order to create the normalization subset, SWAN randomly selects Infinium I and II probes that have one, two and three underlying CpGs; as such, we recommend using set.seed before to ensure that the normalized intensities will be identical, if the normalization is repeated.

The technical differences between Infinium I and II assay designs can result in aberrant beta value distributions (Figure 1, panel “Raw”). Using SWAN corrects for the technical differences between the Infinium I and II assay designs and produces a smoother overall \(\beta\) value distribution (Figure 1, panel “SWAN”).

mSet <- preprocessRaw(rgSet)

mSetSw <- SWAN(mSet,verbose=TRUE)

## [SWAN] Preparing normalization subset
## 450k
## [SWAN] Normalizing methylated channel
## [SWAN] Normalizing array 1 of 6
## [SWAN] Normalizing array 2 of 6
## [SWAN] Normalizing array 3 of 6
## [SWAN] Normalizing array 4 of 6
## [SWAN] Normalizing array 5 of 6
## [SWAN] Normalizing array 6 of 6
## [SWAN] Normalizing unmethylated channel
## [SWAN] Normalizing array 1 of 6
## [SWAN] Normalizing array 2 of 6
## [SWAN] Normalizing array 3 of 6
## [SWAN] Normalizing array 4 of 6
## [SWAN] Normalizing array 5 of 6
## [SWAN] Normalizing array 6 of 6

par(mfrow=c(1,2), cex=1.25)
densityByProbeType(mSet[,1], main = "Raw")
densityByProbeType(mSetSw[,1], main = "SWAN")

Figure 1: Density distributions of \(eta\) values before and after using SWAN

4 Filter out poor quality probes

Poor quality probes can be filtered out based on the detection p-value. For this example, to retain a CpG for further analysis, we require that the detection p-value is less than 0.01 in all samples.

detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSetSw <- mSetSw[keep,]

5 Extracting Beta and M-values

Now that the data has been SWAN normalised we can extract \(\beta\) and M-values from the object. We prefer to add an offset to the methylated and unmethylated intensities when calculating M-values, hence we extract the methylated and unmethylated channels separately and perform our own calculation. For all subsequent analysis we use a random selection of 20000 CpGs to reduce computation time.

mset_reduced <- mSetSw[sample(1:nrow(mSetSw), 20000),]
meth <- getMeth(mset_reduced)
unmeth <- getUnmeth(mset_reduced)
Mval <- log2((meth + 100)/(unmeth + 100))
beta <- getBeta(mset_reduced)
dim(Mval)

## [1] 20000     6

par(mfrow=c(1,1))
plotMDS(Mval, labels=targets$Sample_Name, col=as.integer(factor(targets$status)))
legend("topleft",legend=c("Cancer","Normal"),pch=16,cex=1.2,col=1:2)

Figure 2: A multi-dimensional scaling (MDS) plot of cancer and normal samples

An MDS plot (Figure 2) is a good sanity check to make sure samples cluster together according to the main factor of interest, in this case, cancer and normal.

6 Testing for differential methylation using

To test for differential methylation we use the limma package (Smyth 2005), which employs an empirical Bayes framework based on Guassian model theory. First we need to set up the design matrix. There are a number of ways to do this, the most straightforward is directly from the targets file. There are a number of variables, with the status column indicating cancer/normal samples. From the person column of the targets file, we see that the cancer/normal samples are matched, with 3 individuals each contributing both a cancer and normal sample. Since the limma model framework can handle any experimental design which can be summarised by a design matrix, we can take into account the paired nature of the data in the analysis. For more complicated experimental designs, please refer to the limma User’s Guide.

group <- factor(targets$status,levels=c("normal","cancer"))
id <- factor(targets$person)
design <- model.matrix(~id + group)
design

##   (Intercept) idid2 idid3 groupcancer
## 1           1     0     1           0
## 2           1     1     0           0
## 3           1     0     1           1
## 4           1     0     0           1
## 5           1     0     0           0
## 6           1     1     0           1
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$id
## [1] "contr.treatment"
## 
## attr(,"contrasts")$group
## [1] "contr.treatment"

Now we can test for differential methylation using the lmFit and eBayes functions from limma. As input data we use the matrix of M-values.

fit.reduced <- lmFit(Mval,design)
fit.reduced <- eBayes(fit.reduced)

The numbers of hyper-methylated (1) and hypo-methylated (-1) can be displayed using the decideTests function in limma and the top 10 differentially methylated CpGs for cancer versus normal extracted using topTable.

summary(decideTests(fit.reduced))

##        (Intercept) idid2 idid3 groupcancer
## Down          7056     0   121         568
## NotSig        3422 20000 19870       18996
## Up            9522     0     9         436

top<-topTable(fit.reduced,coef=4)
top

##               logFC    AveExpr        t      P.Value  adj.P.Val        B
## cg14034197 4.575822 -0.3196148 18.25222 7.004140e-06 0.03505176 4.171036
## cg13181745 4.519766 -0.9046829 17.41519 8.891317e-06 0.03505176 4.019872
## cg07952047 4.561251 -0.4189628 16.76257 1.079351e-05 0.03505176 3.892259
## cg15044384 4.251783 -0.3890717 16.12421 1.314170e-05 0.03505176 3.758349
## cg11946336 4.026042 -1.3425857 15.52598 1.591437e-05 0.03505176 3.623979
## cg25537217 4.800006 -0.2738630 15.09534 1.834732e-05 0.03505176 3.521512
## cg04398581 3.755296  0.1218274 15.08863 1.838857e-05 0.03505176 3.519877
## cg21159370 4.070400  0.4561686 14.96422 1.917434e-05 0.03505176 3.489309
## cg01952234 3.760289 -1.0274718 14.67533 2.115849e-05 0.03505176 3.416629
## cg09011883 4.097902 -1.3080656 14.59797 2.173067e-05 0.03505176 3.396755

Note that since we performed our analysis on M-values, the logFC and AveExpr columns are computed on the M-value scale. For interpretability and visualisation we can look at the \(\beta\) values. The beta values for the top 4 differentially methylated CpGs shown in Figure 3.

cpgs <- rownames(top)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(beta[rownames(beta)==cpgs[i],]~design[,4],method="jitter",
group.names=c("Normal","Cancer"),pch=16,cex=1.5,col=c(4,2),ylab="Beta values",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(cpgs[i],cex.main=1.5)
}

Figure 3: The \(eta\) values for the top 4 differentially methylated CpGs

7 Removing unwanted variation when testing for differential methylation

Like other platforms, 450k array studies are subject to unwanted technical variation such as batch effects and other, often unknown, sources of variation. The adverse effects of unwanted variation have been extensively documented in gene expression array studies and have been shown to be able to both reduce power to detect true differences and to increase the number of false discoveries. As such, when it is apparent that data is significantly affected by unwanted variation, it is advisable to perform an adjustment to mitigate its effects.

missMethyl provides a limma-like interface to functions from the CRAN package ruv that enables the removal of unwanted variation when performing a differential analysis (Maksimovic et al. 2015). All of the methods rely on negative control features to accurately estimate the components of unwanted variation. Negative control features are probes/genes/etc. that are known a priori to not truly be associated with the biological factor of interest, but are affected by unwanted variation. For example, in a microarray gene expression study, these could be house-keeping genes or a set of spike-in controls. Negative control features are extensively discussed in Gagnon-Bartsch and Speed (2012) and Gagnon-Bartsch et al. (2013). Once the unwanted factors are accurately estimated from the data, they are adjusted for in the linear model that describes the differential analysis.

If the negative control features are not known a priori, they can be identified empirically. This can be achieved using a 2-stage approach, RUVm, based on RUV-inverse. Stage 1 involves performing a differential methylation analysis using RUV-inverse and the 613 Illumina negative controls (INCs) as negative control features. This will produce a list of CpGs ranked by p-value according to their level of association with the factor of interest. This list can then be used to identify a set of empirical control probes (ECPs), which will capture more of the unwanted variation than using the INCs alone. ECPs are selected by designating a proportion of the CpGs least associated with the factor of interest as negative control features; this can be done based on either an FDR cut-off or by taking a fixed percentage of probes from the bottom of the ranked list. Stage 2 involves performing a second differential methylation analysis on the original data using RUV-inverse and the ECPs. For simplicity, we are ignoring the paired nature of the cancer and normal samples in this example.

# get M-values for ALL probes
meth <- getMeth(mSet)
unmeth <- getUnmeth(mSet)
M <- log2((meth + 100)/(unmeth + 100))

grp <- factor(targets$status,levels=c("normal","cancer"))
des <- model.matrix(~grp)
des

##   (Intercept) grpcancer
## 1           1         0
## 2           1         0
## 3           1         1
## 4           1         1
## 5           1         0
## 6           1         1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$grp
## [1] "contr.treatment"

INCs <- getINCs(rgSet)
head(INCs)

##          5723646052_R02C02 5723646052_R04C01 5723646052_R05C02
## 13792480        -0.3299654        -1.0955482        -0.5266103
## 69649505        -1.0354488        -1.4943396        -1.0067050
## 34772371        -1.1286422        -0.2995603        -0.8192636
## 28715352        -0.5553373        -0.7599489        -0.7186973
## 74737439        -1.1169178        -0.8656399        -0.6429681
## 33730459        -0.7714684        -0.5622424        -0.7724825
##          5723646053_R04C02 5723646053_R05C02 5723646053_R06C02
## 13792480        -0.6374299         -1.116598        -0.4332793
## 69649505        -0.8854881         -1.586679        -0.9217329
## 34772371        -0.6895514         -1.161155        -0.6186795
## 28715352        -1.7903619         -1.348105        -1.0067259
## 74737439        -0.8872082         -1.064986        -0.9841833
## 33730459        -1.5623138         -2.079184        -1.0445246

Mc <- rbind(M,INCs)
ctl <- rownames(Mc) %in% rownames(INCs)
table(ctl)

## ctl
##  FALSE   TRUE 
## 485512    613

rfit1 <- RUVfit(data=Mc, design=des, coef=2, ctl=ctl) # Stage 1 analysis
rfit2 <- RUVadj(rfit1)

Now that we have performed an initial differential methylation analysis to rank the CpGs with respect to their association with the factor of interest, we can designate the CpGs that are least associated with the factor of interest based on FDR-adjusted p-value as ECPs.

top1 <- topRUV(rfit2, num=Inf)
head(top1)

##                  X1       X1           X1      p.BH     p.ebayes p.ebayes.BH
## cg04743961 4.838190 26.74467 3.812882e-05 0.1401969 3.516091e-07  0.01017357
## cg07155336 5.887409 17.62103 1.608653e-04 0.1401969 3.583107e-07  0.01017357
## cg20925841 4.790211 26.69524 3.837354e-05 0.1401969 3.730375e-07  0.01017357
## cg03607359 4.394397 34.74068 1.542013e-05 0.1401969 4.721205e-07  0.01017357
## cg10566121 4.787914 21.80693 7.717708e-05 0.1401969 5.238865e-07  0.01017357
## cg07655636 4.571758 22.99708 6.424216e-05 0.1401969 6.080091e-07  0.01017357

ctl <- rownames(M) %in% rownames(top1[top1$p.ebayes.BH > 0.5,])
table(ctl)

## ctl
##  FALSE   TRUE 
## 172540 312972

We can then use the ECPs to perform a second differential methylation with RUV-inverse, which is adjusted for the unwanted variation estimated from the data.

# Perform RUV adjustment and fit
rfit1 <- RUVfit(data=M, design=des, coef=2, ctl=ctl) # Stage 2 analysis
rfit2 <- RUVadj(rfit1)

# Look at table of top results
topRUV(rfit2)

##                  X1        X1          X1      p.BH     p.ebayes
## cg07155336 5.769286 15.345069 0.002005546 0.3431163 1.434834e-55
## cg06463958 5.733093 15.434797 0.001978272 0.3431163 6.749298e-55
## cg00024472 5.662959 15.946200 0.001832444 0.3431163 1.319390e-53
## cg02040433 5.651399 10.054445 0.005389436 0.3431163 2.146210e-53
## cg13355248 5.595396  9.963702 0.005504213 0.3431163 2.234891e-52
## cg02467990 5.592707  6.859614 0.013008521 0.3431163 2.499534e-52
## cg00817367 5.527501 13.070583 0.002921656 0.3431163 3.710480e-51
## cg11396157 5.487992 10.931263 0.004436178 0.3431163 1.873636e-50
## cg16306898 5.466780  5.573935 0.020790127 0.3431163 4.448085e-50
## cg03735888 5.396242 15.482605 0.001963955 0.3431163 7.700032e-49
##             p.ebayes.BH
## cg07155336 6.966293e-50
## cg06463958 1.638433e-49
## cg00024472 2.135266e-48
## cg02040433 2.605027e-48
## cg13355248 2.022589e-47
## cg02467990 2.022589e-47
## cg00817367 2.573547e-46
## cg11396157 1.137091e-45
## cg16306898 2.399554e-45
## cg03735888 3.738458e-44

Note, at present does not support contrasts, so only one factor of interest can be interrogated at a time using a design matrix with an intercept term.

8 Testing for differential variability (DiffVar)

8.1 Methylation data

Rather than testing for differences in mean methylation, we may be interested in testing for differences between group variances. For example, it has been hypothesised that highly variable CpGs in cancer are important for tumour progression (Hansen et al. 2011). Hence we may be interested in CpG sites that are consistently methylated in the normal samples, but variably methylated in the cancer samples.

In general we recommend at least 10 samples in each group for accurate variance estimation, however for the purpose of this vignette we perform the analysis on 3 vs 3. In this example, we are interested in testing for differential variability in the cancer versus normal group. Note that when we specify the coef parameter, which corresponds to the columns of the design matrix to be used for testing differential variability, we need to specify both the intercept and the fourth column. The ID variable is a nuisance parameter and not used when obtaining the absolute deviations, however it can be included in the linear modelling step. For methylation data, the function will take either a matrix of M-values, \(\beta\) values or a object as input. If \(\beta\) values are supplied, a logit transformation is performed. Note that as a default, varFit uses the robust setting in the limma framework, which requires the use of the statmod package.

fitvar <- varFit(Mval, design = design, coef = c(1,4))

The numbers of hyper-variable (1) and hypo-variable (-1) genes in cancer vs normal can be obtained using decideTests.

summary(decideTests(fitvar))

##        (Intercept) idid2 idid3 groupcancer
## Down             0     2     4           1
## NotSig       19615 19996 19955       19996
## Up             385     2    41           3

topDV <- topVar(fitvar, coef=4)
topDV

##            SampleVar LogVarRatio DiffLevene         t      P.Value
## cg05593887 13.742723  -0.4103645 -0.8161431 -5.255478 1.483965e-07
## cg23071808  5.551705   3.3795056  2.6509302  5.199805 2.004777e-07
## cg17942639  5.095744   4.0407816  2.6686110  4.745393 2.088204e-06
## cg17484671  4.590943   6.0095322  2.9026688  4.730653 2.245616e-06
## cg09770278  4.817211   2.2606740  2.0574707  4.261769 2.032784e-05
## cg21211367  4.499561   5.5939653  2.7003369  4.259559 2.052983e-05
## cg14479889  7.494275   3.2845218  2.3881498  4.242851 2.211936e-05
## cg18222083  5.066301   3.7112902  2.4550876  4.206153 2.603191e-05
## cg04186360  5.383690   2.2348168  2.0435622  4.200815 2.665310e-05
## cg23543123  4.617776   6.1347668  2.6889606  4.182078 2.894711e-05
##            Adj.P.Value
## cg05593887 0.002004777
## cg23071808 0.002004777
## cg17942639 0.011228080
## cg17484671 0.011228080
## cg09770278 0.057894223
## cg21211367 0.057894223
## cg14479889 0.057894223
## cg18222083 0.057894223
## cg04186360 0.057894223
## cg23543123 0.057894223

An alternate parameterisation of the design matrix that does not include an intercept term can also be used, and specific contrasts tested with contrasts.varFit. Here we specify the design matrix such that the first two columns correspond to the normal and cancer groups, respectively.

design2 <- model.matrix(~0+group+id)
fitvar.contr <- varFit(Mval, design=design2, coef=c(1,2))
contr <- makeContrasts(groupcancer-groupnormal,levels=colnames(design2))
fitvar.contr <- contrasts.varFit(fitvar.contr,contrasts=contr)

The results are identical to before.

summary(decideTests(fitvar.contr))

##        groupcancer - groupnormal
## Down                           1
## NotSig                     19996
## Up                             3

topVar(fitvar.contr,coef=1)

##            SampleVar LogVarRatio DiffLevene         t      P.Value
## cg05593887 13.742723  -0.4103645 -0.8161431 -5.255478 1.483965e-07
## cg23071808  5.551705   3.3795056  2.6509302  5.199805 2.004777e-07
## cg17942639  5.095744   4.0407816  2.6686110  4.745393 2.088204e-06
## cg17484671  4.590943   6.0095322  2.9026688  4.730653 2.245616e-06
## cg09770278  4.817211   2.2606740  2.0574707  4.261769 2.032784e-05
## cg21211367  4.499561   5.5939653  2.7003369  4.259559 2.052983e-05
## cg14479889  7.494275   3.2845218  2.3881498  4.242851 2.211936e-05
## cg18222083  5.066301   3.7112902  2.4550876  4.206153 2.603191e-05
## cg04186360  5.383690   2.2348168  2.0435622  4.200815 2.665310e-05
## cg23543123  4.617776   6.1347668  2.6889606  4.182078 2.894711e-05
##            Adj.P.Value
## cg05593887 0.002004777
## cg23071808 0.002004777
## cg17942639 0.011228080
## cg17484671 0.011228080
## cg09770278 0.057894223
## cg21211367 0.057894223
## cg14479889 0.057894223
## cg18222083 0.057894223
## cg04186360 0.057894223
## cg23543123 0.057894223

The \(\beta\) values for the top 4 differentially variable CpGs can be seen in Figure 4.

cpgsDV <- rownames(topDV)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(beta[rownames(beta)==cpgsDV[i],]~design[,4],method="jitter",
group.names=c("Normal","Cancer"),pch=16,cex=1.5,col=c(4,2),ylab="Beta values",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(cpgsDV[i],cex.main=1.5)
}

Figure 4: The \(eta\) values for the top 4 differentially variable CpGs

8.2 RNA-Seq expression data

Testing for differential variability in expression data is straightforward if the technology is gene expression microarrays. The matrix of expression values can be supplied directly to the varFit function. For RNA-Seq data, the mean-variance relationship that occurs in count data needs to be taken into account. In order to deal with this issue, we apply a voom transformation (Law et al. 2014) to obtain observation weights, which are then used in the linear modelling step. For RNA-Seq data, the varFit function will take a DGElist object as input.

To demonstrate this, we use data from the tweeDEseqCountData package. This data is part of the International HapMap project, consisting of RNA-Seq profiles from 69 unrelated Nigerian individuals (Pickrell et al. 2010). The only covariate is gender, so we can look at differentially variable expression between males and females. We follow the code from the limma vignette to read in and process the data before testing for differential variability.

First we load up the data and extract the relevant information.

library(tweeDEseqCountData)
data(pickrell1)
counts<-exprs(pickrell1.eset)
dim(counts)

## [1] 38415    69

gender <- pickrell1.eset$gender
table(gender)

## gender
## female   male 
##     40     29

rm(pickrell1.eset)
data(genderGenes)
data(annotEnsembl63)
annot <- annotEnsembl63[,c("Symbol","Chr")]
rm(annotEnsembl63)

We now have the counts, gender of each sample and annotation (gene symbol and chromosome) for each Ensemble gene. We can form a DGElist object using the edgeR package.

library(edgeR)
y <- DGEList(counts=counts, genes=annot[rownames(counts),])

We filter out lowly expressed genes by keeping genes with at least 1 count per million reads in at least 20 samples, as well as genes that have defined annotation. Finally we perform scaling normalisation.

isexpr <- rowSums(cpm(y)>1) >= 20
hasannot <- rowSums(is.na(y$genes))==0
y <- y[isexpr & hasannot,,keep.lib.sizes=FALSE]
dim(y)

## [1] 17310    69

y <- calcNormFactors(y)

We set up the design matrix and test for differential variability. In this case there are no nuisance parameters, so coef does not need to be explicitly specified.

design.hapmap <- model.matrix(~gender)
fitvar.hapmap <- varFit(y, design = design.hapmap)

## Converting counts to log counts-per-million using voom.

fitvar.hapmap$genes <- y$genes

We can display the results of the test:

summary(decideTests(fitvar.hapmap))

##        (Intercept) gendermale
## Down             0          2
## NotSig           0      17308
## Up           17310          0

topDV.hapmap <- topVar(fitvar.hapmap,coef=ncol(design.hapmap))
topDV.hapmap

##                       Symbol Chr  SampleVar LogVarRatio DiffLevene         t
## ENSG00000213318 RP11-331F4.1  16 5.69839463   -2.562939 -0.9859943 -8.031243
## ENSG00000129824       RPS4Y1   Y 2.32497726   -2.087025 -0.4585620 -4.957005
## ENSG00000233864       TTTY15   Y 6.79004140   -2.245369 -0.6085233 -4.612934
## ENSG00000176171        BNIP3  10 0.41317384    1.199292  0.3632133  4.219404
## ENSG00000197358      BNIP3P1  14 0.39969125    1.149754  0.3353288  4.058147
## ENSG00000025039        RRAGD   6 0.91837213    1.091229  0.4926839  3.977022
## ENSG00000103671        TRIP4  15 0.07456448   -1.457139 -0.1520583 -3.911300
## ENSG00000171100         MTM1   X 0.44049558   -1.133295 -0.3334619 -3.896490
## ENSG00000149476          DAK  11 0.13289523   -1.470460 -0.1919880 -3.785893
## ENSG00000064886       CHI3L2   1 2.70234584    1.468059  0.8449434  3.782010
##                      P.Value  Adj.P.Value
## ENSG00000213318 7.238039e-12 1.252905e-07
## ENSG00000129824 3.960855e-06 3.428120e-02
## ENSG00000233864 1.496237e-05 8.633290e-02
## ENSG00000176171 6.441668e-05 2.787632e-01
## ENSG00000197358 1.147886e-04 3.973982e-01
## ENSG00000025039 1.527695e-04 4.375736e-01
## ENSG00000103671 1.921104e-04 4.375736e-01
## ENSG00000171100 2.022293e-04 4.375736e-01
## ENSG00000149476 2.956364e-04 4.425050e-01
## ENSG00000064886 2.995692e-04 4.425050e-01

The log counts per million for the top 4 differentially variable genes can be seen in Figure 5.

genesDV <- rownames(topDV.hapmap)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(cpm(y,log=TRUE)[rownames(y)==genesDV[i],]~design.hapmap[,ncol(design.hapmap)],method="jitter",
group.names=c("Female","Male"),pch=16,cex=1.5,col=c(4,2),ylab="Log counts per million",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(genesDV[i],cex.main=1.5)
}

Figure 5: The log counts per million for the top 4 differentially variably expressed genes

9 Gene ontology analysis

Once a differential methylation or differential variability analysis has been performed, it may be of interest to know which gene pathways are targeted by the significant CpG sites. It is not entirely clear from the literature how best to perform such an analysis, however Geeleher et al. (Geeleher et al. 2013) showed there is a severe bias when performing gene ontology analysis with methylation data. This is due to the fact that there are differing numbers of probes per gene on several different array technologies. For the Illumina Infinium HumanMethylation450 array the number of probes per gene ranges from 1 to 1299, with a median of 15 probes per gene. For the EPIC array, the range is 1 to 1487, with a median of 20 probes per gene. This means that when mapping CpG sites to genes, a gene is more likely to be selected if there are many CpG sites associated with the gene.

One way to take into account this selection bias is to model the relationship between the number of probes per gene and the probability of being selected. This can be performed by adapting the goseq method of Young et al. (Young et al. 2010). Each gene then has a prior probability associated with it, and a modified version of a hypergeometric test can be performed, testing for over-representation of the selected genes in each gene set.

The gometh function performs gene set testing on GO categories or KEGG pathways (Phipson, Maksimovic, and Oshlack 2016). The gsameth function is a more generalised gene set testing function which can take as input a list of user specified gene sets. Note that for gsameth, the format for the gene ids for each gene in the gene set needs to be Entrez Gene IDs. For example, the entire curated gene set list (C2) from the Broad’s Molecular Signatures Database can be specified as input. The R version of these lists can be downloaded from http://bioinf.wehi.edu.au/software/MSigDB/index.html. Both functions take a vector of significant CpG probe names as input.

To illustrate how to use gometh, consider the results from the differential methylation analysis with RUVm.

topRUV(rfit2)

##                  X1        X1          X1      p.BH     p.ebayes
## cg07155336 5.769286 15.345069 0.002005546 0.3431163 1.434834e-55
## cg06463958 5.733093 15.434797 0.001978272 0.3431163 6.749298e-55
## cg00024472 5.662959 15.946200 0.001832444 0.3431163 1.319390e-53
## cg02040433 5.651399 10.054445 0.005389436 0.3431163 2.146210e-53
## cg13355248 5.595396  9.963702 0.005504213 0.3431163 2.234891e-52
## cg02467990 5.592707  6.859614 0.013008521 0.3431163 2.499534e-52
## cg00817367 5.527501 13.070583 0.002921656 0.3431163 3.710480e-51
## cg11396157 5.487992 10.931263 0.004436178 0.3431163 1.873636e-50
## cg16306898 5.466780  5.573935 0.020790127 0.3431163 4.448085e-50
## cg03735888 5.396242 15.482605 0.001963955 0.3431163 7.700032e-49
##             p.ebayes.BH
## cg07155336 6.966293e-50
## cg06463958 1.638433e-49
## cg00024472 2.135266e-48
## cg02040433 2.605027e-48
## cg13355248 2.022589e-47
## cg02467990 2.022589e-47
## cg00817367 2.573547e-46
## cg11396157 1.137091e-45
## cg16306898 2.399554e-45
## cg03735888 3.738458e-44

table(rfit2$p.ebayes.BH < 0.01)

## 
##  FALSE   TRUE 
## 424168  61344

At a 1% false discovery rate cut-off, there are still tens of thousands of CpG sites differentially methylated. These will undoubtably map to almost all the genes in the genome, making a gene ontology analysis irrelevant. One option for selecting CpGs in this context is to apply not only a false discovery rate cut-off, but also a \(\Delta\beta\) cut-off. However, for this dataset, taking a relatively large \(\Delta\beta\) cut-off of 0.25 still leaves more than 30000 CpGs differentially methylated.

beta <- getBeta(mSet)
beta_norm <- rowMeans(beta[,des[,2]==0])
beta_can <- rowMeans(beta[,des[,2]==1])
Delta_beta <- beta_can - beta_norm
sigDM <- rfit2$p.ebayes.BH < 0.01 & abs(Delta_beta) > 0.25
table(sigDM)

## sigDM
##  FALSE   TRUE 
## 451760  33748

Instead, we take the top 10000 CpG sites as input to gometh.

topCpGs<-topRUV(rfit2,number=10000)
sigCpGs <- rownames(topCpGs)
sigCpGs[1:10]

##  [1] "cg07155336" "cg06463958" "cg00024472" "cg02040433" "cg13355248"
##  [6] "cg02467990" "cg00817367" "cg11396157" "cg16306898" "cg03735888"

The takes as input a character vector of CpG names, and optionally, a character vector of all CpG sites tested. If the all.cpg argument is omitted, all the CpGs on the array are used as background. To change the array type, the array.type argument can be specified as either “450K” or “EPIC”. The default is “450K”.

If the plot.bias argument is TRUE, a figure showing the relationship between the probability of being selected and the number of probes per gene will be displayed.

For testing of GO terms, the collection argument takes the value “GO”, which is the default setting. For KEGG pathway analysis, set collection to “KEGG”. The function topGO shows the top enriched GO categories. For KEGG testing, use the topKEGG function. The functions goana and kegga are called by for GO and KEGG pathway analysis respectively.

For GO testing on our example dataset:

library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
gst <- gometh(sig.cpg=sigCpGs, all.cpg=rownames(rfit2), collection="GO")

## Warning in alias2SymbolTable(flat$symbol): Multiple symbols ignored for one
## or more aliases

topGO(gst)

##                                              Term Ont    N   DE         P.DE
## GO:0007399             nervous system development  BP 2273  596 8.334148e-39
## GO:0048731                     system development  BP 4678  955 7.281790e-38
## GO:0048856       anatomical structure development  BP 5738 1107 9.157918e-38
## GO:0007275     multicellular organism development  BP 5262 1033 1.857692e-36
## GO:0032502                  developmental process  BP 6148 1148 1.389122e-34
## GO:0032501       multicellular organismal process  BP 7286 1260 2.587120e-34
## GO:0005886                        plasma membrane  CC 5145  926 1.070386e-27
## GO:0030154                   cell differentiation  BP 4042  796 6.657713e-27
## GO:0071944                         cell periphery  CC 5253  939 7.994505e-27
## GO:0022008                           neurogenesis  BP 1521  409 1.356500e-26
## GO:0007417     central nervous system development  BP  925  280 1.521804e-26
## GO:0048869         cellular developmental process  BP 4239  823 2.211930e-26
## GO:0048699                  generation of neurons  BP 1422  389 2.576657e-26
## GO:0030182                 neuron differentiation  BP 1289  361 3.975027e-26
## GO:0009653     anatomical structure morphogenesis  BP 2631  589 3.034001e-25
## GO:0031226 intrinsic component of plasma membrane  CC 1707  380 1.545682e-24
## GO:0005887  integral component of plasma membrane  CC 1645  367 4.404428e-24
## GO:0007267                    cell-cell signaling  BP 1648  399 1.833143e-23
## GO:0031224        intrinsic component of membrane  CC 5212  858 8.001879e-23
## GO:0048513               animal organ development  BP 3371  663 7.181310e-22
##                     FDR
## GO:0007399 1.867349e-34
## GO:0048731 6.839744e-34
## GO:0048856 6.839744e-34
## GO:0007275 1.040586e-32
## GO:0032502 6.224932e-31
## GO:0032501 9.661170e-31
## GO:0005886 3.426152e-24
## GO:0030154 1.864659e-23
## GO:0071944 1.990276e-23
## GO:0022008 3.039373e-23
## GO:0007417 3.099776e-23
## GO:0048869 4.130041e-23
## GO:0048699 4.440968e-23
## GO:0030182 6.361748e-23
## GO:0009653 4.531989e-22
## GO:0031226 2.164535e-21
## GO:0005887 5.805036e-21
## GO:0007267 2.281856e-20
## GO:0031224 9.436321e-20
## GO:0048513 8.045222e-19

For a more generalised version of gene set testing in methylation data where the user can specify the gene set to be tested, the function gsameth can be used. To display the top 20 pathways, topGSA can be called. gsameth can take a single gene set, or a list of gene sets. The gene identifiers in the gene set must be Entrez Gene IDs. To demonstrate gsameth, a toy example is shown below, with gene sets made up of randomly selected genes from the org.Hs.eg.db package.

library(org.Hs.eg.db)

## Loading required package: AnnotationDbi

genes <- toTable(org.Hs.egSYMBOL2EG)
set1 <- sample(genes$gene_id,size=80)
set2 <- sample(genes$gene_id,size=100)
set3 <- sample(genes$gene_id,size=30)
genesets <- list(set1,set2,set3)
gsa <- gsameth(sig.cpg=sigCpGs, all.cpg=rownames(rfit2), collection=genesets)

## Warning in alias2SymbolTable(flat$symbol): Multiple symbols ignored for one
## or more aliases

topGSA(gsa)

##        N DE       P.DE       FDR
## [1,]  30  3 0.04968364 0.1388094
## [2,]  80  5 0.09253957 0.1388094
## [3,] 100  2 0.85641320 0.8564132

Note that if it is of interest to obtain the Entrez Gene IDs that the significant CpGs are mapped to, the getMappedEntrezIDs can be called.

10 Session information

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] org.Hs.eg.db_3.6.0
[2] AnnotationDbi_1.42.0
[3] edgeR_3.22.0
[4] tweeDEseqCountData_1.17.0
[5] minfiData_0.25.0
[6] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0 [7] IlluminaHumanMethylation450kmanifest_0.4.0
[8] minfi_1.26.0
[9] bumphunter_1.22.0
[10] locfit_1.5-9.1
[11] iterators_1.0.9
[12] foreach_1.4.4
[13] Biostrings_2.48.0
[14] XVector_0.20.0
[15] SummarizedExperiment_1.10.0
[16] DelayedArray_0.6.0
[17] BiocParallel_1.14.0
[18] matrixStats_0.53.1
[19] Biobase_2.40.0
[20] GenomicRanges_1.32.0
[21] GenomeInfoDb_1.16.0
[22] IRanges_2.14.0
[23] S4Vectors_0.18.0
[24] BiocGenerics_0.26.0
[25] limma_3.36.0
[26] missMethyl_1.14.0
[27] BiocStyle_2.8.0

loaded via a namespace (and not attached): [1] colorspace_1.3-2
[2] siggenes_1.54.0
[3] mclust_5.4
[4] rprojroot_1.3-2
[5] base64_2.0
[6] bit64_0.9-7
[7] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0 [8] xml2_1.2.0
[9] codetools_0.2-15
[10] splines_3.5.0
[11] methylumi_2.26.0
[12] knitr_1.20
[13] Rsamtools_1.32.0
[14] annotate_1.58.0
[15] GO.db_3.6.0
[16] HDF5Array_1.8.0
[17] readr_1.1.1
[18] compiler_3.5.0
[19] httr_1.3.1
[20] backports_1.1.2
[21] assertthat_0.2.0
[22] Matrix_1.2-14
[23] lazyeval_0.2.1
[24] htmltools_0.3.6
[25] prettyunits_1.0.2
[26] tools_3.5.0
[27] bindrcpp_0.2.2
[28] gtable_0.2.0
[29] glue_1.2.0
[30] GenomeInfoDbData_1.1.0
[31] dplyr_0.7.4
[32] doRNG_1.6.6
[33] Rcpp_0.12.16
[34] multtest_2.36.0
[35] preprocessCore_1.42.0
[36] nlme_3.1-137
[37] rtracklayer_1.40.0
[38] DelayedMatrixStats_1.2.0
[39] xfun_0.1
[40] stringr_1.3.0
[41] rngtools_1.2.4
[42] statmod_1.4.30
[43] XML_3.98-1.11
[44] beanplot_1.2
[45] zlibbioc_1.26.0
[46] MASS_7.3-50
[47] scales_0.5.0
[48] hms_0.4.2
[49] rhdf5_2.24.0
[50] GEOquery_2.48.0
[51] RColorBrewer_1.1-2
[52] yaml_2.1.18
[53] memoise_1.1.0
[54] gridExtra_2.3
[55] ggplot2_2.2.1
[56] pkgmaker_0.22
[57] biomaRt_2.36.0
[58] reshape_0.8.7
[59] stringi_1.1.7
[60] RSQLite_2.1.0
[61] highr_0.6
[62] genefilter_1.62.0
[63] GenomicFeatures_1.32.0
[64] rlang_0.2.0
[65] pkgconfig_2.0.1
[66] bitops_1.0-6
[67] nor1mix_1.2-3
[68] evaluate_0.10.1
[69] lattice_0.20-35
[70] purrr_0.2.4
[71] Rhdf5lib_1.2.0
[72] bindr_0.1.1
[73] ruv_0.9.7
[74] GenomicAlignments_1.16.0
[75] bit_1.1-12
[76] plyr_1.8.4
[77] magrittr_1.5
[78] bookdown_0.7
[79] R6_2.2.2
[80] DBI_0.8
[81] pillar_1.2.2
[82] survival_2.42-3
[83] RCurl_1.95-4.10
[84] tibble_1.4.2
[85] rmarkdown_1.9
[86] progress_1.1.2
[87] grid_3.5.0
[88] data.table_1.10.4-3
[89] blob_1.1.1
[90] digest_0.6.15
[91] xtable_1.8-2
[92] tidyr_0.8.0
[93] illuminaio_0.22.0
[94] openssl_1.0.1
[95] munsell_0.4.3
[96] registry_0.5
[97] BiasedUrn_1.07
[98] quadprog_1.5-5