Bio Workflows Timecourse Pipeline

Expression matrix + time metadata
    |
    v
[1. Temporal DE] ---------> limma splines / DESeq2 LRT / statsmodels spline F-test
    |                            (selects temporally variable genes; FDR <0.05)
    v
[2. Filter] --------------> Significant temporal genes only (clustering input)
    |
    v
[3. Soft Clustering] -----> Mfuzz (R) / tslearn TimeSeriesKMeans (Python) on z-scored profiles
    |                            +---> QC: membership >0.5, no empty clusters, sweep k
    |
    v
[4b. GAM Trajectory] -----> mgcv / pygam GAM on standardized cluster-mean profiles
    |
    v
[5. Pathway Enrichment] --> clusterProfiler / gseapy per cluster
    |                            background = temporal genes (NOT the genome)
    v
Temporal gene modules + enriched pathways + trajectory fits

    OPTIONAL branch, gated (NOT on the default path):
    IF design covers >=2 full cycles AND >=6-8 samples/cycle (even spacing)
       AND collection order was randomized:
        [4a. Rhythm Detection] --> MetaCycle meta2d / CosinorPy fit_group
                                   period/phase/amplitude + BH q across genes
    ELSE: SKIP (design does not license a rhythm test)

library(limma)
library(splines)

expr <- as.matrix(read.csv('counts_normalized.csv', row.names = 1))
meta <- read.csv('metadata.csv')

# Natural cubic spline on time; df=3 is enough for most courses, raise to 4-5 for >10 timepoints
design <- model.matrix(~ ns(meta$time, df = 3))
fit <- lmFit(expr, design)
fit <- eBayes(fit)

# Joint F-test on all spline coefficients = "expression changes over time"
temporal_results <- topTable(fit, coef = 2:ncol(design), number = Inf, sort.by = 'F')
# topTable already returns adj.P.Val (BH-corrected); use it directly (not $FDR, which does not exist)

library(DESeq2)

counts <- as.matrix(read.csv('raw_counts.csv', row.names = 1))
meta <- read.csv('metadata.csv')
meta$time <- factor(meta$time)

# LRT: full model (time as factor) vs reduced (intercept) = any between-timepoint change
dds <- DESeqDataSetFromMatrix(counts, colData = meta, design = ~ time)
dds <- DESeq(dds, test = 'LRT', reduced = ~ 1)
res <- results(dds)   # BH-adjusted padj

import pandas as pd
import numpy as np
from statsmodels.stats.multitest import multipletests
from patsy import dmatrix
from scipy import stats

expr = pd.read_csv('counts_normalized.csv', index_col=0)
meta = pd.read_csv('metadata.csv')

spline_basis = dmatrix('bs(time, df=3)', data=meta, return_type='dataframe')
# patsy already includes an Intercept column; do NOT prepend another np.ones (that duplicates the
# intercept, inflates model df, and miscalibrates the F-test -> wrong temporal-DE FDR). Use it directly.
design_full = spline_basis.values                          # [Intercept, bs1, bs2, bs3] = 4 cols
design_reduced = np.ones((len(meta), 1))
df_diff = design_full.shape[1] - design_reduced.shape[1]   # 4 - 1 = 3
df_resid = len(meta) - design_full.shape[1]                # n - 4

pvals = []
for gene in expr.index:
    y = expr.loc[gene].values
    ss_full = np.sum((y - design_full @ np.linalg.lstsq(design_full, y, rcond=None)[0]) ** 2)
    ss_red = np.sum((y - design_reduced @ np.linalg.lstsq(design_reduced, y, rcond=None)[0]) ** 2)
    f_stat = ((ss_red - ss_full) / df_diff) / (ss_full / df_resid)
    pvals.append(1 - stats.f.cdf(f_stat, df_diff, df_resid))

# multipletests default is Holm-Sidak; force BH explicitly
_, fdr, _, _ = multipletests(pvals, method='fdr_bh')
temporal_genes = expr.index[fdr < 0.05].tolist()

sig_genes <- temporal_results[temporal_results$adj.P.Val < 0.05, ]
n_sig <- nrow(sig_genes)
message(sprintf('Significant temporal genes: %d', n_sig))
# <100: underpowered clustering; >10000: check batch/normalization before proceeding
if (n_sig < 100) message('WARNING: Few temporal genes. Check timepoint spacing or relax FDR.')
if (n_sig > 10000) message('WARNING: Many temporal genes. Inspect batch effects / normalization.')

# FDR <0.05 standard; 0.1 acceptable for exploratory clustering only
sig_genes <- rownames(temporal_results[temporal_results$adj.P.Val < 0.05, ])
expr_sig <- expr[sig_genes, ]
message(sprintf('Genes passing FDR <0.05: %d', length(sig_genes)))

library(Mfuzz)

eset <- ExpressionSet(assayData = as.matrix(expr_sig))
eset <- standardise(eset)   # per-gene mean 0, sd 1: makes distance shape-based, not magnitude-based

# mestimate() implements Schwaemmle & Jensen (2010): the smallest m that stops fuzzy c-means
# from clustering RANDOMIZED data. It is dominated by the number of timepoints; inspect the
# returned value and the membership distribution rather than trusting it blindly (or hardcoding m=2).
m <- mestimate(eset)
message(sprintf('Estimated fuzzifier m = %.2f', m))

# k is a resolution CHOICE, not a result: start ~sqrt(n_genes/2), then sweep and validate (below)
n_clusters <- 8
cl <- mfuzz(eset, c = n_clusters, m = m)

# Membership >0.5 = core (confident) genes; lower to 0.3 only for exploratory overlap
core_genes <- acore(eset, cl, min.acore = 0.5)

from tslearn.clustering import TimeSeriesKMeans

# Row-wise z-score: normalize each gene across its own timepoints (shape, not level)
expr_scaled = (expr_sig.values - expr_sig.values.mean(axis=1, keepdims=True)) / expr_sig.values.std(axis=1, keepdims=True)

# soft-DTW tolerates phase-shifted profiles Euclidean would split; gamma smooths the DTW geometry.
# Use plain 'euclidean' when absolute phase is biologically meaningful (morning vs evening genes).
model = TimeSeriesKMeans(n_clusters=8, metric='softdtw', metric_params={'gamma': 0.01},
                         max_iter=50, random_state=42)
labels = model.fit_predict(expr_scaled.reshape(expr_scaled.shape[0], expr_scaled.shape[1], 1))

library(cluster)
cluster_sizes <- table(cl$cluster)
print(cluster_sizes)
if (any(cluster_sizes == 0)) message('WARNING: Empty clusters. Reduce n_clusters.')

for (i in seq_along(core_genes)) {
    message(sprintf('Cluster %d: %d core genes (membership >0.5)', i, nrow(core_genes[[i]])))
}

# Silhouette on the z-scored profiles; triangulate k with a sweep, do not crown one index
sil <- silhouette(cl$cluster, dist(exprs(eset)))
message(sprintf('Mean silhouette: %.3f', mean(sil[, 3])))

CIRCADIAN_DESIGN = False  # the un-computable precondition (circadian design + randomized order); set True only if it holds

n_cycles = (meta['time'].max() - meta['time'].min()) / 24.0
samples_per_cycle = meta['time'].nunique() / max(n_cycles, 1e-9)
gate_design = n_cycles >= 2 and samples_per_cycle >= 6   # the COMPUTABLE part of the gate
if CIRCADIAN_DESIGN and gate_design:
    pass  # run rhythm detection (CosinorPy/MetaCycle below)
elif not gate_design:
    print(f'Rhythm detection SKIPPED: inadequate design ({n_cycles:.1f} cycles, {samples_per_cycle:.1f}/cycle; need >=2 and >=6-8).')
else:
    print('Rhythm detection SKIPPED: design meets the cycle/sampling floor but CIRCADIAN_DESIGN is not set (randomized-order / circadian precondition unconfirmed).')

library(MetaCycle)
expr_for_meta <- expr_sig
colnames(expr_for_meta) <- meta$time
write.csv(expr_for_meta, 'expr_for_metacycle.csv')

# Circadian search window 20-28h; ARS/JTK require EVEN integer sampling and drop out silently
# (analysisStrategy='auto') on uneven/replicated data, leaving LS only.
meta2d('expr_for_metacycle.csv', filestyle = 'csv',
       minper = 20, maxper = 28,
       timepoints = sort(unique(meta$time)),
       outdir = 'metacycle_results')
# Filter on meta2d_BH.Q (BH FDR) AND meta2d_rAMP (relative amplitude); significance alone over-detects.

from CosinorPy import cosinor   # note: import name is capitalized CosinorPy, not cosinorpy

# fit_group expects long-format columns 'x' (time), 'y' (expression), 'test' (gene id).
# period=24 for circadian; n_components=2 adds the 12h harmonic for non-sinusoidal shapes.
results = cosinor.fit_group(expr_long, period=24, n_components=1)

# fit_group ALREADY returns a BH-adjusted 'q' column across the fitted group; use it, do NOT
# threshold the raw per-gene 'p'. Add a RELATIVE-amplitude filter (fit_group has no rAMP column,
# so compute rAMP = amplitude/mesor): significance alone over-detects rhythms. rAMP>0.1 = >=10% of baseline.
results['rAMP'] = results['amplitude'] / results['mesor']
rhythmic = results[(results['q'] < 0.05) & (results['rAMP'] > 0.1)]

library(mgcv)

cluster_trajectories <- list()
for (cl_id in 1:n_clusters) {
    cl_genes <- names(cl$cluster[cl$cluster == cl_id])
    mean_profile <- colMeans(expr_sig[cl_genes, ])
    df_gam <- data.frame(time = meta$time, expr = mean_profile)

    # k is a flexibility CEILING (max basis dimension), NOT the number of knots/bends to fit.
    # REML (not GCV) then picks the wiggliness penalty; realized complexity is reported as edf.
    # Keep k < number of unique timepoints (identifiability); k=5 suits >=6-8 timepoints.
    gam_fit <- gam(expr ~ s(time, k = 5), data = df_gam, method = 'REML')

    cluster_trajectories[[cl_id]] <- list(fit = gam_fit,
                                          r_squared = summary(gam_fit)$r.sq,
                                          edf = summary(gam_fit)$edf)
    message(sprintf('Cluster %d: R^2 = %.3f, EDF = %.2f (edf~1 => linear; edf~k-1 => highly non-linear)',
                    cl_id, summary(gam_fit)$r.sq, summary(gam_fit)$edf))
}

from pygam import LinearGAM, s

for cl_id in range(n_clusters):
    mean_profile = expr_scaled[labels == cl_id].mean(axis=0)
    # n_splines is the basis-dimension ceiling (like mgcv k); the penalty picks realized wiggliness.
    gam = LinearGAM(s(0, n_splines=5)).fit(meta['time'].values.reshape(-1, 1), mean_profile)
    print(f'Cluster {cl_id}: GCV = {gam.statistics_["GCV"]:.4f}, edof = {gam.statistics_["edof"]:.2f}')

library(clusterProfiler)
library(org.Hs.eg.db)

all_temporal_entrez <- bitr(rownames(expr_sig), fromType = 'SYMBOL', toType = 'ENTREZID',
                            OrgDb = org.Hs.eg.db)

enrichment_results <- list()
for (i in seq_along(core_genes)) {
    entrez <- bitr(core_genes[[i]]$NAME, fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = org.Hs.eg.db)
    ego <- enrichGO(gene = entrez$ENTREZID,
                    universe = all_temporal_entrez$ENTREZID,   # background = temporal genes
                    OrgDb = org.Hs.eg.db, ont = 'BP', pAdjustMethod = 'BH',
                    pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE)
    if (nrow(as.data.frame(ego)) > 0) {
        ego <- simplify(ego, cutoff = 0.7, by = 'p.adjust')   # collapse redundant parent-child GO terms
    }
    enrichment_results[[i]] <- ego
    message(sprintf('Cluster %d: %d significant GO terms', i, nrow(as.data.frame(ego))))
}

import gseapy as gp

all_temporal_genes = list(expr_sig.index)   # background = temporal genes, not the genome

for cl_id in range(n_clusters):
    cl_genes = [g for g, l in zip(expr_sig.index, labels) if l == cl_id]
    # enrichr hits the Enrichr web API; pass background=temporal genes and outdir=None (no files).
    # For a strictly offline hypergeometric test with a custom background, gp.enrich(gene_sets=<gmt/dict>,
    # background=all_temporal_genes) computes it locally instead.
    enr = gp.enrichr(gene_list=cl_genes, gene_sets='GO_Biological_Process_2023',
                     organism='human', background=all_temporal_genes, outdir=None)
    sig_terms = enr.results[enr.results['Adjusted P-value'] < 0.05]
    print(f'Cluster {cl_id}: {len(sig_terms)} significant GO terms')

clusters_with_terms <- sum(sapply(enrichment_results, function(x) nrow(as.data.frame(x)) > 0))
message(sprintf('Clusters with significant GO terms: %d / %d', clusters_with_terms, length(enrichment_results)))
if (clusters_with_terms < 3) message('WARNING: Few clusters enriched. Check gene ID mapping or thresholds.')