bio-temporal-genomics-periodicity-detection

name	bio-temporal-genomics-periodicity-detection
description	Discovers a periodic signal of UNKNOWN period in time-series omics data and puts a defensible significance on it, especially when sampling is IRREGULAR (dropped timepoints, pooled harvests) so FFT/Welch/JTK are invalid. Estimates the dominant period with Lomb-Scargle / generalized Lomb-Scargle (scipy, astropy), corroborates with autocorrelation, resolves transient/time-varying periodicity with the wavelet CWT (pywt), and screens genome-wide with false-alarm probabilities under BH FDR. Use when finding an oscillation whose period is not known a priori, analyzing cell-cycle or ultradian rhythms, or handling unevenly sampled time courses. Not for testing a KNOWN 24-hour rhythm (see temporal-genomics/circadian-rhythms).
tool_type	python
primary_tool	scipy

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

Method	Sampling	Estimates	Significance model	Fails when
Generalized Lomb-Scargle (astropy default, or scipy floating_mean=True)	Uneven OK	Global dominant period(s), amplitude, phase	Baluev analytic FAP (screen); bootstrap (few hits)	<2 cycles; strong harmonics; red noise mistaken for a peak
Classic Lomb-Scargle (scipy, pre-centered)	Uneven OK	Global dominant period	Permutation / Baluev via astropy	Nonzero-mean data if not centered; sparse points
Autocorrelation (statsmodels)	Even only	Fundamental period (coarse)	Bartlett bands (weak)	Uneven sampling; trend; poor resolution -> use as CHECK only
Wavelet CWT (pywt Morlet)	Even (interp small gaps)	Time-varying / transient period	Torrence-Compo AR(1) chi-square + COI mask	Edge claims (COI); needs a grid -> interpolation caveat
Welch PSD (scipy)	Even required	Smoothed global PSD	Segment-averaging variance reduction	Any gaps; long/close periods vs nperseg
Fisher's g-test (Wichert 2004)	Even only	Single dominant frequency, exact p	Exact analytic (white null)	Uneven sampling; red noise (white-null limitation)

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

Ridge extraction

Track argmax over scale at each time for the instantaneous dominant period, but only AFTER masking the COI and only for points above the red-noise contour; enforce continuity (a real ridge does not teleport between distant periods sample-to-sample). The global wavelet spectrum (time-average of power) is the wavelet analogue of a Fourier/LS spectrum, with its own chi-square test at reduced dof.

Welch PSD (evenly sampled only)

from scipy.signal import welch

# nperseg is the bias/variance knob: longer segments -> finer frequency resolution but
# fewer segments -> noisier PSD; shorter -> smoother but cannot resolve close/long periods.
# n//2 is a middling, defensible-but-arbitrary compromise; you cannot see a period longer
# than one segment. Welch CANNOT handle gaps -> feeding it interpolated data reintroduces
# interpolation bias. Even sampling only; otherwise use Lomb-Scargle.
freqs_w, psd = welch(values_even, fs=1/dt, nperseg=len(values_even)//2, detrend='constant')
periods_w = 1 / freqs_w[1:]

Genome-Wide Screening

Goal: Turn per-gene FAPs into a genome-wide error rate without inflating the hit list with harmonics or trends.

Approach: Compute one FAP per gene, control FDR across genes, and guard against harmonic contamination and non-white noise.

from statsmodels.stats.multitest import multipletests

# FAP is per-gene: thresholding 15,000 genes at FAP<0.01 yields ~150 false positives.
# Convert per-gene FAPs to q-values with Benjamini-Hochberg (assumes independence /
# positive dependence; gene-gene correlation makes it mildly conservative).
reject, qvals, _, _ = multipletests(fap_per_gene, method='fdr_bh')
n_periodic = int((qvals < 0.05).sum())

Screening subtleties:

• Harmonic contamination: a non-sinusoidal 24 h oscillation has real power at 12 h, 8 h, 6 h. The 12 h peak is a genuine harmonic, NOT an independent 12 h rhythm. A naive screen reports a phantom cohort of 12 h genes. Guard: check whether a putative P/2 peak co-occurs with a stronger peak at P; fit the fundamental plus harmonics jointly with astropy LombScargle(t, y, nterms=k) and attribute power correctly; distrust exact 2:1 period ratios.
• Permutation null validity: shuffling values across the fixed times destroys ALL temporal structure -> the null becomes WHITE noise. That is correct only if the alternative is "periodicity vs i.i.d. noise." Real omics noise is autocorrelated / red (1/f): a gene with a smooth trend or slow drift beats a white null and is falsely flagged periodic. On un-detrended, autocorrelated data a shuffle test is ANTI-CONSERVATIVE. Fixes: detrend first; use an AR(1) surrogate / block-bootstrap null that preserves short-range autocorrelation; or use the analytic Baluev FAP. State which null was used and what alternative it implies.
• Fisher's g-test (Wichert, Fokianos & Strimmer 2004): for EVENLY sampled series, g = (max periodogram ordinate)/(sum of ordinates) has a known exact distribution under the white-noise null, giving an exact per-gene p with no simulation - the canonical microarray cell-cycle screen. Even sampling only; shares the white-null limitation.
• Interpolate-then-FFT biases the spectrum red; never interpolate uneven data to run an even-sampling screen - use LS/GLS directly.
• Genuine oscillation vs 1/f vs trend: a real oscillation is a NARROW peak riding above the smooth 1/f background and recurring at harmonics; 1/f humps are broad and non-harmonic; a trend is indistinguishable from a very long period over <2 cycles. Require >=2-3 cycles, detrend, and test against a colored (AR(1)) null, not white.

Common Errors

Trap	Why it is wrong	Fix
Raw (nonzero-mean) values into scipy.signal.lombscargle	Classic LS assumes zero mean; the offset leaks into the sinusoid -> inflated power	floating_mean=True (GLS) OR pass y - y.mean(); precenter deprecated in 1.17
Ordinary frequency (1/period) in scipy	scipy takes ANGULAR omega=2pi/period; period off by 2pi	angular = 2*np.pi/periods; use astropy for ordinary-frequency grids
Pre-centering then handing to astropy	astropy already fits the mean (fit_mean=True) -> redundant / conceptual muddle	Pass RAW values to astropy; it is GLS by default
Ad-hoc np.linspace(f_min,f_max,1000) grid	Too-coarse spacing steps over the finite-width peak	autopower(samples_per_peak=5..10) or space finer than ~1/T_span
Reporting FAP without the grid	FAP grows with grid width / oversampling; numbers not comparable	Always report [f_min,f_max] and samples_per_peak
Interpolate-then-FFT/Welch on uneven data	Interpolation is a low-pass filter -> spurious red power, killed high-freq signal	Analyze the uneven series directly with LS/GLS
Wavelet scalogram with no COI mask	Edge coefficients are artifacts against padding; worse at long periods	Compute the COI (Morlet half-width sqrt(2)*scale), grey it out
mean + 2*SD wavelet significance	Wavelet power is not Gaussian; variance varies with scale	Torrence-Compo AR(1) background x chi2(0.95,2)/2 contour; report alpha
Shuffle-time permutation on un-detrended data	White null; a trend / red-noise gene beats it -> false "periodic"	Detrend first, or AR(1)/block-bootstrap surrogate, or Baluev FAP
ACF as the period estimator	Poor resolution; assumes even sampling; trend-fooled	Use ACF only as a corroborating check after detrending
Treating the 12 h (P/2) peak as an independent rhythm	It is the harmonic of a non-sinusoidal 24 h signal	Check co-occurrence with a stronger peak at P; fit

References

• Lomb 1976. Least-squares frequency analysis of unequally spaced data. Astrophys Space Sci 39(2):447-462.
• Scargle 1982. Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data. Astrophys J 263:835-853.
• VanderPlas 2018. Understanding the Lomb-Scargle Periodogram. Astrophys J Suppl Ser 236(1):16.
• Baluev 2008. Assessing the statistical significance of periodogram peaks. Mon Not R Astron Soc 385(3):1279-1285.
• Torrence & Compo 1998. A Practical Guide to Wavelet Analysis. Bull Amer Meteor Soc 79(1):61-78.
• Wichert, Fokianos & Strimmer 2004. Identifying periodically expressed transcripts in microarray time series data. Bioinformatics 20(1):5-20.

Related Skills

temporal-genomics/circadian-rhythms - Known-period (24 h) rhythm testing with cosinor and JTK_CYCLE temporal-genomics/temporal-clustering - Group genes by periodicity characteristics temporal-genomics/trajectory-modeling - Non-periodic trajectory fitting with GAMs