In silico inference of gene regulatory networks
2026-10-01
Source:vignettes/theory-network-inference.qmd
Note 1: Notation (aligned with the DeCovarT manuscript)
Symbol Role g=1,\ldots,G genes j=1,\ldots,J cell types i=1,\ldots,N purified samples / bootstrap replicates x_{gji} latent expression of gene g in cell type j, replicate i \mathcal{X}=(x_{gji}) tensor \mathcal{M}_{G\times J\times N} of single-cell (purified) profiles \boldsymbol{\mu}=(\mu_{gj}) mean signature \mathcal{M}_{G\times J}; column \boldsymbol{\mu}_{\cdot j} is the aggregated reference for type j \boldsymbol{\Sigma}_j cell-type covariance \mathcal{M}_{G\times G} \boldsymbol{\Omega}_j\equiv\boldsymbol{\Sigma}_j^{-1} sparse precision of the type-j GGM (not a generic parameter \theta) DeCovarT treats purified draws \boldsymbol{x}_{\cdot j}^{(i)}\sim\mathcal{N}_G(\boldsymbol{\mu}_{\cdot j}, \boldsymbol{\Sigma}_j) as latent and forms bulk mixtures \boldsymbol{y}_{\cdot i}=\sum_j p_j\,\boldsymbol{x}_{\cdot j}^{(i)}. GRN inference instead targets the conditional dependence structure among genes within each purified profile—or across a pooled design matrix. When \boldsymbol{p} is already known,
CSNetestimates type-level co-expression from the same convolution; the mapping onto DeCovarT’s precision \boldsymbol{\Omega}_j is in the perspectives vignette (Su et al. 2024).
General reviews
Marginal co-expression is a poor proxy for regulation. Shared upstream TFs, cell-type programmes, and common stimuli induce correlation without a direct edge; time lags and nonlinear kinetics can hide true edges. Conditional graphical models target direct dependence via the precision \boldsymbol{\Omega}, while causal / interventional designs are needed for direction and mechanism (Gene regulatory network inference in the era of single-cell multi-omics; Svinin and Glaab (2025); Network inference in systems biology).
Badia-i-Mompel et al. (2023) review single-cell multi-omics GRN tools and prior databases. Svinin and Glaab (2025) place association, prior-informed, and causal network analyses in one taxonomy. Together they motivate separating undirected skeletons (Sec. 2) from DAG / perturbation work (Sec. 3).
Undirected probabilistic strategies
General inference
Gaussian graphical models encode conditional independence in \boldsymbol{\Omega}. The MLE exists only under sample-size / chordality conditions formalised by Uhler et al. (2012): if n is too small relative to clique size, the covariance MLE fails. Sparse \ell_1 estimation (graphical lasso; Friedman et al. (2008)) is the workhorse when p\gg n (Chiquet 2015; Menéndez et al. 2010).
Two practical workflows appear repeatedly:
- Seeds then expand — select a small discriminating gene set (e.g. low condition number / high centroid separation, in the spirit of Aliee and Theis (2021)), then grow a prior-informed network with tools such as NeKo.
- Structure then parameters — estimate a sparse skeleton (glasso, NeighbourNet, CeSpGRN), then fit continuous weights, possibly under hard zeros on \Omega_{ij} for (i,j)\notin E.
Count-aware extensions replace the Gaussian likelihood by Poisson–log-normal (PLN) or copula models (Chiquet et al. 2021, 2018) and CeSpGRN, which matters when DeCovarT moves from Gaussian convolutions to sequencing counts.
Regularisation and topological constraints
Three standard devices:
\max_{\boldsymbol{\Omega} \succ 0} \log\det\boldsymbol{\Omega} - \operatorname{tr}(S\boldsymbol{\Omega}) \quad\text{s.t.}\quad \boldsymbol{\Omega}_{ij}=0\;\forall\,(i,j)\notin E \tag{1}
(hard covariance selection);
\max_{\boldsymbol{\Omega} \succ 0} \log\det\boldsymbol{\Omega} - \operatorname{tr}(S\boldsymbol{\Omega}) - \sum_{i\neq j}\lambda_{ij}\,|\boldsymbol{\Omega}_{ij}| \tag{2}
(edge-specific / hybrid penalties); and structured group or fused graphical lasso for hubs, modules, or shared biological constraints (Chiquet 2015). Trees and chordal graphs remain PD-friendly (Uhler et al. (2012)). DeCovarT’s affine spectral shift of a binary adjacency is a constructive PD completion compatible with hard zeros off E.
Algorithm unrolling and deep inference
Algorithm unrolling turns iterative sparse-precision solvers into trainable layers (Monga et al. (2021); GLAD). For hierarchical count ecosystems, Chaussard et al. (2024) introduce PLN-Tree: a tree-structured extension of PLN with backward Markov variational inference and identifiability results—an attractive route when taxonomy or gene modules should inform \boldsymbol{\Omega}.
Prior knowledge and PPI
Priors from TF–target, chromatin, and PPI databases can enter as \lambda_{ij} masks or as Bayesian edge probabilities (Badia-i-Mompel et al. (2023); Svinin and Glaab (2025); NeKo). This is the undirected analogue of pathway contextualisation in Sec. 3: association graphs become biologically constrained without claiming full causality.
Directed approaches
Gaussian Bayesian networks orient edges under Markov and faithfulness assumptions; interventions (Perturb-seq and related CRISPR screens) supply the environments needed for stronger causal claims (Peters et al. (2016); Svinin and Glaab (2025)). Pathway tools (NPA, SignalingProfiler, LIONESS/PANDA, NLBayes) contextualise expression on signed graphs but often stop short of a full weighted GRN MLE—treat them as priors or downstream scorers rather than drop-in replacements for \boldsymbol{\Omega}.
Note 2: From factorised Gaussian BNs to a joint multivariate normal
Consider a directed acyclic graph on G genes with topological order 1,\ldots,G. A linear Gaussian Bayesian network specifies, for each gene i,
X_i \mid \mathrm{Pa}(i) \sim \mathcal{N}\!\left( \mu_i + \sum_{j\in\mathrm{Pa}(i)} \beta_{ij} X_j,\, \sigma_i^{2} \right), \tag{3}
where \mathrm{Pa}(i) is the parent set in the DAG and \beta_{ij}=0 when j\not\in\mathrm{Pa}(i). Equivalently,
X_i = \mu_i + \sum_{j=1}^{G} b_{ij} X_j + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_i^{2}), \qquad b_{ij}=0 \text{ unless } j\in\mathrm{Pa}(i). \tag{4}
Collect coefficients in a strictly upper-triangular matrix B=(b_{ij}) (after ordering genes along the DAG) and noise variances D=\mathrm{diag}(\sigma_1^{2},\ldots,\sigma_G^{2}). Then Eq. 4 is \boldsymbol{X}=B\boldsymbol{X}+\boldsymbol{\varepsilon} with \boldsymbol{\varepsilon}\sim\mathcal{N}_G(0,D), hence
(\boldsymbol{I}_G-B)\boldsymbol{X}=\boldsymbol{\mu}+\boldsymbol{\varepsilon}. \tag{5}
When (\boldsymbol{I}_G-B) is invertible (always true for a DAG with consistent ordering),
\boldsymbol{X}\sim\mathcal{N}_G(\boldsymbol{\mu},\boldsymbol{\Sigma}), \qquad \boldsymbol{\Sigma} =(\boldsymbol{I}_G-B)^{-1} D (\boldsymbol{I}_G-B)^{-\mathsf{T}}. \tag{6}
The joint precision is
\boldsymbol{\Omega}=\boldsymbol{\Sigma}^{-1} =(\boldsymbol{I}_G-B)^{\mathsf{T}} D^{-1}(\boldsymbol{I}_G-B), \tag{7}
which is generally not sparse even when the DAG skeleton is sparse: moralisation and conditioning introduce fill-in. That is why undirected GGMs (Sec. 2) and directed BNs answer different questions: Eq. 6 is the directed generative model; Eq. 1 targets conditional independences in an undirected graph.
In R,
bnlearn::gbn2mvnorm()(Scutari 2026) returns (\boldsymbol{\mu}, \boldsymbol{\Sigma}) from a fitted Gaussian BN—the directed counterpart of taking \boldsymbol{\Sigma}=\boldsymbol{\Omega}^{-1} in an undirected GGM. For DeCovarT simulation, the same construction can populate cell-type-specific \boldsymbol{\Sigma}_j when a DAG is the ground truth for purified profiles \boldsymbol{x}_{\cdot j}^{(i)}, before bulk convolution with \boldsymbol{p}.
Benchmarks
Comparative reverse-engineering studies and DREAM-style challenges remain the default yardsticks for undirected and directed GRN callers (Penfold & Wild (2011); Network inference in systems biology; Menéndez et al. (2010)). Emulation tests for TF–target function complement topological accuracy when gold-standard DAGs are incomplete. Any DeCovarT-related network estimator should report against these protocols before biological claims.
Perspectives
Multi-omics
Collaborative / multi-layer glasso and tensor-fusion GNNs couple transcriptomic layers with metabolomic, methylation, or miRNA views. They generalise the single-precision story to coupled \boldsymbol{\Omega} blocks—relevant when references cease to be expression-only.
Boolean meets Bayesian
Noisy logic models (e.g. NLBayes) combine Boolean pathway structure with Bayesian TF-activity inference. They are powerful for regulon-level questions but are not substitutes for a dense weighted precision used inside DeCovarT’s Gaussian convolution.
Dynamics
Static snapshots average over regulatory delays. Time-series reviews (Natalia et al. (2023)), Dictys, and related dynamic biomarker methods target \boldsymbol{\Omega}_t along trajectories. Real-time or densely time-stamped assays would further reduce confounding between co-expression and causation.
Graph neural networks
OOD-robust GNNs and invariant causal prediction (Peters et al. (2016)) address batch / cell-type shift when message-passing replaces classical glasso. Unrolling (Sec. 2.3) remains preferable when interpretability of each layer as an optimisation step matters for statistical GRN recovery.
Link to DeCovarT simulation
simulate_hierarchical_grn_moments() currently draws a random-graph adjacency and completes it to a PD precision. The literature above suggests three near-term upgrades, in increasing ambition:
- replace the random skeleton by a prior-informed or glasso skeleton on real purified profiles;
- use edge-specific \lambda_{ij} (PPI / TF) instead of a uniform
precision_scale; - move from a shared Gaussian \boldsymbol{\Omega} to PLN / PLN-Tree or cell-type-specific copula models when count structure dominates.