Bio Temporal Genomics Temporal Grn

Version Compatibility

Reference examples tested with: statsmodels 0.14+, numpy 1.26+, pandas 2.2+, dynGENIE3 (GitHub vahuynh/dynGENIE3), bnlearn 4.9+, R 4.x

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: bulk time-series GRN inference is low-precision and assumption-heavy. Every edge is a HYPOTHESIS. Results are dominated by the sampling design (interval vs the minutes-scale of transcription, number of timepoints, replicate count), not by the algorithm. A tiny p-value from 6-12 timepoints is not evidence of regulation.

Temporal Gene Regulatory Network Inference

"Infer causal regulatory relationships from my time-series expression data" -> Rank directed, time-delayed TF->target edges from bulk temporal expression, to prioritize perturbation experiments.

• Python: statsmodels.tsa.stattools.grangercausalitytests() (VAR F-test on predictive precedence)
• R: dynGENIE3::dynGENIE3() (tree ensembles on ODE-derived derivatives); bnlearn::hc() + boot.strength() (dynamic Bayesian network)

The governing principle: inference produces ranked HYPOTHESES, not validated causal edges

Bulk temporal GRN inference turns a time course into a ranked list of candidate directed edges whose only honest downstream use is prioritizing perturbation experiments (knockdown / overexpression + re-measure). Two hard facts set the ceiling and must be stated up front, not buried.

1. Granger is PREDICTIVE precedence, not mechanism. It tests whether past X improves prediction of future Y, which is neither necessary nor sufficient for regulation. It collapses in three routine biological situations:
   • Unobserved common driver (confounding). An unmeasured TF, or a shared circadian/cell-cycle oscillation driving hundreds of genes, makes X "Granger-cause" Y with zero direct regulation. Pairwise methods are structurally blind to this; a shared sinusoid manufactures dense, entirely spurious directed structure whose lags are just phase offsets.
   • Sampling coarser than the regulation timescale (aliasing). Transcription acts in minutes; bulk courses are sampled every 1-6 h. When the interval exceeds the regulatory delay, cause and effect land in the same sampled timepoint and directionality becomes unidentifiable. No statistic recovers information the sampling threw away.
   • Non-stationarity. VAR-Granger assumes weak stationarity, but the interesting biology (a stimulus response, a developmental transient, a monotone induction) IS the non-stationary trend, and differencing it away removes the signal (see the differencing dilemma below).
2. Community benchmarks put a LOW ceiling on precision and no single method wins. DREAM5 (Marbach 2012 Nat Methods 9:796) evaluated 30+ methods and found time-series network inference is low-precision, no method is best across datasets, and the robust win is the "wisdom of crowds": integrating independent methods beats any one. Prior information (restricting regulators to known TFs) is the other reliable lever.

Operational consequence: restrict regulators to annotated TFs, run more than one method, keep edges recovered by >=2 methods and stable across replicate series, match density before comparing conditions, and hand the top edges to perturbation. This skill is bounded to BULK real-clock-time data; single-cell pseudotime GRN is a different problem (gene-regulatory-networks/scenic-regulons).

Method selection

Methodology evolves; verify current best practice against each tool's latest documentation before committing to one. The defensible default is to run more than one and intersect.

Granger causality (Python / statsmodels)

Goal: Rank TF->target pairs by whether past TF expression improves prediction of future target expression, with honest multiple-testing control.

Approach: Difference all genes uniformly to approach stationarity, select a single lag per pair by BIC (so the reported p-value is not the best-of-several), run ONE F-test at that lag, then BH-correct across pairs. Test only TF->target pairs to shrink the family and encode the TF prior.

The F-test compares an unrestricted VAR (Y on its own lags AND X's lags) to a restricted model (Y on its own lags only); statsmodels reports it as ssr_ftest, matching R's lmtest::grangertest. Two constraints dominate:

• Degrees-of-freedom floor. After lagging, n_eff = n - maxlag rows fit 2*maxlag+1 parameters, so the test is only defined for n > 3*maxlag + 1, and barely-defined means no power. With n=8 and maxlag=2 the F-test has ~1 residual DoF: a coin flip. This, not compute, is why genome-wide pairwise Granger fails. Prefer maxlag=1 on short courses.
• Lag selection is itself a multiple test. Taking the minimum p-value over lags 1..maxlag and reporting it as a single test inflates significance. Fix by selecting one lag a priori, or by BIC (below), or by Bonferroni across lags before the across-pairs BH.

import numpy as np
import pandas as pd
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import grangercausalitytests
from statsmodels.stats.multitest import multipletests

# expr_df: genes x timepoints DataFrame; columns MUST be in temporal order.
# Difference uniformly to approach stationarity. Uniform (not per-gene) differencing
# keeps every series on the same footing: mixing I(0) and differenced I(1) series in one
# VAR corrupts the F-test reference distribution. Cost: over-differencing already-stationary
# genes. The deeper tradeoff: differencing removes the trend that CARRIES the regulatory
# signal, so on short courses prefer maxlag=1 over aggressive differencing.
expr_diff = expr_df.diff(axis=1).iloc[:, 1:]

tf_genes = ['TF1', 'TF2', 'TF3']
target_genes = ['geneA', 'geneB', 'geneC']
maxlag = 1  # short courses have ~no DoF beyond lag 1 (need n > 3*maxlag+1)

def granger_pvalue(pair_data, maxlag):
    # column 0 = response Y (target), column 1 = predictor X (TF): tests X -> Y.
    # Select ONE lag by BIC, then run a SINGLE test at it -> avoids the min-p-over-lags
    # multiple test. Guard BIC=0 (no lag structure) up to 1.
    lag = max(1, int(VAR(pair_data).select_order(maxlag).bic))
    res = grangercausalitytests(pair_data, maxlag=[lag])  # list -> tests only this lag
    return res[lag][0]['ssr_ftest'][1], lag  # (p_value, lag)

records = []
for tf in tf_genes:
    for target in target_genes:
        if tf == target:
            continue
        pair = np.column_stack([expr_diff.loc[target].values, expr_diff.loc[tf].values])
        p, lag = granger_pvalue(pair, maxlag)
        records.append({'tf': tf, 'target': target, 'p_value': p, 'lag': lag})

results_df = pd.DataFrame(records)
# multipletests default is Holm-Sidak, NOT BH; force fdr_bh explicitly.
results_df['q_value'] = multipletests(results_df['p_value'], method='fdr_bh')[1]
significant = results_df[results_df['q_value'] < 0.05].sort_values('q_value')

Pairwise Granger cannot separate direct regulation from a chain X->Z->Y or a fork Z->{X,Y}. The correct fix is conditional (multivariate) Granger, conditioning on all other regulators' lags, but that explodes the parameter count and is infeasible at transcriptomic sample sizes. Label pairwise output as a CONFOUNDED candidate set, not direct interactions.

dynGENIE3 (R)

Goal: Rank regulator->target edges non-linearly by how much a regulator's current expression predicts a target's temporal derivative.

Approach: dynGENIE3 models each gene as dx_i/dt = f_i(x) - alpha_i * x_i, estimates dx_i/dt by finite differences between consecutive timepoints, and trains a tree ensemble to regress that derivative-plus-decay target on candidate-regulator expression; summed variable importance becomes the edge weight.

library(dynGENIE3)

# TS.data: list of genes x timepoints matrices (one per replicate/series).
# time.points: matching list of time vectors (real deltas -> handles uneven spacing).
expr_list <- list(as.matrix(expr_series1), as.matrix(expr_series2), as.matrix(expr_series3))
time_list <- list(c(0, 4, 8, 12, 24, 48), c(0, 4, 8, 12, 24, 48), c(0, 4, 8, 12, 24, 48))

# Restrict regulators to known TFs (AnimalTFDB / PlantTFDB). This helps TWICE: fewer
# features searched per split (faster) AND a non-TF can never be reported as a regulator
# (higher precision). Single highest-yield precision lever.
tf_indices <- which(rownames(expr_list[[1]]) %in% tf_names)

# tree.method DEFAULTS to 'RF' (Random Forests). Extra-Trees is opt-in: tree.method='ET'
# (the config GENIE3 used to win DREAM4). alpha='from.data' (default) estimates per-gene
# mRNA decay from the data; pass a numeric vector to inject measured half-lives (4sU/BRIC-seq).
res <- dynGENIE3(TS.data = expr_list, time.points = time_list, regulators = tf_indices)

# get.link.list (DOT form) is the dynGENIE3 function. The camelCase getLinkList belongs to
# the separate Bioconductor GENIE3 package -- do not swap them.
link_list <- get.link.list(res$weight.matrix, report.max = 1000)

Two properties gate interpretation:

• The weight matrix is a RANKING with no null, no p-value, no calibrated threshold. A "top edge" is top only relative to the others in this run; thresholding by rank (top-K) is unavoidably arbitrary. This is why cross-method agreement and stability matter more here than anywhere.
• Finite-difference derivatives amplify noise. With few, unevenly-spaced timepoints (0,4,8,12,24,48 h is typical) each dx/dt rests on one noisy pair and late wide intervals blur short-timescale regulation into a single slope. More REPLICATES (independent derivative samples averaging the noise down) help far more than adding one or two timepoints.