Example data

Start by loading the package and have the methylation data ready in the workspace. The minimal data requirements for this package to run, is a list of numerics, with sequentially ordered test statistics, on which the sliding window is used. In this example, we will be using an example data set to generate the needed input.

library(DMRScan)
data(DMRScan.methylationData) ## Load methylation data from chromosome 22, with 52018 CpGs measured
data(DMRScan.phenotypes) ## Load phenotype (end-point for methylation data)

Test statistics

Note that the DMRScan do not require raw data to be used, only test statistics from a CpG wise model. To generate test statistics for use in DMRScan, we use logistic regression on the example DNA methylation data, and the phenotype data as endpoint, but any other model could as well be used.

#observations <- apply(DMRScan.methylationData,1,function(x,y){
#                        summary(glm(y ~ x, 
#                        family = binomial(link = "logit")))$coefficients[2,3]}, 
#                                            y = DMRScan.phenotypes)

observations <- apply(DMRScan.methylationData,1,function(x,y){
                        summary(lm(x ~ y))$coefficients[2,3]}, 
                                           y = DMRScan.phenotypes)
head(observations)

## chr22.17061793 chr22.17062497 chr22.17062614 chr22.17063057 chr22.17063105 
##     0.02894454     0.65418440     0.85931406    -0.39210304    -0.73461053 
## chr22.17063107 
##     0.28545561

This will return a sequence of test statistics; which is the data that is used further. The output which will be used further is given below:

The sliding window function in DMRScan is not dependent on raw methylation data, but rather the test statistic for each genomic position. This allows for much grater flexibility when applying different models to the raw data, and can also be used with meta-analysis values, and not only primary studies.

Clustering of the test statics into “Chunks”

To apply the sliding window on a set of test statistics, we cluster them into chunks with not to much space between the probes. The maximum allowed distance between two probes in the same region is given by max.gap, and the minimum number of probes within a cluster is given by min.cpg.

The clustering is done automatically by make.cpg.cluster() and requires four inputs; - A sequence of test statistics - The corresponding genomic position for the test statistics - The maximum allowed gap (in base pairs) inside a sliding window, i.e. the maximum gap allowed in a cluster - The minimum number of probes in a cluster. This is set to 2 at default.

To identify the coordinate to each test statistic given in the previous table, we split the names inherited from the methylation values into its chromosomal part, and the base pair location.

pos <- matrix(as.integer(unlist(strsplit(names(observations),split="chr|[.]"))), 
                                                ncol = 3, byrow = TRUE)[,-1] 
head(pos)

##      [,1]     [,2]
## [1,]   22 17061793
## [2,]   22 17062497
## [3,]   22 17062614
## [4,]   22 17063057
## [5,]   22 17063105
## [6,]   22 17063107

The chromosome number and genomic position is here stored in the same object, a two column matrix. This is not needed, and if more convinient, the two peaces of information can be stored in separate objects.

To generate the clusters, we also need to set the additional clustering parameters:

## Minimum number of CpGs in a tested cluster 
min.cpg <- 3  

## Maxium distance (in base-pairs) within a cluster 
## before it is broken up into two seperate cluster
max.gap <- 750  

This, together with the test statistic vector from earlier, allows us to generate a list of clusters where the sliding window is operated.

regions <- makeCpGregions(observations  = observations, chr = pos[,1], pos = pos[,2], 
                                              maxGap = 750, minCpG = 3)

Calculating window thresholds

Three different approaches to estimate the window thresholds is given in the package. The different options are: - A full MCMC model with simulation of null model model = "mcmc", with given correlation structure for the null data, provided by arima.sim(). - An faster MCMC called “important sampling”, where only a subset of the data is sampled. This is implemented in model = "sampling", and is up to 700 times faster than the MCMC model, and simulation studies indicates that these to options are comparable in performance. - An analytical solution to calculate the window thresholds, based on Siegmund et.al (2012). Given by the model model = "siegmund". Since this is a closed form expression, it is much faster to calculate compared to the two other thresholds. However, this option tends to give more conservative thresholds, with somewhat lower power, but fewer false positive windows. An additional consideration with this method, is that it assumes that the test statistics follows an Ornstein-Uhlenbeck process, which may not always be the case.

To calculate the threshold(s), a few parameters needs to be set; - The method by which to estimate the threshold (see section above). - The window sizes. Can be either a single window size, or a sequence of window sizes. - The total number of observations in the study.

Additionally, for the Monte Carlo options, the number of iterations need to be specified, as well as the correlation structure. We will illustrate using only important sampling (which has a fixed correlation structure on the form of AR(2)) and the closed form expression from Siegmund et.al.

window.sizes <- 3:7 ## Number of CpGs in the sliding windows 
## (can be either a single number or a sequence)
n.CpG        <- sum(sapply(regions, length)) ## Number of CpGs to be tested
## Estimate the window threshold, based on the number of CpGs and window sizes 
## using important sampling
window.thresholds.importancSampling <- estimateWindowThreshold(nProbe = n.CpG, windowSize = window.sizes, 
                                                                    method = "sampling", mcmc = 10000)

## Estimating the window threshold using the closed form expression
window.thresholds.siegmund <- estimateWindowThreshold(nProbe = n.CpG, windowSize = window.sizes, 
                                                          method = "siegmund")

Identifying Differentially Methylated Regions

We now have all the data and parameters needed to run the dmrscan() function. First using a window threshold estimated with importanc sampling, we have

window.thresholds.importancSampling <- estimateWindowThreshold(nProbe = n.CpG, windowSize = window.sizes,
method = "sampling", mcmc = 10000)
dmrscan.results   <- dmrscan(observations = regions, windowSize = window.sizes, 
                                     windowThreshold = window.thresholds.importancSampling)
## Print the result
print(dmrscan.results)

## GRanges object with 6 ranges and 3 metadata columns:
##       seqnames            ranges strand |    no.cpg        pVal      tVal
##          <Rle>         <IRanges>  <Rle> | <integer>   <numeric> <numeric>
##   [1]       22 22874575-22874588      * |         3 1.21294e-04   2.73826
##   [2]       22 23801059-23801103      * |         5 5.69147e-04   2.34921
##   [3]       22 23801271-23801316      * |         6 8.67715e-04   2.20588
##   [4]       22 23801377-23801402      * |         3 4.66514e-05   2.99540
##   [5]       22 23801581-23801591      * |         3 3.73211e-05   3.32117
##   [6]       22 29704113-29704131      * |         5 4.85174e-04   2.39349
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

This will give the following result, with two genome wide significant regions.

When using a more stringent cut-off, generated by using the “siegmund” option in dmrscan(), we get no significant regions on the same data set. Exemplified by the syntax below:

dmrscan.results   <- dmrscan(observations = regions, windowSize = window.sizes, 
                                windowThreshold = window.thresholds.siegmund)
## Print the result
print(dmrscan.results)

## GRanges object with 1 range and 3 metadata columns:
##       seqnames            ranges strand |    no.cpg        pVal      tVal
##          <Rle>         <IRanges>  <Rle> | <integer>   <numeric> <numeric>
##   [1]       22 23801581-23801591      * |         3 3.73211e-05   3.32117
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths