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.

Installation

Install the package:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("RhyME")

Simulation example

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     R004

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())

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            TRUE

The 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 and period, 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] 146

The 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] 2000

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$meta

Apply 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).

Known simulation truth versus calls at BH < 0.05.
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

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()

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          28

Alternatively, 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).

Other options

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] 1200

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.

Session information

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 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] 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] sass_0.4.10         SparseArray_1.13.4  lattice_0.23-1     
#>  [4] digest_0.6.39       magrittr_2.0.5      evaluate_1.0.5     
#>  [7] grid_4.6.1          RColorBrewer_1.1-3  fastmap_1.2.0      
#> [10] jsonlite_2.0.0      Matrix_1.7-6        gnm_1.1-5          
#> [13] nnet_7.3-21         BiocManager_1.30.27 scales_1.4.0       
#> [16] codetools_0.2-20    jquerylib_0.1.4     abind_1.4-8        
#> [19] cli_3.6.6           rlang_1.3.0         XVector_0.53.0     
#> [22] withr_3.0.3         cachem_1.1.0        DelayedArray_0.39.8
#> [25] yaml_2.3.12         MetaCycle_1.2.1     S4Arrays_1.13.2    
#> [28] parallel_4.6.1      tools_4.6.1         dplyr_1.2.1        
#> [31] ggplot2_4.0.3       buildtools_1.0.0    vctrs_0.7.3        
#> [34] R6_2.6.1            lifecycle_1.0.5     relimp_1.0-5       
#> [37] MASS_7.3-66         pkgconfig_2.0.3     pillar_1.11.1      
#> [40] bslib_0.12.0        gtable_0.3.6        glue_1.8.1         
#> [43] Rcpp_1.1.2          tidyselect_1.2.1    tibble_3.3.1       
#> [46] xfun_0.62           sys_3.4.3           knitr_1.52         
#> [49] qvcalc_1.0.4        farver_2.1.2        htmltools_0.5.9    
#> [52] labeling_0.4.3      rmarkdown_2.32      maketools_1.3.2    
#> [55] compiler_4.6.1      S7_0.2.2

Using RhyME for circadian methylation analysis