GPTomics/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)

| Class | What changed | Interaction signature | |-------|--------------|-----------------------| | Gain of rhythm | arrhythmic in reference, rhythmic in test | interaction significant; reference amplitude ~0 | | Loss of rhythm | rhythmic in reference, arrhythmic in test | interaction significant; test amplitude ~0 | | Phase change | rhythmic in both, peak time shifted | interaction significant; amplitudes similar, acrophases differ | | Amplitude change | rhythmic in both, oscillation damped/amplified | interaction significant; same phase, amplitudes differ | | Unchanged-rhythmic (same) | rhythmic in both, same amplitude+phase | interaction NOT significant; both rhythmic | | Arrhythmic-in-both | flat in both | neither main effect nor interaction; excluded by amplitude filter |

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)

| Method | Pick when | Mechanism | Fails / caveat | |--------|-----------|-----------|----------------| | LimoRhyde + limma/edgeR/DESeq2 | 2 conditions; want to fold DR into a standard DE pipeline with covariates/batch; count or microarray data | limorhyde() adds sin/cos time columns; condition:time interaction = DR, condition main effect = DE | Two-condition framing; the interaction test flags THAT a rhythm changed, not which class - read per-condition amplitude/phase to classify | | dryR (dryseq) | >=2 conditions; want a parsimonious per-gene MODEL assignment (shared vs independent rhythm parameters) | BIC model selection over a family of shared/independent-parameter models, tailored to RNA-seq noise | Model-selection categories depend on the model family and BIC penalty; needs enough timepoints for BIC to discriminate models | | compareRhythms | want gain/loss/change/same DIRECTLY, built to replace the Venn approach; microarray or RNA-seq | wraps model-selection (mod_sel) or hypothesis tests (dodr/limma/voom/deseq2/edger/cosinor); classifies vs reference | Two groups only; mod_sel needs no DE package but deseq2/edger/voom do; a feature must clear amp_cutoff in >=1 group to be reported | | DODR | direct two-condition differential-rhythmicity test on already-detected rhythmic features | robust/rank comparison of rhythm shape (amplitude, phase, signal-to-noise) between conditions | Tests differential rhythmicity given rhythmicity; pre-filter to features rhythmic in >=1 group first | | CircaCompare | a FEW targeted genes; want explicit estimates + p-values for the mesor/amplitude/phase DIFFERENCE | non-linear regression fitting both curves jointly with difference parameters | Compares two groups only if BOTH are rhythmic (amplitude non-zero); not a genome-scale screen | | diffCircadian | a few genes; want likelihood-ratio tests separating differential amplitude vs phase vs basal vs fit | likelihood-based tests (LR_diff) on two conditions | Two conditions; targeted rather than transcriptome-wide throughput |

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)

| Trap | Fix | |------|-----| | Defining "lost rhythm" as (rhythmic in WT) minus (rhythmic in KO) via two thresholded lists | Fit one model across both conditions and test the condition:time INTERACTION (LimoRhyde/dryR/compareRhythms); the Venn approach overestimates reprogramming (Pelikan 2022) | | Calling a condition MAIN-effect (mean-shift) hit "differential rhythmicity" | Main effect = differential EXPRESSION; differential RHYTHMICITY is the condition:time INTERACTION. Report them separately | | Unmatched or unevenly-spaced timepoints across conditions | Use a matched, evenly-sampled grid in every condition; unmatched times make the interaction terms non-estimable or confound them with condition | | Interpreting a "phase shift" for features arrhythmic in both conditions | Amplitude-filter FIRST; classify phase/amplitude change only for features confidently rhythmic in >=1 condition | | Calling reduced BULK amplitude "loss of rhythm" | Ensemble amplitude damps from cell DESYNCHRONY too; report "reduced ensemble amplitude" and separate with single-cell or imaging assays | | Claiming the clock "rewired" from LD (entrained) data | An LD rhythm change can be a driven/masking change (light/feeding); endogenous rewiring needs free-running (DD) conditions | | Running DR on a design too sparse to even DETECT a rhythm | DR needs MORE power than detection (it estimates an interaction); ensure >=2 cycles, >=6/cycle, replicates in EACH condition before comparing | | Reading the interaction F-test as the CLASS | The interaction says a rhythm changed, not which class; fit per-condition amplitude/phase (or LimoRhyde2 posteriors) to assign gain/loss/phase/amplitude | | Using compareRhythms deseq2/edger/voom on already-normalized log data | Those methods expect RAW counts; use mod_sel/limma/cosinor for normalized or microarray data |

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