Bio Temporal Genomics Circadian Rhythms

Version Compatibility

Reference examples tested with: CosinorPy 3.1 (requires numpy<2.0 - v3.1 calls the removed np.round_), pandas 2.2+, statsmodels 0.14+, MetaCycle 1.2+, RAIN 1.x (Bioconductor), DiscoRhythm 1.x (Bioconductor).

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

1. 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)
2. Prepare the time-series matrix (features x timepoints), decide log-vs-linear scale and any detrending BEFORE testing
3. Fit cosinor models or apply rhythmicity tests at the specified period
4. Extract parameters: amplitude, relative amplitude (rAMP), phase (acrophase), MESOR, p-value
5. Control FDR (BH), then apply an EFFECT-SIZE filter (rAMP / fold-change) - significance alone over-detects
6. For between-condition questions, fit a differential-rhythmicity model (never intersect two separate rhythm lists)

Method Selection (which to pick and why)

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)]

Population-mean cosinor (replicated / multi-subject)

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, ]

RAIN (R/Bioconductor)

Goal: Detect ASYMMETRIC waveforms (fast induction, slow decay) that cosinor/JTK miss.

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.