A basic task in the analysis of count data from Whole Genome Bisulfite-Sequencing is the detection of differentially methylated regions. The count data consist of, for each sample, the number of methylated reads and the total number of reads covering CpG. An important analysis question is to detect regions (collections of neighboring CpGs) with systematic differences between conditions, as compared to within-condition variability. These so-called Differentially Methylated Regions (DMRs) are thought to be more informative than single CpGs in terms of of biological function. Although several methods exist to quantify and perform statistical inference on changes at the individual CpG level, detection of DMRs is still limited to aggregating signifiant CpGs without proper inference at the region level. The package **dmrseq** addresses this gap by providing a rigorous permutation-based approach to detect and perform inference for differential methylation by use of generalized least squares models that account for inter-individual and inter-CpG variability to generate region-level statistics that can be comparable across the genome. This allows the framework to perform well even on samples as small as two per group. This vignette explains the use of the package and demonstrates typical workflows. This vignette was generated with dmrseq package version 1.20.0

dmrseq 1.20.0

**If you use dmrseq in published research, please cite:**

Korthauer, K., Chakraborty, S., Benjamini, Y., and Irizarry, R.A. Detection and accurate False Discovery Rate control of differentially methylated regions from Whole Genome Bisulfite Sequencing

Biostatistics, 2018 (in press).

This package builds upon the
bsseq package (Hansen, Langmead, and Irizarry 2012),
which provides efficient storage and manipulation of bisulfite
sequencing data and inference for differentially methylated CpGs.
The main goal of **dmrseq** (Korthauer et al. 2018)
is to provide inference for differentially methylated *regions*, or
groups of CpGs.

Here we show the most basic steps for a differential methylation
analysis. There are a variety of steps upstream of **dmrseq** that result
in the generation of counts of methylated reads and total reads covering each
CpG for each sample, including mapping of sequencing reads to a reference
genome with and without bisulfite conversion. You can use the software
of your preference for this step (one option is
Bismark), as
long as you are able to obtain counts of methylation and coverage (as
opposed to solely methylation proportions, as discussed below).

This package uses a specific data structure to store and manipulate
bisulfite sequencing data introduced by the **bsseq** package. This data
structure is a *class* called `BSseq`

. Objects of the class `BSseq`

contain
all pertinent information for a bisulfite sequencing experiment, including
the number of reads corresponding to methylation, and the total number
of reads at each
CpG site, the location of each CpG site, and experimental metadata on the
samples. Note that here we focus on CpG methylation, since this is the
most common form of methylation in humans and many other organisms; take
care when applying this method to other types of methylation and make sure
that it will
be able to scale to the number of methylation sites, and that similar
assumptions can be made regarding spatial correlation. Also note that
the default settings for smoothing parameters and spacing/gap parameters
are set to values that we found useful, but may need to be altered for
datasets for other organisms.

To store your data in a `BSseq`

object, make sure you have the following
neccessary components:

genomic positions, including chromosome and location, for methylation loci.

a (matrix) of M (Methylation) values, describing the number of reads supporting methylation covering a single loci. Each row in this matrix is a methylation loci and each column is a sample.

a (matrix) of Cov (Coverage) values, describing the total number of reads covering a single loci. Each row in this matrix is a methylation loci and each column is a sample.

The following code chunk assumes that `chr`

and `pos`

are vectors of
chromosome names and positions, respectively, for each CpG in the dataset.
You can also provide a `GRanges`

object instead of `chr`

and `pos`

. It
also assumes that the matrices of methylation and coverage values (described
above) are named `M`

and `Cov`

, respectively. Note, `M`

and `Cov`

can also
be data stored on-disk (not in memory) using HDF5 files with the `HDF5Array`

package or `DelayedMatrix`

with the `DelayedArray`

package.

The `sampleNames`

and `trt`

objects are
vectors with sample labels and condition labels for each sample. A condition
label could be something like
treatment or control, a tissue type, or a continous measurement.
This is the covariate for which you wish to test for differences in
methylation. Once the `BSseq`

object is constructed and the sample covariate
information is added, DMRs are obtained by running the `dmrseq`

function.
A continuous covariate is assumed if the data type of the `testCovariate`

arugment in `dmrseq`

is
continuous, with the exception of if there are only two unique values
(then a two group comparison is carried out).

```
bs <- BSseq(chr = chr, pos = pos,
M = M, Cov = Cov,
sampleNames = sampleNames)
pData(bs)$Condition <- trt
regions <- dmrseq(bs=bs, testCovariate="Condition")
```

For more information on constructing and manipulating `BSseq`

objects,
see the bsseq vignettes.

- If you used
*Bismark*to align your bisulfite sequencing data, you can use the`read.bismark`

function to read bismark files into`BSseq`

objects. See below for more details.

Please post **dmrseq** questions to the
**Bioconductor support site**, which serves as a searchable knowledge
base of questions and answers:

https://support.bioconductor.org

Posting a question and tagging with “dmrseq” will automatically send an alert to the package authors to respond on the support site. See the first question in the list of Frequently Asked Questions (FAQ) for information about how to construct an informative post.

As input, the **dmrseq** package expects count data as obtained, e.g.,
from Bisulfite-sequencing. The value in the *i*-th row and the *j*-th column of
the `M`

matrix tells how many methylated reads can be assigned to CpG *i*
in sample *j*. Likewise, the value in the *i*-th row and the *j*-th column of
the `Cov`

matrix tells how many total reads can be assigned to CpG *i*
in sample *j*. Although we might be tempted to combine these matrices into
one matrix that contains the methylation *proportion* (`M`

/`Cov`

) at each CpG
site, it is critical to notice that this would be throwing away a lot of
information. For example, some sites have much higher coverage than others,
and naturally, we have more confidence in those with many reads mapping to them.
If we only kept the proportions, a CpG with 2 out of 2 reads methylated would
be treated the same as a CpG with 30 out of 30 reads methylated.

To use **dmrseq**, you need to have at least 2 samples in each condition.
Without this replicates, it is impossible to distinguish between biological
variability due to condition/covariate of interest, and inter-individual
variability within condition.

If your experiment contains additional samples, perhaps from other conditions
that are not of interest in the current test, these should be filtered out
prior to running **dmrseq**. Rather than creating a new filtered object,
the filtering step can be included in the call to the main function `dmrseq`

.
For more details, see the
Filtering CpGs and samples Section.

If you used Bismark for mapping and methylation level extraction, you can
use the `read.bismark`

function from the **bsseq** package to read the
data directly into
a `BSeq`

object.

The following example is from the help page of the function. After running
Bismark’s methylation extractor, you should have output files with names
that end in `.bismark.cov.gz`

. You can specify a vector of file names with
the `file`

argument, and a corresponding vector of `sampleNames`

. It is
recommended that you set `rmZeroCov`

to TRUE in order to remove CpGs with
no coverage in any of the samples, and set `strandCollapse`

to TRUE in order
to combine CpGs on opposite strands into one observation (since CpG methylation)
is symmetric.

```
## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationDbi'
```

```
infile <- system.file("extdata/test_data.fastq_bismark.bismark.cov.gz",
package = 'bsseq')
bismarkBSseq <- read.bismark(files = infile,
rmZeroCov = TRUE,
strandCollapse = FALSE,
verbose = TRUE)
bismarkBSseq
```

```
## An object of type 'BSseq' with
## 2013 methylation loci
## 1 samples
## has not been smoothed
## All assays are in-memory
```

See the bsseq help pages for more information on using this function.

If you haven’t used Bismark, but you have count data for number of methylated
reads and total coverage for each CpG, along with their corresponding chromosome
and position information, you can construct a `BSseq`

object from scratch,
like below. Notice that the `M`

and `Cov`

matrices have the same dimension, and
`chr`

and `pos`

have the same number of elements as rows in the count matrices
(which corresponds to the number of CpGs). Also note that the number of columns
in the count matrices matches the number of elements in `sampleNames`

and the
condition variable ’celltype`.

```
## <6 x 4> matrix of class DelayedMatrix and type "double":
## imr90.r1 imr90.r2 h1.r1 h1.r2
## [1,] 18 23 47 63
## [2,] 66 63 91 78
## [3,] 23 27 54 58
## [4,] 12 17 34 46
## [5,] 46 38 46 42
## [6,] 3 3 20 22
```

```
## <6 x 4> matrix of class DelayedMatrix and type "double":
## imr90.r1 imr90.r2 h1.r1 h1.r2
## [1,] 115 137 93 112
## [2,] 140 137 108 91
## [3,] 168 175 138 152
## [4,] 153 165 116 108
## [5,] 174 166 136 123
## [6,] 43 42 79 69
```

`## [1] "chr21" "chr21" "chr21" "chr21" "chr21" "chr21"`

`## [1] 9719962 9720054 9720224 9720323 9720909 9721175`

`## [1] 300684 4`

`## [1] 300684 4`

`## [1] 300684`

`## [1] 300684`

`## [1] "imr90.r1" "imr90.r2" "h1.r1" "h1.r2"`

`## [1] "imr90" "imr90" "h1" "h1"`

```
## An object of type 'BSseq' with
## 300684 methylation loci
## 4 samples
## has not been smoothed
## All assays are in-memory
```

The example data contains CpGs from chromosome 21 for four samples
from Lister et al. (2009). To load this data directly (already in the `BSseq`

format),
simply type `data(BS.chr21)`

.
Two of the samples are replicates of the cell type ‘imr90’
and the other two are replicates of the cell type ‘h1’. Now that we have the
data loaded into a `BSseq`

object, we can use **dmrseq**
to find regions of the genome where these two cell types have significantly
different methylation levels. But first, we need to add the sample metadata
that indicates which samples are from which cell type (the `celltype`

varialbe above). This information, which we call ‘metadata’,
will be used by the `dmrseq`

function to decide
which samples to compare to one another. The next section shows how to add
this information to the `BSseq`

object.

To add sample metadata, including the covariate of interest, you can add it
to the
`BSseq`

object by adding columns to the `pData`

slot. You must have at least
one column of `pData`

, which contains the covariate of interest. Additional
columns are optional.

```
pData(bs)$CellType <- celltype
pData(bs)$Replicate <- substr(sampleNames,
nchar(sampleNames), nchar(sampleNames))
pData(bs)
```

```
## DataFrame with 4 rows and 2 columns
## CellType Replicate
## <character> <character>
## imr90.r1 imr90 1
## imr90.r2 imr90 2
## h1.r1 h1 1
## h1.r2 h1 2
```

We will then tell the `dmrseq`

function which metadata variable to use
for testing for methylation differences by setting the `testCovariate`

parameter equal to its column name.

Note that unlike in **bsseq**, you do not need to carry out the smoothing step
with a separate function. In addition, you should not use **bsseq**’s `BSmooth()`

function to smooth the methylation levels, since **dmrseq** smooths in a very
different way. Briefly, **dmrseq** smooths methylation *differences*, so it
carries out the smoothing step once. This is automatically done with the main
`dmrseq`

function. **bsseq** on the other hand, smooths each sample
independently, so smoothing needs to be carried out once per sample.

For pairwise comparisons, **dmrseq** analyzes all CpGs that have at least one
read in at least one sample per group.
Thus, if your dataset contains CpGs with zero reads in every sample within a
group, you should filter them out prior to running `dmrseq`

. Likewise,
if your `bsseq`

object contains extraneous samples that are part of the
experiment but not the differential methylation testing of interest, these
should be filtered out as well.

Filtering `bsseq`

objects is straightforward:

- Subset rows to filter CpG loci
- Subset columns to filter samples

If we wish to remove all CpGs that have no coverage in at least one sample and only keep samples with a CellType of “imr90” or “h1”, we would do so with:

```
# which loci and sample indices to keep
loci.idx <- which(DelayedMatrixStats::rowSums2(getCoverage(bs, type="Cov")==0) == 0)
sample.idx <- which(pData(bs)$CellType %in% c("imr90", "h1"))
bs.filtered <- bs[loci.idx, sample.idx]
```

Note that this is a trivial example, since our toy example object `BS.chr21`

already contains only loci with coverage at least one read in all samples as well
as only samples from the “imr90” and “h1” conditions.

Also note that instead of creating a separate object, the filtering step
can be combined with the call to `dmrseq`

by replacing the `bs`

input with a
filtered version `bs[loci.idx, sample.idx]`

.

There are two ways to adjust for covariates in the dmrseq model. The first way
is to specify the `adjustCovariate`

parameter of the `dmrseq()`

function as
a column of the `pData()`

slot that contains the covariate you
would like to adjust for. This will include that covariate directly in the
model. This is ideal if the adjustment covariate is continuous or has more
than two groups.

The second way is to specify the `matchCovariate`

parameter of the `dmrseq`

function as a column of the `pData()`

slot that contains the covariate you
would like to match on. This will restrict the permutations considered to only
those where the `matchCovariate`

is balanced. For example, the `matchCovariate`

could represent the sex of each sample. In that case, a permutation that
includes all males in one group and all females in another would not be
considered (since there is a plausible biological difference that may induce
the null distribution to resemble non-null). This matching adjustment is ideal
for two-group comparisons.

The standard differential expression analysis steps are wrapped
into a single function, `dmrseq`

. The estimation steps performed
by this function are described briefly below, as well as in
more detail in the **dmrseq** paper. Here we run the results for a subset
of 20,000 CpGs in the interest of computation time.

```
testCovariate <- "CellType"
regions <- dmrseq(bs=bs[240001:260000,],
cutoff = 0.05,
testCovariate=testCovariate)
```

`## Assuming the test covariate CellType is a factor.`

`## Condition: imr90 vs h1`

`## Parallelizing using 4 workers/cores (backend: BiocParallel:MulticoreParam).`

`## Computing on 1 chromosome(s) at a time.`

`## Detecting candidate regions with coefficient larger than 0.05 in magnitude.`

```
## ...Chromosome chr21: Smoothed (0.07 min). 380 regions scored (0.15 min).
## * 380 candidates detected
## Performing balanced permutations of condition across samples to generate a null distribution of region test statistics
##
## Beginning permutation 1
## ...Chromosome chr21: Smoothed (0.08 min). 173 regions scored (0.08 min).
## * 1 out of 2 permutations completed (173 null candidates)
##
## Beginning permutation 2
## ...Chromosome chr21: Smoothed (0.07 min). 151 regions scored (0.05 min).
## * 2 out of 2 permutations completed (151 null candidates)
```

Progress messages are printed to the console if `verbose`

is TRUE.
The text, `condition h1 vs imr90`

, tells you that positive methylation
differences mean h1 has higher methylation than imr90 (see below for
more details).

The results object is a `GRanges`

object with the coordiates
of each candidate region, and contains the following metadata columns (which
can be extracted with the `$`

operator:

`L`

= the number of CpGs contained in the region,`area`

= the sum of the smoothed beta values`beta`

= the coefficient value for the condition difference (Note: if the test covariate is categorical with more than 2 groups, there will be more than one beta column),`stat`

= the test statistic for the condition difference,`pval`

= the permutation*p*-value for the significance of the test statistic, and`qval`

= the*q*-value for the test statistic (adjustment for multiple comparisons to control false discovery rate).`index = an`

IRanges` containing the indices of the region’s first CpG to last CpG.

```
## GRanges object with 380 ranges and 7 metadata columns:
## seqnames ranges strand | L area beta
## <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric>
## [1] chr21 43605625-43606688 * | 24 12.73094 -1.243455
## [2] chr21 44356429-44357250 * | 18 9.17973 -1.201346
## [3] chr21 44516290-44525100 * | 154 26.03485 -0.438070
## [4] chr21 44079752-44084490 * | 89 14.32149 -0.451037
## [5] chr21 44477228-44484542 * | 183 22.96151 -0.354231
## ... ... ... ... . ... ... ...
## [376] chr21 43628482-43628903 * | 8 0.536424 -0.0754552
## [377] chr21 43951797-43952020 * | 5 0.277526 0.0713785
## [378] chr21 43769879-43770694 * | 9 0.905185 0.0447750
## [379] chr21 44403584-44403679 * | 6 0.340733 -0.0329621
## [380] chr21 43584709-43584897 * | 5 0.302315 0.0292136
## stat pval qval index
## <numeric> <numeric> <numeric> <IRanges>
## [1] -18.5650 0.00307692 0.0259829 845-868
## [2] -12.5613 0.00307692 0.0259829 15147-15164
## [3] -12.2668 0.00307692 0.0259829 18812-18965
## [4] -11.7191 0.00307692 0.0259829 9742-9830
## [5] -11.0596 0.00307692 0.0259829 17954-18136
## ... ... ... ... ...
## [376] -0.431015 0.904615 0.914239 1356-1363
## [377] 0.348837 0.923077 0.930422 6801-6805
## [378] 0.287736 0.938462 0.943427 4490-4498
## [379] -0.234484 0.944615 0.947108 16341-16346
## [380] 0.132277 0.972308 0.972308 351-355
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
```

The above steps are carried out on a very small subset of data (20,000 CpGs).
This package loads data into memory one chromosome at a
time. For on human data, this means objects with a few million
entries per sample (since there are roughly 28.2 million total CpGs in the human
genome, and the largest chromosomes will have more than 2 million CpGs).
This means that whole-genome `BSseq`

objects for several samples can use up
several GB of RAM. In order to improve speed, the package allows for easy
parallel processing of chromosomes, but be aware that using more cores will
also require the use of more RAM.

To use more cores, use the `register`

function of
BiocParallel. For example,
the following chunk (not evaluated here), would register 4 cores, and
then the functions above would
split computation over these cores.

**dmrseq** is a two-stage approach that first detects candidate regions and then
explicitly evaluates statistical significance at the region level while
accounting for known sources of variability.
Candidate DMRs are defined by segmenting the genome into groups of CpGs
that show consistent evidence of differential methylation.
Because the methylation levels of neighboring CpGs are highly correlated,
we first smooth the signal to combat loss of power due to low coverage as done
in **bsseq**.

In the second stage, we compute a statistic for each candidate DMR that takes into account variability between biological replicates and spatial correlation among neighboring loci. Significance of each region is assessed via a permutation procedure which uses a pooled null distribution that can be generated from as few as two biological replicates, and false discovery rate is controlled using the Benjamini-Hochberg procedure.

For more details, refer to the **dmrseq** paper (Korthauer et al. 2018).

The default smoothing parameters (`bpSpan`

, `minInSpan`

, and `maxGapSmooth`

)
are designed to focus on local DMRs, generally in the range of hundreds to
thousands of bases. In some applications, such as cancer, it is of interest
to effectively ‘zoom out’ in order to detect larger (lower-resolution)
methylation blocks on the order of hundreds of thousands to millions of bases.
To do so, you can
set the `block`

argument to true, which will only include candidate regions with
at least `blockSize`

basepairs (default = 5000). This setting will also merge
candidate regions that (1) are in the same direction and (2) are less than 1kb
apart with no covered CpGs separating them. The region-level model used is also
slightly modified - instead of a loci-specific intercept for each CpG in the
region, the intercept term is modeled as a natural spline with one interior
knot per each 10kb of length (up to 10 interior knots).

In addition, detecting large-scale blocks requires that
the smoothing window be increased to minimize the impact of noisy local
methylation measurements. To do so, the values of the
smoothing parameters should be increased. For example, to use a smoothing window
that captures at least 500 CpGs or 50,000 basepairs that are spaced apart by no
more than 1e6 bases, use `minInSpan=500`

, `bpSpan=5e4`

, and `maxGapSmooth=1e6`

.
In addition, to avoid a block being broken up simply due to a gap with no
covered CpGs, you can increase the `maxGap`

parameter.

```
testCovariate <- "CellType"
blocks <- dmrseq(bs=bs[120001:125000,],
cutoff = 0.05,
testCovariate=testCovariate,
block = TRUE,
minInSpan = 500,
bpSpan = 5e4,
maxGapSmooth = 1e6,
maxGap = 5e3)
```

`## Searching for large scale blocks with at least 5000 basepairs.`

`## Assuming the test covariate CellType is a factor.`

`## Condition: imr90 vs h1`

`## Parallelizing using 4 workers/cores (backend: BiocParallel:MulticoreParam).`

`## Computing on 1 chromosome(s) at a time.`

`## Detecting candidate regions with coefficient larger than 0.05 in magnitude.`

```
## ...Chromosome chr21: Smoothed (0 min). 3 regions scored (0.35 min).
## * 3 candidates detected
## Performing balanced permutations of condition across samples to generate a null distribution of region test statistics
##
## Beginning permutation 1
## ...Chromosome chr21: Smoothed (0 min). No candidates found.
## No regions found.
## * 1 out of 2 permutations completed ( null candidates)
##
## Beginning permutation 2
## ...Chromosome chr21: Smoothed (0 min). No candidates found.
## No regions found.
## * 2 out of 2 permutations completed ( null candidates)
```

```
## Warning in dmrseq(bs = bs[120001:125000, ], cutoff = 0.05, testCovariate =
## testCovariate, : No candidate regions found in permutation, so inference can't
## be carried out. Try decreasing the cutoff, or running on a larger dataset if
## you are currently using a subset.
```

```
## GRanges object with 3 ranges and 7 metadata columns:
## seqnames ranges strand | L area beta
## <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric>
## [1] chr21 32345307-32520539 * | 1458 345.9647 -0.582317
## [2] chr21 32700803-32766102 * | 775 54.9953 -0.187589
## [3] chr21 32789489-32811767 * | 260 18.6190 -0.272892
## stat pval qval index
## <numeric> <logical> <logical> <IRanges>
## [1] -38.37178 <NA> <NA> 1-1458
## [2] -8.55934 <NA> <NA> 3670-4444
## [3] -7.08400 <NA> <NA> 4736-4995
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
```

The top hit is 175 thousand basepairs wide. In general, it also might be advised to decrease the cutoff when detecting blocks, since a smaller methylation difference might be biologically significant if it is maintained over a large genomic region. Note that block-finding can be more computationally intensive since we are fitting region-level models to large numbers of CpGs at a time. In the toy example above we are only searching over 5,000 CpGs (which span 467 thousand basepairs), so we do not find enough null candidate regions to carry out inference and obtain significance levels.

How many regions were significant at the FDR (*q*-value) cutoff of 0.05? We
can find this by counting how many values in the `qval`

column of the `regions`

object were less than 0.05.
You can also subset the regions by an FDR cutoff.

`## [1] 144`

You can determine the proportion of regions with hyper-methylation by counting how many had a positive direction of effect (positive statistic).

`## [1] 0.25`

To interpret the direction of effect, note that for a two-group comparison
**dmrseq** uses alphabetical order of the covariate of interest.
The condition with a higher alphabetical rank will become the reference category.
For example, if
the two conditions are “A” and “B”, the “A” group will be the reference category,
so a positive direction of effect means that
“B” is hyper-methylated relative to “A”. Conversely, a negative direction of
effect means that “B” is hypo-methylated relative to “A”.

It can be useful to visualize individual DMRs, so we provide a plotting
function that is based off of **bsseq**’s plotting functions. There is also
functionality to add annotations using the
annotatr package to
see the nearby CpG categories (island, shore, shelf, open sea) and nearby
coding sequences.

To retrieve annotations for genomes supported by **annotatr**, use the
helper function `getAnnot`

, and pass this annotation object to the `plotDMRs`

function as the `annoTrack`

parameter.

```
# get annotations for hg18
annoTrack <- getAnnot("hg18")
plotDMRs(bs, regions=regions[1,], testCovariate="CellType",
annoTrack=annoTrack)
```

Here we also plot the top methylation block from the block analysis: