GPTomics/temporal-genomics/trajectory-modeling

Version Compatibility

Reference examples tested with: mgcv 1.9+, tradeSeq 1.16+, segmented 2.0+, ruptures 1.1+, numpy 1.26+, pandas 2.2+

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.

Note: the modeled quantity must be on a scale the family assumes. A Gaussian GAM is valid only on variance-stabilized/log-transformed expression (vst, rlog, log-CPM); raw RNA-seq counts require family=nb() with a offset(log(library_size)). Successive timepoints are correlated, so a plain gam() (which assumes independent residuals) inflates smooth-term significance unless the AR structure is modeled or the residual ACF is checked.

Temporal Trajectory Modeling

"Fit smooth curves to my gene expression over real time, compare trajectories, and find abrupt shifts" -> model a continuous function of MEASURED time f(time), test whether it changes / differs between conditions, and locate discrete regime changes.

• R: mgcv::gam()/gamm()/bam() for penalized-spline GAMs; segmented::segmented() for slope breaks
• Python: ruptures for level/distribution changepoints

The governing principle: measured time is not pseudotime, and the model class must match the mechanism

This skill models trajectories where the x-axis is ACTUAL experimental time (hours, days, developmental stage) or a pseudobulk value aggregated over real time. Time is measured and shared across every sample at a timepoint, so replicates are exchangeable, timepoints are few and fixed, and residual autocorrelation across ordered timepoints is real. This is categorically different from single-cell pseudotime, which is a latent per-cell ordering estimated with error and belongs to single-cell/trajectory-inference (and to tradeSeq's native use case). Conflating the two is the deepest error in the area.

Three decisions dominate correctness before any p-value is read:

1. Distribution. A Gaussian GAM on raw counts gives wrong SEs, wrong p-values, and can predict negatives. Model counts with nb() and a library-size offset, or fit Gaussian on a variance-stabilized scale.
2. Autocorrelation. Positive residual correlation across timepoints shrinks the effective sample size, so a plain gam() under-estimates SEs and over-calls temporal trends. Name the assumption and model it (corAR1/bam(rho=)) or at least inspect the residual ACF.
3. Mechanism. A smooth GAM assumes a gradually curving process; a changepoint model assumes a genuinely abrupt regime shift. Imposing changepoints on smooth data invents regime shifts; smoothing a true step Gibbs-rings over it. Match the model to the biology.

Smooth GAM vs changepoint: choosing the model class

| Question | Model | Use when | Do NOT use when | |----------|-------|----------|-----------------| | Does expression change / curve over time? | GAM s(time) (mgcv) | process is gradual (induction/decay kinetics, developmental ramps) | the process is a discrete switch -> a smooth smears the discontinuity | | Do two conditions' trajectories diverge? | ordered-factor difference smooth (mgcv) | testing whether treated shape departs from control | groups have no shared reference / only a constant offset differs (use a parametric term) | | When does the regime shift (slope break)? | segmented (broken-line) | continuous piecewise-LINEAR change in slope | the shift is a level jump, not a slope change (use ruptures l2) | | When does the regime shift (level/distribution)? | ruptures Pelt/Binseg | step change in mean (l2) or distribution (rbf) | the curve is smooth -> any liberal penalty fabricates breaks | | Is a curve even warranted? | AIC/edf of s(time) vs linear | deciding non-linear vs linear is enough | over-interpreting edf near 1 as a real curve |

Methodology evolves; before committing verify current defaults and recommendations against the latest mgcv/ruptures/segmented documentation.

mgcv GAM (R)

Goal: Fit a smooth non-linear curve to expression over measured time and test whether it changes, on the correct distributional scale.

Approach: Use a penalized regression spline s(time); set the basis-dimension ceiling k generously but below the number of unique timepoints and let REML choose the realized wiggliness; use nb() + a library-size offset for raw counts, or Gaussian on a variance-stabilized scale.

library(mgcv)

# k is a CEILING (max basis dimension), NOT the number of knots the curve will use.
# The REML-chosen penalty picks the realized wiggliness; edf (below) reports it.
# k must be < number of unique timepoints; realistically k <= (#timepoints - 1).
# method='REML': better-behaved objective than GCV, resists under/over-smoothing (Wood 2011).
fit <- gam(expression ~ s(time, k = 6, bs = 'tp'), data = gene_df, method = 'REML')

summary(fit)
# s.table columns: edf, Ref.df, F (Gaussian/unknown scale), p-value.
# edf ~ 1 => penalty shrank the smooth to linear; edf near k-1 => nearly full flexibility.
# The smooth-term p-value is APPROXIMATE (conditional on estimated lambda; Wood 2013):
# treat it as categorical significant/not, do not compare tiny magnitudes.

Raw counts: NB family with a library-size offset

Goal: Model overdispersed RNA-seq counts on the count scale without violating the Gaussian assumption.

Approach: Fit family=nb() (mgcv estimates theta by REML) with offset(log(library_size)) so the smooth describes rate, not depth.

# Counts are mean-variance coupled and overdispersed: Gaussian is wrong on raw counts.
# nb() estimates theta jointly with the smoothing parameters under REML.
# offset(log(libsize)) absorbs sequencing depth so s(time) models expression rate.
fit_nb <- gam(counts ~ s(time, k = 6) + offset(log(library_size)),
              data = gene_df, family = nb(), method = 'REML')
# With a known-scale family the test column becomes Chi.sq, not F.
summary(fit_nb)$s.table

Residual autocorrelation across timepoints

Goal: Prevent inflated smooth-term significance caused by correlation between successive timepoints.

Approach: Model a lag-1 AR structure grouped by the replication unit with gamm(correlation=corAR1()), or fix rho in bam() for genome-wide fits after reading the lag-1 residual ACF.

# Plain gam() assumes independent residuals; positive AR(1) shrinks effective n,
# under-estimates SEs, and lets the smooth get too wiggly -> false temporal trends.
# corAR1 models e_t = rho * e_{t-1} + noise; group by subject/animal/plate.
fit_ar <- gamm(expression ~ s(time, k = 6),
               correlation = corAR1(form = ~ time | subject),
               data = gene_df, method = 'REML')
# fit_ar$gam holds the smooth; fit_ar$lme holds the correlation estimate.

# Genome-wide alternative: fit once without AR, read lag-1 residual ACF, set rho, refit.
# AR.start flags the first observation of each independent series.
fit_bam <- bam(expression ~ s(time, k = 6), data = gene_df,
               rho = 0.4, AR.start = series_start, method = 'fREML')

With independent biological replicates AT EACH timepoint the correlation is often weak or unidentifiable and plain gam() is defensible; a single series sampled repeatedly over many timepoints is where AR bites hardest. Always inspect the residual ACF before trusting the smooth p-value.

Comparing conditions with an ordered-factor difference smooth

Goal: Directly test whether the treated trajectory's shape diverges from control, with its own p-value.

Approach: Make the grouping an ORDERED factor so s(time, by=grp) becomes a difference smooth (level minus reference); keep the reference global smooth AND the parametric main effect.

# Unordered by= gives each group's curve vs zero -- NOT a divergence test.
# Ordered factor: s(time, by=grp) is the difference smooth; its single p-value
# directly tests whether the trajectories diverge. Extends cleanly to >2 groups.
gene_df$condition <- as.ordered(gene_df$condition)
fit_diff <- gam(expression ~ condition + s(time, k = 6) + s(time, k = 6, by = condition),
                data = gene_df, method = 'REML')
# The parametric 'condition' term is REQUIRED: centered smooths cannot carry the
# group's overall level, so without it a constant offset is misattributed to the smooth.
summary(fit_diff)

A numeric 0/1 by=is_treated indicator is a valid shortcut for a single 2-level contrast (the second smooth is the treatment deviation), but the ordered-factor form is the general, canonical idiom.

Diagnostics: gam.check, k-index, concurvity

Goal: Decide whether the basis is adequate and whether smooth terms are mutually identifiable.

Approach: Read gam.check()/k.check(); respond to a low k-index by doubling k and refitting, not by reflexively cranking k; use concurvity() only for multi-smooth models.

gam.check(fit)   # 4 residual plots + k.check(): reports k', edf, k-index, p-value

# k-index < 1 with a small p means residual pattern the basis is too rigid to capture.
# CORRECT diagnostic = double k and refit: if edf rises substantially, k was too low;
# if edf barely moves, k was fine and the low k-index reflects autocorrelation or
# a distributional problem -- do not just raise k.

# Concurvity = the smooth analog of collinearity; only meaningful for multi-term models.
# Near 1 => partial attribution between smooths is unstable; > 0.8 is a worry, not a cutoff.
concurvity(fit_diff, full = TRUE)

Prediction and pointwise intervals

Goal: Visualize the fitted trajectory with an uncertainty band, within the sampled range only.

Approach: Predict on a fine grid with se.fit=TRUE; band = fit +/- 1.96*SE (pointwise, not simultaneous); never extrapolate.

grid <- data.frame(time = seq(min(gene_df$time), max(gene_df$time), length.out = 200))
pred <- predict(fit, newdata = grid, se.fit = TRUE)
grid$fitted <- pred$fit
# 1.96*SE is a POINTWISE 95% band; whole-curve (simultaneous) coverage is < 95%,
# so overlapping condition bands are NOT a formal test -- use the difference smooth p-value.
grid$lower <- pred$fit - 1.96 * pred$se.fit
grid$upper <- pred$fit + 1.96 * pred$se.fit
# Beyond [min(time), max(time)] the spline and its SE diverge -- do not predict outside range.

Genome-wide GAM fitting + FDR

Goal: Rank genes by temporal significance across the transcriptome.

Approach: Fit s(time) per gene, collect the smooth p-value, apply BH across genes (the per-gene p-values are approximate, so the FDR is approximate; permutation calibration is the gold standard for strong claims).

# expr_mat here must be variance-stabilized (vst / log-CPM); for raw counts use the family=nb() + offset
# form above instead of this default-Gaussian fit, or the per-gene SEs and p-values are invalid.
results <- data.frame()
for (gene in rownames(expr_mat)) {
    df <- data.frame(expression = as.numeric(expr_mat[gene, ]), time = timepoints)
    fit <- gam(expression ~ s(time, k = 6), data = df, method = 'REML')
    s_tab <- summary(fit)$s.table
    results <- rbind(results, data.frame(gene = gene, edf = s_tab[, 'edf'],
                                         p_value = s_tab[, 'p-value']))
}
results$q_value <- p.adjust(results$p_value, method = 'BH')  # q<0.05: standard FDR floor
temporal_genes <- results[results$q_value < 0.05, ]

tradeSeq (R/Bioconductor) -- off-label for bulk

tradeSeq is BUILT for single-cell pseudotime lineages, not bulk real-time. fitGAM(counts, pseudotime, cellWeights, nknots) expects a gene x cell count matrix, a cell x lineage pseudotime matrix, and cell x lineage soft-assignment weights, and fits an NB GAM per gene per lineage. Its tests (associationTest, startVsEndTest, conditionTest, patternTest) are keyed to pseudotime lineages.

# Only reasonable when the design is genuinely a pseudobulk-over-a-lineage.
# For a standard bulk time-course with replicates at fixed timepoints, use mgcv directly:
# the same NB GAM, with full control over by= contrasts, offsets, and AR structure,
# and without single-cell scaffolding. nknots plays the same ceiling role as k (choose it with evaluateK()).
library(tradeSeq)
sce <- fitGAM(counts = count_mat, pseudotime = pt_mat, cellWeights = cw_mat, nknots = 6)
assoc_res <- associationTest(sce)   # matrix of Wald stat + df + p-value, NOT an SCE

segmented (R) -- broken-line slope break

Goal: Locate a continuous change in SLOPE and test whether a break exists at all.

Approach: Pre-test with davies.test before fitting a break; estimate the breakpoint with segmented() from a starting value; limit to one break unless the data are dense.

library(segmented)
lm_fit <- lm(expression ~ time, data = gene_df)

# davies.test H0: difference-in-slopes = 0 (no breakpoint). It searches candidate
# locations and corrects the minimum p for the search (slightly conservative).
# Decide WHETHER a break exists before interpreting one. pscore.test() is more powerful
# for a single break.
davies.test(lm_fit, seg.Z = ~time)

# segmented models a CONTINUOUS slope change (not a level jump); psi = starting value(s),
# NA auto-initializes. The estimator is a local search: sensitive to psi, fragile with
# multiple breaks (closely spaced breaks are non-identifiable).
seg_fit <- segmented(lm_fit, seg.Z = ~time, psi = NA)
summary(seg_fit)$psi   # breakpoint estimate + SE (break CI is approximate/often too narrow)

ruptures (Python) -- level/distribution changepoints

Goal: Detect discrete times where the mean (or whole distribution) shifts.

Approach: Factorize as (search method) x (cost model) x (penalty); the penalty choice IS the number-of-changepoints choice; match the cost model to the shift type and estimate the noise variance, not the total variance.

import numpy as np
import ruptures as rpt

signal = np.asarray(expression_values)

# model='l2' detects changes in MEAN (piecewise-constant level -- the natural choice for
#   step-like expression regimes). 'rbf' detects changes in the whole distribution
#   (mean AND variance) -- more general, hungrier for data, less interpretable.
# min_size=2: minimum segment length. Pelt is exact (O(n) via pruning); it takes a penalty
#   and RETURNS the number+location of breaks -- there is no separate 'how many' knob.
n = len(signal)
# Noise variance, NOT total variance: np.var(signal) includes between-regime variation, so
# a BIC-style penalty log(n)*np.var(signal) is too large and UNDER-detects. Estimate noise
# from lag-1 differences instead. BIC is derived for the l2/Gaussian-mean cost -- pairing it
# with model='rbf' is theoretically mismatched; use l2 with BIC, or calibrate rbf empirically.
sigma2 = np.var(np.diff(signal)) / 2.0
penalty = np.log(n) * sigma2
bkps = rpt.Pelt(model='l2', min_size=2).fit(signal).predict(pen=penalty)
# predict returns break indices INCLUDING the terminal index n:
n_changepoints = len(bkps) - 1   # bkps[:-1] are the actual break locations

# Binseg: greedy, approximate, fast; takes a KNOWN number of breaks (sanity check vs Pelt).
bkps_binseg = rpt.Binseg(model='l2', min_size=2).fit(signal).predict(n_bkps=2)

Guard against fabricated breaks: a liberal penalty always "finds" changepoints in a smooth ramp. Require a pre-test (a break exists) or compare a piecewise fit against a smooth-GAM fit by AIC -- if the smooth wins, the "changepoint" is a sampling-noise artifact. With few timepoints, be extremely skeptical of more than one break.

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | p-values wrong / fitted curve predicts negative expression | Gaussian GAM on raw overdispersed counts | family=nb() + offset(log(library_size)), or fit Gaussian on vst/log-CPM | | Many genes "significantly change over time" implausibly | residual autocorrelation inflates smooth-term significance | gamm(..., correlation=corAR1(form=~time|subject)) or bam(..., rho=, AR.start=); check residual ACF | | Treating k as "the number of bends I want" | k is the flexibility CEILING, not realized complexity | set k generously (< #timepoints), let REML pick lambda; read edf, not k | | Cranking k whenever k-index < 1 | low k-index can mean autocorrelation/heteroscedasticity, not low basis | double k and refit -- if edf jumps, raise k; if not, look at correlation/distribution | | Reading p=1e-30 as thirty orders of certainty | smooth p-values are approximate (ignore full lambda uncertainty) | treat as categorical significant/not; apply BH FDR across genes | | Unordered by= used "to test if curves differ" | it gives each group vs zero, not a divergence test | as.ordered(condition) -> the difference smooth's p-value IS the divergence test | | by= smooth without the parametric main effect | centered smooths cannot carry the group level | include condition + alongside s(time, by=condition) | | tradeSeq::fitGAM on a plain bulk time-course | tradeSeq is single-cell pseudotime machinery (needs cellWeights) | use mgcv directly for bulk real-time; reserve tradeSeq for pseudobulk lineages | | ruptures under-detects real changepoints | pen=log(n)*np.var(signal) uses TOTAL variance -> penalty too large | estimate noise from np.var(np.diff(signal))/2; sweep the penalty | | model='rbf' with a BIC (log n * var) penalty | BIC penalty is derived for the l2/Gaussian-mean cost | use model='l2' with BIC, or calibrate the rbf penalty empirically | | Changepoints "found" in a clearly smooth ramp | a liberal penalty fabricates breaks in gradual data | pre-test with davies.test; compare piecewise vs smooth-GAM AIC; the smooth often wins | | Fitted curve behaves wildly past the last timepoint | extrapolating a penalized spline beyond the data | predict only within [min(time), max(time)] | | "The condition bands overlap, so no difference" | 1.96*SE bands are pointwise, not simultaneous | use the difference-smooth p-value or simultaneous (posterior-simulation) intervals | | segmented with several psi gives unstable breaks | multiple breakpoints are weakly identifiable with few noisy points | limit to 1 break unless data are dense; supply good psi starts; check convergence |

Related Skills

• temporal-clustering - group genes by trajectory shape after fitting
• circadian-rhythms - periodic (known-period) trajectory models rather than smooth trends
• periodicity-detection - discover unknown-period oscillation instead of a smooth trend
• differential-expression/timeseries-de - linear/spline model alternatives for temporal DE
• single-cell/trajectory-inference - single-cell pseudotime (latent inferred ordering), the case tradeSeq is built for

References

• Wood SN. 2011. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J R Stat Soc B 73(1):3-36. doi:10.1111/j.1467-9868.2010.00749.x. (REML smoothing-parameter selection, better-behaved than GCV.)
• Wood SN. 2013. On p-values for smooth components of an extended generalized additive model. Biometrika 100(1):221-228. doi:10.1093/biomet/ass048. (Smooth-term p-values are approximate; pointwise interval coverage.)
• Wood SN. 2017. Generalized Additive Models: An Introduction with R, 2nd ed. Chapman & Hall/CRC. ISBN 9781498728331. (Basis-penalty framework, gam.check, concurvity.)
• Pedersen EJ, Miller DL, Simpson GL, Ross N. 2019. Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7:e6876. doi:10.7717/peerj.6876. (Global-plus-difference-smooth and factor-smooth condition-comparison structure.)
• Van den Berge K, Roux de Bezieux H, Street K, Saelens W, Cannoodt R, Saeys Y, Dudoit S, Clement L. 2020. Trajectory-based differential expression analysis for single-cell sequencing data. Nat Commun 11(1):1201. doi:10.1038/s41467-020-14766-3. (tradeSeq: NB-GAM DE along pseudotime lineages, hence off-label for bulk.)
• Muggeo VMR. 2003. Estimating regression models with unknown break-points. Stat Med 22(19):3055-3071. doi:10.1002/sim.1545. (Broken-line estimator behind segmented and davies.test.)
• Killick R, Fearnhead P, Eckley IA. 2012. Optimal detection of changepoints with a linear computational cost. J Am Stat Assoc 107(500):1590-1598. doi:10.1080/01621459.2012.737745. (PELT exact O(n) penalized algorithm behind rpt.Pelt.)
• Truong C, Oudre L, Vayatis N. 2020. Selective review of offline change point detection methods. Signal Processing 167:107299. doi:10.1016/j.sigpro.2019.107299. (The cost x search x constraint taxonomy; the ruptures reference paper.)