Bio Temporal Genomics Periodicity Detection

Version Compatibility

Reference examples tested with: numpy 2.2+, scipy 1.15+ (lombscargle floating_mean present), astropy 8.0+, PyWavelets 1.8+, statsmodels 0.14+, matplotlib 3.8+

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

• Python: pip show <package> then help(module.function) to check signatures
• scipy: scipy.signal.lombscargle changed in 1.17 (precenter deprecated -> removed 1.19; use floating_mean=True or pre-center). Check scipy.__version__.

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

Periodicity Detection

"Find a periodic pattern of unknown period in my time-series data" -> Estimate the dominant period, attach a false-alarm probability, and (for transient signals) localize it in time. This is the complement of temporal-genomics/circadian-rhythms, which TESTS a known 24 h period; here the period is unknown and sampling is often uneven.

Governing Principle

Unknown-period discovery is an ESTIMATION problem stacked on a DETECTION problem, and uneven sampling corrupts both. (1) Estimation: the period is a continuous quantity with a confidence region, not a yes/no on 24 h. (2) Detection: a periodogram peak is a random variable even under pure noise, and the MAX over a frequency grid is extreme-value-distributed, so significance is the False Alarm Probability (FAP) of that max, NOT the raw peak height. (3) Sampling design dominates: regular-grid methods (FFT, Welch, JTK, RAIN) assume equal spacing; genomics rarely delivers it (a timepoint fails QC, harvests are pooled, sampling is denser early). For uneven sampling the Lomb-Scargle least-squares periodogram replaces the FFT.

The single most important consequence: do NOT interpolate-then-FFT. Interpolation is a low-pass filter that injects spurious low-frequency power, suppresses power near the average Nyquist, and biases the whole spectrum red. Analyze the uneven series directly with Lomb-Scargle; interpolate only when a wavelet transform forces a grid, and flag it as an assumption.

Method Selection

Decision spine: uneven times -> GLS + Baluev FAP -> BH FDR for a genome-wide screen; suspected transient -> CWT with COI + red-noise chi-square; even and want a quick classical p -> Fisher's g; ACF/Welch are even-sampling sanity checks, never the primary estimate. Re-verify the astropy FAP methods and scipy floating_mean/normalize semantics against the installed versions before relying on them.

Lomb-Scargle Periodogram

Goal: Estimate the dominant period of an unevenly sampled series and attach a false-alarm probability.

Approach: Fit a sinusoid (plus a floating offset) at each trial frequency over a grid whose min/max and density are set from the record length, then read the peak period and its FAP. Prefer the generalized LS (floating mean) so a poorly-determined offset cannot leak into the amplitude.

Lomb-Scargle is NOT an FFT with jitter tolerance; it is a least-squares sinusoid fit at each frequency, invariant to time-origin, and (for even sampling) equivalent to the classical periodogram. Classical LS assumes zero-mean data; expression is always positive, so a nonzero mean leaks into the sinusoid and inflates power. Two correct fixes: pre-center (y - y.mean(), plug-in mean, old-school) or fit a per-frequency floating mean (generalized LS, the modern default).

scipy: classic (pre-centered) vs generalized

import numpy as np
from scipy.signal import lombscargle

# scipy uses ANGULAR frequency (rad/time): omega = 2*pi/period. Mixing ordinary
# frequency (1/period) is the #1 units bug -> periods off by 2*pi.
periods = np.linspace(6.0, 72.0, 2000)   # biologically plausible periods
angular = 2 * np.pi / periods

# Classic LS assumes zero mean: pre-center or power is inflated/distorted.
power_classic = lombscargle(times, values - values.mean(), angular, normalize=True)

# Generalized LS: fit a floating offset per frequency (Zechmeister & Kuerster).
# Preferred with sparse/uneven data; pass RAW values. floating_mean added in scipy
# 1.15; before that, pre-center. `precenter` is deprecated in 1.17, removed in 1.19.
power_gls = lombscargle(times, values, angular, normalize=True, floating_mean=True)

dominant_period = periods[np.argmax(power_gls)]

astropy: generalized LS by default, with real FAP

Goal: Get a peak period plus a Baluev false-alarm probability, using a grid astropy sizes automatically.

Approach: astropy's LombScargle fits a floating mean by default; pass RAW values, let autopower build a correctly oversampled grid, then convert the peak power to a FAP.

from astropy.timeseries import LombScargle

# fit_mean=True and center_data=True by DEFAULT -> generalized LS out of the box.
# Pass RAW values; pre-centering here is redundant. Uses ORDINARY frequency
# (cycles/time), NOT angular -> do not reuse a scipy omega-grid.
ls = LombScargle(times, values)          # add dy=sigma for heteroskedastic weighting
freq, power = ls.autopower(
    minimum_frequency=1 / 72.0,          # <= 1/record-span; need >=2 cycles to trust
    maximum_frequency=1 / 6.0,           # pseudo-Nyquist; uneven times can exceed 1/(2*dt_mean)
    samples_per_peak=10)                 # oversample: a peak has finite width ~1/T_span

best_period = 1 / freq[np.argmax(power)]

# FAP = P(noise produces a peak this high anywhere on the grid). It is GRID-DEPENDENT
# (rises with wider range / more oversampling), so always report the grid.
fap = ls.false_alarm_probability(power.max(), method='baluev')   # analytic, fast, screen-safe
levels = ls.false_alarm_level([0.1, 0.05, 0.01], method='baluev')
# method='bootstrap' (method_kwds={'n_bootstraps': 1000}) is most faithful but ~1000x cost

Frequency-grid design (where peaks get missed)

The grid is not cosmetic: a bad grid silently loses signals.

• Minimum frequency comes from the record length: a period longer than ~the span cannot be claimed. Require >=2 full cycles (f_min ~ 1/(T_span/2)); 1 cycle cannot be told from a trend.
• Maximum frequency is a pseudo-Nyquist. There is no single Nyquist for irregular times; the average 1/(2*dt_mean) is only a heuristic, and genuine unevenness permits probing far above it. Set nyquist_factor to a few, understanding it is a convention.
• Oversampling: a peak has finite width ~1/T_span; a grid coarser than a fraction of that steps over the true peak and under-reports its height and location. Use samples_per_peak 5-10 (astropy's autopower does this) rather than an ad-hoc linspace(...,1000).
• Leakage / aliasing: gapped sampling has a window function that convolves the true spectrum, creating sidelobes and aliases at 1/(1/P +/- 1/P_sampling). Eyeball the window-function periodogram (LS of a constant signal at the same times) to see where sampling itself manufactures peaks.

Autocorrelation (a sanity check, not an estimator)

Goal: Corroborate an LS-estimated period, not measure it.

Approach: Detrend, then look for a harmonic comb of ACF peaks at multiples of the period; treat a single band-crossing as weak evidence.

from statsmodels.tsa.stattools import acf

# statsmodels assumes EVEN sampling and takes NO time vector -> invalid on the
# uneven data that motivates this skill. Detrend first: a trend gives a slowly
# decaying ACF that mimics periodicity.
vals_dt = values_even - np.polyval(np.polyfit(np.arange(len(values_even)), values_even, 1), np.arange(len(values_even)))
acf_vals, confint = acf(vals_dt, nlags=len(vals_dt) // 2, alpha=0.05)   # index 0 is lag-0 = 1.0

ACF is a WEAK period estimator: its resolution is quantized to the sampling step (period only to +/- one step), it requires even sampling, and trends fool it. A periodic signal repeats at lags P, 2P, 3P; the partial ACF (PACF) helps isolate the fundamental from harmonic echoes. Use it as a robustness box after detrending, never as the primary estimate.

Wavelet CWT (transient / time-varying periodicity)

Global LS/Welch report one spectrum for the whole record and blur a signal that oscillates only early (cell-cycle synchrony decaying as cells desynchronize) or shifts period after a stimulus. The CWT resolves power in the time x period plane, at the cost of lower frequency precision and edge artifacts.

Goal: Map how the dominant period changes over time and mark which of that map is trustworthy and significant.

Approach: Transform with a complex Morlet, mask the cone of influence (edge artifacts), and test power against a Torrence-Compo AR(1) red-noise background rather than an ad-hoc mean + 2*SD.

import pywt
import numpy as np
from scipy.stats import chi2

# Complex Morlet 'cmorB-C' (bandwidth B, center freq C): complex -> amplitude AND phase.
wavelet = 'cmor1.5-1.0'
C = pywt.central_frequency(wavelet)        # 1.0 for cmor1.5-1.0
dt = times[1] - times[0]                    # requires even sampling
periods_to_test = np.arange(6, 49, 0.5)
scales = C * periods_to_test / dt           # scale = C * period / dt ; verify below
# sanity: pywt.scale2frequency(wavelet, scales) / dt  ->  1/periods_to_test

# Pass sampling_period, else returned freqs are in per-sample units (period axis off by dt).
coeffs, freqs = pywt.cwt(signal, scales, wavelet, sampling_period=dt)
power = np.abs(coeffs) ** 2                  # power = |coefficients|^2

Cone of influence (COI) - mask edge artifacts BEFORE reading ridges

Near each end of a finite record a wavelet overlaps the edge and its coefficients are computed against padding - edge artifacts, not signal. The COI widens for LONGER periods (larger scales reach farther from the edge). For the Morlet the e-folding half-width is ~sqrt(2)*scale in time; since period = scale*dt/C and C=1, a period is inside the COI (unreliable) where sqrt(2)*period > distance-to-nearest-edge. Grey out / mask the two triangular corners before ridge extraction or peak reading. Claiming a long-period oscillation that lives only in the first/last fraction of the record, inside the COI, is a classic false positive.

n = len(signal)
edge_dist = np.minimum(np.arange(n), np.arange(n)[::-1]) * dt   # time to nearest edge
coi_max_period = edge_dist / np.sqrt(2)                          # longest reliable period per time
inside_coi = periods_to_test[:, None] > coi_max_period[None, :]
power_masked = np.where(inside_coi, np.nan, power)               # ignore masked cells downstream

Red-noise significance (Torrence & Compo 1998), NOT mean + 2*SD

mean(power) + 2*SD(power) is indefensible: wavelet power under noise is not Gaussian and its variance changes with scale. Model the null as red noise - an AR(1) process with lag-1 autocorrelation alpha estimated from the data (white noise, alpha=0, is too permissive for intrinsically red / 1/f omics). The theoretical background is P_k = (1-alpha^2)/(1 - 2*alpha*cos(2*pi*dt/period) + alpha^2); local power normalized by this background is chi-square with 2 dof (complex wavelet), so the 95% contour is (noise level) * P_k * chi2(0.95, 2)/2. Peaks poking through it are significant. Report the alpha assumed. pywt's Morlet power is not in variance units, so calibrate one wavelet constant k from the ROBUST (median) noise level of the scalogram - the median is dominated by noise cells, so it estimates the white-noise power per unit variance without the strong signal cells inflating it.

alpha = np.corrcoef(signal[:-1], signal[1:])[0, 1]   # lag-1 autocorrelation
variance = signal.var(ddof=1)
Pk = (1 - alpha**2) / (1 - 2*alpha*np.cos(2*np.pi*dt/periods_to_test) + alpha**2)
k = np.nanmedian(power_masked / (variance * Pk[:, None]))   # robust wavelet power constant
sig95 = variance * Pk * k * chi2.ppf(0.95, df=2) / 2 # per-period 95% level, 2 dof
significant = power_masked > sig95[:, None]          # inside COI is NaN -> False