Detecting Circadian Rhythmicity in DNA Methylation with RhyME

Introduction

RhyME is an R package for detecting circadian rhythmicity in DNA methylation data. DNA methylation beta values are bounded between 0 and 1 and often exhibit strong mean-variance dependence, making standard Gaussian cosinor models less well suited for modeling methylation dynamics. In addition, circadian methylation studies commonly involve repeated measurements from the same subjects.

RhyME addresses these issues by fitting cosinor-based beta regression models for methylation beta values. The model represents rhythmic methylation patterns using sine and cosine terms, while optional subject-specific random intercepts account for repeated measurements. To improve stability in genome-wide analyses, RhyME also incorporates empirical Bayesian shrinkage of the beta-regression precision parameter.

The modeling framework uses beta regression (Ferrari and Cribari-Neto, 2004), cosinor rhythmometry, and Template Model Builder’s automatic differentiation and Laplace approximation (Kristensen et al., 2016).

This vignette demonstrates how to use RhyME for:

  1. Fitting fixed- and mixed-effect beta-regression cosinor models to methylation data.
  2. Detecting rhythmic CpG sites and ranking results.
  3. Visualizing fitted circadian methylation patterns.

Installation

After RhyME becomes available through Bioconductor, install it with:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("RhyME")

Data Preparation

RhyME expects a SummarizedExperiment with: 1. An assay containing methylation beta values (CpGs x samples); 2.colData with at minimum a numeric time column (in hours) 3.Optionally, a subject identifier column for repeated measures.

The package includes example data (3000 CpGs, 8 subjects, 6 timepoints). To keep this vignette quick to build, the evaluated example uses 50 CpGs:

library(RhyME)
library(SummarizedExperiment)

# Load bundled example data
example <- readRDS(system.file("extdata", "example_data.rds",
    package = "RhyME"))
se <- example$se[seq_len(50), ]
se
#> class: SummarizedExperiment 
#> dim: 50 48 
#> metadata(0):
#> assays(1): beta
#> rownames(50): cg00001 cg00002 ... cg00049 cg00050
#> rowData names(0):
#> colnames(48): S1_T0 S1_T4 ... S8_T16 S8_T20
#> colData names(3): time subject Sample_Name

Statistical Model

For CpG \(g\), subject \(i\), and observation \(j\), RhyME assumes

\[ Y_{gij} \sim \operatorname{Beta}(\mu_{gij}, \phi_g), \]

with

\[ \operatorname{logit}(\mu_{gij}) = \beta_{0g} + \beta_{cg}\cos(2\pi t_{ij}/P) + \beta_{sg}\sin(2\pi t_{ij}/P) + \mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\gamma}_g + b_{gi}. \]

For the mixed-effects model, the subject-specific random intercept follows

\[ b_{gi} \sim \mathcal{N}(0,\tau_g^2). \]

For longitudinal data, \(b_{gi}\) is a subject-specific random intercept. It is omitted when fitting the fixed-effects model. By default, the period \(P\) is 24 hours. Optional covariates \(\mathbf{x}_{ij}\), such as sex or age, are included in both the rhythmic and null models as adjustment terms. By default, rhythmicity is assessed using a joint Wald test of \(H_0:\beta_{cg}=\beta_{sg}=0\). The test concerns only the sine and cosine coefficients rather than the adjustment covariates.

Empirical-Bayes Precision Shrinkage

RhyME optionally stabilizes CpG-specific precision estimates by placing an empirical-Bayes prior on the log-precision parameter,

\[ \log(\phi_g) \sim \mathcal{N}(\mu_\phi, \sigma_\phi^2). \]

With prior = "auto", the prior center and scale are estimated robustly from a pilot subset of CpGs. The estimated prior is then used to shrink noisy site-specific precision estimates toward the shared prior distribution. Setting prior = "none" disables precision shrinkage.

Fixed-Effects Analysis

For cross-sectional studies or datasets in which observations are independent, RhyME fits a fixed-effects beta-regression cosinor model when subject_col is not specified.

se_fixed <- fitCosinor(se,
    time_col = "time",
    prior = "none"
)
se_fixed

Optional sample- or subject-level confounders can be added with, for example, covariates = c("sex", "age"). They are included as adjustment terms, while the reported pvalue tests rhythmicity.

Longitudinal Analysis (Mixed-Effects)

For longitudinal studies, repeated observations from the same subject are generally correlated. Specifying subject_col fits a subject-specific random intercept to account for this within-subject correlation.

set.seed(1)
se_demo <- fitCosinor(se,
    time_col = "time",
    subject_col = "subject",
    prior = "none",
    verbose = FALSE
)

For a genome-wide analysis, enable automatic prior estimation with prior = "auto".

fit_auto <- fitCosinor(
    example$se,
    time_col = "time",
    subject_col = "subject",
    prior = "auto"
)

Parameter Estimates and Result Interpretation

The output of fitCosinor() is the input SummarizedExperiment object with model results appended to rowData() and analysis-level information stored in metadata().

Specifically:

  • rowData(se_demo) contains per-CpG model estimates, circadian summary hypothesis-test p-values, adjusted p-values, and convergence indicators.
  • metadata(se_demo)$RhyME contains analysis-level information, including the model design, prior and test settings, and fitting summary.

All per-CpG model results are available from rowData(se_demo):

res <- rowData(se_demo)
head(res)
#> DataFrame with 6 rows and 12 columns
#>                b0     b_sin     b_cos amplitude_link amplitude_response
#>         <numeric> <numeric> <numeric>      <numeric>          <numeric>
#> cg00001  0.872055  0.213472  0.387246       0.442188          0.0910976
#> cg00002  0.497807  0.397548  0.184976       0.438475          0.1009806
#> cg00003  0.692676  0.153464  0.332504       0.366210          0.0807892
#> cg00004 -1.106705  0.262553  0.296635       0.396140          0.0735327
#> cg00005  0.568513  0.284908  0.118302       0.308493          0.0707801
#> cg00006 -0.482277  0.407751  0.415961       0.582482          0.1335751
#>             phase mesor_link mesor_response         tau      pvalue        padj
#>         <numeric>  <numeric>      <numeric>   <numeric>   <numeric>   <numeric>
#> cg00001   1.92440   0.872055       0.700758 9.00344e-02 1.83984e-06 3.40710e-06
#> cg00002   4.33652   0.497807       0.618238 1.97527e-01 6.56190e-07 1.49134e-06
#> cg00003   1.65168   0.692676       0.664114 1.31115e-05 4.35966e-06 7.26610e-06
#> cg00004   2.76748  -1.106705       0.253939 2.02582e-01 9.14374e-08 2.28593e-07
#> cg00005   4.49668   0.568513       0.636915 6.95681e-06 9.32588e-02 9.32588e-02
#> cg00006   2.96193  -0.482277       0.386943 1.68350e-01 4.06039e-09 1.19423e-08
#>         converged
#>         <logical>
#> cg00001      TRUE
#> cg00002      TRUE
#> cg00003      TRUE
#> cg00004      TRUE
#> cg00005      TRUE
#> cg00006      TRUE

The main circadian summaries are:

  • amplitude_link: fitted oscillation amplitude on the logit-link scale.
  • amplitude_response: estimated oscillation amplitude on the methylation beta-value scale.
  • phase: estimated timing of peak methylation, reported in hours.
  • mesor_link: fitted MESOR on the logit-link scale.
  • mesor_response: rhythm-adjusted mean methylation level on the beta-value scale.

In mixed-effects mode, response-scale summaries are reported as population-average quantities. The fitted methylation curve is averaged over the subject-specific random-intercept distribution, and plotCosinorFit() displays the same population-average curve by default.

Top-ranked CpG sites

The function topCpGs() returns the top-ranked CpG sites according to the rhythmicity test p-value or adjusted p-value.

topCpGs(se_demo, n = 10)
#> DataFrame with 10 rows and 6 columns
#>            cpg amplitude_response     phase mesor_response      pvalue
#>    <character>          <numeric> <numeric>      <numeric>   <numeric>
#> 1      cg00012           0.126640   2.33704       0.459666 5.40677e-18
#> 2      cg00031           0.153276   2.99069       0.539040 4.31393e-13
#> 3      cg00038           0.159600   3.20391       0.603375 1.03131e-12
#> 4      cg00043           0.104492   4.12467       0.376758 2.48179e-12
#> 5      cg00032           0.114898   1.21763       0.531620 2.73400e-12
#> 6      cg00046           0.162601   3.22724       0.475082 3.40613e-12
#> 7      cg00047           0.134797   2.83950       0.606229 3.78471e-11
#> 8      cg00028           0.122757   3.52351       0.471980 4.74209e-11
#> 9      cg00017           0.128277   3.17616       0.652890 5.67394e-11
#> 10     cg00013           0.120354   4.07705       0.513148 8.26264e-11
#>           padj
#>      <numeric>
#> 1  2.70338e-16
#> 2  1.07848e-11
#> 3  1.71885e-11
#> 4  2.73400e-11
#> 5  2.73400e-11
#> 6  2.83844e-11
#> 7  2.70336e-10
#> 8  2.96380e-10
#> 9  3.15219e-10
#> 10 4.13132e-10

Visualization

The fitted circadian methylation pattern can be visualized using plotCosinorFit(). By default, the function displays fitted 24-hour cosinor curves for the top-ranked CpG sites.

plotCosinorFit(se_demo)

Session Info

sessionInfo()
#> 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] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [3] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [5] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [7] BiocGenerics_0.59.12        generics_0.1.4             
#>  [9] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [11] RhyME_0.99.0                BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10         SparseArray_1.13.2  lattice_0.23-1     
#>  [4] digest_0.6.39       evaluate_1.0.5      grid_4.6.1         
#>  [7] RColorBrewer_1.1-3  fastmap_1.2.0       jsonlite_2.0.0     
#> [10] Matrix_1.7-6        BiocManager_1.30.27 scales_1.4.0       
#> [13] codetools_0.2-20    jquerylib_0.1.4     abind_1.4-8        
#> [16] cli_3.6.6           rlang_1.3.0         XVector_0.53.0     
#> [19] withr_3.0.3         cachem_1.1.0        DelayedArray_0.39.6
#> [22] yaml_2.3.12         S4Arrays_1.13.0     tools_4.6.1        
#> [25] parallel_4.6.1      BiocParallel_1.47.0 ggplot2_4.0.3      
#> [28] buildtools_1.0.0    vctrs_0.7.3         R6_2.6.1           
#> [31] lifecycle_1.0.5     bslib_0.12.0        gtable_0.3.6       
#> [34] glue_1.8.1          xfun_0.61           sys_3.4.3          
#> [37] knitr_1.52          farver_2.1.2        htmltools_0.5.9    
#> [40] labeling_0.4.3      rmarkdown_2.32      maketools_1.3.2    
#> [43] TMB_1.9.25          compiler_4.6.1      S7_0.2.2