Tests and estimates rhythmicity at a PRE-SPECIFIED period (canonically 24h) in time-series omics using cosinor regression (CosinorPy), JTK_CYCLE/ARSER/Lomb-Scargle meta-analysis (MetaCycle meta2d), and non-parametric tests for asymmetric waveforms (RAIN, DiscoRhythm); estimates phase (acrophase), amplitude, and MESOR, and controls FDR with an effect-size (rAMP) filter against over-detection. Use when testing for 24-hour or other known-period oscillations in a single condition (circadian, feeding-fasting, or light-dark experiments) and estimating their phase/amplitude. Not for unknown-period discovery (see temporal-genomics/periodicity-detection) or comparing rhythms between conditions (see temporal-genomics/differential-rhythmicity).
Tests and estimates rhythmicity at a PRE-SPECIFIED period (canonically 24h) in time-series omics using cosinor regression (CosinorPy), JTK_CYCLE/ARSER/Lomb-Scargle meta-analysis (MetaCycle meta2d), and non-parametric tests for asymmetric waveforms (RAIN, DiscoRhythm); estimates phase (acrophase), amplitude, and MESOR, and controls FDR with an effect-size (rAMP) filter against over-detection. Use when testing for 24-hour or other known-period oscillations in a single condition (circadian, feeding-fasting, or light-dark experiments) and estimating their phase/amplitude. Not for unknown-period discovery (see temporal-genomics/periodicity-detection) or comparing rhythms between conditions (see temporal-genomics/differential-rhythmicity).
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.
CosinorPy 3.1 imports as from CosinorPy import cosinor, cosinor1, file_parser (capitalized package, lowercase submodules); older 0.x/1.x releases used the lowercase cosinorpy package name.
Known-Period Rhythm Testing
Governing principle: known-period TESTING, not unknown-period DISCOVERY
This skill answers "at THIS period (usually 24h), is a feature rhythmic, and what are its phase, amplitude, and MESOR?" - a hypothesis test plus parameter estimation at a period the analyst specifies. That is categorically different from asking "what period does this feature have?", which is a spectral DISCOVERY search (Lomb-Scargle periodogram, wavelets, FFT) handled by temporal-genomics/periodicity-detection. Conflating them is the field's most common conceptual error: opening a wide period window turns a test into a search and inflates false positives, because structured noise can always be "fit" better at SOME period in a broad window.
The load-bearing consequence: temporal conclusions are dominated by SAMPLING DESIGN, not by the algorithm. Nyquist (>=2 samples/cycle) is a mathematical FLOOR that only prevents aliasing - it gives zero robustness to noise and no ability to estimate phase/amplitude. Real detection needs >=6 (ideally 8-12) samples/cycle AND >=2 full cycles. Resolving that a period exists < estimating its phase < estimating its amplitude, in ascending sampling demand. A design good enough to say "yes, 24h" is usually too thin to trust its phase and far too thin to trust its amplitude.
Core Workflow
Declare the light regime (LD=entrained, use ZT; DD=free-running, use CT) and the period window (entrained: fix minper=maxper=24; free-running: allow ~22-26h for tau != 24h)
Prepare the time-series matrix (features x timepoints), decide log-vs-linear scale and any detrending BEFORE testing
Fit cosinor models or apply rhythmicity tests at the specified period
Each harmonic costs 2 df; with 6-8 pts/cycle a 3-component model is near-saturated and fits noise - use AIC/BIC or automatic model selection
JTK_CYCLE
Evenly sampled at integer-hour intervals; robust rank-based test; genome-scale speed
Correlates the series against reference cosines of all phases (Jonckheere-Terpstra + Kendall tau); best phase = matched reference
Requires EVEN integer sampling, no gaps; with few timepoints/1 replicate the tau null is DISCRETE so p-values are quantized and ANTI-conservative (source of "everything is rhythmic")
eJTK / BooteJTK
Short/sparse or few-replicate series where JTK p-values are untrustworthy; asymmetric/spiky waveforms
Umbrella/Mack-Wolfe (Jonckheere-Terpstra) test with SEPARATE rising and falling limbs
LOWER power than cosinor/JTK for genuinely symmetric sinusoids; gives a coarse phase/peak-shape, not clean amplitude CIs
MetaCycle meta2d
A robust consensus RANK across methods is wanted on a standard even design
Runs a subset of {ARS,JTK,LS}, combines p by Fisher's method -> meta2d_pvalue (BH -> meta2d_BH.Q), averages period, circular-averages phase
Fisher assumes INDEPENDENT p; ARS/JTK/LS on the same data are correlated, so meta2d_BH.Q is NOT a literal FDR (read it as a rank aid); analysisStrategy='auto' SILENTLY drops ARS/JTK on uneven/replicated data (may run LS only); averaged period is meaningless when methods disagree
For between-condition comparison, see temporal-genomics/differential-rhythmicity - none of the single-condition tests above answer it correctly.
Design constraints (upstream of any method; non-negotiable)
=2 full cycles (48h circadian minimum; 3 cycles / 72h improves power and reveals damping). One cycle cannot distinguish an oscillation from a monotone trend or a single transient.
=6, ideally 8-12+, samples/cycle. 2-4h spacing is standard; 1-2h is needed to resolve waveform shape or fast harmonics. Sparse designs are exactly where JTK's calibration fails.
=2-3 biological replicates/timepoint. Single-replicate designs cripple FDR calibration (no within-timepoint variance; empirical-null/bootstrap corrections cannot work). Replication in TIME and AT a timepoint buy different things - do not trade all of one for the other.
Harvest-ORDER confound (the silent killer): collecting/extracting/sequencing timepoints in temporal order aliases any drift (reagent lots, RIN, lane position) PERFECTLY onto ZT and manufactures spurious 24h rhythms. No rhythmicity test detects this. Fix by DESIGN: randomize processing order, balance replicates across batches, model batch as a covariate (trivial in a limma/DESeq2 design). It cannot be repaired analytically because batch and the rhythm are the same axis.
CosinorPy (Python)
Goal: Test each feature for rhythmicity at a known period and estimate amplitude, relative amplitude, acrophase, and MESOR with FDR control.
Approach: Fit cosine curves per feature with fit_group (batch), use its built-in BH q column (or recompute BH over a chosen correction set), then filter on both q and relative amplitude.
Single- and multi-component fit
Fits y = M + A*cos(2*pi*t/T + phi) where M = MESOR (rhythm-adjusted midline, NOT the arithmetic mean unless sampling is balanced), A = amplitude, phi = acrophase stored as atan2(-gamma, beta) (usually negative).
from CosinorPy import cosinor, cosinor1, file_parser
df = file_parser.read_csv('expression_timecourse.csv') # long format: columns x (time), y (value), test (feature id)# Single-component (sinusoidal). period=24: standard circadian period in hours.# fit_me returns a 5-tuple: (results, statistics, rhythm_params, X_test, Y_fit_test).
single = cosinor.fit_me(df[df['test'] == 'Arntl']['x'].values,
df[df['test'] == 'Arntl']['y'].values,
period=24, n_components=1)
# Multi-component adds harmonics for non-sinusoidal shape; add ONLY with dense sampling + a biological reason.# fit_me takes a SINGLE n_components (an int); n_components=2 adds one 12h harmonic to the 24h fundamental.
two_comp = cosinor.fit_me(df[df['test'] == 'Dbp']['x'].values,
df[df['test'] == 'Dbp']['y'].values,
period=24, n_components=2)
# To let CosinorPy PICK the harmonic order by information criterion, fit a range with fit_group over a# candidate list, then select per feature with get_best_models (do not pass a list to fit_me).
group_multi = cosinor.fit_group(df, period=24, n_components=[1, 2, 3], plot=False)
best_models = cosinor.get_best_models(df, group_multi, n_components=[1, 2, 3])
Batch analysis with built-in q-values
Goal: Score every feature genome-wide and keep confident, high-amplitude oscillators.
Approach:fit_group returns per-feature statistics INCLUDING a BH-adjusted q column; add an rAMP effect-size filter on top of q.
import numpy as np
from statsmodels.stats.multitest import multipletests
# fit_group returns columns: test, period, n_components, p, q, p_reject, q_reject, RSS, R2, R2_adj,# log-likelihood, amplitude, acrophase, mesor, peaks, heights, troughs, heights2, ME, resid_SE.
results = cosinor.fit_group(df, period=24, n_components=1, plot=False)
# 'q' is already BH-adjusted across the fitted group. Recompute BH only if the correction SET should differ# (e.g. exclude non-expressed features first). Default multipletests method is Holm-Sidak, so pass fdr_bh explicitly.
valid = results['p'].notna()
results.loc[valid, 'q_bh'] = multipletests(results.loc[valid, 'p'], method='fdr_bh')[1]
# rAMP = amplitude / MESOR normalizes out expression level so calls are comparable across features.# rAMP > 0.1 (>=10% of baseline) is a conventional biological-relevance floor - sweep it, do not treat as law.
results['rAMP'] = results['amplitude'] / results['mesor']
rhythmic = results[(results['q'] < 0.05) & (results['rAMP'] > 0.1)]
Goal: Get group-level amplitude/phase with CIs that propagate BETWEEN-subject variance, instead of pseudoreplicating.
Approach: Fit one cosinor per subject and combine the estimates - pooling all subjects' points into one fit understates uncertainty.
# cosinor1.population_fit_cosinor returns a DICT with keys: test, names, values, means, confint (nested amp/acr/mesor CIs),# p_value, p_amp, p_acr, p_mesor (all underscore; e.g. pop['confint']['amp'], pop['p_amp']).
pop = cosinor1.population_fit_cosinor(subject_df, period=24, plot_on=False)
# cosinor1.population_fit_group(df, period=24) batches this across groups; cosinor1.population_test_cosinor_pairs# compares two populations' rhythms (a differential-rhythmicity test on replicated data).
Convert acrophase to peak-hour with peak_h = (-acrophase) * T / (2*pi) % T; sanity-check against a known clock gene (mouse liver Arntl/Bmal1 peaks ~CT22-0, Nr1d1 ~CT4-6, Dbp ~CT8-10).
MetaCycle meta2d (R)
Goal: Produce a robust consensus rhythmicity rank on an evenly sampled design.
Approach: Run meta2d over {JTK,ARS,LS}; read meta2d_BH.Q as a ranking aid (Fisher over correlated nulls, not a literal FDR), and distrust the averaged period/phase when constituents disagree.
library(MetaCycle)# minper=maxper=24 for entrained (LD) data; 22-26 for free-running (DD) where tau != 24h.# timepoints must match column order. ARS/JTK need EVEN integer sampling with no missing values / no replicates;# analysisStrategy='auto' silently drops ineligible methods (may leave LS only) - check the per-method columns.# timepoints span 0-68h at 4h resolution: >=2 full 24h cycles at ~6 samples/cycle (the design floor this skill sets).
meta2d(infile ='expression_matrix.csv', filestyle ='csv', outdir ='metaout',
timepoints = seq(0,68, by =4), cycMethod =c('JTK','ARS','LS'),
minper =24, maxper =24, outputFile =TRUE, outRawData =FALSE)
res <- read.csv('metaout/meta2d_expression_matrix.csv')# meta2d_pvalue (Fisher-combined), meta2d_BH.Q, meta2d_period, meta2d_phase (hours from ZT0, peak time),# meta2d_Base (baseline/MESOR), meta2d_AMP, meta2d_rAMP (= AMP/Base). Filter on rank AND relative amplitude.
rhythmic <- res[res$meta2d_BH.Q <0.05& res$meta2d_rAMP >0.1,]
Approach: Transpose to one-row-per-timepoint, declare replicate count, adjust p for multiple testing.
library(rain)# x needs ONE ROW PER TIMEPOINT (transpose a features x timepoints matrix). deltat = sampling interval (h).# nr.series = replicates per timepoint (interleaved r1t1,r2t1,r1t2,...). method='independent' vs 'longitudinal'# sets replicate handling. peak.border controls the allowed rising-fraction (asymmetry) window.
res <- rain(t(expression_mat), period =24, deltat =4, nr.series =2, method ='independent')
res$q <- p.adjust(res$pVal, method ='BH')# output columns: pVal, phase, peak.shape, period
rhythmic <- res[res$q <0.05,]
DiscoRhythm (R/Bioconductor)
Goal: Run Cosinor/JTK/LS/ARS under one interface with built-in QC/PCA (scripted or Shiny).
library(DiscoRhythm)
se <- discoGetSimu(TRUE)# bundled demo SummarizedExperiment
disco <- discoBatch(se, osc_method ='CS', report =NULL, osc_period =24)# osc_method='CS'=Cosinor; report=NULL skips the HTML report
Comparing rhythms BETWEEN conditions (differential rhythmicity: gain/loss/phase-shift/amplitude-change with LimoRhyde/dryR/compareRhythms, and the detect-then-Venn anti-pattern) is a distinct analysis - see temporal-genomics/differential-rhythmicity. Do NOT infer "genes that lost rhythm in the KO" by subtracting two independently thresholded single-condition rhythm lists.
Common Errors (trap -> fix)
Trap
Fix
Opening a wide period window on a known-period test
Fix minper=maxper=24 (entrained) or 22-26h (free-running); a wide window is discovery, not testing, and inflates false positives
Trusting JTK p-values / BH-Q from a single-replicate sparse design
Expect anti-conservative, quantized p-values; use eJTK/BooteJTK (empirical/bootstrap null) and inspect the genome-wide p-value HISTOGRAM before believing FDR
Reading meta2d_BH.Q as a literal FDR
Fisher integration over correlated ARS/JTK/LS p-values is not calibrated; use it as a consensus RANK and distrust averaged period/phase when methods disagree
Calling reduced BULK amplitude "arrhythmic"
Ensemble amplitude damps from cell DESYNCHRONY too; report "reduced ensemble amplitude" and use single-cell or imaging assays to separate loss-of-rhythm vs loss-of-synchrony
Claiming an "endogenous circadian rhythm" from LD data
LD rhythms can be light/feeding-DRIVEN (masking); endogeneity requires free-running (DD/constant) conditions. Diurnal != circadian. Use ZT for entrained, CT for free-running
Claiming a rhythm is "clock-CONTROLLED" from wild-type data alone
Persistence in DD proves endogeneity, not clock control; genetic dependence needs a clock-gene perturbation (compare WT vs clock-mutant, see temporal-genomics/differential-rhythmicity)
Mixing phase units/conventions (radians vs hours, +phi vs -phi, ZT vs CT)
State the convention; convert CosinorPy acrophase via peak_h = (-acrophase)*T/(2*pi) % T; sanity-check against a known clock gene's phase
Ranking features by RAW amplitude across the genome
Raw amplitude scales with expression and normalization; use relative amplitude (AMP/MESOR) or peak-to-trough fold-change for cross-feature comparison and the amplitude filter
Reporting significant rhythms with NO effect-size filter
Significance alone over-detects (Laloum 2020); add an rAMP/fold-change cutoff and report the amplitude DISTRIBUTION of the hit list, not just the count
Trusting phase/amplitude POINT estimates for near-threshold features
Estimation is unreliable where detection is marginal; interpret parameters only for confidently rhythmic features
Overfitting with n_components=3 on 6-8 points/cycle
Harmonics cost 2 df each; use AIC/BIC or automatic model selection; add harmonics only with dense sampling and a biological reason
Harvest-order drift confounded with ZT
Randomize PROCESSING order, balance replicates across batches, model batch as a covariate; no rhythmicity test detects this
Feeding replicates to cosinor as one pooled single-fit
Use population-mean cosinor (subject = replication unit) so between-subject variance enters the CI and the test
The over-detection controversy (state it as live)
Laloum & Robinson-Rechavi (2020) showed that across seven popular methods (ARS, LS, RAIN, JTK, eJTK, GeneCycle, meta2d) rhythm calls are consistent and biologically meaningful ONLY for strong-amplitude signals; weak-signal calls are method-dependent and largely non-functional. There is no consensus "correct" method. The pragmatic (not full) response: (1) require an amplitude/rAMP effect-size filter IN ADDITION to FDR; (2) prefer methods with calibrated empirical nulls (eJTK, BooteJTK) over raw JTK on sparse data; (3) verify the genome-wide p-value histogram is roughly uniform with a spike near 0 before trusting any q. Report the amplitude distribution of the hit list, not just "N% of the transcriptome is rhythmic."
Parameter Guide
Parameter
Typical value
Rationale
Period
24h (12h for ultradian)
Specified a priori; this is a test, not a search
Period window
24 (LD) / 22-26 (DD)
Entrained locks to 24h; free-running tau != 24h. Wide windows inflate false positives
Sampling interval
2-4h
Nyquist (<=12h) is a floor, not a target; shape resolution needs 1-2h
Cycles
>=2 (>=3 better)
One cycle cannot separate rhythm from trend/transient
Samples/cycle
>=6 (8-12+ better)
Six gives stable fit df; more resolves waveform and calibrates FDR
Replicates/timepoint
>=2-3
Single replicate has no within-timepoint variance; FDR miscalibrates
FDR threshold
q < 0.05
Necessary but not sufficient; always pair with an amplitude filter
Relative amplitude
rAMP > 0.1
>=10% of baseline as a biological-relevance floor; a convention to sweep, not a law
Related Skills
temporal-genomics/differential-rhythmicity - Comparing rhythms between conditions (gain/loss/phase/amplitude change)
temporal-genomics/periodicity-detection - Unknown-period discovery with Lomb-Scargle and wavelets
temporal-genomics/temporal-clustering - Group rhythmic genes by phase/shape
differential-expression/timeseries-de - Temporal differential expression (a monotone trend, not rhythmicity)
data-visualization/heatmaps-clustering - Circular phase heatmaps and phase-ordered maps
References
Hughes ME, Hogenesch JB, Kornacker K. 2010. JTK_CYCLE: an efficient nonparametric algorithm for detecting rhythmic components in genome-scale data sets. J Biol Rhythms 25(5):372-380. doi:10.1177/0748730410379711
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
Thaben PF, Westermark PO. 2014. Detecting rhythms in time series with RAIN. J Biol Rhythms 29(6):391-400. doi:10.1177/0748730414553029
Wu G, Anafi RC, Hughes ME, Kornacker K, Hogenesch JB. 2016. MetaCycle: an integrated R package to evaluate periodicity in large scale data. Bioinformatics 32(21):3351-3353. doi:10.1093/bioinformatics/btw405
Yang R, Su Z. 2010. Analyzing circadian expression data by harmonic regression based on autoregressive spectral estimation (ARSER). Bioinformatics 26(12):i168-i174. doi:10.1093/bioinformatics/btq189
Hutchison AL, Maienschein-Cline M, Chiang AH, et al. 2015. Improved statistical methods enable greater sensitivity in rhythm detection for genome-wide data (eJTK). PLoS Comput Biol 11(3):e1004094. doi:10.1371/journal.pcbi.1004094
Cornelissen G. 2014. Cosinor-based rhythmometry. Theor Biol Med Model 11:16. doi:10.1186/1742-4682-11-16
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
Mei W, Jiang Z, Chen Y, Chen L, Sancar A, Jiang Y. 2021. Genome-wide circadian rhythm detection methods: systematic evaluations and practical guidelines. Brief Bioinform 22(3):bbaa135. doi:10.1093/bib/bbaa135
Moškon M. 2020. CosinorPy: a python package for cosinor-based rhythmometry. BMC Bioinformatics 21:485. doi:10.1186/s12859-020-03830-w
