Paper deep dive
CoDiffGRN: Rethinking Gene Regulatory Network Inference via the BEELINE-KGC Benchmark and Co-evolutionary Discrete Diffusion
Jiaze Song, Runhao Zhao, Minghao Xu, Bin Cui, Wentao Zhang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 7/16/2026, 3:29:41 AM
Summary
The paper addresses limitations in gene regulatory network (GRN) inference from single-cell transcriptomic data, specifically the misalignment of existing benchmarks with real-world needs for inductive generalization and top-ranked predictions. It introduces BEELINE-KGC, a benchmark that reformulates GRN inference as an inductive knowledge graph completion task using ranking metrics like Hits@K and MRR. To solve this, the authors propose CoDiffGRN, a co-evolutionary discrete diffusion framework that jointly models discretized gene expression states and regulatory interactions. CoDiffGRN incorporates cell-cluster discretization, a joint node-edge diffusion process, and TF-ALL Subgraph Sampling (TASS) for scalable training, achieving state-of-the-art performance in novel regulatory discovery.
Entities (10)
Relation Signals (10)
CoDiffGRN → addresses → Gene Regulatory Network
confidence 98% · Inferring gene regulatory networks (GRNs) from single-cell transcriptomic data is crucial for biological discovery
BEELINE-KGC → reformulatesas → Knowledge Graph Completion
confidence 96% · By recasting evaluation as a Knowledge Graph Completion (KGC) task [4], BEELINE-KGC uses Hits@K and Mean Reciprocal Rank (MRR)
BEELINE-KGC → evaluateswith → Hits@K
confidence 95% · BEELINE-KGC uses Hits@K and Mean Reciprocal Rank (MRR) to directly measure the top-KK ranking consistency
BEELINE-KGC → evaluateswith → MRR
confidence 95% · BEELINE-KGC uses Hits@K and Mean Reciprocal Rank (MRR) to directly measure the top-KK ranking consistency
CoDiffGRN → uses → TASS
confidence 95% · We further introduce TF-ALL Subgraph Sampling (TASS) for scalable training.
scRNA-seq → feeds → Gene Regulatory Network
confidence 94% · Inferring gene regulatory networks (GRNs) from single-cell transcriptomic data is crucial for biological discovery
CoDiffGRN → uses → Cell-cluster discretization
confidence 94% · we introduce a Cell-cluster-based Discretization strategy that transforms scRNA-seq profiles into biologically coherent discrete states
BEELINE-KGC → →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Inferring gene regulatory networks (GRNs) from single-cell transcriptomic data is crucial for biological discovery, yet existing approaches suffer from a fundamental misalignment with real-world needs. Researchers typically seek a small set of high-confidence regulatory interactions for experimental validation, often involving previously unseen genes. However, current benchmarks rely on transductive splits with global classification metrics, while prevailing models struggle to generalize under inductive settings. To bridge this gap, we reformulate GRN inference as an inductive, ranking-centric graph completion problem and introduce \textbf{\benchmark}, a new benchmark that incorporates an inductive gene-holdout split together with knowledge graph completion metrics to better evaluate top-ranked predictions. Building on this, we propose \textbf{\method}, the first co-evolutionary discrete diffusion framework that jointly models biologically coherent discretized gene expression states and regulatory interactions for robust inductive generalization and improved top-ranked regulatory discovery. We further introduce TF-ALL Subgraph Sampling (TASS) for scalable training. Extensive experiments on {\benchmark} show that {\method} establishes new state-of-the-art performance, significantly outperforming existing methods in novel regulatory discovery, and ablation studies further verify the effectiveness of our design.
Tags
Links
- Source: https://arxiv.org/abs/2607.13120v1
- Canonical: https://arxiv.org/abs/2607.13120v1
Trouble viewing inline? Open PDF directly →
Full Text
62,848 characters extracted from source content.
Expand or collapse full text
CoDiffGRN: Rethinking Gene Regulatory Network Inference via the BEELINE-KGC Benchmark and Co-evolutionary Discrete Diffusion Jiaze Song Peking University jzsong25@stu.pku.edu.cn &Runhao Zhao National University of Defense Technology runhaozhao@nudt.edu.cn &Minghao Xu Peking University minghao.xu@stu.pku.edu.cn &Bin Cui Peking University bin.cui@pku.edu.cn &Wentao Zhang Peking University wentao.zhang@pku.edu.cn Corresponding author Abstract Inferring gene regulatory networks (GRNs) from single-cell transcriptomic data is crucial for biological discovery, yet existing approaches suffer from a fundamental misalignment with real-world needs. Researchers typically seek a small set of high-confidence regulatory interactions for experimental validation, often involving previously unseen genes. However, current benchmarks rely on transductive splits with global classification metrics, while prevailing models struggle to generalize under inductive settings. To bridge this gap, we reformulate GRN inference as an inductive, ranking-centric graph completion problem and introduce BEELINE-KGC, a new benchmark that incorporates an inductive gene-holdout split together with knowledge graph completion metrics to better evaluate top-ranked predictions. Building on this, we propose CoDiffGRN, the first co-evolutionary discrete diffusion framework that jointly models biologically coherent discretized gene expression states and regulatory interactions for robust inductive generalization and improved top-ranked regulatory discovery. We further introduce TF-ALL Subgraph Sampling (TASS) for scalable training. Extensive experiments on BEELINE-KGC show that CoDiffGRN establishes new state-of-the-art performance, significantly outperforming existing methods in novel regulatory discovery, and ablation studies further verify the effectiveness of our design. 1 Introduction Gene regulatory networks (GRNs) dictate cellular identity and physiological function by encoding the complex interactions between transcription factors (TFs) and target genes [37, 8, 14, 39]. While the advent of single-cell RNA sequencing (scRNA-seq) has fundamentally advanced GRN inference by resolving bulk-level signal averaging [21], the resulting data is notoriously heterogeneous and sparse [34, 20]. To tackle these challenges, recent progress has shifted from classical statistical methods to graph neural networks (GNNs), particularly graph autoencoder (GAE) architectures [7, 25, 46, 15, 42, 36]. Despite these architectural advances, we argue that computational GRN inference suffers from a fundamental misalignment with practical biological discovery in two critical dimensions: First, there is an evaluation misalignment. Standard benchmarks such as BEELINE [29] evaluate GRN inference under transductive settings, where all genes are observed during training. However, discovering regulatory mechanisms for novel modules or rare cell types inherently requires inductive generalization to unseen genes [9, 11]. Moreover, these benchmarks formalize GRN inference as a binary classification task measured by AUROC and AUPRC. In practice, however, biologists rely on AI models to prioritize a small, high-confidence set of novel interactions for costly experimental validation. Global metrics such as AUROC are easily dominated by low-ranked edges, making them poor proxies for the top-K ranking quality that ultimately drives scientific discovery. Figure 1: Overview of current limitations and our proposed solutions. (a) Current benchmarks mismatch with realistic demand of inductive generalization and top-K ranking quality, while existing methods degrade with unseen genes and fixed edge transitions. (b) Our BEELINE-KGC introduces an inductive, ranking-based evaluation protocol, and CoDiffGRN leverages cell-cluster discretization, co-evolutionary discrete diffusion with TASS. Abbr., G.T.: Ground-truth, N.S.: Negative Sample. Second, there exists a gene-level inductive generalization gap. Under realistic inductive settings, existing GAE-based methods rely heavily on predefined prior graphs. When holdout genes become disconnected during inference, they lose critical contextual signals and suffer severe performance degradation (Figure 1a). Recent diffusion-based methods partially alleviate this topological isolation, but still fail to generalize effectively to unseen genes, as the highly heterogeneous and sparse nature of scRNA-seq data prevents reliable feature correspondence and can lead to feature collapse. More fundamentally, current approaches fail to capture the conditional nature of gene regulations. Regulatory interactions are not isolated events, but are conditioned on the paired states of regulator and target genes. Applying a uniform transition to all edges during diffusion, regardless of the activity states of their endpoint genes, violates this fundamental principle of GRNs. To address these limitations, we introduce a new evaluation paradigm alongside a novel generative framework. First, we propose BEELINE-KGC, a systematic extension of BEELINE that enforces an Inductive Gene-Holdout Split. By recasting evaluation as a Knowledge Graph Completion (KGC) task [4], BEELINE-KGC uses Hits@K and Mean Reciprocal Rank (MRR) to directly measure the top-K ranking consistency of regulatory predictions. Under this rigorous protocol, existing state-of-the-art models exhibit severe generalization bottlenecks. On the modeling side, we propose CoDiffGRN, the first gene-regulation co-evolutionary discrete diffusion model for GRN inference. Built upon the D3PM framework [1], CoDiffGRN departs from prior-dependent paradigms and is designed for robust inductive generalization. To mitigate feature collapse on unseen genes, we introduce a Cell-cluster-based Discretization strategy that transforms scRNA-seq profiles into biologically coherent discrete states, allowing unseen genes to inherit prior distributional structure during diffusion. To capture the conditional nature of gene regulation, we introduce a joint node-edge diffusion process that directly models the co-evolution of discretized gene expression states and network topology at the level of diffusion dynamics. By devising a gene-state-dependent transition matrix, edge transitions are explicitly conditioned on the states of their endpoint genes, promoting the coherent neighborhoods required for top-K regulatory discovery. Furthermore, to enable scalable generative training under data scarcity, we introduce TF-ALL Subgraph Sampling (TASS), a data augmentation strategy that preserves TF-target integrity while substantially improving scalability and generalization. Comprehensive experimental results on BEELINE-KGC show that CoDiffGRN achieves state-of-the-art performance and consistently outperforms existing models. Extensive ablation studies examine the contributions of the proposed gene-regulation co-evolutionary discrete diffusion and TASS framework. These results collectively support the effectiveness of our approach in improving KGC-based ranking performance, enabling more accurate and reliable biological discovery. 2 Related Work Evaluation protocols for GRN inference. BEELINE [29] established the standard for GRN inference, primarily utilizing transductive edge-centric splits and global metrics such as AUROC and AUPRC [18]. To improve discriminative sensitivity, negative sampling has evolved from random selection [26, 43] to heuristic strategies [36, 46]. However, these protocols exhibit a fundamental misalignment with practical biological discovery: transductive settings neglect the imperative for generalization to unseen genes, while global metrics fail to prioritize the top-K candidates required for experimental validation. In this work, we introduce BEELINE-KGC to bridge this gap by recasting evaluation as an inductive KGC task, utilizing ranking-based metrics (Hits@K, MRR) to assess top-K regulatory predictions. Gene regulatory network inference methods. Early GRN inference relied on statistical associations like GENIE3 [17] and GRNBoost2 [26], or feature-based deep learning models such as CNNC [43] and DeepSEM [30]. Recently, graph representation learning has become the dominant paradigm, with GAE-based frameworks like GENELink [7], GENELink+ [46] and GNNLink [25] utilizing prior regulatory graphs, while methods like GCLink [42] and GRNFormer [15] introduce contrastive learning and VGAE-based frameworks. More recently, generative models such as RegDiffusion [47], DigNet [35], and Planet [41] have applied continuous and discrete diffusion processes [16, 1] to model the evolution nature of gene regulations. However, existing methods face critical limitations: traditional approaches struggle with noisy single-cell data, while graph-based models often hinge on prior-graph connectivity. Furthermore, current generative frameworks struggle with feature collapse in inductive scenarios and fail to capture the essential co-evolution between gene and regulation with severe data scarcity. In this work, we propose CoDiffGRN, which employs cell-cluster discretization and joint discrete diffusion to model regulatory dynamics, alongside TASS specifically designed to address the challenges of data scarcity. 3 Rethinking for Practical Discovery: The BEELINE-KGC Benchmark Table 1: Statistics of the BEELINE-KGC benchmark. For each cell type and reference network, we report the number of cells, the number of source TFs, and the edge counts under inductive splits. Values are shown as TFs+500 (TFs+1000). Abbr., S. TFs: Source TFs. Cell Types #Cells Specific Non-Specific #S. TFs #Genes #Train #Valid #Test #S. TFs #Genes #Train #Valid #Test hESC 758 34 (34) 815 (1260) 3535 (5995) 505 (544) 505 (545) 283 (292) 753 (1138) 2546 (3531) 447 (543) 448 (543) hHEP 425 30 (31) 874 (1331) 7724 (12466) 1107 (1546) 1108 (1546) 322 (332) 825 (1217) 2993 (3490) 568 (930) 568 (931) mDC 383 20 (21) 443 (684) 536 (903) 110 (145) 110 (145) 250 (254) 634 (969) 2115 (2704) 476 (607) 476 (607) mESC 421 88 (89) 977 (1385) 22694 (33645) 3459 (4575) 3460 (4575) 516 (522) 890 (1214) 4677 (5738) 1108 (1146) 1108 (1146) mHSC-E 1071 29 (33) 691 (1177) 9286 (17825) 1135 (2075) 1136 (2075) 144 (147) 442 (674) 1002 (1370) 211 (295) 212 (295) mHSC-GM 889 22 (23) 618 (1089) 5796 (10901) 784 (1617) 784 (1617) 82 (88) 297 (526) 496 (925) 123 (216) 124 (217) mHSC-L 847 16 (16) 525 (640) 3354 (3887) 522 (646) 522 (647) 35 (37) 164 (192) 175 (217) 52 (50) 52 (50) Cell Types #Cells STRING LOF/GOF #S. TFs #Genes #Train #Valid #Test #S. TFs #Genes #Train #Valid #Test hESC 758 343 (351) 511 (695) 2862 (3433) 697 (858) 698 (858) – – – – – hHEP 425 409 (414) 646 (874) 4960 (6269) 1281 (1367) 1282 (1367) – – – – – mDC 383 264 (273) 479 (664) 3277 (4066) 769 (916) 769 (916) – – – – – mESC 421 495 (499) 638 (785) 5099 (5830) 1331 (1324) 1332 (1325) 34 (34) 775 (1099) 3322 (4851) 423 (445) 424 (446) mHSC-E 1071 156 (161) 291 (413) 950 (1330) 210 (248) 211 (248) – – – – – mHSC-GM 889 92 (100) 201 (344) 505 (1108) 121 (101) 122 (102) – – – – – mHSC-L 847 39 (40) 70 (81) 89 (101) 24 (26) 24 (27) – – – – – We introduce BEELINE-KGC, a comprehensive benchmark designed to bridge the gap between computational GRN inference and realistic biological discovery. Built upon the BEELINE [29] (C BY-NC 4.0) collection of experimental scRNA-seq datasets, BEELINE-KGC establishes a rigorous evaluation framework tailored to real-world laboratory constraints. 3.1 Formalization and Problem Setting We formalize a gene regulatory network as a directed graph =(,ℰ)G=(V,E), where nodes v∈v represent genes and edges (u,v)∈ℰ(u,v) denote regulatory interactions. We partition V into TFs (TFV_TF) and target genes (tgtV_tgt). Within the TFs, we further distinguish source TFs (src⊆TFV_src _TF), defined as regulators with at least one observed outgoing edge in the ground truth. The corresponding scRNA-seq data is represented by an expression matrix ∈ℝN×FXX ^N× F_X, where i,jX_i,j denotes the expression level of gene vi∈v_i in cell j∈1,…,FXj∈\1,…,F_X\. Given the expression matrix X and a subset of observed interactions ℰobsE_obs, the primary objective of GRN inference is to recover the unobserved regulatory structure by predicting potential links in ℰ∖ℰobsE _obs. 3.2 A Unified Framework for Practical GRN Discovery Standard GRN inference benchmarks adopt a transductive setting and rely on global metrics such as AUPRC. However, this paradigm fundamentally diverges from the practical demands of real biological discovery, which requires inductive generalization to novel regulators and high top-K ranking quality under limited experimental budgets. Pillar 1: Inductive gene-holdout split. In realistic discovery scenarios, researchers predominantly focus on characterizing novel or unannotated TFs by identifying their downstream regulatory targets, necessitating an inductive setting where models must generalize to unseen regulators [31, 38, 45, 44, 28, 2, 9, 11]. To simulate this, we extract a subset of srcV_src as holdout TFs (holdV_hold), introducing an inductive protocol: (1) Training: The model observes only the subgraph induced by ∖holdV _hold. Expression profiles and all edges involving holdout TFs are withheld. (2) Evaluation: The model receives the expression profiles for the full gene set V and predicts regulations across the entire network. Performance is evaluated solely on edges involving at least one holdout TF. Pillar 2: KGC reformulation and multi-scale ranking protocol. Practical biological discovery is governed by a discrete experimental budget K, where true interactions are useful only if ranked within the top-K and lower ranks are effectively equivalent to being undiscovered. However, global metrics like AUPRC are insensitive to top-ranked precision, and are further distorted by the severe class imbalance in GRNs leaving many candidate pairs unevaluated under standard negative sampling. To resolve this, we reformulate GRN inference as a Knowledge Graph Completion (KGC) problem. Each regulation is modeled as a structural triple (vi,regulates,vj)(v_i, regulates,v_j), where viv_i is the source TF (head) and vjv_j is the target gene (tail). The task of finding targets for a holdout TF naturally translates to tail entity prediction: given a query (vi,regulates,?)(v_i, regulates,?), the model must rank all candidate tails t∈t such that the ground-truth targets are placed as high as possible. Aligned with this KGC formulation, we evaluate GRN inference as a resource-constrained ranking problem, adopting Hits@K for top-K recovery under a fixed validation budget and Mean Reciprocal Rank (MRR) for capturing overall ranking quality across heterogeneous budgets: Hits@K=1||∑q∈[rank(tq)≤K],MRR=1||∑q∈1rank(tq),Hits@K= 1|Q| _q 1[rank(t_q)≤ K], = 1|Q| _q 1rank(t_q), (1) where Q is the set of evaluation queries targeting the holdout TFs, and rank(tq)rank(t_q) is the rank of the ground-truth tail. Following standard KGC practice, we employ the filtered setting, which excludes other known valid targets for the same query when calculating the rank. This reformulation guarantees that our benchmark directly quantifies a model’s practical capacity to accelerate discovery within realistic experimental budgets. 4 Modeling Gene-Regulation Co‑evolution: The CoDiffGRN Framework 4.1 Preliminaries Diffusion models consist of a forward noising process and a reverse denoising process [16]. The forward process gradually corrupts data 0x_0 into noise Tx_T through a Markov chain q(1:T∣0)=∏t=1Tq(t∣t−1)q(x_1:T _0)= _t=1^Tq(x_t _t-1). The reverse Markov process is parameterized by neural networks and progressively denoises the latent states from Tx_T to 0x_0. The reverse trajectory is modeled as pθ(0:T)=q(T)∏t=1Tpθ(t−1∣t).p_θ(x_0:T)=q(x_T) _t=1^Tp_θ(x_t-1 _t). (2) Diffusion process. We base our formulation on discrete diffusion spaces. Following D3PM [1], we consider one-hot discrete random variables with K categories, where 0,1,…,t∈0,1Kx_0,x_1,…,x_t∈\0,1\^K. The forward transition probabilities are encoded by matrices Qt∈ℝK×KQ_t ^K× K, where [Qt]ij=q(xt=j∣xt−1=i)[Q_t]_ij=q(x_t=j x_t-1=i). The transition distribution is defined as: q(t∣t−1)=Cat(t;=t−1Qt),q(x_t _t-1)=Cat(x_t;p=x_t-1Q_t), (3) where Cat(;)Cat(x;p) denotes a categorical distribution over the one-hot vector x with probabilities p. The marginal at timestep t follows directly from the Markov property: q(t∣0)=Cat(t;=0Q¯t),whereQ¯t=Q1Q2⋯Qt.q(x_t _0)=Cat(x_t;p=x_0 Q_t), Q_t=Q_1Q_2·s Q_t. (4) Reverse process. The reverse process initializes T∼π()x_T π(x) and generates 0x_0 by reversing the timesteps t=T,T−1,…,0t=T,T-1,…,0. Following prior work [16, 1], a neural network nθ(t)n_θ(x_t) predicts p~θ(~0∣t) p_θ( x_0 _t). The reverse transition is then computed as: pθ(t−1∣t)∝∑~0q(t−1,t∣~0)p~θ(~0∣t).p_θ(x_t-1 _t) _ x_0q(x_t-1,x_t x_0) p_θ( x_0 _t). (5) By Bayes’ theorem and the Markov property, the joint posterior admits a tractable form: q(t−1,t∣0)=q(t−1∣t,0)q(t−1∣0)q(t∣0)=Cat(t−1;=tQt⊤⊙0Q¯t−10Q¯tt⊤).q(x_t-1,x_t _0)= q(x_t-1 _t,x_0)q(x_t-1 _0)q(x_t _0)=Cat (x_t-1;p= x_tQ_t _0 Q_t-1x_0 Q_tx_t ). (6) Choice of transition matrix and conditioning. Different designs of QtQ_t induce distinct inductive biases [1]. A simple yet effective parameterization is Qt=αt+(1−αt)⊤Q_t= _tI+(1- _t)1m , where m denotes the marginal distribution and αt _t controls the noise scale. The cumulative transition admits a closed form Q¯t=α¯t+(1−α¯t)⊤ Q_t= α_tI+(1- α_t)1m , where α¯t=∏τ=1tατ α_t= _τ=1^t _τ. The noise schedule is governed by a cosine function [27]. Furthermore, to enable conditional generation, the reverse process learns a conditional data distribution pθ(t−1∣t,)p_θ(x_t-1 _t,c), where auxiliary condition c is injected into the denoiser to guide the reverse trajectory. Figure 2: Overview of the CoDiffGRN framework. The framework combines TF-ALL Subgraph Sampling (TASS) with a co-evolutionary discrete diffusion model. TASS enables efficient training of a multi-conditioned TF-Aware Graph Denoising Module. At inference, predictions from multiple subgraphs are aggregated via the reverse diffusion process. Abbr., Aggr.: Aggregation. 4.2 Joint Discrete Diffusion on Gene Regulatory Network As illustrated in Figure 2, we propose a joint discrete diffusion framework for GRNs. In contrast to prior works that decouple gene features from edge diffusion [47, 35], our framework unifies discretized gene expression and regulatory interactions into a single co-evolutionary state space. Cell-cluster discretization. To mitigate feature collapse on unseen genes, we transform continuous scRNA-seq profiles into biologically coherent discrete states. Prior studies show that discretized gene activation patterns capture stable functional structures across diverse biological settings [12, 5]. We cluster cells into k groups and binarize gene expression within each cluster via thresholding, producing a k-bit activation pattern per gene. Each binary pattern is mapped to a categorical variable and represented as a one-hot vector in ℝFDR^F_D, where FD=2kF_D=2^k. Joint modeling. We follow a joint diffusion formulation [23] to unify gene expression and regulatory edges. Edges are represented as ∈ℝN×N×2E ^N× N× 2 (indicating absence/presence), and node features as ∈ℝN×FDD ^N× F_D. The global graph token for each node i is constructed via concatenation: [i]=[i‖i1‖⋯∥iN]G[i]= [d_i\;\|\;e_i1\;\|\;·s\;\|\;e_iN ], where i∈ℝFDd_i ^F_D denotes the discretized node features, and ij∈ℝ2e_ij ^2 is the one-hot encoding of the edge from node i to node j. Design of the transition matrix. We formulate a joint transition matrix QG∈ℝFG×FGQ_G ^F_G× F_G, constructed from four component transition probability matrices: QDQ_D (gene → gene), QEQ_E (regulation → regulation), QEDQ_ED (regulation → gene), and QDEQ_DE (gene → regulation): QG=[QDN⊤⊗QDEN⊗QEDN×N⊗QE],Q_G= bmatrixQ_D&1_N Q_DE\\ 1_N Q_ED&1_N× N Q_E bmatrix, (7) where ⊗ denotes the Kronecker product, and N1_N, N×N1_N× N are all-ones vectors and matrices of appropriate dimensions. Following Vignac et al. [33], QDQ_D and QEQ_E are derived from the marginal distributions Dm_D and Em_E computed from the training data. For the cross-dependencies, EDm_ED and DEm_DE capture the conditional co-occurrence probabilities between gene expression types and regulation types. This explicitly ensures that edge transitions are biologically conditioned on the states of their endpoint genes. Forcing true D technique. During training, the reverse process approximates pθ(Gt−1∣Gt)∝∑G~0q(Gt−1,Gt∣G~0)p~θ(G~0∣Gt)p_θ(G_t-1 G_t) _ G_0q(G_t-1,G_t G_0) p_θ( G_0 G_t). To mitigate error accumulation and sharpen the learning signal for unknown edge structures, we apply a teacher-forcing-style technique that replaces the predicted node component ~0 D_0 within the generated G~0 G_0 with the ground-truth 0D_0. 4.3 Conditional Graph Denoising Architecture Conditioning mechanism. We incorporate multiple sources of conditioning information into the denoiser to guide the reverse process. Combining sinusoidal timestep embeddings it_i, continuous gene expression features ix_i, noisy discrete node states t,id_t,i, flattened incoming/outgoing edge features t,i:e_t,i:, and a boolean TF indicator iTFz_i^TF, the conditioned input feature for node i is constructed as: i(0)=[i‖i‖t,i‖t,i:‖iTF].h_i^(0)= [t_i\;\|\;x_i\;\|\;d_t,i\;\|\;e_t,i:\;\|\;z_i^TF ]. (8) Denoiser backbone. We parameterize the reverse process using a TF-aware graph denoiser built upon GATv2 [6]. Each layer updates node representations dynamically: i(l+1)=σ(∑j∈(i)αij(l)(l)j(l)),h_i^(l+1)=σ ( _j (i) _ij^(l)W^(l)h_j^(l) ), (9) where the attention coefficients αij(l) _ij^(l) are computed via: αij(l)=softmaxj(⊤LeakyReLU((l)[i(l)∥j(l)])). _ij^(l)=softmax_j (a \,LeakyReLU (W^(l)[h_i^(l)\|h_j^(l)] ) ). (10) The node representations are then projected into TF-specific (itfh^tf_i) and target-specific (jtgh^tg_j) embedding spaces. A score head yields a scalar compatibility score sij=fscore([itf∥jtg])s_ij=f_score([h^tf_i\|h^tg_j]). Binary edge logits are formed as [−sij,sij][-s_ij,\,s_ij] and added residually to the noisy edge input tE_t. Node logits are obtained via a linear projection of averaged embeddings with a residual connection to tD_t. Training objective. The denoiser is trained using a weighted combination of discrete edge reconstruction loss (ℒedgeL_edge) and a node-level auxiliary regression objective (ℒnodeL_node). Let e^i,j e_i,j denote the predicted probability of edge existence and ei,j∈0,1e_i,j∈\0,1\ the ground-truth. The total loss is defined as: ℒ=−∑i,j(ei,jloge^i,j+(1−ei,j)log(1−e^i,j))+λ1N∑i=1N‖^0,i−i‖22,L=- _i,j (e_i,j e_i,j+(1-e_i,j) (1- e_i,j) )+λ 1N _i=1^N \| X_0,i-X_i \|_2^2, (11) where the auxiliary MSE loss preserves continuous biological signals within the discretized space, and λ≪1λ 1 controls its strength. 4.4 TF-ALL Subgraph Sampling (TASS) Diffusion models typically require large training datasets, but GRN inference is inherently data-scarce since each dataset constitutes a single large graph. To address this, we propose TF-ALL Subgraph Sampling (TASS), which constructs a training set by sampling subgraphs from the input graph (Figure 1(c)). The model is trained on these subgraphs and, at inference time, repeatedly samples subgraphs, predicts edges, and aggregates results to reconstruct the full graph. This expands effective supervision by exposing diverse gene–topology patterns, improving generalization under limited data. Sampling. As shown in Figure 2, we adopt node-based sampling due to the lack of reliable edge priors, following uniform node sampling [22]. Exploiting the bipartite-like structure of GRNs (TFs → targets), each subgraph of size k is formed by sampling a fraction t (TF-to-ALL ratio) of nodes from the TF set, and the remainder from all genes, with resampling for TF duplicates. Let p=k/np=k/n. Following theoretical bounds for graph recovery [22], we set the number of subgraphs to m=⌈p−2lognlog(1/δ)⌉m= p^-2 n (1/δ) . In GRNs, this bipartite sampling strategy inherently yields higher edge coverage due to the sparse, hierarchical nature of regulatory interactions. Reconstruction. At inference, we apply the same sampling strategy, generate predictions across subgraphs, and map them back to the global edge indices. Outputs are converted to probabilities via softmax, and overlapping predictions are averaged to form the final consensus. Uncovered edges are assigned a default baseline prior. 5 Experiments 5.1 Experimental Setup Evaluation protocols. We primarily evaluate performance using Hits@K and MRR on Inductive Gene-Holdout Split introduced in Section 3. To reflect practical experimental settings, we report Hits@10 and Hits@50, corresponding to low- and high-throughput validation regimes. Baseline models. We compare against a comprehensive set of baselines, including two traditional methods (GENIE3 [17], GRNBoost2 [26]), two feature-based deep learning models (CNNC [43], DeepSEM [30]), four graph-based models (GNNLink [25], GENELink+ [46], GRNFormer [15], GCLink [42]), and two diffusion-based methods (RegDiffusion [47], DigNet [35]). Model setups. For fair comparison, all baselines are implemented following their original configurations. Our model adopts the TF-Aware Graph Denoiser, built upon a multi-layer GATv2 architecture with three attention heads per layer and progressively reduced hidden dimensions (128, 64, 64, 32), producing a 16-dimensional output embedding. We use a LeakyReLU with negative slope 0.2 and concatenate multi-head outputs. Score head fscoref_score is implemented as a two-layer MLP. Training setups. All baselines are trained according to their original protocols. Our model is optimized using AdamW with a learning rate of 1×10−41× 10^-4, weight decay 1×10−31× 10^-3, batch size 64, and trained for 500 epochs. We use a discrete state space with k=4k=4 clusters (yielding 16 classes) and 500 diffusion steps. Subgraphs are sampled with 100 nodes each. All experiments are conducted on a local server with 100 CPU cores and 4 NVIDIA RTX 4090 GPUs (24 GB), and are run over three random seeds (0, 1, and 2), with both mean and standard deviation reported. Table 2: Benchmark results on BEELINE-KGC (Specific network). We report mean over runs. The best, second-best, and third-best are denoted by bold, underline, and italic, respectively. Δ indicates change relative to the baseline: green up, red down. Leading zeros are omitted. Abbr., Trad.: Traditional methods, DL.: Feature-based deep learning models, GNN.: GNN-based models, Diff.: Diffusion-based models. Model hESC hHEP mDC mESC mHSC-E mHSC-GM mHSC-L H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR Setting: TFs+500 Trad. GENIE3 .000 .002 .003 .001 .004 .007 .000 .009 .003 .004 .013 .007 .001 .005 .004 .003 .006 .005 .004 .004 .005 GRNBoost2 .000 .002 .003 .000 .004 .006 .000 .009 .003 .001 .013 .006 .000 .005 .004 .003 .006 .006 .000 .004 .004 DL. CNNC .006 .006 .005 .009 .009 .013 .009 .009 .004 .026 .036 .019 .000 .005 .004 .000 .000 .005 .000 .000 .004 DeepSEM .000 .004 .003 .001 .002 .007 .000 .000 .003 .005 .017 .007 .000 .002 .004 .000 .003 .005 .000 .000 .004 GNN. GNNLink .026 .110 .015 .059 .200 .045 .012 .037 .004 .037 .157 .020 .049 .141 .032 .044 .119 .037 .055 .152 .024 GENELink+ .052 .218 .030 .141 .385 .100 .025 .077 .009 .082 .343 .050 .124 .355 .073 .115 .290 .073 .122 .304 .067 GRNFormer .036 .150 .020 .093 .248 .072 .017 .047 .006 .061 .250 .030 .076 .189 .052 .069 .172 .044 .077 .216 .038 GCLink .033 .131 .016 .079 .224 .060 .014 .040 .005 .046 .202 .024 .066 .180 .043 .057 .148 .040 .064 .170 .030 Diff. RegDiffusion .033 .151 .021 .077 .230 .058 .023 .071 .009 .051 .195 .033 .133 .362 .074 .094 .249 .065 .066 .162 .036 DigNet .063 .242 .029 .138 .355 .097 .019 .055 .006 .051 .236 .035 .168 .361 .088 .079 .216 .056 .139 .314 .085 CoDiffGRN .079 .264 .034 .173 .499 .126 .028 .122 .014 .123 .412 .060 .175 .407 .099 .167 .463 .108 .162 .361 .111 Δ(%) (\%) +24.9% +9.1% +13.0% +22.1% +29.6% +25.1% +11.3% +57.1% +52.7% +49.3% +20.1% +20.2% +4.1% +12.6% +11.8% +44.7% +59.5% +48.9% +16.6% +14.8% +30.4% Setting: TFs+1000 Trad. GENIE3 .000 .000 .002 .001 .003 .005 .000 .007 .002 .001 .004 .004 .001 .280 .018 .001 .643 .017 .002 .006 .007 GRNBoost2 .000 .000 .002 .000 .003 .004 .000 .007 .002 .001 .004 .004 .000 .280 .017 .000 .643 .017 .000 .006 .006 DL. CNNC .005 .015 .004 .006 .008 .008 .000 .000 .002 .011 .022 .009 .017 .295 .030 .010 .652 .024 .000 .003 .006 DeepSEM .000 .002 .002 .001 .003 .005 .000 .000 .002 .003 .008 .004 .001 .287 .018 .003 .647 .018 .000 .006 .006 GNN. GNNLink .010 .037 .006 .029 .142 .016 .007 .028 .006 .038 .081 .029 .050 .160 .020 .112 .301 .052 .086 .303 .035 GENELink+ .019 .070 .010 .067 .340 .038 .014 .061 .009 .086 .211 .058 .073 .301 .046 .249 .648 .101 .166 .531 .073 GRNFormer .015 .051 .007 .045 .216 .024 .008 .034 .007 .058 .114 .040 .058 .196 .031 .144 .367 .066 .105 .420 .047 GCLink .012 .043 .006 .035 .163 .021 .007 .031 .006 .044 .093 .035 .055 .186 .025 .129 .333 .062 .098 .361 .038 Diff. RegDiffusion .022 .089 .014 .068 .347 .043 .009 .046 .007 .083 .191 .055 .060 .230 .041 .277 .602 .106 .100 .306 .038 DigNet .021 .071 .011 .067 .335 .044 .013 .065 .009 .117 .266 .077 .063 .249 .039 .217 .609 .089 .135 .441 .055 CoDiffGRN .027 .108 .015 .079 .398 .053 .016 .069 .010 .139 .299 .093 .091 .429 .051 .334 .711 .140 .276 .587 .117 Δ(%) (\%) +23.5% +22.3% +12.3% +15.3% +14.7% +20.2% +13.9% +5.8% +14.1% +19.4% +12.1% +21.1% +25.4% +42.5% +11.7% +20.7% +9.0% +31.4% +66.1% +10.5% +61.6% 5.2 Results We report the performance of representative models on the Specific network under the proposed BEELINE-KGC in Table 2. Based on these results, we highlight the following key findings: • The proposed BEELINE-KGC benchmark is substantially more challenging than existing protocols. Traditional and feature-based methods largely collapse under inductive gene holdout, yielding near-zero Hits@10 and MRR. Even several GNN-based models degrade substantially, indicating that BEELINE-KGC reveals inductive generalization failures obscured by transductive benchmarks and global metrics. • GNN-based models show limited inductive generalization. While GENELink+ remains competitive, most GNN methods suffer from low performance. Their heavy reliance on predefined prior graphs leads to disconnected components for unseen genes, causing error accumulation during message passing and widening the gap against diffusion-based approaches. • Diffusion-based models offer partial robustness yet remain insufficient. According to Table 2, diffusion-based methods generally outperform GNN baselines and match GENELink+ in some cases, but still fall behind CoDiffGRN. This suggests that existing diffusion models alleviate reliance on explicit graph structure, yet remain vulnerable to feature collapse on unseen genes. Moreover, their single-graph optimization paradigm appears suboptimal for diffusion. • CoDiffGRN achieves consistent SOTA performance. CoDiffGRN ranks first across all cell types and metrics in both TFs+500 and TFs+1000 settings, with relative gains up to +184.8% (and +24.5% on average) over the strongest baselines on the most competitive Specific network. These results validate the effectiveness of our co-evolutionary discrete diffusion framework with gene-state-dependent transitions and the TASS strategy for addressing inductive GRN inference and top-K ranking quality. Figure 3: Ablation study results on BEELINE-KGC (Specific network). We report mean performance of Hits@10 and Hits@50 over runs, under the TFs+500 (top) and TFs+1000 (bottom) settings. Abbr., w/o: without. 5.3 Ablation Study To verify the effectiveness of our designs, we conduct extensive ablation experiments on BEELINE-KGC, focusing on the joint diffusion mechanism and the TASS augmentation strategy. The variants are: (i) Edge Only: a variant that only using edge diffusion; (i) w/o TASS: the model trained on the original graph without TASS; and (i) RandSS: replacing TASS with uniform random subgraph sampling. Effect of joint modeling. As shown in Figure 3, removing joint modeling (Edge Only) consistently degrades performance across datasets and metrics. This indicates that decoupling gene expression states from regulatory edges limits the model’s ability to capture conditional regulatory interactions and reduces top-K ranking quality, especially when generalizing to unseen genes in inductive settings. Effect of subgraph sampling. Removing TASS leads to a substantial performance drop across all experimental settings. In contrast, uniform random sampling (RandSS) partially recovers the performance and achieves superior results in a few isolated cases, but remains consistently less effective overall. These results highlight the importance of informative subgraph construction and demonstrate that the task-aware sampling strategy of TASS provides more reliable improvements than generic random sampling under inductive and data-scarce GRN inference settings. 6 Conclusion In this work, we address the gap between computational GRN inference and real-world biological discovery. We introduce BEELINE-KGC, an application-aligned framework that enforces an inductive, KGC-based evaluation protocol to reflect the top-K constraints of wet-lab validation. To overcome the exposed gene-level generalization gap, we propose CoDiffGRN, the first joint node-edge discrete diffusion model for GRNs. By capturing the co-evolution of gene expression and network topology, and further enhanced by our TASS strategy, CoDiffGRN achieves robust inductive generalization. Together, our benchmark and model establish a unified paradigm for reliable and budget-aware biological discovery in systems biology. Limitations and Future Work. Despite its performance, our approach has several limitations that suggest directions for future research. The current evaluation is restricted to established scRNA-seq systems, which may not fully represent the complexity of rare tissues or cross-species regulatory dynamics. We plan to extend our benchmark to include more diverse biological contexts and validate the model’s transferability across different species. Moreover, the joint diffusion process and TASS strategy introduce non-trivial computational overhead during training and inference. Future work will focus on developing optimized sampling kernels and parallelized subgraph processing to enhance scalability. References [1] J. Austin, D. D. Johnson, J. Ho, D. Tarlow, and R. Van Den Berg (2021) Structured denoising diffusion models in discrete state-spaces. Advances in neural information processing systems 34, p. 17981–17993. Cited by: §1, §2, §4.1, §4.1, §4.1. [2] E. Bastakis, J. Gerke, S. Özkan, R. Harting, T. Lienard, C. Sasse, E. S. Xylakis, M. Aden, A. Strohdiek, G. Heinrich, et al. (2025) Molecular circuit between aspergillus nidulans transcription factors msna and velb to coordinate fungal stress and developmental responses. PLoS Genetics 21 (7), p. e1011578. Cited by: §3.2. [3] E. Becht, L. McInnes, J. Healy, C. Dutertre, I. W. Kwok, L. G. Ng, F. Ginhoux, and E. W. Newell (2019) Dimensionality reduction for visualizing single-cell data using umap. Nature biotechnology 37 (1), p. 38–44. Cited by: Figure 4, §B.2, §B.2. [4] A. Bordes, N. Usunier, A. Garcia-Duran, J. Weston, and O. Yakhnenko (2013) Translating embeddings for modeling multi-relational data. Advances in neural information processing systems 26. Cited by: §1. [5] G. A. Bouland, A. Mahfouz, and M. J. Reinders (2021) Differential analysis of binarized single-cell rna sequencing data captures biological variation. NAR genomics and bioinformatics 3 (4), p. lqab118. Cited by: §4.2. [6] S. Brody, U. Alon, and E. Yahav (2021) How attentive are graph attention networks?. arXiv preprint arXiv:2105.14491. Cited by: §4.3. [7] G. Chen and Z. Liu (2022) Graph attention network for link prediction of gene regulations from single-cell rna-sequencing data. Bioinformatics 38 (19), p. 4522–4529. Cited by: §1, §2. [8] C. Chronis, P. Fiziev, B. Papp, S. Butz, G. Bonora, S. Sabri, J. Ernst, and K. Plath (2017) Cooperative binding of transcription factors orchestrates reprogramming. Cell 168 (3), p. 442–459. Cited by: §1. [9] J. W. Ding, C. N. Kim, M. S. Ostrowski, Y. Abeykoon, B. J. Pavlovic, J. L. Wallace, N. K. Schaefer, T. J. Nowakowski, and A. A. Pollen (2026) Dissecting gene regulatory networks governing human cortical cell fate. Nature, p. 1–11. Cited by: §1, §3.2. [10] L. Garcia-Alonso, C. H. Holland, M. M. Ibrahim, D. Turei, and J. Saez-Rodriguez (2019) Benchmark and integration of resources for the estimation of human transcription factor activities. Genome research 29 (8), p. 1363–1375. Cited by: Appendix D. [11] R. B. Gilmore, Y. Liu, C. E. Stoddard, M. S. Chung, G. G. Carmichael, and J. Cotney (2024) Identifying key underlying regulatory networks and predicting targets of orphan c/d box snord116 snornas in prader–willi syndrome. Nucleic Acids Research 52 (22), p. 13757–13774. Cited by: §1, §3.2. [12] D. Grün, A. Lyubimova, L. Kester, K. Wiebrands, O. Basak, N. Sasaki, H. Clevers, and A. Van Oudenaarden (2015) Single-cell messenger rna sequencing reveals rare intestinal cell types. Nature 525 (7568), p. 251–255. Cited by: §4.2. [13] H. Han, J. Cho, S. Lee, A. Yun, H. Kim, D. Bae, S. Yang, C. Y. Kim, M. Lee, E. Kim, et al. (2018) TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic acids research 46 (D1), p. D380–D386. Cited by: Appendix D. [14] M. Hecker, S. Lambeck, S. Toepfer, E. Van Someren, and R. Guthke (2009) Gene regulatory network inference: data integration in dynamic models—a review. Biosystems 96 (1), p. 86–103. Cited by: §1. [15] A. Hegde and J. Cheng (2026) GRNFormer: accurate gene regulatory network inference using graph transformer. Bioinformatics 42 (4), p. btag144. Cited by: §1, §2, §5.1. [16] J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. Advances in neural information processing systems 33, p. 6840–6851. Cited by: §2, §4.1, §4.1. [17] V. A. Huynh-Thu, A. Irrthum, L. Wehenkel, and P. Geurts (2010) Inferring regulatory networks from expression data using tree-based methods. PloS one 5 (9), p. e12776. Cited by: §2, §5.1. [18] I. Isewon, M. Kanonte, E. Alagbe, and J. Oyelade (2025) Benchmarking frameworks for gene regulatory network inference methods: a review. NIPES JSTR SPECIAL ISSUE 7 (1), p. 2224–2228. Cited by: §2. [19] M. S. Kagda, B. Lam, C. Litton, C. Small, C. A. Sloan, E. Spragins, F. Tanaka, I. Whaling, I. Gabdank, I. Youngworth, et al. (2025) Data navigation on the encode portal. Nature Communications 16 (1), p. 9592. Cited by: Appendix D. [20] P. V. Kharchenko, L. Silberstein, and D. T. Scadden (2014) Bayesian approach to single-cell differential expression analysis. Nature methods 11 (7), p. 740–742. Cited by: §1. [21] A. A. Kolodziejczyk, J. K. Kim, V. Svensson, J. C. Marioni, and S. A. Teichmann (2015) The technology and biology of single-cell rna sequencing. Molecular cell 58 (4), p. 610–620. Cited by: §1. [22] S. Limnios, P. Selvaraj, M. Cucuringu, C. Maple, G. Reinert, and A. Elliott (2023) Sagess: sampling graph denoising diffusion model for scalable graph generation. arXiv preprint arXiv:2306.16827. Cited by: §4.4. [23] G. Liu, J. Xu, T. Luo, and M. Jiang (2024) Graph diffusion transformers for multi-conditional molecular generation. Advances in Neural Information Processing Systems 37, p. 8065–8092. Cited by: §4.2. [24] Z. Liu, C. Wu, H. Miao, and H. Wu (2015) RegNetwork: an integrated database of transcriptional and post-transcriptional regulatory networks in human and mouse. Database 2015, p. bav095. Cited by: Appendix D. [25] G. Mao, Z. Pang, K. Zuo, Q. Wang, X. Pei, X. Chen, and J. Liu (2023) Predicting gene regulatory links from single-cell rna-seq data using graph neural networks. Briefings in Bioinformatics 24 (6), p. bbad414. Cited by: §1, §2, §5.1. [26] T. Moerman, S. Aibar Santos, C. Bravo González-Blas, J. Simm, Y. Moreau, J. Aerts, and S. Aerts (2019) GRNBoost2 and arboreto: efficient and scalable inference of gene regulatory networks. Bioinformatics 35 (12), p. 2159–2161. Cited by: §2, §2, §5.1. [27] A. Q. Nichol and P. Dhariwal (2021) Improved denoising diffusion probabilistic models. In International conference on machine learning, p. 8162–8171. Cited by: §4.1. [28] W. Niu, Z. J. Lu, M. Zhong, M. Sarov, J. I. Murray, C. M. Brdlik, J. Janette, C. Chen, P. Alves, E. Preston, et al. (2011) Diverse transcription factor binding features revealed by genome-wide chip-seq in c. elegans. Genome research 21 (2), p. 245–254. Cited by: §3.2. [29] A. Pratapa, A. P. Jalihal, J. N. Law, A. Bharadwaj, and a. M. Murali (2020) Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature methods 17 (2), p. 147–154. Cited by: Appendix D, §1, §2, §3. [30] H. Shu, J. Zhou, Q. Lian, H. Li, D. Zhao, J. Zeng, and J. Ma (2021) Modeling gene regulatory networks using neural network architectures. Nature Computational Science 1 (7), p. 491–501. Cited by: §2, §5.1. [31] S. Subramanian, J. A. Thoms, Y. Huang, P. Cornejo-Páramo, F. C. Koch, S. Jacquelin, S. Shen, E. Song, S. Joshi, C. Brownlee, et al. (2023) Genome-wide transcription factor–binding maps reveal cell-specific changes in the regulatory architecture of human hspcs. Blood 142 (17), p. 1448–1462. Cited by: §3.2. [32] D. Szklarczyk, R. Kirsch, M. Koutrouli, K. Nastou, F. Mehryary, R. Hachilif, A. L. Gable, T. Fang, N. T. Doncheva, S. Pyysalo, et al. (2023) The string database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic acids research 51 (D1), p. D638–D646. Cited by: Appendix D. [33] C. Vignac, I. Krawczuk, A. Siraudin, B. Wang, V. Cevher, and P. Frossard (2022) Digress: discrete denoising diffusion for graph generation. arXiv preprint arXiv:2209.14734. Cited by: §4.2. [34] A. Wagner, A. Regev, and N. Yosef (2016) Revealing the vectors of cellular identity with single-cell genomics. Nature biotechnology 34 (11), p. 1145–1160. Cited by: §1. [35] C. Wang and Z. Liu (2025) Diffusion-based generation of gene regulatory networks from scrna-seq data with dignet. Genome Research 35 (2), p. 340–354. Cited by: §2, §4.2, §5.1. [36] J. Wang, Y. Chen, and Q. Zou (2023) Inferring gene regulatory network from single-cell transcriptomes with graph autoencoder model. PLoS Genetics 19 (9), p. e1010942. Cited by: §1, §2. [37] Y. Wang, S. E. Peterson, and J. F. Loring (2014) Protein post-translational modifications and regulation of pluripotency in human stem cells. Cell research 24 (2), p. 143–160. Cited by: §1. [38] C. R. Wiggers, B. Yüzügüldü, N. G. Tadros, T. B. Heavican-Foral, E. Y. Cho, Z. C. Eisenbies, M. Ozdemir, S. B. Kulp, Y. Chae, A. Gutierrez, et al. (2025) Genome-wide crispr screen identifies irf1 and tfap4 as transcriptional regulators of galectin-9 in t cell acute lymphoblastic leukemia. Science advances 11 (12), p. eads8351. Cited by: §3.2. [39] P. J. Wittkopp and G. Kalay (2012) Cis-regulatory elements: molecular mechanisms and evolutionary processes underlying divergence. Nature Reviews Genetics 13 (1), p. 59–69. Cited by: §1. [40] H. Xu, C. Baroukh, R. Dannenfelser, E. Y. Chen, C. M. Tan, Y. Kou, Y. E. Kim, I. R. Lemischka, and A. Ma’ayan (2013) ESCAPE: database for integrating high-content published data collected from human and mouse embryonic stem cells. Database 2013, p. bat045. Cited by: Appendix D. [41] S. Xu, N. Yu, D. Zhang, and C. Wang (2025) Attention-guided probabilistic diffusion model for generating cell-type-specific gene regulatory networks from gene expression profiles. Genes 16 (11), p. 1255. Cited by: §2. [42] W. Yu, Z. Lin, M. Lan, and L. Ou-Yang (2025) GCLink: a graph contrastive link prediction framework for gene regulatory network inference. Bioinformatics 41 (3), p. btaf074. Cited by: §1, §2, §5.1. [43] Y. Yuan and Z. Bar-Joseph (2019) Deep learning for inferring gene relationships from single-cell expression data. Proceedings of the National Academy of Sciences 116 (52), p. 27151–27158. Cited by: §2, §2, §5.1. [44] J. Zhang, B. Liu, Y. Song, Y. Chen, J. Fu, J. Liu, T. Ma, Z. Xi, and H. Liu (2021) Genome-wide (chip-seq) identification of target genes regulated by wrky33 during submergence stress in arabidopsis. BMC genomic data 22 (1), p. 16. Cited by: §3.2. [45] Q. Zhang, W. Liu, H. Zhang, G. Xie, Y. Miao, M. Xia, and A. Guo (2020) HTFtarget: a comprehensive database for regulations of human transcription factors and their targets. Genomics, proteomics & bioinformatics 18 (2), p. 120–128. Cited by: §3.2. [46] W. Zhang, B. Shao, W. Li, W. Guo, J. Lyu, G. Chen, C. Wang, and Z. Liu (2025) Inferring cell-type-specific gene regulatory network from cellular transcriptomics data with genelink+. Briefings in Bioinformatics 26 (4), p. bbaf359. Cited by: §1, §2, §2, §5.1. [47] H. Zhu and D. Slonim (2024) From noise to knowledge: diffusion probabilistic model-based neural inference of gene regulatory networks. Journal of Computational Biology 31 (11), p. 1087–1103. Cited by: §2, §4.2, §5.1. [48] Z. Zou, T. Ohta, and S. Oki (2024) ChIP-atlas 3.0: a data-mining suite to explore chromosome architecture together with large-scale regulome data. Nucleic Acids Research 52 (W1), p. W45–W53. Cited by: Appendix D. Appendix A Results on Previous Transductive Setting Table 3: Results on transductive setting (Specific network). Caption settings are the same as Table 2. Model hESC hHEP mDC mESC mHSC-E mHSC-GM mHSC-L H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR H@10 H@50 MRR Setting: TFs+500 Trad. GENIE3 .012 .012 .008 .023 .023 .020 .009 .018 .005 .034 .036 .024 .008 .008 .010 .010 .010 .014 .013 .013 .014 GRNBoost2 .006 .006 .007 .001 .001 .007 .000 .000 .003 .005 .005 .010 .002 .002 .005 .000 .000 .005 .002 .002 .006 DL. CNNC .028 .140 .016 .047 .198 .028 .011 .109 .013 .058 .237 .033 .021 .103 .016 .032 .149 .018 .016 .095 .015 DeepSEM .024 .071 .012 .054 .059 .030 .045 .073 .022 .094 .149 .042 .078 .086 .049 .074 .075 .053 .065 .067 .036 GNN. GNNLink .015 .051 .013 .001 .006 .007 .027 .130 .015 .000 .020 .006 .000 .002 .003 .010 .022 .006 .017 .025 .015 GENELink+ .040 .169 .022 .137 .377 .057 .016 .071 .010 .103 .358 .053 .329 .598 .167 .048 .213 .029 .130 .330 .057 GRNFormer .039 .187 .026 .106 .323 .049 .033 .125 .015 .118 .382 .058 .034 .119 .020 .040 .209 .024 .035 .137 .024 GCLink .098 .296 .048 .146 .394 .072 .049 .141 .031 .188 .461 .092 – – – – – – .058 .213 .032 Diff. RegDiffusion .031 .148 .017 .035 .140 .019 .022 .060 .018 .053 .176 .025 .040 .208 .024 .087 .311 .044 .081 .257 .036 DigNet .043 .168 .023 .055 .198 .029 .049 .120 .019 .064 .281 .031 .049 .184 .025 .132 .434 .066 .048 .235 .028 CoDiffGRN .102 .318 .056 .247 .539 .128 .049 .179 .033 .259 .524 .137 .470 .769 .259 .435 .755 .227 .349 .669 .189 Δ(%) (\%) +4.1% +7.3% +15.9% +68.7% +36.9% +76.3% +0.0% +26.9% +4.8% +37.6% +13.9% +48.7% +42.8% +28.7% +54.9% +228.7% +74.0% +244.1% +167.6% +102.7% +231.2% Setting: TFs+1000 Trad. GENIE3 .004 .004 .004 .010 .010 .010 .014 .014 .007 .021 .022 .015 .003 .283 .020 .006 .646 .022 .011 .011 .016 GRNBoost2 .000 .000 .002 .001 .001 .006 .000 .000 .002 .005 .005 .008 .001 .281 .017 .002 .643 .018 .000 .000 .006 DL. CNNC .014 .096 .013 .032 .151 .018 .004 .042 .006 .043 .162 .023 .153 .520 .066 .129 .464 .056 .022 .142 .017 DeepSEM .024 .046 .008 .044 .049 .023 .014 .035 .008 .064 .112 .032 .046 .322 .043 .043 .668 .045 .043 .043 .037 GNN. GNNLink .004 .021 .005 .000 .005 .005 .014 .046 .007 .000 .000 .003 .000 .516 .019 .005 .340 .016 .026 .059 .020 GENELink+ .041 .144 .022 .085 .295 .042 .025 .099 .010 .130 .356 .066 .445 .698 .251 .390 .670 .225 .084 .307 .036 GRNFormer .029 .112 .019 .062 .254 .035 .021 .053 .010 .062 .268 .033 .449 .660 .258 .396 .676 .207 .039 .134 .024 GCLink .042 .203 .024 .091 .274 .055 .042 .141 .021 .126 .351 .064 – – – – – – .103 .319 .044 Diff. RegDiffusion .022 .080 .013 .024 .089 .013 .021 .057 .012 .034 .119 .017 .064 .354 .036 .090 .334 .041 .154 .429 .082 DigNet .027 .107 .013 .026 .116 .016 .042 .064 .008 .042 .191 .023 .063 .346 .038 .055 .288 .034 .121 .407 .060 CoDiffGRN .079 .269 .035 .203 .439 .115 .053 .159 .025 .245 .502 .128 .464 .708 .267 .410 .702 .235 .311 .609 .144 Δ(%) (\%) +88.6% +32.7% +47.6% +122.8% +49.1% +110.3% +25.0% +12.5% +20.3% +88.9% +41.0% +94.1% +3.3% +1.5% +3.5% +3.5% +3.8% +4.4% +101.1% +42.0% +76.2% Table 4: Results of AUPRC on transductive setting (Specific network). Model TFs+500 TFs+1000 hESC hHEP mDC mESC mHSC-E mHSC-GM mHSC-L hESC hHEP mDC mESC mHSC-E mHSC-GM mHSC-L Trad. GENIE3 .156 .395 .055 .314 .566 .532 .525 .158 .384 .056 .312 .547 .532 .484 GRNBoost2 .154 .387 .069 .327 .578 .522 .597 .150 .380 .059 .324 .547 .527 .481 DL. CNNC .255 .469 .066 .484 .746 .682 .648 .271 .499 .058 .597 .774 .737 .565 DeepSEM .196 .461 .054 .314 .567 .520 .537 .192 .419 .056 .320 .565 .543 .523 GNN. GNNLink .521 .751 .251 .766 .881 .898 .855 .515 .787 .216 .785 .930 .931 .864 GENELink+ .519 .813 .153 .796 .943 .937 .870 .521 .817 .122 .776 .954 .951 .853 GRNFormer .503 .713 .121 .802 .874 .869 .764 .522 .744 .125 .797 .917 .915 .743 GCLink .465 .707 .207 .739 .856 .849 .814 .476 .740 .159 .722 .898 .870 .840 Diff. RegDiffusion .447 .808 .117 .761 .932 .925 .863 .483 .814 .132 .784 .944 .944 .859 DigNet .495 .799 .112 .775 .717 .703 .871 .527 .813 .139 .798 .625 .612 .869 CoDiffGRN .582 .817 .155 .836 .923 .909 .840 .598 .822 .154 .856 .931 .926 .839 Δ(%) (\%) +11.8% +0.5% -38.4% +4.3% -2.0% -3.0% -3.6% +13.3% +0.5% -28.4% +7.2% -2.4% -2.6% -3.5% Appendix B Cell-cluster Discretization This section provides additional details and supporting evidence for the proposed Cell-cluster-based Discretization strategy. We first present the algorithmic procedure, followed by visualizations that demonstrate the biological coherence of the induced discrete gene state space. B.1 Algorithm Algorithm 1 Cell-cluster Discretization 1:Gene expression matrix ∈ℝN×FXX ^N× F_X, number of clusters k 2:Discrete one-hot representations ∈ℝN×FDD ^N× F_D, where FD=2kF_D=2^k 3:Construct X_vis using only visible data 4:Apply k-means on cell-wise expression profiles of X_vis to obtain clusters 1,…,k\C_1,…,C_k\ 5:For each cluster jC_j, compute τj=1|j|⋅N∑c∈j∑g=1Ng,c _j= 1|C_j|· N _c _j _g=1^NX_g,c 6:for each gene g do 7: for each cluster j do 8: Compute mean expression μg,j=1|j|∑c∈jg,c _g,j= 1|C_j| _c _jX_g,c 9: Set binary activation bg,j=(μg,j≥τj)b_g,j=I( _g,j≥ _j) 10: end for 11: Form binary vector g∈0,1kb_g∈\0,1\^k 12: Map gb_g to category index cg∈0,…,FD−1c_g∈\0,…,F_D-1\ 13: Convert cgc_g to one-hot vector ∈ℝFDd ^F_D 14:end for 15:return =g=1ND=\d\_g=1^N B.2 Design Motivation and Biological Evidence We provide the design motivation and biological evidence supporting the proposed Cell-cluster-based Discretization. The goal of this design is to construct a biologically coherent discrete state space that enables stable inductive diffusion over unseen genes. Figure 4: Cell-state structure induced by clustering. Visualization of cell clusters obtained from scRNA-seq profiles using UMAP [3]. Each cluster corresponds to a distinct cellular state, supporting the use of cell-level structure as the basis for discretization. Motivation. A natural alternative is to directly cluster gene expression profiles using standard unsupervised methods (e.g., k-means). However, such approaches may produce clusters dominated by statistical artifacts such as sparsity patterns or zero inflation in scRNA-seq data, rather than biologically meaningful structure. In contrast, scRNA-seq measurements are inherently organized by cellular contexts (e.g., conditions or states), suggesting that discretization should be grounded in cell-level structure to preserve functional coherence. Figure 5: Gene activation consistency across cell clusters. Representative genes show stable and coherent activation patterns within each cell cluster, indicating that discretization preserves biologically meaningful gene activity structure. Design Principle. Our method explicitly leverages this structure by first clustering cells into k biologically meaningful groups, and then constructing gene-level discrete representations based on expression consistency across these cell states. This ensures that discretization reflects stable functional activity rather than superficial similarity in gene expression vectors. Figure Organization and Evidence. We use UMAP [3] to visualize the resulting cell clusters obtained from scRNA-seq data in Figure 4. The clusters correspond to distinct cellular states, indicating that the clustering captures meaningful biological organization rather than arbitrary partitions. Figure 5 presents representative genes across different cell clusters. We observe that genes exhibit consistent activation patterns within clusters, suggesting that discretization captures stable gene-level functional behavior conditioned on cellular context. This is further proved by Figure 6, which visualizes the induced gene representation manifold in the discrete state space using UMAP [3]. The resulting structure exhibits smooth organization and local coherence, suggesting that the discretization preserves biologically meaningful continuity while providing a structured state space suitable for diffusion modeling. Figure 6: Discrete gene manifold structure. The induced gene representation space exhibits smooth and coherent structure, indicating that the discretization preserves biologically meaningful organization while enabling stable diffusion over discrete states. Together, these results support the role of Cell-cluster-based Discretization as a biologically grounded representation scheme that enables stable inductive generalization under sparse and heterogeneous scRNA-seq data. Appendix C TF-ALL Subgraph Sampling (TASS) In this section, we present the complete pipeline of TF-ALL Subgraph Sampling (TASS), which consists of a stochastic subgraph sampling strategy followed by a consensus-based reconstruction procedure for global edge inference. C.1 Sampling Algorithm 2 TF-ALL Subgraph Sampling (TASS) 1:Graph G=(V,E)G=(V,E), TF set VTFV_TF, gene set VGV_G, subgraph size k 2:Subgraph collection =G1,…,GmS=\G_1,…,G_m\ 3:Compute n←|V|n←|V|, p←k/np← k/n, TF-ALL ratio t←|VTF|/|VG|t←|V_TF|/|V_G| 4:Calculate theoretical bound mbase←⌈p−2lognlog(1/δ)⌉m_base← p^-2 n (1/δ) 5:Set m=100⋅⌈mbase/100⌉m=100· m_base/100 ⊳ Round up to nearest multiple of 100 6:Initialize ←∅S← 7:for i=1i=1 to m do 8: Sample TF nodes: VTF(i)∼Uniform(VTF,⌊|VTF|⋅t⌋)V^(i)_TF (V_TF, |V_TF|· t ) 9: Sample gene nodes: VG(i)∼Uniform(VG,k−|VTF(i)|)V^(i)_G (V_G,k-|V^(i)_TF|) 10: Duplicates Vdup←VTF(i)∩VG(i)V_dup← V^(i)_TF∩ V^(i)_G, ndup←|Vdup|n_dup←|V_dup| 11: while ndup≠0n_dup≠ 0 do ⊳ Resolve duplicates 12: Resample a new gene g∼Uniform(VG)g (V_G) 13: if g∉VTF(i)∪VG(i)g ∈ V^(i)_TF∪ V^(i)_G then 14: Add g into the VG(i)V^(i)_G, ndup←ndup−1n_dup← n_dup-1 15: end if 16: end while 17: Construct induced subgraph: Gi=G[VTF(i)∪VG(i)]G_i=G[V^(i)_TF∪ V^(i)_G] 18: ←∪GiS ∪\G_i\ 19:end for 20:return S C.2 Reconstruction Algorithm 3 TF-ALL Subgraph Sampling (TASS): Reconstruction 1:Global graph G=(V,E)G=(V,E), sampled node sets =Vii=1mV=\V_i\_i=1^m, trained model Φ , placeholder λ 2:Global edge score matrix ∈ℝ|V|×|V|S ^|V|×|V| 3:Initialize accumulation matrix ←A 0, count matrix ←C 0 4:for each sampled node set Vi∈V_i do 5: Construct induced subgraph Gi=G[Vi]G_i=G[V_i] 6: Predict raw edge logits: i=Φ(Gi)E_i= (G_i) 7: Normalize edge scores: i=Softmax(i)P_i=Softmax(E_i) ⊳ within-subgraph normalization 8: Map local edges to global index space: (U,V)←π(Vi×Vi)(U,V)←π(V_i× V_i) 9: Accumulate predictions: U,V←U,V+i,U,V←U,V+1A_U,V _U,V+P_i,C_U,V _U,V+1 10:end for 11:Consensus aggregation: 12:for each potential edge (u,v)∈V×V(u,v)∈ V× V do 13: if u,v>0C_u,v>0 then 14: u,v=u,v/u,vS_u,v=A_u,v/C_u,v ⊳ mean across subgraphs 15: else 16: u,v=λS_u,v=λ ⊳ unobserved edge prior 17: end if 18:end for 19:return S Appendix D Dataset Information of BEELINE-KGC The original BEELINE benchmark is constructed from five experimental single-cell RNA-seq datasets covering seven cell types (hESC, hHEP, mDC, mESC, mHSC-E, mHSC-GM, mHSC-L) across human and mouse systems [29]. These datasets span diverse biological contexts, including differentiation and developmental processes, and provide measurements that enable pseudotime inference. Each dataset is paired with up to four types of ground-truth networks:: cell-type-specific (Specific) network derived from ChIP-seq experiments and curated databases such as ENCODE [19], ChIP-Atlas [48], and ESCAPE [40]; non-specific (Non-Specific) network aggregated from large-scale resources such as DoRothEA [10], RegNetwork [24], and TRRUST [13]; functional interaction networks from STRING [32], capturing broader biological associations that may include indirect regulatory relationships; loss-/gain-of-function (LOF/GOF) evidence from ESCAPE [40] (available only for mESC). For each dataset, two gene subsets are constructed by selecting the top 500 or 1,000 most variably expressed genes with transcription factors (Abbr.: TFs+500 and TFs+1000, respectively), culminating in total 44 distinct evaluation tasks. Among these settings, future work may place particular emphasis on the Specific networks, as their cell-type-matched regulatory evidence provides the closest alignment between the expression measurements and the evaluation targets, making them a more representative testbed for context-dependent GRN discovery. The other network types remain valuable as complementary settings for assessing robustness to broader, heterogeneous, or functional interaction evidence.