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).
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).
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)
<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.
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.
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 nterms>1
Claiming a period longer than the record
<2 cycles cannot separate oscillation from trend
Require >=2 (ideally 3) cycles; cap max_period <= T_span/2
Assuming 24 h for cell cycle
Cell-cycle period is cell-type / condition dependent, usually != 24 h
Estimate the period; use circadian-rhythms only when 24 h is the hypothesis
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
