Bio 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

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)