Importing and filtering long-read modification data with bsseq

Introduction

This guide outlines how to import and filter Oxford Nanopore sequencing data using the bsseq package in R. Designed for DNA methylation analysis, bsseq provides tools for importing, filtering, analyzing, and visualizing modification data.

The first part focuses on preprocessing. It begins by processing raw POD5 sequencing files into basecalled and modification-called BAM files using dorado. These reads are then mapped to a reference genome while preserving modification information using SAMtools and Minimap2. Finally, modification calls are summarized into read-based bedMethyl files using Modkit.

The second part demonstrates how to import bedMethyl files as MethylCounts objects and filter loci based on coverage and the likelihood of representing homozygous and/or heterozygous CpG loci. The workflow is illustrated using both single-sample and multi-sample MethylCounts objects.

Terminology

The following terms are used throughout this document:

Basecalling: Determining the sequenced nucleotide sequence (A,C,G or T) from the raw sequencing signals.
Modification calling: Identifying the modification state from the raw sequencing signal. In this guide, limited to C in CpG-context modifications: C, 5hmC or 5mC.
CpG site: A cytosine followed by a guanine in a DNA strand (5’ to 3’).
CpG locus: The combined term for the CpG site the forward and reverse strand in double stranded DNA.
CpG state: The inferred status of a locus based on the sequencing data. A locus may be classified as homozygous CpG, heterozygous CpG, homozygous or heterozygous CpG (“allCpG”), or non-CpG.
.99 homozygous CpG filtering: Filtering for loci with scaled likelihood above 99% of being a homozygous CpG locus given the data.
.99 heterozygous CpG filtering: Filtering for loci with scaled likelihood of above 99% of being a heterozygous CpG locus given the data.
.99 ‘allCpG’ filtering: Filtering for loci with scaled likelihood of above 99% of being a homozygous CpG locus or heterozygous CpG locus given the data.
Reference-guided filtering: Restricting analysis to loci that overlap CpG loci in the reference genome (using the –cpg or –preset traditional in Modkit).

Citation

If you use the likelihood-based filtering methods described in this vignette, please cite our preprint (Hansen et al. 2025). General use of the package can be cited using the BSmooth paper (Hansen et al. 2012).

Dependencies

library(bsseq)
library(tidyverse)
#Additional software modules needed for preprocessing: 
#Dorado https://github.com/nanoporetech/dorado
#SAMtools https://github.com/samtools
#Minimap2 https://github.com/lh3/minimap2
#Modkit https://github.com/nanoporetech/modkit

Preprocessing

Before importing modification data from Oxford Nanopore sequencing into the bsseq package, the raw POD5 data must be processed to generate base and modification called reads, mapped to a reference genome, and summarized in the pileup format, bedMethyl.

Base and modification calling with Dorado

To obtain the base and modification calls from the raw signal data, we use the dorado basecaller. The following command will base and modification call the POD5 files in the input directory using a CpG-context model to call 5-methylcytosine and 5-hydroxymethylcytosine and output an unaligned BAM file. The unaligned BAM file will contain basecalls and modification calls for each read.

# Set input and output directories
input_directory= #/insert/input/directory/with/POD5/files/here
output_directory= #/insert/output/directory/for/unaligned/bam/files/here

# Run dorado basecaller with modification calling
dorado basecaller sup,5mCG_5hmCG \ 
$input_directory/ > $output_directory/unaligned.bam

Mapping with Minimap2 and SAMtools

To map the unaligned BAM file to a reference genome, we use minimap2 and samtools. The following command converts the BAM to FASTQ, maps the reads, and output a sorted and indexed BAM file with modification tags.

# Set output directory and reference genome
output_directory= #/insert/output/directory/for/unaligned/bam/files/here
reference_genome= #/insert/reference/genome/here
cd $output_directory

# Map reads to the reference genome
samtools fastq -TMM,ML unaligned.bam | \
  minimap2 -ax map-ont -y $reference_genome - | \
  samtools view -bS -| \
  samtools sort - > aligned.bam 
  
# Index the aligned BAM file  
samtools index aligned.bam

Pileup modifications with Modkit

To generate a modification pileup from the aligned BAM file, we use modkit. The pileup can be read-based approach to include all the CpG loci observed in the reads (recommended), or reference-guided and restricted to the CpG loci in the reference genome (not recommended).

Read-based pileup

For read-based pileup we do not set any flags related to the reference genome CpG loci (–CpG or —-motif CG 0). However, since we used a CpG-context model for modification calling, only the reference positions with at least one mapped CpG site are included in the pileup.

# Set directory
output_directory= #/insert/output/directory/for/unaligned/bam/files/here
cd $output_directory

# Pileup modification from all CpG loci observed in the reads
modkit pileup aligned.bam all_GpG.bedMethyl

Reference-guided pileup

For reference-guided pileup we use the --cpg flag to restrict the analysis to the CpG loci in the reference genome.

# Set directory and reference genome
output_directory= #/insert/output/directory/for/unaligned/bam/files/here
reference_genome= #/insert/reference/genome/here
cd $output_directory

# Pileup modification from the CpG loci present in the reference genome
modkit pileup --$reference_genome --cpg aligned.bam ref_GpG.bedMethyl

The commands above represent the minimum requirements needed to process the raw Oxford Nanopore sequencing data into a bedMethyl file. We strongly recommend reading the documentation for each program to adjust the parameters for your project and data.

Importing bedMethyl files

The function read.bedMethyl() reads one or more bedMethyl file(s) and returns a MethylCounts object which can be used for likelihood filtering, when setting output = "MethylCounts". Setting strandCollapse = TRUE merges data from the forward and reverse strand into a single representation.

Importing read-based pileup

The read-based bedMethyl files include all CpG loci observed in the reads should be read as MethylCounts objects as follows:

files <- list.files(system.file("extdata/longread/silversides/all_cpgs", 
                                package = "bsseq"),
                    full.names = TRUE)
mc_all<-read.bedMethyl(files, strandCollapse = T, output = "MethylCounts")
## Validating bedMethyl files and collecting metadata ...
mc_all
## An object of type 'MethylCounts' with
##   7675 loci
##   9 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

The MethylCounts object includes all loci that are a CpG sites in at least one read mapped to the position. This includes non-reference CpG loci and false positive CpG loci introduced by sequencing and mapping errors.

Importing reference CpG loci only

The reference-guided bedMethyl files include only the CpG loci present in the reference genome, and can be read as BSseq objects or as MethylCounts objects as follows:

files <- list.files(system.file("extdata/longread/silversides/ref_cpgs", 
                                package = "bsseq"),
                    full.names = TRUE)
mc_cpg<-read.bedMethyl(files=files, strandCollapse = T, output = "MethylCounts")
mc_cpg
## An object of type 'MethylCounts' with
##   2176 loci
##   9 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

This object includes only loci that are a CpG site in at least one read, and are a CpG loci in the reference genome.

Filter a single sample

Modification analysis is often restricted to a subset of the loci in a sample. Filtering can be based on coverage thresholds or based on the likelihood of a locus being a homozygous or heterozygous CpG.

Coverage filtering

A MethylCounts object can be filtered using getMethylCounts(), where the coverage represent the number of times a CpG site is mapped at a specific locus.

All CpG loci

#get the first sample
mc_all_1 <- mc_all[,1]
#get the indices of loci with coverage >= 5
loci.idx <- which(getMethylCounts(mc_all_1, type="Cov")>= 5)
#filter the object to retain these loci
mc_all_1_filtered <- mc_all_1[loci.idx,]
mc_all_1_filtered
## An object of type 'MethylCounts' with
##   2245 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

Reference CpG loci

mc_cpg_1 <- mc_cpg[,1]
loci.idx <- which(getMethylCounts(mc_cpg_1, type="Cov")>= 5)
mc_cpg_1_filtered<-mc_cpg_1[loci.idx,]
mc_cpg_1_filtered
## An object of type 'MethylCounts' with
##   1929 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

In addition to the sample and the sequencing depth of a sample, the number of coverage filtered loci depend on whether non-reference CpG loci are included, and which coverage threshold is applied.

Likelihood filtering

A MethylCounts object imported using read.bedMethyl() includes both the CpG coverage and the non-CpG coverage all loci, which we use to estimate the error rate and call the CpG status i.e. determine if a locus is a homozygous CpG, a heterozygous CpG or not a CpG at all.

.99 “allCpG” filtering

To get the total (homozygous and heterozygous) CpG loci in a sample, one can use the getCpGs() function with type set to "allCpG".

This can be for all the CpG loci observed in the reads:

#get the first sample
mc_all_1 <- mc_all[,1]
#get the indices of loci with scaled likelihood above 0.99 of being a "allCpG"
loci.idx <- getCpGs(mc_all_1, type = "allCpG", threshold = 0.99)
#filter the object to retain these loci
mc_all_1_filtered <- mc_all_1[loci.idx,]
mc_all_1_filtered
## An object of type 'MethylCounts' with
##   2263 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

Or for all the reference CpG loci observed in the reads:

mc_cpg_1 <- mc_cpg[,1]
loci.idx <- getCpGs(mc_cpg_1, type = "allCpG", threshold = 0.99)
mc_cpg_1_filtered <- mc_cpg_1[loci.idx,]
mc_cpg_1_filtered
## An object of type 'MethylCounts' with
##   1947 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

The number of “AllCpG” loci depends on the likelihood threshold and whether non-reference CpG loci are included.

.99 homozygous CpG

To get the homozygous CpG loci in a sample we use getCpGs() with type = "homozygous".

This can be for all the CpG loci observed in the reads:

#get the first sample
mc_all_1 <- mc_all[,1]
#get the indices of loci with scaled likelihood above 0.99 of being homozygous
loci.idx <- getCpGs(mc_all_1, type = "homozygous", threshold = 0.99)
#filter the object to retain these loci
mc_all_1_filtered <- mc_all_1[loci.idx,]
mc_all_1_filtered
## An object of type 'MethylCounts' with
##   1978 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

Or for all the reference CpG loci observed in the reads:

mc_cpg_1 <- mc_cpg[,1]
loci.idx <- getCpGs(mc_cpg_1, type = "homozygous", threshold = 0.99)
mc_cpg_1_filtered <- mc_cpg_1[loci.idx,]
mc_cpg_1_filtered
## An object of type 'MethylCounts' with
##   1764 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

The number of homozygous loci depends on the likelihood threshold and whether non-reference CpG loci are included.

.99 heterozygous CpG

To get the heterozygous CpG loci in a sample, we use the getCpGs() function with type = "heterozygous".

This can be for all the CpG loci observed in the reads:

#get the first sample
mc_all_1 <- mc_all[,1]
#get the indices of loci with scaled likelihood above 0.99 of being heterozygous
loci.idx <- getCpGs(mc_all_1, type = "heterozygous", threshold = 0.99)
#filter the object to retain these loci
mc_all_1_filtered <- mc_all_1[loci.idx,]
mc_all_1_filtered
## An object of type 'MethylCounts' with
##   152 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

Or for all the reference CpG loci observed in the reads:

mc_cpg_1 <- mc_cpg[,1]
loci.idx <- getCpGs(mc_cpg_1, type = "heterozygous", threshold = 0.99)
mc_cpg_1_filtered <- mc_cpg_1[loci.idx,]
mc_cpg_1_filtered
## An object of type 'MethylCounts' with
##   74 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

The number of heterozygous CpG loci depends on the likelihood threshold applied and whether non-reference CpG loci are retained. Reference-guided pileup generally removes ~half of the heterozygous loci in a sample.

Filter multiple samples

A project often includes multiple samples, and filtering can be applied to all samples in a MethylCounts object. Samples often share CpG loci and it is therefore advantageous to filter the loci in the multi-sample MethylCounts object to avoid false positive CpG loci introduced by sequencing and mapping errors.

Coverage filtering

A multi-sample MethylCounts object can be coverage filtered using getMethylCounts().

Coverage threshold

The coverage threshold can restrict the analysis to the loci with a coverage above a threshold in all samples (here, 9 samples):

Samples threshold

The coverage filtering can be less conservative by including all loci passing the threshold of 5X in some of the 9 samples:

Example

Filtering a MethylCounts object to retain only loci with a coverage of at least 5X in at least six of the nine samples can be done using this command:

loci.idx <- which(
  DelayedMatrixStats::rowSums2(getMethylCounts(mc_all, type="Cov")>= 5) >= 6)
mc_coverage_filtered <- mc_all[loci.idx,]
mc_coverage_filtered
## An object of type 'MethylCounts' with
##   2339 loci
##   9 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

Likelihood filtering

For likelihood filtering of multiple samples, we recommend utilizing the functions getCpGMatrix() and getMaxLikelihoodMatrix() to obtain the CpG matrix and the maximum likelihood matrix for the samples.

The function getCpGMatrix() returns a matrix with the most likely CpG call for the loci and sample with the same dimensions as the MethylCounts object. In the the default setting homozygous CpG loci are represented by 0, heterozygous CpG loci by 1 and non-CpG loci by 2.

G_all <- getCpGMatrix(mc_all)
head(G_all)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    2    2    2    2    2    2    2    2    2
## [2,]    0    0    0    0    0    0    0    0    0
## [3,]    2    2    2    2    2    2    2    2    2
## [4,]    2    2    2    2    2    2    2    2    2
## [5,]    2    2    2    2    2    2    2    2    2
## [6,]    2    2    2    2    2    2    2    2    2

The function getMaxLikelihoodMatrix returns a matrix with the corresponding scaled likelihood of most likely CpG call for the loci and sample with the same dimensions.

Q_all <- getMaxLikelihoodMatrix(mc_all)
head(round(Q_all,3))
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
## [1,] 0.333 1.000 0.333 0.333 0.333 0.333 0.333 0.333 1.000
## [2,] 0.987 1.000 0.872 1.000 1.000 1.000 1.000 1.000 1.000
## [3,] 0.333 0.333 0.333 1.000 0.333 0.333 0.333 0.333 0.333
## [4,] 0.333 0.999 0.333 0.333 0.333 0.333 0.333 0.333 0.333
## [5,] 0.333 0.999 0.333 0.333 0.333 0.333 0.333 0.333 0.333
## [6,] 0.333 0.333 0.653 0.333 0.333 0.333 0.333 0.333 0.333

Both functions can be run using allCpG = TRUE, where 0 in the CpG matrix represent homozygous or heterozygous CpG and the MaxLikelihoodMatrix represent the joint probability of homozygous or heterozygous CpG. This is useful for filtering the data based on likelihood thresholds.

G_all <- getCpGMatrix(mc_all, allCpG = TRUE)
head(G_all)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    2    2    2    2    2    2    2    2    2
## [2,]    0    0    0    0    0    0    0    0    0
## [3,]    2    2    2    2    2    2    2    2    2
## [4,]    2    2    2    2    2    2    2    2    2
## [5,]    2    2    2    2    2    2    2    2    2
## [6,]    2    2    2    2    2    2    2    2    2
Q_all <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE)
head(round(Q_all,3))
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
## [1,] 0.333 1.000 0.333 0.333 0.333 0.333 0.333 0.333 1.000
## [2,] 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## [3,] 0.333 0.333 0.333 1.000 0.333 0.333 0.333 0.333 0.333
## [4,] 0.333 0.999 0.333 0.333 0.333 0.333 0.333 0.333 0.333
## [5,] 0.333 0.999 0.333 0.333 0.333 0.333 0.333 0.333 0.333
## [6,] 0.333 0.333 0.653 0.333 0.333 0.333 0.333 0.333 0.333

In both settings, locus with a coverage of 0 are represented as non-CpGs with likelihood of 1/3.

Likelihood threshold

The multi-sample MethylCounts object can be filtered for loci scaled likelihood of being a homozygous or heterozygous CpG loci in all nine samples above a certain likelihood threshold.

Samples threshold

The likelihood threshold can be combined with the samples threshold for more relaxed filtering for loci scaled likelihood of being a homozygous or heterozygous CpG loci in X out of all nine samples.

Example

The likelihood and samples thresholds can be combined to filter for e.g. loci with a likelihood above 0.99 of being a homozygous or heterozygous CpG loci in at least six of the nine samples:

G_all <- getCpGMatrix(mc_all, allCpG = TRUE)
Q_all <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE)

loci.idx <- which(
  DelayedMatrixStats::rowSums2(Q_all >= .99 & G_all==0) >= 6)
mc_filtered <- mc_all[loci.idx,]
mc_filtered
## An object of type 'MethylCounts' with
##   2331 loci
##   9 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

HDF5 storage

The raw or filtered MethylCounts object can be saved as an HDF5 file using the saveHDF5SummarizedExperiment() function from the HDF5Array package. This allows for efficient storage and retrieval of large datasets.

#hdf5_path <- 
#    "~/Documents/bsseq/inst/extdata/longread/silversides/mc_filtered.hdf5"
#HDF5Array::saveHDF5SummarizedExperiment(mc_filtered, hdf5_path)

# Load pre-saved HDF5 file
files <- system.file("extdata/longread/silversides/mc_filtered.hdf5", 
                                package = "bsseq")

mc_filtered <- HDF5Array::loadHDF5SummarizedExperiment(files)
mc_filtered
## An object of type 'MethylCounts' with
##   2331 loci
##   9 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## Some assays are HDF5Array-backed

Converting to BSseq

For downstream analysis the filtered MethylCounts object can be converted to a BSseq object using the BSseq() constructor function.

# Filter and process a sample from the MethylCounts object
mc_sample <- mc_filtered[, 1]
mc_sample_filtered <- mc_sample[
  getCpGs(mc_sample, type = "homozygous", threshold = 0.99)]

# Display the filtered MethylCounts object
mc_sample_filtered
## An object of type 'MethylCounts' with
##   1952 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## Some assays are HDF5Array-backed

Depending on the analysis, you can specify the modification type (mods) as “5mC+5hmC” (default).

# Convert to BSseq with both 5mC and 5hmC data
bs <- BSseq(mc = mc_sample_filtered, mods = "5mC+5hmC")
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
bs
## An object of type 'BSseq' with
##   1952 loci
##   1 samples
##   5mC+5hmC values are stored in 'M'
## has not been smoothed
## Some assays are HDF5Array-backed
# Calculate the mean methylation using raw data
mean(getMeth(bs, type = "raw"))
## [1] 0.5310431

Or only one modification type, e.g., “5mC”:

# Convert to BSseq with only 5mC data
bs_M <- BSseq(mc = mc_sample_filtered, mods = "5mC")
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
bs_M
## An object of type 'BSseq' with
##   1952 loci
##   1 samples
##   5mC values are stored in 'M'
## has not been smoothed
## Some assays are HDF5Array-backed
# Calculate the mean methylation for 5mC
mean(getMeth(bs_M, type = "raw"))
## [1] 0.4927917

or “5hmC”:

# Convert to BSseq with only 5hmC data
bs_H <- BSseq(mc = mc_sample_filtered, mods = "5hmC")
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
## using unknown matrix fallback for 'HDF5ArraySeed'
bs_H
## An object of type 'BSseq' with
##   1952 loci
##   1 samples
##   5hmC values are stored in 'M'
## has not been smoothed
## Some assays are HDF5Array-backed
# Calculate the mean methylation for 5hmC
mean(getMeth(bs_H, type = "raw"))
## [1] 0.03825143

Summary

In summary bedMethyl files from modkit can be imported in to bsseq and filtered for loci with a specific coverage or a specific likelihood of being a homozygous, heterozygous or homozygous/heterozygous CpG loci in all or a subset of samples. Multiple samples can be imported and filtered for loci with a scaled likelihood above 0.99 of being a homozygous orheterozygous CpG loci in at least six samples, and a coverage of at least 5 in at least six samples:

files <- list.files(system.file("extdata/longread/silversides/all_cpgs", 
                                package = "bsseq"),
                    full.names = TRUE)
mc_all<-read.bedMethyl(files, strandCollapse = T, output = "MethylCounts")

G <- getCpGMatrix(mc_all, allCpG = TRUE)
Q <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE)

loci.idx <- which(
  DelayedMatrixStats::rowSums2(Q >= .99 & G==0) >= 6 &
  DelayedMatrixStats::rowSums2(getMethylCounts(mc_all, type="Cov")>= 5) >= 6)
mc_filtered <- mc_all[loci.idx,]
mc_filtered
## An object of type 'MethylCounts' with
##   2317 loci
##   9 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

Specific samples in the MethylCounts object can be filtered for loci with a high probability of being e.g. homozygous CpG loci:

mc_sample<- mc_filtered[,1]
mc_sample_filtered<-mc_sample[
  getCpGs(mc_sample, type = "homozygous", threshold = 0.99)]
mc_sample_filtered
## An object of type 'MethylCounts' with
##   1950 loci
##   1 samples
##   5mC and 5hmC data for all samples
##   CG-context modification model detected for all samples
## All assays are in-memory

And can be converted to a BSseq object for downstream modification analyses:

bs_sample_filtered<-BSseq(mc = mc_sample_filtered)
bs_sample_filtered
## An object of type 'BSseq' with
##   1950 loci
##   1 samples
##   5mC+5hmC values are stored in 'M'
## has not been smoothed
## All assays are in-memory

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [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       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] lubridate_1.9.5             forcats_1.0.1              
##  [3] stringr_1.6.0               dplyr_1.2.1                
##  [5] purrr_1.2.2                 readr_2.2.0                
##  [7] tidyr_1.3.2                 tibble_3.3.1               
##  [9] ggplot2_4.0.3               tidyverse_2.0.0            
## [11] bsseqData_0.51.0            bsseq_1.49.2               
## [13] SummarizedExperiment_1.43.0 Biobase_2.73.2             
## [15] MatrixGenerics_1.25.0       matrixStats_1.5.0          
## [17] GenomicRanges_1.65.1        Seqinfo_1.3.0              
## [19] IRanges_2.47.2              S4Vectors_0.51.6           
## [21] BiocGenerics_0.59.12        generics_0.1.4             
## [23] BiocStyle_2.41.0           
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.1-0              permute_0.9-10           
##  [3] rlang_1.3.0               magrittr_2.0.5           
##  [5] otel_0.2.0                compiler_4.6.1           
##  [7] DelayedMatrixStats_1.35.0 vctrs_0.7.3              
##  [9] pkgconfig_2.0.3           crayon_1.5.3             
## [11] fastmap_1.2.0             XVector_0.53.0           
## [13] labeling_0.4.3            Rsamtools_2.29.0         
## [15] rmarkdown_2.31            tzdb_0.5.0               
## [17] xfun_0.60                 cachem_1.1.0             
## [19] beachmat_2.29.0           cigarillo_1.3.1          
## [21] jsonlite_2.0.0            rhdf5filters_1.25.4      
## [23] DelayedArray_0.39.6       Rhdf5lib_2.1.0           
## [25] BiocParallel_1.47.0       parallel_4.6.1           
## [27] R6_2.6.1                  bslib_0.12.0             
## [29] stringi_1.8.9             RColorBrewer_1.1-3       
## [31] limma_3.69.4              rtracklayer_1.73.0       
## [33] jquerylib_0.1.4           Rcpp_1.1.2               
## [35] knitr_1.51                R.utils_2.13.0           
## [37] BiocBaseUtils_1.15.1      timechange_0.4.0         
## [39] Matrix_1.7-6              tidyselect_1.2.1         
## [41] abind_1.4-8               yaml_2.3.12              
## [43] codetools_0.2-20          curl_7.1.0               
## [45] lattice_0.23-1            withr_3.0.3              
## [47] S7_0.2.2                  evaluate_1.0.5           
## [49] Biostrings_2.81.6         pillar_1.11.1            
## [51] BiocManager_1.30.27       RCurl_1.98-1.19          
## [53] hms_1.1.4                 sparseMatrixStats_1.25.0 
## [55] scales_1.4.0              gtools_3.9.5             
## [57] glue_1.8.1                maketools_1.3.2          
## [59] tools_4.6.1               BiocIO_1.23.3            
## [61] sys_3.4.3                 data.table_1.18.4        
## [63] BSgenome_1.81.1           locfit_1.5-9.12          
## [65] GenomicAlignments_1.49.1  buildtools_1.0.0         
## [67] XML_3.99-0.23             rhdf5_2.57.11            
## [69] grid_4.6.1                HDF5Array_1.41.2         
## [71] restfulr_0.0.17           cli_3.6.6                
## [73] S4Arrays_1.13.0           gtable_0.3.6             
## [75] R.methodsS3_1.8.2         sass_0.4.10              
## [77] digest_0.6.39             SparseArray_1.13.2       
## [79] rjson_0.2.23              farver_2.1.2             
## [81] htmltools_0.5.9           R.oo_1.27.1              
## [83] lifecycle_1.0.5           h5mread_1.5.1            
## [85] httr_1.4.8                statmod_1.5.2

References

Hansen, Kasper D, Benjamin Langmead, and Rafael A Irizarry. 2012. “BSmooth: from whole genome bisulfite sequencing reads to differentially methylated regions.” Genome Biology 13 (10): R83. https://doi.org/10.1186/gb-2012-13-10-r83.
Hansen, Soeren B., Kasper D. Hansen, and Morten T. Limborg. 2025. “Sample-specific CpG loci are important for accurate long-read methylation analysis.” bioRxiv, ahead of print. https://doi.org/10.1101/2025.04.30.651558.