Using RhyME for circadian methylation analysis
8 November 2026
Contents
1 Introduction
RhyME is an R package for detecting circadian rhythms in DNA methylation data. It fits beta-regression cosinor models to methylation beta values. The package supports independent samples and repeated measurements from the same subjects, with optional covariates such as age and sex. An empirical Bayes (EB) prior can improve the stability of precision estimates across CpGs.
RhyME identifies rhythmic CpG sites and estimates their amplitude, peak time, and mean methylation level. This vignette walks through installation, data preparation, model fitting, viewing results, and plotting fitted curves.
2 Installation
Install the package:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("RhyME")3 Simulation example
3.1 Simulation settings
This example uses simulated data with 15 subjects,
112 observed sampling records across nine planned time points,
and 2,000 CpGs comprising 200 rhythmic and 1,800 non-rhythmic sites.
Baseline methylation, precision, subject-level random-intercept variability,
and phase were jointly sampled from rows of CpG-specific parameter
estimates obtained from the real dataset. The logit-scale amplitude was
fixed at 0.09 for rhythmic sites and zero for non-rhythmic sites.
The parameter subset and sampling design are bundled in
inst/extdata/example_parameters.rds. Run inst/scripts/make_example_data.R
to regenerate the example; inst/scripts/README describes the inputs and usage.
library(RhyME)
library(SummarizedExperiment)
library(BiocParallel)
example <- readRDS(system.file("extdata", "example_data.rds",
package = "RhyME"))
se <- example$se
dim(se)
#> [1] 2000 112
head(colData(se))
#> DataFrame with 6 rows and 2 columns
#> time subject
#> <numeric> <factor>
#> R004_07 7 R004
#> R004_10 10 R004
#> R004_16 16 R004
#> R004_19 19 R004
#> R004_25 25 R004
#> R004_28 28 R0043.2 Fit the model
For these repeated measurements, specify subject_col to include a random
intercept per subject. The default prior = "auto" estimates an empirical Bayes
prior from the data. Set period to the rhythm length in hours and the default is 24.
Sampling times can be irregular or span several days.
fit <- fitCosinor(se, subject_col = "subject", prior = "auto", BPPARAM = SerialParam())3.3 View results
Results are stored in rowData(fit) with the RhyME_ prefix.
Here we view the first two rhythmic and first two non-rhythmic CpGs,
including their detection results and estimated rhythm parameters.
result <- rowData(fit)
truth <- result$true_status == "Rhythmic"
selected_cpgs <- rownames(fit)[c(head(which(truth), 2),
head(which(!truth), 2))]
columns <- c("true_status", "RhyME_pvalue", "RhyME_padj",
"RhyME_amplitude", "RhyME_phase", "RhyME_mesor",
"RhyME_converged")
result[selected_cpgs, columns]
#> DataFrame with 4 rows and 7 columns
#> true_status RhyME_pvalue RhyME_padj RhyME_amplitude RhyME_phase
#> <character> <numeric> <numeric> <numeric> <numeric>
#> cg09882730 Rhythmic 5.51092e-05 0.00134413 0.00470347 17.67165
#> cg03203646 Rhythmic 1.14089e-04 0.00232834 0.00671324 14.91218
#> cg13545212 Non-rhythmic 3.90288e-01 0.85108845 0.00303331 7.77422
#> cg02609848 Non-rhythmic 5.52659e-01 0.90897806 0.00055133 12.61443
#> RhyME_mesor RhyME_converged
#> <numeric> <logical>
#> cg09882730 0.0570941 TRUE
#> cg03203646 0.9102037 TRUE
#> cg13545212 0.0866832 TRUE
#> cg02609848 0.0241537 TRUEThe rhythm summaries describe the fitted mean curve on the beta-value scale. For models with a subject random intercept, they are population-averaged quantities obtained by marginalizing over the fitted random-intercept distribution.
RhyME_amplitude: half the fitted peak-to-trough difference on the beta-value scale.RhyME_phase: peak time in hours, between 0 andperiod, relative to the input time origin.RhyME_mesor: fitted mean methylation level averaged over one period.
Select CpGs with BH-adjusted p-values below 0.05:
detected <- result$RhyME_padj < 0.05
significant <- result[detected, , drop = FALSE]
nrow(significant)
#> [1] 146The simulation labels in rowData() are annotations and are not used to fit
the model. Below, we compare the calls from each method with these labels.
View the fitting summary and estimated EB prior:
metadata(fit)$RhyME$fitSummary
#> $total
#> [1] 2000
#>
#> $usable_tests
#> [1] 2000
#>
#> $usable_parameters
#> [1] 2000
metadata(fit)$RhyME$ebPrior
#> $mean
#> [1] 6.162728
#>
#> $sd
#> [1] 1.03046
#>
#> $n
#> [1] 20003.4 Run MetaCycle and compare calls with simulation truth
MetaCycle is an R package for detecting rhythmic signals in time-series data using ARSER, JTK_CYCLE, and Lomb–Scargle (LS). For the replicated observations at regularly spaced time points in this example, its default meta2d() strategy selects JTK_CYCLE and LS and combines their p-values using Fisher’s method. We specify a 24-hour period to match RhyME. Here, meta2d() uses the sampling times without modeling within-subject correlation, whereas RhyME includes a subject-specific random intercept. We compare the resulting rhythmicity calls with the known simulation truth.
Y <- data.frame(CycID = rownames(se), assay(se))
meta_fit <- MetaCycle::meta2d(
infile = "simulated_beta", inDF = Y, filestyle = "csv",
timepoints = se$time, minper = 24, maxper = 24,
outputFile = FALSE)
meta_result <- meta_fit$metaApply BH adjustment separately for each method across all 2,000 CpGs and use the same threshold of 0.05. The table compares known truth with detection: TP and FN are the detected and missed rhythmic CpGs; FP and TN are the non-rhythmic CpGs called rhythmic and not called rhythmic. Power is TP / 200, and the observed FDP is FP / (TP + FP).
| Method | TP | FN | FP | TN | Power (%) | Observed FDP (%) |
|---|---|---|---|---|---|---|
| RhyME | 140 | 60 | 6 | 1794 | 70.0 | 4.11 |
| JTK | 103 | 97 | 3 | 1797 | 51.5 | 2.83 |
| LS | 36 | 164 | 0 | 1800 | 18.0 | 0.00 |
| MetaCycle | 103 | 97 | 2 | 1798 | 51.5 | 1.90 |
3.5 Plot fitted curves
Plot the same four loci shown in the result table. Points are observed beta values, including between-subject variation, and red lines are fitted population-averaged mean curves. The first two panels are rhythmic loci; the last two are non-rhythmic. Each panel shows the locus ID from the table and uses its own beta-value range to show small changes.
plotCosinorFit(fit, cpgs = selected_cpgs, ncol = 2) +
ggplot2::theme_gray()
4 Prepare real data for RhyME
The study by Ammala et al. (2025)
provides longitudinal blood DNA methylation data through
OSF. Here we use 2,000 randomly sampled CpGs from
the processed baseline sleep/wake data, retaining all 112 available samples
from 15 participants. The subset is supplied as two CSV files: a beta-value
matrix and sample metadata. Sampling times are hours since the first planned
baseline sample, spanning 0–24 hours. Estimated phases use this same origin.
The original time values are kept in source_time.
Use a numeric matrix of beta values after quality control and normalization, with CpGs in rows and samples in columns. Sample metadata should contain collection time in hours and, for repeated measurements, a subject identifier. CpG and sample names should be unique, and sample names must match between the matrix and metadata.
Read the bundled files and assemble a SummarizedExperiment. The first
column of each file contains row names; read.csv() reads the compressed
beta-value file directly.
beta_file <- system.file("extdata", "real_beta.csv.gz", package = "RhyME")
sample_file <- system.file("extdata", "real_samples.csv", package = "RhyME")
beta <- as.matrix(read.csv(beta_file, row.names = 1,
check.names = FALSE))
sample_info <- read.csv(sample_file, row.names = 1)
sample_info <- sample_info[colnames(beta), , drop = FALSE]
real_se <- SummarizedExperiment(
assays = list(beta = beta),
colData = DataFrame(sample_info))
dim(real_se)
#> [1] 2000 112
head(colData(real_se))
#> DataFrame with 6 rows and 3 columns
#> time subject source_time
#> <integer> <character> <integer>
#> R004_07 3 R004 7
#> R004_10 6 R004 10
#> R004_16 12 R004 16
#> R004_19 15 R004 19
#> R004_25 21 R004 25
#> R004_28 24 R004 28Alternatively, extract beta values from a preprocessed minfi
GenomicRatioSet, retaining genomic coordinates and sample metadata:
se <- SummarizedExperiment(
assays = list(beta = minfi::getBeta(grSet)),
rowRanges = rowRanges(grSet),
colData = colData(grSet))For either input, use assay_name = "beta" in fitCosinor().
Missing beta values are allowed, but collection time, subject identifiers
when required, and any adjustment covariates must be present and complete
in colData(se).
5 Other options
5.1 Independent samples and the default fit
For independent samples, omit subject_col. The bundled example also contains
independent observations for 1,200 of the same CpGs, with 48 samples at six
collection times and no subject effects.
se_independent <- example$se_independent
fit_fixed <- fitCosinor(se_independent, BPPARAM = SerialParam())
metadata(fit_fixed)$RhyME$ebPrior
#> $mean
#> [1] 6.216663
#>
#> $sd
#> [1] 1.028647
#>
#> $n
#> [1] 12005.2 Adjust for covariates
Add covariates to colData() and pass their column names. Numeric and categorical variables are supported. Here we assign illustrative ages to subjects and analyze 25 rhythmic and 25 non-rhythmic CpGs. We use prior = “none” because this small subset provides limited information for estimating the empirical Bayes prior.
subject_age <- setNames(seq(28, 56, by = 2), levels(se$subject))
se$age <- unname(subject_age[as.character(se$subject)])
i <- c(head(which(truth), 25), head(which(!truth), 25))
adjusted <- fitCosinor(se[i, ],
time_col = "time", subject_col = "subject", covariates = "age",
prior = "none", BPPARAM = SerialParam(), verbose = FALSE)For multiple covariates, use covariates = c("age", "sex") when both
columns are present. Do not include the subject identifier as a covariate.
6 Session information
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB 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
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] BiocParallel_1.47.0 SummarizedExperiment_1.43.0
#> [3] Biobase_2.73.2 GenomicRanges_1.65.4
#> [5] Seqinfo_1.3.2 IRanges_2.47.5
#> [7] S4Vectors_0.51.10 BiocGenerics_0.59.12
#> [9] generics_0.1.4 MatrixGenerics_1.25.0
#> [11] matrixStats_1.5.0 RhyME_0.99.3
#> [13] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.62 bslib_0.12.0
#> [4] ggplot2_4.0.3 lattice_0.23-1 vctrs_0.7.3
#> [7] tools_4.6.1 parallel_4.6.1 tibble_3.3.1
#> [10] pkgconfig_2.0.3 Matrix_1.7-6 RColorBrewer_1.1-3
#> [13] S7_0.2.2 lifecycle_1.0.5 compiler_4.6.1
#> [16] farver_2.1.2 tinytex_0.61 codetools_0.2-20
#> [19] htmltools_0.5.9 sass_0.4.10 yaml_2.3.12
#> [22] pillar_1.11.1 jquerylib_0.1.4 MASS_7.3-66
#> [25] DelayedArray_0.39.8 cachem_1.1.0 magick_2.9.1
#> [28] abind_1.4-8 tidyselect_1.2.1 digest_0.6.39
#> [31] gnm_1.1-5 dplyr_1.2.1 bookdown_0.48
#> [34] labeling_0.4.3 qvcalc_1.0.4 fastmap_1.2.0
#> [37] grid_4.6.1 cli_3.6.6 SparseArray_1.13.4
#> [40] magrittr_2.0.5 S4Arrays_1.13.2 dichromat_2.0-1
#> [43] relimp_1.0-5 withr_3.0.3 scales_1.4.0
#> [46] rmarkdown_2.32 XVector_0.53.0 otel_0.3.0
#> [49] nnet_7.3-21 evaluate_1.0.5 knitr_1.52
#> [52] MetaCycle_1.2.1 rlang_1.3.0 Rcpp_1.1.2
#> [55] glue_1.8.1 BiocManager_1.30.27 jsonlite_2.0.0
#> [58] R6_2.6.1