Bio Temporal Genomics Temporal Clustering

Version Compatibility

Reference examples tested with: Mfuzz 2.64+, TCseq 1.14+, DEGreport 1.30+ (R/Bioconductor); tslearn 0.8+, scikit-learn 1.4+ (Python).

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

• Python: pip show tslearn scikit-learn then help(module.function) to check signatures
• R: packageVersion('Mfuzz') 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.

Temporal Gene Clustering

"Group my time-course genes by expression pattern shape" -> Partition PRE-SELECTED temporally variable genes into co-expression modules by trajectory shape (fuzzy c-means, hierarchical, or DTW), producing candidate temporal programs.

• R: Mfuzz::mfuzz() (fuzzy/soft), TCseq::timeclust(), DEGreport::degPatterns()
• Python: tslearn.clustering.TimeSeriesKMeans (Euclidean / DTW / soft-DTW)

Governing Principle - read before clustering anything

Clustering answers only "which genes share a temporal SHAPE." It is DESCRIPTIVE and UNSUPERVISED: it has no null model, no p-value, and no notion of a "true" cluster count, so it ALWAYS returns clusters from whatever it is handed. It is strictly DOWNSTREAM of gene selection.

• It does NOT answer "which genes are rhythmic" (that is temporal-genomics/circadian-rhythms) and does NOT answer "is this gene significantly changing" (that is differential-expression/timeseries-de: LRT, spline-DE, maSigPro). Clustering adds description, not inference.
• The input MUST already be the temporally variable genes - the output of timeseries-DE or, at minimum, a variance filter. Never the full expression matrix.
• Feeding in flat/all genes is the #1 error. Per-gene z-scoring (mandatory, below) rescales a flat gene's pure noise to unit variance, so it lands in a "cluster" of noise that mimics a real program. Z-scoring erases the one signal (near-zero variance) that flagged the gene as flat, which is exactly why prefiltering is a gate, not optional hygiene.
• Clusters are HYPOTHESES. A centroid is a candidate program; membership is not evidence a gene is regulated - that evidence came (or did not) from the upstream DE step.

If a user asks "cluster my RNA-seq time course," the first question is always: have these genes already been selected for temporal change, and how? If the answer is "no, it is all 20,000 genes," stop and prefilter.

Core Workflow

1. Confirm the input is pre-selected temporally variable genes (DE hits or top-variance); if not, prefilter
2. Standardize each gene's profile (z-score across timepoints) - mandatory
3. Choose a distance metric (Euclidean-on-zscore / correlation / DTW), then an algorithm and k
4. Assign genes to clusters (soft membership or hard labels); filter by membership if fuzzy
5. Validate by stability (bootstrap/consensus), then interpret centroids and run per-cluster enrichment with the correct background

Soft vs Hard, and Why Standardization Is Mandatory

Soft (fuzzy) clustering is preferred for expression. Genes participate in multiple regulatory programs, so forcing one gene into one cluster (hard k-means) is biologically false at boundaries and brittle: a gene between two centroids flips clusters under trivial noise. Futschik & Carlisle (2005) established fuzzy c-means as noise-ROBUST for expression time courses - low-membership (ambiguous, likely-noise) genes are down-weighted in centroid estimation, so centroids track the high-confidence core of each program, and ambiguity is exposed as a continuous membership score to threshold rather than hidden inside a hard label.

Z-score per gene is mandatory (Mfuzz standardise(), TCseq standardize=TRUE, tslearn TimeSeriesScalerMeanVariance()). Without it, MAGNITUDE dominates SHAPE: a high-abundance housekeeping gene sits far (Euclidean) from a low-abundance gene of identical shape, while two high-abundance genes co-cluster on abundance alone. Clustering-by-shape requires removing each gene's mean and scaling to unit variance across timepoints.

Mfuzz (R/Bioconductor)

Goal: Group temporally variable genes into soft co-expression clusters by trajectory shape.

Approach: Build an ExpressionSet, gate out flat genes (filter.std), z-score (standardise), estimate then VALIDATE the fuzzifier, run fuzzy c-means, and filter genes by membership. Mfuzz wraps e1071::cmeans (it does not implement its own optimizer); distance is Euclidean on z-scored profiles.

Setup and Preprocessing

library(Mfuzz)
library(Biobase)

# Rows = genes (already selected as temporally variable), columns = timepoints (mean across replicates)
expr_mat <- as.matrix(read.csv('temporal_expression.csv', row.names = 1))
eset <- ExpressionSet(assayData = expr_mat)

# filter.std: flat-gene GATE (keeps the governing principle true). min.std=0.5 is a starting
# point; inspect the SD distribution and set it above the flat-gene noise floor for your data.
eset <- filter.std(eset, min.std = 0.5)

# Per-gene mean 0, sd 1 across timepoints (British spelling; no 'standardize' alias)
eset <- standardise(eset)

Fuzzifier Estimation - inspect, do not trust blindly

Goal: Pick a fuzzifier m that keeps clusters informative for THIS number of timepoints.

Approach: mestimate() implements Schwaemmle & Jensen (2010): it returns the smallest m that stops fuzzy c-means from finding tight clusters in RANDOMIZED data. The estimate is dominated by D (number of timepoints) via a D^-2 term, so it can go degenerate at the extremes - inspect the returned m AND the membership distribution rather than trusting either the estimate or the historical m=2 default.

# With FEW timepoints (small D), mestimate pushes m HIGH -> over-fuzzy: memberships flatten
# toward 1/c and an acore(0.5) filter can discard nearly everything.
# With MANY timepoints (large D), m falls toward ~1.05-1.2 -> near-hard, soft advantage evaporates.
m <- mestimate(eset)
cat(sprintf('Estimated fuzzifier m: %.2f\n', m))

cl <- mfuzz(eset, c = 8, m = m)  # c=8: starting point for 6-12 timepoints; refine below

# VALIDATE m: what fraction of genes clears the alpha-core cutoff? If very few do, m is too high.
max_mem <- apply(cl$membership, 1, max)
cat(sprintf('Genes with max membership >= 0.5: %.0f%%\n', 100 * mean(max_mem >= 0.5)))
# Sanity check the estimate's own criterion: cluster a permuted copy; it should NOT form tight clusters.

Membership Filtering and Cluster Selection

# acore returns, per cluster, genes with MAX membership >= min.acore ("alpha cores").
# 0.5 is a convention; it discards a data-dependent fraction (larger m -> more discarded).
# Relaxing to 0.3 is legitimate for exploratory work but admits more noise. Always report the retained fraction.
core_genes <- acore(eset, cl, min.acore = 0.5)

# Minimum centroid distance vs k: as k grows the closest centroid pair collapses; a knee hints at
# over-splitting. This is a WEAK, monotone-ish signal, not an oracle -- triangulate with stability (below).
min_dist <- sapply(4:20, function(k) {
    d <- as.matrix(dist(mfuzz(eset, c = k, m = m)$centers))
    diag(d) <- Inf
    min(d)
})
plot(4:20, min_dist, type = 'b', xlab = 'k', ylab = 'Min centroid distance')

Visualization

mfuzz.plot2(eset, cl, mfrow = c(2, 4), time.labels = colnames(expr_mat), centre = TRUE, x11 = FALSE)
overlap.plot(cl, over = overlap(cl), thres = 0.05)  # centroid-overlap view; merges hint at over-clustering

TCseq (R/Bioconductor)

TCseq was built for time-course SEQUENCING (RNA-seq/ATAC-seq); upstream DE/peak steps live in the same package, and timeclust clusters the summarized (per-gene, per-timepoint) matrix.

library(TCseq)

# algo='cm': fuzzy c-means (soft, Mfuzz-like). Also 'km' (hard k-means), 'pam', 'hc' (hierarchical).
# standardize=TRUE does the mandatory per-gene z-score.
tc <- timeclust(expr_mat, algo = 'cm', k = 6, standardize = TRUE)
timeclustplot(tc, value = 'z-score', cols = 3)

tc_km <- timeclust(expr_mat, algo = 'km', k = 6, standardize = TRUE)  # hard alternative

DEGreport degPatterns (R)

Goal: Hierarchical clustering with automatic k and design-aware grouping.

Approach: degPatterns takes replicate-level data plus metadata, collapses samples within each (time, col) group to a MEAN internally, then clusters on correlation distance and cuts the tree. Convenient, but "auto k" is really "cut + merge under minc," a heuristic - not an optimum.

library(DEGreport)

# time, col: COLUMN NAMES in metadata (col defaults to NULL). minc=15: minimum cluster size;
# clusters smaller than minc are DROPPED -- this both blocks singletons AND silently discards genes,
# so it can yield fewer clusters than the tree suggested. Set deliberately.
patterns <- degPatterns(expr_mat, metadata = sample_info, time = 'timepoint', col = 'condition', minc = 15)

cluster_df <- patterns$df                     # gene -> cluster assignments
degPlotCluster(patterns$normalized, time = 'timepoint', color = 'condition')  # note: 'color', not 'col'

tslearn (Python) - Euclidean / DTW / soft-DTW

Goal: Cluster time-series profiles, optionally warping the time axis for phase-shifted genes.

Approach: Z-score, then TimeSeriesKMeans. The DISTANCE METRIC matters more than the algorithm - default to Euclidean-on-zscore (which, after standardization, is monotone in Pearson correlation and captures "same shape, different amplitude"). Escalate to DTW ONLY for real, expected phase shifts, and ALWAYS constrain it.

import numpy as np
from tslearn.clustering import TimeSeriesKMeans, silhouette_score
from tslearn.preprocessing import TimeSeriesScalerMeanVariance

# expr_mat: (n_genes, n_timepoints) of PRE-SELECTED temporally variable genes
expr_scaled = TimeSeriesScalerMeanVariance().fit_transform(expr_mat[:, :, np.newaxis])

# Default, safe choice: Euclidean on z-scored profiles (phase-SENSITIVE, cheap, no fabricated structure)
model = TimeSeriesKMeans(n_clusters=8, metric='euclidean', max_iter=50, random_state=42)
labels = model.fit_predict(expr_scaled)

DTW - powerful for phase shifts, but constrain the band or it invents structure

DTW (Sakoe & Chiba 1978) warps the time axis so a profile peaking one timepoint later can still match - the ONLY reason to reach for it (signaling cascades, developmental heterochrony, unequal sampling). Its default failure mode is the SINGULARITY: unconstrained DTW maps one point of series A onto a long run of points of series B, manufacturing apparent co-regulation from noise. tslearn's default global_constraint=None is exactly this singularity-prone configuration.

# The Sakoe-Chiba BAND caps how far in time a point may be matched -- kills most singularities AND
# cuts cost. This constraint is mandatory, not optional, for DTW clustering.
# sakoe_chiba_radius: warping-window half-width in timepoints; small (1-2) for tight sampling.
model = TimeSeriesKMeans(
    n_clusters=8, metric='dtw',
    metric_params={'global_constraint': 'sakoe_chiba', 'sakoe_chiba_radius': 2},
    max_iter=50, random_state=42)
labels = model.fit_predict(expr_scaled)

# Soft-DTW: replaces DTW's hard min with a soft-min -> DIFFERENTIABLE loss, enabling proper
# soft-DTW barycenters (cluster centers). It is NOT "faster" -- still quadratic; use it for smooth,
# well-defined averaging, not speed. gamma via metric_params (NOT the deprecated gamma_sdtw kwarg).
soft = TimeSeriesKMeans(n_clusters=8, metric='softdtw', metric_params={'gamma': 0.5},
                        max_iter=50, random_state=42)

When DTW is worth it: only when phase shift is real and expected, the band is set, AND DTW has been checked against fabricating structure. On data with NO phase shifts, DTW should not beat Euclidean - if it "finds more clusters" there, that is invented structure, not signal.

Selecting k - score under the SAME geometry that formed the clusters

# Scoring DTW clusters with a EUCLIDEAN silhouette is geometrically inconsistent: clusters were
# formed under DTW geometry but ranked under Euclidean, which can pick a DIFFERENT (wrong) k.
# tslearn.clustering.silhouette_score takes metric='dtw'/'softdtw' and precomputes the matching
# distances internally -- score under the SAME geometry that formed the clusters.
dtw_params = {'global_constraint': 'sakoe_chiba', 'sakoe_chiba_radius': 2}
scores = {}
for k in range(3, 11):
    km = TimeSeriesKMeans(n_clusters=k, metric='dtw', metric_params=dtw_params, max_iter=30, random_state=42)
    labels_k = km.fit_predict(expr_scaled)
    scores[k] = silhouette_score(expr_scaled, labels_k, metric='dtw', metric_params=dtw_params)
best_k = max(scores, key=scores.get)

Under a pure-Euclidean pipeline, sklearn.metrics.silhouette_score(expr_scaled.squeeze(), labels) is consistent and fast. It is only the DTW/Euclidean MISMATCH that mis-ranks k.

Choosing k - the honest story

No index is authoritative; triangulate and let biology and stability decide.

Over-clustering FRAGMENTS one real program across centroids (the same GO terms then reappear in three clusters); under-clustering MERGES distinct programs into an averaged centroid matching no gene. Report a stable partition, not a single silhouette peak.