--- title: "Using RhyME for circadian methylation analysis" output: BiocStyle::html_document: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Using RhyME for circadian methylation analysis} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE, fig.crop = FALSE) ``` # 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: ```{r install, eval=FALSE} 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. ```{r data} library(RhyME) library(SummarizedExperiment) library(BiocParallel) example <- readRDS(system.file("extdata", "example_data.rds", package = "RhyME")) se <- example$se dim(se) head(colData(se)) ``` ## 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. ```{r fit} 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. ```{r results} 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] ``` 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: ```{r significant} detected <- result$RhyME_padj < 0.05 significant <- result[detected, , drop = FALSE] nrow(significant) ``` 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: ```{r summary} metadata(fit)$RhyME$fitSummary metadata(fit)$RhyME$ebPrior ``` ## 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. ```{r metacycle, results='hide'} 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 ``` ```{r align-metacycle, include=FALSE} meta_result <- meta_result[match(rownames(se), meta_result$CycID), ] ``` 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). ```{r method-comparison, echo=FALSE} pvalues <- data.frame( RhyME = result$RhyME_pvalue, JTK = meta_result$JTK_pvalue, LS = meta_result$LS_pvalue, MetaCycle = meta_result$meta2d_pvalue) padj <- sapply(pvalues, p.adjust, method = "BH") rownames(padj) <- rownames(se) comparison <- do.call(rbind, lapply(colnames(padj), function(method) { called <- !is.na(padj[, method]) & padj[, method] < 0.05 tp <- sum(called & truth) fp <- sum(called & !truth) data.frame(Method = method, TP = tp, FN = sum(truth) - tp, FP = fp, TN = sum(!truth) - fp, `Power (%)` = 100 * tp / sum(truth), `Observed FDP (%)` = if (sum(called)) 100 * fp / sum(called) else NA_real_, check.names = FALSE) })) knitr::kable(comparison, digits = c(0, 0, 0, 0, 0, 1, 2), row.names = FALSE, caption = "Known simulation truth versus calls at BH < 0.05.") ``` ## 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. ```{r plot, fig.width=8, fig.height=5.5} plotCosinorFit(fit, cpgs = selected_cpgs, ncol = 2) + ggplot2::theme_gray() ``` # Prepare real data for RhyME The study by [Ammala et al. (2025)](https://doi.org/10.1111/jsr.14438) provides longitudinal blood DNA methylation data through [OSF](https://osf.io/qx7w3/). 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. ```{r real-data} 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) head(colData(real_se)) ``` Alternatively, extract beta values from a preprocessed `minfi` `GenomicRatioSet`, retaining genomic coordinates and sample metadata: ```{r minfi-input, eval=FALSE} 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. ```{r fixed-default} se_independent <- example$se_independent fit_fixed <- fitCosinor(se_independent, BPPARAM = SerialParam()) metadata(fit_fixed)$RhyME$ebPrior ``` ## 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. ```{r covariates} 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 session, echo=FALSE} sessionInfo() ```