Bio Temporal Genomics Differential Rhythmicity

Version Compatibility

Reference examples tested with: limorhyde 1.0+, limma 3.50+, compareRhythms 1.0+, dryR (GitHub naef-lab), DODR 0.99+, CircaCompare 0.2+, R 4.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: differential rhythmicity estimates an INTERACTION (condition x time), so it needs MORE power than single-condition detection - the same matched, evenly-sampled, replicated grid must exist in BOTH conditions, and unmatched timepoints across conditions break the interaction model.

Differential Rhythmicity (comparing rhythms between conditions)

"How does this gene's rhythm change in the knockout / high-fat diet / aged tissue?" -> classify each feature's rhythm CHANGE relative to a reference condition into gain, loss, phase, or amplitude change, using a model that borrows strength across conditions.

• R: limorhyde() (sin/cos basis for a limma/edgeR/DESeq2 interaction model); compareRhythms() (direct gain/loss/change/same classification); dryseq() (dryR, BIC model selection); DODR, circacompare(), diffCircadian (targeted comparators)

This is the sibling of temporal-genomics/circadian-rhythms, which asks the single-condition question "is this feature rhythmic at 24h?" (cosinor/JTK/RAIN detection). Do not re-run detection here; this skill answers the categorically different question of how a rhythm DIFFERS between groups.

The governing principle: differential rhythmicity is a distinct question, and detect-then-Venn overestimates it

The #1 error is the detect-then-Venn anti-pattern: defining "genes that lost rhythm in the KO" as "rhythmic in WT MINUS rhythmic in KO" from two independently thresholded rhythmicity lists. This systematically OVERESTIMATES reprogramming, because near p~0.05 the two lists differ mostly from threshold noise, not biology - a gene at q=0.04 in WT and q=0.06 in KO is called "lost" when nothing changed (Pelikan 2022 FEBS J 289:6605). The fix is to never intersect two lists: fit ONE model spanning both conditions and test the condition x time INTERACTION directly, so a single calibrated test asks "did the rhythm change?" with strength borrowed across conditions.

The second load-bearing distinction is differential EXPRESSION vs differential RHYTHMICITY. A gene whose mean level shifts between conditions but whose oscillation is unchanged is differentially EXPRESSED, not differentially rhythmic. In the sin/cos framing these are orthogonal: the condition MAIN effect = differential expression (mean shift, adjusting for time), the condition:time INTERACTION = differential rhythmicity (amplitude/phase change). Reporting a main-effect hit as "rhythm reprogramming" is a common and wrong conflation.

The four canonical outcome classes (relative to a reference group)

Interpret phase and amplitude CHANGE only for features confidently rhythmic in AT LEAST ONE condition; a "phase shift" between two genes that are arrhythmic in both is noise. Amplitude-filter first (peak-to-trough or relative amplitude), then classify.

Method selection (which to pick and why)

Methodology here is evolving (LimoRhyde2 reframes around effect-size/posterior shrinkage; benchmarks disagree on the best classifier) - verify current best practice against the latest tool docs before committing to one method, and prefer a screen (LimoRhyde/dryR/compareRhythms) followed by targeted confirmation (CircaCompare/diffCircadian) on hits.

LimoRhyde + limma interaction test

Goal: Rank features by differential rhythmicity between two conditions while separately quantifying differential expression, in a single linear model.

Approach: Decompose measured time into a sin/cos basis with limorhyde(), fit condition*(time_cos+time_sin), moderated-F-test the two interaction coefficients for DR and the condition main effect for DE.

library(limorhyde); library(limma)

# limorhyde() decomposes measured time into a cosinor basis; prefix 'time_' names them time_cos, time_sin.
# period = 24h circadian; time is MEASURED (ZT/CT), not inferred pseudotime.
meta <- cbind(meta, limorhyde(meta$time, 'time_', period = 24))

# Differential RHYTHMICITY = condition:time interaction; differential EXPRESSION = condition main effect.
design <- model.matrix(~ condition * (time_cos + time_sin), data = meta)
fit <- eBayes(lmFit(expr, design))   # expr: features x samples, columns aligned to meta rows

dr_cols <- grep('conditionKO:time_', colnames(design), value = TRUE)  # the two interaction coefficients
dr <- topTable(fit, coef = dr_cols, number = Inf, sort.by = 'F')      # BH-adjusted adj.P.Val ranks DR
de <- topTable(fit, coef = 'conditionKO', number = Inf, sort.by = 'p')  # condition main effect = DE

The interaction F-test says a rhythm changed; it does not name the class. To assign gain/loss/phase/amplitude, fit a per-condition cosinor (or read limorhyde2 posterior estimates) and compare amplitudes and acrophases between conditions for the significant features. For count data, run the same design through voom+limma, edgeR, or DESeq2 with the sin/cos and interaction columns.

compareRhythms direct classification

Goal: Assign every feature directly to gain / loss / change / same relative to a reference condition, without intersecting two detection lists.

Approach: Pass a features x samples matrix plus an exp_design data.frame (numeric time, 2-level factor group), choose a method, and read the returned category per feature.

library(compareRhythms)

# exp_design: one row per sample; numeric 'time', factor 'group' with EXACTLY 2 levels (reference first).
# data: numeric matrix, rows = features (rownames = ids), columns = samples matching exp_design rows.
# method='mod_sel' = BIC model selection (no DE package); 'deseq2'/'edger'/'voom' for RNA-seq counts;
# 'limma' for log-microarray; 'dodr'/'cosinor' also available. amp_cutoff = peak-to-trough floor (>=1 group).
res <- compareRhythms(data, exp_design = exp_design, period = 24,
                      method = 'mod_sel', amp_cutoff = 0.5, criterion = 'bic')
# res: data.frame with id + category (gain / loss / change / same, relative to the reference group).

For >2 conditions, dryR does BIC model selection across all conditions at once: dryseq(counts, group, time) assigns each gene a rhythm-parameter-sharing model and returns per-condition amplitude/phase/mean. Prefer it over pairwise interaction tests when the design has three or more groups and a parsimonious classification is wanted.

The confounds a differential-rhythmicity claim must address

Reduced bulk amplitude can be loss of SYNCHRONY, not loss of per-cell rhythm. A bulk/tissue readout is the sum over many single-cell oscillators; if cells DESYNCHRONIZE (dephase) between conditions, the ensemble amplitude damps toward zero even though every cell still oscillates. Bulk "amplitude change" or "loss of rhythm" therefore has three indistinguishable causes: true loss of cell-autonomous rhythmicity, loss of inter-cell synchrony, or reduced single-cell amplitude. Report it as "reduced ensemble amplitude" and use single-cell or live-imaging assays to separate the causes; see single-cell/preprocessing for cell-level analysis.

A rhythm change under light-dark may be a driven change, not a clock change. Under an entraining LD cycle, an apparent rhythm can be masked (driven directly by light/feeding/temperature). A between-condition difference (e.g. a feeding-time or lighting manipulation) can shift the DRIVEN component without touching the endogenous clock. Only free-running (DD/constant) conditions license "the clock rewired"; under LD, a differential-rhythmicity hit may reflect a change in the environmental drive. Declare the light regime and use ZT (entrained) vs CT (free-running) accordingly.

Design constraints (shared with detection, stricter here)

• Matched sampling GRID across conditions: the SAME timepoints, evenly spaced, in every condition. Unmatched timepoints break the interaction model (the condition:time terms become non-estimable or confounded) - this is the constraint most often violated when two datasets are compared post hoc.

• =2 full cycles and >=6 (ideally 8-12) samples/cycle in EACH condition (the genome-scale rhythm-analysis design guidelines of Hughes et al. 2017 apply per condition). One cycle cannot separate a rhythm from a trend, so it certainly cannot compare rhythms.

• Replicates per timepoint in each condition. DR estimates an interaction (a difference of differences), so it is hungrier for power than detection; a design just adequate to detect a rhythm is usually too thin to confidently call a rhythm CHANGE.
• The amplitude filter is not optional. Classify phase/amplitude change only for features confidently rhythmic in at least one condition; significance without an effect-size floor over-detects (Laloum 2020 PLoS Comput Biol 16:e1007666), and the effect is worse for a difference test.

Common Errors (trap -> fix)

Related Skills

temporal-genomics/circadian-rhythms - Single-condition rhythm DETECTION and parameter estimation (cosinor/JTK/RAIN); run it first, this skill compares its results across conditions differential-expression/timeseries-de - Temporal differential expression (a monotone trend or between-timepoint change), which is differential EXPRESSION over time, not differential rhythmicity temporal-genomics/temporal-clustering - Group differentially-rhythmic genes by the shape of their change single-cell/preprocessing - Entry to cell-level analysis, to separate reduced ensemble amplitude (desynchrony) from true loss of per-cell rhythm

References

• Singer JM, Hughey JJ. 2019. LimoRhyde: a flexible approach for differential analysis of rhythmic transcriptome data. J Biol Rhythms 34(1):5-18. doi:10.1177/0748730418813785
• Weger BD, Gobet C, David FPA, et al. 2021. Systematic analysis of differential rhythmic liver gene expression mediated by the circadian clock and feeding rhythms (dryR). PNAS 118(3):e2015803118. doi:10.1073/pnas.2015803118
• Pelikan A, Herzel H, Kramer A, Ananthasubramaniam B. 2022. Venn diagram analysis overestimates the extent of circadian rhythm reprogramming (compareRhythms). FEBS J 289(21):6605-6621. doi:10.1111/febs.16095
• Thaben PF, Westermark PO. 2016. Differential rhythmicity: detecting altered rhythmicity in biological data (DODR). Bioinformatics 32(18):2800-2808. doi:10.1093/bioinformatics/btw309
• Parsons R, Parsons R, Garner N, Oster H, Rawashdeh O. 2020. CircaCompare: a method to estimate and statistically support differences in mesor, amplitude and phase, between circadian rhythms. Bioinformatics 36(4):1208-1212. doi:10.1093/bioinformatics/btz730
• Ding H, Meng L, Liu AC, et al. 2021. Likelihood-based tests for detecting circadian rhythmicity and differential circadian patterns in transcriptomic applications (diffCircadian). Brief Bioinform 22(6):bbab224. doi:10.1093/bib/bbab224
• Hughes ME, Abruzzi KC, Allada R, et al. 2017. Guidelines for genome-scale analysis of biological rhythms. J Biol Rhythms 32(5):380-393. doi:10.1177/0748730417728663
• Laloum D, Robinson-Rechavi M. 2020. Methods detecting rhythmic gene expression are biologically relevant only for strong signal. PLoS Comput Biol 16(3):e1007666. doi:10.1371/journal.pcbi.1007666