Bio Differential Expression Timeseries De

Version Compatibility

Reference examples tested with: DESeq2 1.42+, edgeR 4.0+, limma 3.58+, splines (base R), maSigPro 1.74+, ImpulseDE2 1.10+ (Bioconductor archive; verify availability), variancePartition / dream 1.32+, Mfuzz 2.62+, TCseq 1.26+, ggplot2 3.5+, pheatmap 1.0+

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

• 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.

Time-Series Differential Expression

"Find genes that change over time" -> Define what "change" means -- any non-zero time effect (LRT), a smooth nonlinear trend (splines), a transient impulse (ImpulseDE2), or differing trajectories between groups (interaction LRT) -- and choose the model that asks that specific question while handling repeated-measures correctly.

The Single Most Important Modern Insight -- Most dedicated time-course tools UNDERPERFORM pairwise + LRT on short series

Spies, Renz, Beyer, Ciaudo 2019 Brief Bioinform 20:288 benchmarked dedicated time-course tools (ImpulseDE2, splineTC, maSigPro, EBSeqHMM, TimeReg) against naive DESeq2/edgeR pairwise comparisons + LRT for omnibus, on simulated and real time-courses. Finding: on short series (<8 time points), naive pairwise pattern-of-significance OUTPERFORMS dedicated TC tools because of high false-positive rates in the latter. The exception is ImpulseDE2, which holds up better than the others -- IF its impulse assumption (rise-then-plateau or fall-then-plateau) actually fits the biology.

For most experimental time courses (3-6 time points, common in pharmacology and developmental biology), the right tool is DESeq2 with test='LRT' and a sensible reduced model. Reserve splines for >5 evenly-spaced time points; reserve ImpulseDE2 for monotonic-then-saturating dynamics; reserve DREAM for repeated measures.

A second insight that is constantly violated: pseudoreplication. If 3 subjects each contribute 4 time points (12 samples), the effective sample size for testing TIME effects is closer to 3, not 12 -- the within-subject observations are not independent. Treating them as independent inflates type-I error dramatically. Either include subject as a fixed effect, use DREAM (mixed model), or collapse to per-subject means (loses time info).

Algorithmic Taxonomy

Decision Tree by Scenario

DESeq2 LRT for Time -- The Canonical Pattern

Goal: Test whether time has ANY effect (omnibus), or whether time trajectories differ between groups (interaction).

Approach: Specify a full design including the time terms; specify a reduced design dropping those terms; LRT compares.

library(DESeq2)

dds <- DESeqDataSetFromMatrix(counts, colData, design = ~ time)
dds <- DESeq(dds, test = 'LRT', reduced = ~ 1)
res <- results(dds)

The LRT p-value tests "any difference among time points". The log2FoldChange column reports the LAST coefficient in resultsNames(dds) -- NOT the omnibus effect. Use padj only from LRT results; extract individual Wald per time point for effect sizes.

Interaction (treatment-specific time response):

dds <- DESeqDataSetFromMatrix(counts, colData,
    design = ~ treatment + time + treatment:time)
dds <- DESeq(dds, test = 'LRT', reduced = ~ treatment + time)
res_interaction <- results(dds)

This tests "does the time trajectory differ between treatment groups?" -- a different question than "is there a time effect" or "is there a treatment effect".

Time as Factor vs Continuous vs Spline

Rule of thumb: df <= unique_time_points / 2.

library(splines)

dds <- DESeqDataSetFromMatrix(counts, colData,
    design = ~ treatment * ns(time, df = 3))
dds <- DESeq(dds, test = 'LRT', reduced = ~ treatment + ns(time, df = 3))

limma-voom + Splines (Modern Alternative)

library(limma)
library(edgeR)
library(splines)

y <- DGEList(counts = counts)
y <- normLibSizes(y)
keep <- filterByExpr(y, group = metadata$treatment)
y <- y[keep, , keep.lib.sizes = FALSE]

design <- model.matrix(~ treatment * ns(time, df = 3), data = metadata)
v <- voom(y, design, plot = TRUE)
fit <- lmFit(v, design)
fit <- eBayes(fit, robust = TRUE)

interaction_cols <- grep(':ns\\(time', colnames(design))
tt <- topTable(fit, coef = interaction_cols, number = Inf)

Tests the joint significance of all interaction spline coefficients.

DREAM for Repeated Measures

Goal: Properly model longitudinal data where the same subject is sampled multiple times.

Approach: Linear mixed model per gene with subject as a random intercept (or slope), via variancePartition::dream. Uses voom-weighted linear regression internally. Default ddf='adaptive' uses Kenward-Roger for n <= 20 samples (the typical longitudinal-DE regime) and Satterthwaite otherwise; force either explicitly with ddf='Kenward-Roger' or ddf='Satterthwaite'.

library(variancePartition)
library(edgeR)
library(BiocParallel)

y <- DGEList(counts = counts)
y <- normLibSizes(y)
keep <- filterByExpr(y, group = metadata$treatment)
y <- y[keep, , keep.lib.sizes = FALSE]

formula <- ~ treatment + time + treatment:time + (1 | subject)

vobj <- voomWithDreamWeights(y, formula, metadata)
fitmm <- dream(vobj, formula, metadata)
fitmm <- eBayes(fitmm)
tt <- topTable(fitmm, coef = grep(':time', colnames(coefficients(fitmm))),
               number = Inf)

When to use DREAM over ~ subject + treatment + time (subject as fixed effect):

• More than 2 time points per subject (more random-effect signal to estimate)
• Many subjects (random effects more parsimonious than fixed)
• Random slopes needed (e.g., (1 + time | subject))
• Multiple random effects (e.g., (1 | donor) + (1 | batch))

duplicateCorrelation (limma) is the older approximate alternative -- assumes a SINGLE within-subject correlation across all genes. Adequate for technical replicates within biological replicates; less so for proper longitudinal data with multiple time points.

maSigPro (Nueda 2014 for RNA-seq)

library(maSigPro)

edesign <- data.frame(
    Time      = metadata$time,
    Replicate = metadata$replicate,
    Control   = as.numeric(metadata$treatment == 'Control'),
    Treatment = as.numeric(metadata$treatment == 'Treatment')
)
rownames(edesign) <- metadata$sample

design <- make.design.matrix(edesign, degree = 3)

fit <- p.vector(norm_counts, design, Q = 0.05, MT.adjust = 'BH')
tstep <- T.fit(fit, step.method = 'backward', alfa = 0.05)
sigs <- get.siggenes(tstep, rsq = 0.6, vars = 'groups')
see.genes(sigs$sig.genes, show.fit = TRUE, dis = design$dis,
          cluster.method = 'hclust', k = 9)

CITATION: Nueda MJ, Tarazona S, Conesa A (2014) "Next maSigPro: updating maSigPro bioconductor package for RNA-seq time series." Bioinformatics 30(18):2598-2602. The earlier Conesa et al. 2006 Bioinformatics 22:1096 is the microarray-era original. Cite 2014 for RNA-seq use; many secondary refs miscite 2006.

ImpulseDE2

Goal: Detect transient impulse-like expression patterns (rise then decay, or constant-then-monotonic).

Approach: Fits constant, monotonic, and impulse (6-parameter sigmoid: baseline, peak, post-peak baseline, rise rate, decay rate) models per gene; selects the best-fitting and tests for differential dynamics.

library(ImpulseDE2)

dfAnnotation <- data.frame(
    Sample    = colnames(counts),
    Time      = metadata$time,
    Condition = metadata$condition,
    Batch     = metadata$batch
)

imp <- runImpulseDE2(
    matCountData    = as.matrix(counts),
    dfAnnotation    = dfAnnotation,
    boolCaseCtrl    = TRUE,
    vecConfounders  = c('Batch'),
    scaNProc        = 4
)

sig <- imp$dfImpulseDE2Results[imp$dfImpulseDE2Results$padj < 0.05, ]

The impulse model FAILS on:

• Oscillatory expression (circadian, cell cycle) -- model can't represent periodicity
• Monotonic-but-non-asymptotic responses (linear increases that don't plateau)
• Multi-phase responses (rise-fall-rise)

For oscillatory biology, use Fourier basis or dedicated tools (JTK_CYCLE, MetaCycle). For non-asymptotic monotonic, splines fit better.

NOTE: ImpulseDE2 was removed from Bioconductor at the 3.13 release (May 2021); last hosted version was 3.10 -- install from the BiocArchive (pin to Bioc 3.10) or the YosefLab GitHub mirror. The Spies 2019 benchmark showed ImpulseDE2 holds up on longer time courses (~8+ time points) but is sensitive to noise on short series; below ~8 time points, DESeq2 pairwise + LRT typically outperforms.

Trajectory Clustering of DE Genes

Goal: After identifying DE-over-time genes, group them by trajectory shape.

Approach: Standardize per gene (subtract mean, divide by SD per gene) -- otherwise clusters reflect mean level, not shape. Then apply DPGP (nonparametric), Mfuzz (fuzzy c-means), or splines + k-means.

library(Mfuzz)

eset <- ExpressionSet(assayData = as.matrix(norm_counts[sig_genes, ]))
eset_std <- standardise(eset)

m <- mestimate(eset_std)
cl <- mfuzz(eset_std, c = 9, m = m)

mfuzz.plot(eset_std, cl, mfrow = c(3, 3))

The Mfuzz fuzzifier m is critical -- too low gives crisp clusters (loses fuzzy advantage); too high collapses everything. mestimate() implements Schwämmle & Jensen 2010 Bioinformatics 26:2841 to estimate m from the data.

CITATION CARE: the Mfuzz PACKAGE paper is Kumar L, Futschik ME (2007) Bioinformation 2(1):5-7 (the journal is Bioinformation, NOT Bioinformatics). The Schwämmle 2010 paper is the fuzzifier-estimation methodology paper, Bioinformatics. Many references confuse the two.

For DPGP (Dirichlet Process Gaussian Process; nonparametric in cluster number AND trajectory shape):

CITATION: McDowell IC, Manandhar D, Vockley CM, Schmid AK, Reddy TE, Engelhardt BE (2018) "Clustering gene expression time series data using an infinite Gaussian process mixture model." PLoS Comput Biol 14(1):e1005896. CITATION CARE: this is PLoS Comp Biol, NOT Genome Research (a common miscitation).

Splines + k-means (fast alternative):

library(splines)
spline_coefs <- t(apply(norm_counts[sig_genes, ], 1, function(x) {
    fit <- lm(x ~ ns(metadata$time, df = 4))
    coef(fit)
}))
km <- kmeans(scale(spline_coefs), centers = 6)