1. Introduction

Understanding gene regulatory networks (GRNs) has been a long-standing pursuit in biology. The discovery of the "on-off" regulation mechanism in the 1960s lac operon experiment established gene expression as a switchable process. Consequently, the reconstruction of large-scale GRNs has emerged as a major goal in systems biology^[1]. Traditional GRN inference methods based on bulk omics data have enabled the characterization of genome-wide regulatory events, but they fail to capture cell type or state-specific GRNs in mixed samples. With the advancement of single-cell sequencing technologies, especially single-cell transcriptomics, it has become possible to infer specific GRNs and to track GRN changes across different developmental stages of cells. However, significant expression variation of certain genes across these stages poses considerable challenges for inferring GRNs from single-cell transcriptomic data.

To address this challenge, numerous computational methods have been proposed for inferring GRNs from single-cell RNA sequencing (scRNA-seq) data. A common strategy is to infer regulatory associations between transcription factors (TFs) and their candidate target genes from gene expression profiles, including methods capable of capturing non-linear dependencies. Some methods use dynamic modeling frameworks, such as SCODE^[2], to model temporal changes in gene expression. Since pseudotime and time-labeled single-cell transcriptomes are considered to provide dynamic information reflecting cell differentiation, methods such as LEAP^[3] apply time-delay correlation analysis to measure dynamic changes in gene expression. More recently, perturbation-aware GRN modeling studies have suggested that network topology, including hierarchical organization, modular structure, and sparsity, strongly influences the distribution of perturbation effects^[4]. These findings indicate that the structural properties of GRNs are important for understanding how regulatory systems respond to molecular perturbations, and further highlight the necessity of considering the global organization of GRNs rather than only individual regulatory interactions. In addition, recent cell-specific causal network approaches have been developed to depict dynamic biological processes at single-cell resolution, enabling the characterization of regulatory heterogeneity across individual cells^[5]. In recent years, deep learning-based approaches have further advanced GRN inference from single-cell transcriptomics. For example, DeepSEM^[6] and DeepRIG^[7] have been proposed to capture complex regulatory relationships without requiring strong assumptions about data distributions. In addition, scGREAT^[8] employs gene embeddings and a transformer-based architecture to learn representations of TF–target gene pairs by integrating scRNA-seq expression profiles with gene-related textual information, thereby capturing complex regulatory interaction patterns. More recently, TransGRN^[9] introduced a transfer learning-based framework for inferring GRNs across cell lines. It adopts a cross-cell-line pre-training strategy that combines scRNA-seq data from multiple source cell lines with biological knowledge derived from large language models, enabling effective GRN inference in few-shot settings where prior regulatory information is limited or unavailable. Overall, existing methods have improved the performance and scalability of GRN inference from high-throughput single-cell data. However, many approaches still focus mainly on pairwise regulatory associations and insufficiently exploit the global topological organization of GRNs, such as hierarchy, modularity, sparsity, and coordinated regulatory programs^[1,4]. Recent studies further show that GRN topology strongly influences perturbation effects, highlighting the need to model global regulatory structures beyond individual TF-target interactions^[4]. This motivates graph-based models that explicitly incorporate network topology for GRN inference.

Therefore, we propose a novel deep graph model for single-cell GRN inference, named ZINB-GRAN. This framework effectively infers GRNs through adversarial networks. Specifically, we first use correlation-based measurement methods to construct a weighted gene co-expression network (WGCN) from scRNA-seq data, thereby generating an initial regulatory graph. Subsequently, we transform the prediction of GRNs into a link prediction task, using graph convolutional networks (GCNs) as an encoder, combined with a decoder aimed at reconstructing the network topology structure. Meanwhile, ZINB-GRAN introduces an adversarial discriminator, which, through the construction of a generator-discriminator mechanism, forces the distribution of the latent expression space learned by the encoder to be consistent with the zero-inflated negative binomial (ZINB)-derived prior, thereby explicitly introducing biological statistical distributions as regularization terms. This mechanism promotes the distributional robustness and generalization capability of expression modeling.

In order to comprehensively evaluate the performance of the ZINB-GRAN model in GRN inference, it was tested on six real-world datasets and compared with other models. Benchmarking results show that ZINB-GRAN can accurately reconstruct GRNs and outperform existing methods.

2. Materials and Methods

2.1 Overview

Based only on the preprocessed scRNA-seq expression matrix, we construct an initial prior graph G_prior = (V,E_prior,X), where V = {v₁,...,v_n} denotes the set of genes, E_prior represents the co-expression-based weighted edge set, and X∈R^n×m denotes the gene expression feature matrix with n genes and m cells. For each gene pair (i,j), the Spearman rank correlation coefficient r_ij is calculated using their expression profiles across cells. The weighted adjacency matrix of the initial graph is defined as:

(1)$A_{ij}=\{\begin{array}{cc}{r}_{ij},& i\ne j\\ 0,& i=j\end{array}$

Thus, A is constructed solely from scRNA-seq expression data and does not use any ground-truth regulatory edges or test-set labels. For graph convolution, self-loops are added to A and symmetric normalization is applied before inputting the adjacency matrix into the GCN. For specific symbol definitions, see Table 1.

ZINB-GRAN takes single-cell gene expression profiles as input and generates a regulatory score matrix using a graph autoencoder (GAE) framework, as illustrated in Figure 1. First, the scRNA-seq dataset is preprocessed to filter out low-expression genes and low-quality cells. Then, Spearman rank correlation coefficients are calculated between gene pairs to construct a WGCN, which is used as the initial prior regulatory graph. Based on this prior graph, the GAE model, consisting of a GCN encoder and a scoring-function decoder, learns latent gene representations and reconstructs the GRN. Adversarial training is further introduced to regularize the latent gene representations by aligning them with a continuous ZINB-derived prior distribution, thereby improving the robustness of the inferred GRN topology. Finally, TF–gene pairs are ranked according to their predicted regulatory scores to reconstruct the GRNs.

Display Full Size

Download

Figure 1. Overview of the ZINB-GRAN: Starting with the count matrix from scRNA-seq data as input, ZINB-GRAN first constructs a WGCN from gene expression data. Based on this WGCN, it builds an initial regulatory graph for the genes. The initial regulatory graph and gene expression profiles are then input into a GAE. The GAE model consists of a GCN and a scoring function: the GCNs serve as the encoder, learning the global regulatory structure and embedding it into gene representations, while the scoring function acts as the decoder, scoring the gene pairs’ representations and reconstructing the GRN. The GAE aligns the latent representation Z with the prior distribution. scRNA-seq: single-cell RNA sequencing; WGCN: weighted gene co-expression network; GAE: graph autoencoder; GCNs: graph convolutional networks; GRN: gene regulatory network.

2.2 Constructing the initial regulatory graph

ZINB-GRAN constructs a WGCN^[10] based on the Spearman rank correlation coefficient of single-cell simulated expression data. All correlation coefficients are retained to build a prior regulatory graph for the genes, with the goal of reconstructing the GRN from the global regulatory structure. Specifically, given the log-transformed expression values log(TPM+1) or log(CPM+1) for genes i and j, (where transcripts per million (TPM) denotes transcripts per million, and counts per million (CPM) denotes counts per million), as X_i and X_j, their Spearman rank correlation coefficients score is calculated as follows:

(2)$\begin{array}{c}{r}_{ij}=\frac{\mathrm{cov}(R_{X_i},R_{X_j})}{{\sigma }_{R_{X_i}}\cdot {\sigma }_{R_{X_j}}}\\ =\frac{\sum _{n=1}^N(R_{X_i},n-{\overline{R}}_{X_i})(R_{X_j,n}-{\overline{R}}_{X_j})}{\sqrt{\sum _{n=1}^N{(R_{X_i,n}-{\overline{R}}_{X_i})}^2}\sqrt{\sum _{n=1}^N{(R_{X_j,n}-{\overline{R}}_{X_j})}^2}}\end{array}$

In the formula, n∈{1,2,..., N} indexes single cells; cov(R_Xi,R_Xj) represents the covariance between the ranks of the expression values of genes i and j across all cells; R_Xi denotes the vector of ranks of gene i’s expression values across cells, and R_Xi,n is the rank of gene i’s expression value in cell n; ${\overline{R}}_{X_i}$ is the mean rank of gene i across all cells, and σR_X,i is the standard deviation of the rank sequence of gene i. It is worth noting that r_ij∈[-1,1], where a value of indicates a perfect positive correlation between genes i and j, -1 indicates a perfect negative correlation, and 0 indicates no correlation (independent relationship). By calculating the Spearman rank correlation coefficients for all gene pairs, the edge weights between genes in the gene co-expression network can be obtained.

2.3 Graph autoencoder for reconstructing the GRN

After obtaining the initial regulatory graph and preprocessed gene expression data (see Section 3 for full preprocessing details), we use a GAE model^[11] to learn from the gene expression profiles and embed the global regulatory structure into the latent representations of genes, thereby reconstructing the regulatory network. The encoder in this study is implemented based on a GCN to process both gene expression data and the adjacency matrix of the prior regulatory graph. It combines the input features with the adjacency matrix through convolution operations, updating the latent representation of each gene. Specifically, the encoder firstly maps the input data to a new space using a weight matrix, then performs graph convolution operations to propagate the information from the adjacency matrix to each node (gene) and applies a nonlinear transformation through an activation function. The encoder also includes a normalization step for the adjacency matrix to ensure that the nodes and edges in the graph are properly handled. In this way, the encoder captures the complex interactions between genes and provides latent representations for subsequent network reconstruction.

To ensure that nodes with similar domains have similar latent representations, the latent representations of the nodes are learned through a hierarchical transformation using a GCN based on node v_i and the adjacency matrix A, resulting in the latent representations of the nodes z_i:

(3)$Z^{(l+1)}=f(Z^{\left(l\right)},A\mid W^{\left(l\right)})$

Here, W^(l) is the filter parameter matrix that we need to learn in the neural network, Z^(l) is the input to the convolution, and Z^(l+1) is the output after the convolution. For our problem, we have Z⁰ = X∈R^n×m (n nodes and m features).

Each layer of the GCN can be expressed using the function f(Z^(l),A|W^(l)), as shown below:

(4)$f(Z^{\left(l\right)},A\mid W^{\left(l\right)})=\phi \left({\tilde{D}}^{-\frac{1}{2}}\tilde{A}{\tilde{D}}^{-\frac{1}{2}}{Z}^{\left(l\right)}{W}^{\left(l\right)}\right)$

Here, $\tilde{A}=A+I$ and ${\tilde{D}}_{ii}=\sum _j{\tilde{A}}_{ij}$, I is the identity matrix, and φ is the activation function. Overall, the graph encoder G(X,A) is constructed using a two-layer GCN.

The construction of the graph encoder is as follows:

(5)$Z^{\left(1\right)}=f_{\mathrm{ReLU}}(X,A\mid W^{\left(0\right)})$

The first and second layers use the ReLU(•) function and a linear activation function, respectively. Our graph convolutional encoder G(Z,A) = q(Z|X,A) encodes the graph structure and node features into a representation Z = q(Z|X,A) = Z⁽²⁾.

The ZINB-GRAN model uses a decoder module with DistMult factorization to score regulatory relationships between gene pairs from convolutional representations e∈R^c. The model evaluates gene pair relationships (i,j) through single-layer DistMult operations^[12] as the core scoring mechanism to reconstruct regulatory network topology. The single DistMult function is formulated as:

(6)$\mathrm{score}(i,j)=e_i^{T}W{e}_j$

Where e_i^T represents the transpose of the representation e_i, and W∈R^c×c is a c-dimensional trainable weight matrix. Similarly, the scoring operation for all gene pairs can be efficiently implemented as ZWZ^T.

The supervised reconstruction loss for regulatory association prediction is computed only on the labeled gene pairs in the training mask. Let Ω_train denote the set of labeled TF–gene pairs used for training. The mask-based binary cross-entropy loss is defined as

(7)$L_{\mathrm{BCE}}=-\frac{1}{\left|{\Omega }_{\text{train}}\right|}\sum _{(i,j)\in {\Omega }_{\text{train}}}[y_{ij}\log \sigma \left({\text{score}}_{ij}\right)+(1-y_{ij})\log (1-\sigma \left({\text{score}}_{ij}\right))]$

Where σ(∙) denotes the logistic sigmoid function, y_ij∈{0,1} is the ground-truth edge label for gene pair (i,j) in the training mask, and score_ij is the predicted regulatory score for the pair (i,j). The held-out test labels are not used during model training and are used only for final evaluation.

2.4 Discriminator and generator network

The core idea of the model is to match the latent representation Z with the prior distribution through adversarial training. In this framework, the model uses a standard multi-layer perceptron structure, where the output layer has only one dimension and the result is transformed into a probability value through a sigmoid activation function^[11]. This adversarial model acts as the discriminator, with the primary task of distinguishing whether the latent representation Z comes from the predefined prior distribution p_z (i.e., positive samples) or from the latent representation G(X,A) generated by the graph encoder (i.e., negative samples). To achieve this, the model is trained by minimizing the cross-entropy loss of a binary classifier. As its prior distribution, the model employs a continuous ZINB-derived prior. The three ZINB parameters are fixed hyperparameters: π = 0.6, r = 2.0, and p = 0.4, where π denotes the zero-inflation probability, r denotes the dispersion parameter, and denotes the success probability. We first sample a discrete random variable U from a ZINB distribution:

(8)$U\sim \mathrm{ZINB}(\pi ,r,p)$

The probability mass function of U is defined as:

(9)$P(U=k)=\pi 1_{\{k=0\}}+(1-\pi )\frac{\Gamma (k+r)}{\Gamma \left(r\right)k!}{p}^r(1-p{)}^k,k=0,1,2,\dots .$

Because the encoder-derived latent representation Z is continuous, we do not directly align Z with the discrete variable U. Instead, U is transformed into a continuous ZINB-derived prior through logarithmic transformation, normalization, and Gaussian perturbation:

(10)$Z_{\text{prior}}=\frac{\log (1+U)-\mu }{\sqrt{{\sigma }^2+ϵ}}+\eta ,{\eta }^{\sim }N(0,{\sigma }_{\eta }^2)$

Here, μ and σ² denote the mean and variance of log(1+U) estimated from sampled prior values, ε is a small constant for numerical stability, and η is Gaussian perturbation used to obtain a continuous prior. The resulting Z_prior is used as the positive prior sample in adversarial training.

In the adversarial training process, the discriminator D is trained to distinguish samples drawn from the continuous ZINB-derived prior Z_prior from the encoder-derived latent representations Z = G_θ(X,A). The discriminator loss is defined as:

(11)$L_D=-\frac{1}{2}{E}_{Z_{\text{prior}}}[\log D\left(Z_{\text{prior}}\right)]-\frac{1}{2}{E}_X[\log (1-D\left(G_{\theta }\right(X,A\left)\right))]$

The graph encoder, which acts as the generator, is optimized to make the latent representations indistinguishable from the ZINB-derived prior samples. The generator adversarial loss is defined as:

(12)$L_{G,\mathrm{adv}}=-\frac{1}{2}{E}_X[\log D\left(G_{\theta }\right(X,A\left)\right)]$

Finally, ZINB-GRAN jointly optimizes the mask-based supervised reconstruction loss and the adversarial regularization loss. The total training objective is:

(13)$L_{\text{total}}=L_{\mathrm{BCE}}+{\lambda }_{\mathrm{adv}}{L}_{G,\mathrm{adv}}$

Where λ_adv is a hyperparameter controlling the contribution of adversarial regularization.

Cross-entropy loss measures the discriminator’s accuracy in classifying latent representations. By optimizing this loss, the generator (graph encoder) and discriminator compete with each other, enabling the generator to produce latent representations that better align with the prior distribution. Adversarial training is key, continuously improving the latent representation while avoiding overfitting.

Adversarial training forces the generator to learn not only from the data’s surface features but also from the latent space structure, optimizing the underlying patterns of the regulatory network. Through adversarial training, the generator continuously optimizes its latent representations, ensuring they not only match the statistical features of gene expression data but also effectively capture complex regulatory relationships between genes.

2.5 Evaluation metrics

We employed four numerical metrics used in BEELINE to quantitatively assess ZINB-GRAN’s performance in GRN inference. In addition to the common metrics of Area Under the Receiver Operating Characteristic curve (AUROC) and Area Under the Precision-Recall Curve (AUPRC), we also calculated the Early Precision Ratio (EPR) and AUPRC ratio (AUPR-R) values as defined by BEELINE^[13]. These additional metrics evaluate the true positive inference capability and the AUPRC value obtained by considering network density. EPR and AUPR-R represent the ratios of true positives and AUPRC values, respectively, between the model and a random predictor for the top k inferred edges, where k represents the number of edges in the true network.

3. Data Processing and Experimental Settings

Following the BEELINE^[13] benchmark framework, we evaluated the performance of the ZINB-GRAN model for the GRN inference task. This framework employs two categories of datasets for validation: the first comprises computationally simulated datasets, which include single-cell gene expression data generated from six synthetic regulatory networks and four manually curated regulatory networks; the second consists of datasets from real biological systems, covering six real scRNA-seq datasets derived from human and mouse cell lines. The six real scRNA-seq datasets include one human cell line and five mouse cell lines: human embryonic stem cells (hESC), mouse dendritic cells (mDC), mouse embryonic stem cells (mESC), and three lineages derived from mouse hematopoietic stem cells, erythroid lineage (mHSC-E), lymphoid lineage (mHSC-L), and granulocyte-macrophage lineage (mHSC-GM)^[9]. Detailed characteristics of these datasets are summarized in Table S1. To assess GRN inference performance, we used ChIP-seq data from ChIP-Atlas as the benchmark. The six scRNA-seq datasets were preprocessed using Scanpy^[14], retaining only high-quality cells and genes expressed in more than 9% of cells. After TPM and CPM normalization, the data were log-transformed. In accordance with the BEELINE standard^[13], we filtered genes using Bonferroni correction (p-value < 0.01) and fully incorporated transcription factor genes from the RegNetwork^[15] and TRRUST^[16] databases. Following deduplication, we selected cell type-specific TFs and the top 500 highly variable non-transcription factor genes for GRN inference in each cell type.

ZINB-GRAN is a GRN inference model that integrates prior regulatory networks with gene expression data. During evaluation, the datasets were randomly partitioned into non-overlapping training and test sets, and performance was validated using three-fold cross-validation. Detailed hyperparameter settings are provided in Table S2.

3.1 Compared approaches

We benchmarked ZINB-GRAN against eleven state-of-the-art GRN inference methods: PIDC^[17], PPCOR^[18], LEAP^[3], GRNBOOST^[19], SCODE^[2], SCRIBE^[20], DeepSEM^[6], CNNC^[21], 3DCEMA^[22], scGREAT^[8], and DeepRIG^[7]. For CNNC, a deep learning approach with limited TF-gene pair output, we adjusted its results to include a complete ranked edge list for fair comparison. While all methods provide confidence-based edge rankings, several (PIDC, GRNBOOST, DeepSEM) produce only unsigned networks, unable to distinguish activation from repression. Since real ChIP-seq data also lacks sign information, we evaluated all baselines on three unsigned GRNs with experimental validation. ZINB-GRAN offers the advantage of generating signed regulatory networks that indicate both activation (positive weights) and repression (negative weights), better capturing biological network complexity.

To ensure fairness at the evaluation level, all methods were assessed on the same benchmark datasets, the same reference GRNs, and using the same evaluation metrics, including AUROC, AUPRC, EPR, and AUPRC-R. For each baseline method, we used its original implementation or recommended configuration whenever possible, and did not use the held-out test set for hyperparameter tuning.

4. Results

4.1 Benchmarking Methods

Through systematic comparative analysis of computationally simulated gene expression data and real scRNA-seq data, this study comprehensively evaluated the performance of ZINB-GRAN compared to existing methods in GRN inference. The predictive capability of the model was also systematically validated across multiple types of datasets.

4.1.1 Evaluation of ZINB-GRAN using computationally simulated scRNA-Seq data

This study systematically evaluated ZINB-GRAN and other comparative methods using six categories of synthetic regulatory networks and their corresponding simulated single-cell gene expression profiles. The synthetic networks were constructed to represent typical regulatory topologies, including linear (LI), long-chain linear (LL), circular (CY), branched (BF), branch-fusion (BFC), and three-branch (TF) structures (Figure 2a). These networks provide benchmark datasets with fully known regulatory relationships for algorithm evaluation. Within this framework, the synthetic networks serve as reference standards for ground truth regulatory relationships, and performance was quantified by measuring the consistency between the networks inferred by each method and the standard networks. Among the methods compared, ZINB-GRAN, DeepRIG, CNNC, and 3DCEMA are supervised learning approaches, which rely on training data with explicit regulatory labels and optimize model parameters by minimizing the discrepancy between predicted outputs and the ground network structures.

Display Full Size

Download

Figure 2. Performance evaluation of the ZINB-GRAN model and other comparative methods on six synthetic regulatory networks. (a) Six types of synthetic regulatory networks used; (b) Performance comparison of ZINB-GRAN and other GRN inference algorithms based on four metrics. Each row presents the results of different models on various datasets under the same evaluation metric. If the EPR value of a model equals 1, it indicates that its predictive performance is equivalent to that of a random predictor; if the EPR value is less than 1, it signifies that its predictive performance is inferior to that of a random predictor. GRN: gene regulatory network; EPR: Early Precision Ratio.

Figure 2b presents the GRN inference performance of each model across six synthetic regulatory networks, and the corresponding quantitative results for AUROC, AUPRC, EPR, and AUPRC-R are summarized in Table S3.

The median AUROC values of ZINB-GRAN exceed 0.75 on all synthetic networks, with particularly outstanding performance observed in the LI, BF, BFC, and TF networks, where median AUROC values surpass 0.95, significantly outperforming all comparative methods. In terms of the AUPRC metric, ZINB-GRAN achieves the best performance in four of the six network types (LI, BF, BFC, and TF), with corresponding AUPRC values of 0.91, 0.98, 0.95, and 0.97, demonstrating a significant advantage over other algorithms. However, in the LL network, the median AUPRC of ZINB-GRAN is 0.4605, which is lower than that of the 3DCEMA method.

To further evaluate the model’s ability to identify true regulatory associations and validate its performance on the AUPRC-R metric, we extracted the “Top-k” regulatory relationships from the ZINB-GRAN inferred network, where k corresponds to the number of regulatory edges in the corresponding true reference network. The results show that ZINB-GRAN achieved the highest EPR values across the four synthetic network types: LI, BF, BFC, and TF. However, in the CY network, its EPR value was slightly lower than that of DeepRIG. Regarding the AUPRC-R metric, ZINB-GRAN significantly outperformed all comparison algorithms on the LI and CY networks, though it performed marginally worse than DeepRIG on the BFC and TF networks. Comprehensive analysis indicates that in synthetic networks with dense regulatory relationships, such as LI, CY, BF, BFC, and TF, ZINB-GRAN demonstrates superior inference performance, significantly surpassing current state-of-the-art algorithms. This confirms the model’s ability to effectively identify complex regulatory associations from computationally simulated data.

Although synthetic regulatory networks provide comprehensive regulatory information, they often fail to accurately simulate the biological regulatory mechanisms within cells. Therefore, this study manually selected four biologically meaningful regulatory networks from the literature to evaluate the practical inference capability of ZINB-GRAN (Figure 3a). These networks reflect the complex regulatory relationships involved in cell development and differentiation, including mCAD (5 genes), VSC (8 genes), HSC (12 genes), and GSD (19 genes). For each network, 30 sets of single-cell gene expression data were generated through computational simulation. To more realistically simulate the zero-inflation characteristic of single-cell data, three levels of expression dropout were systematically introduced: 10 sets with no dropout, 10 sets with a 50% dropout rate, and 10 sets with a 70% dropout rate. Subsequently, using the manually curated networks as the gold standard, ZINB-GRAN and other comparative algorithms were tested on all simulated datasets, and their performance in reconstructing GRNs was evaluated.

Display Full Size

Download

Figure 3. Performance of ZINB-GRAN on four manually selected real regulatory networks. (a) Four manually curated real regulatory networks serving as experimental benchmarks; (b) Performance evaluation of the ZINB-GRAN model and other GRN inference algorithms based on four metrics: AUROC, AUPRC, EPR, and AUPRC-R. Each row presents comparative results of different models across various datasets under a single evaluation metric. If a model’s EPR value is exactly 1, it indicates predictive performance equivalent to that of a random predictor; if the EPR value is less than 1, it signifies worse predictive performance than a random predictor. GRN: gene regulatory network; AUROC: Area Under the Receiver Operating Characteristic curve; AUPRC: Area Under the Precision-Recall Curve; EPR: Early Precision Ratio; AUPRC-R: AUPRC ratio; mCAD: mammalian cortical area development; HSC: hematopoietic stem cells development; VSC: ventral spinal cord; GSD: gonadal sex development.

Figure 3b shows the performance comparison across four manually curated regulatory networks, and the corresponding quantitative results for AUROC, AUPRC, EPR, and AUPRC-R are summarized in Table S4. ZINB-GRAN achieved the highest median AUROC and AUPRC values across all four networks, significantly outperforming all baseline methods. Notably, in the mCAD network, most baseline models exhibited median AUROC values below 0.5, except for SCODE, indicating that their inference performance was no better than that of a random predictor. In contrast, ZINB-GRAN attained a median AUROC of 0.996 on this network. Furthermore, in the HSC and GSD networks, all baseline methods showed relatively low median AUPRC values, while ZINB-GRAN maintained a median AUPRC above 0.85. In terms of the EPR metric, ZINB-GRAN clearly outperformed other methods in the HSC and GSD networks, demonstrating its ability to effectively recover true regulatory interactions among the top-k high-confidence predictions. For the AUPRC-R metric, ZINB-GRAN also performed best across all four manually curated networks. Taken together, these results indicate that ZINB-GRAN can accurately reconstruct the multi-layered regulatory relationships within complex, manually selected regulatory network models.

To evaluate the robustness of model performance under different missing rate conditions, we conducted paired Wilcoxon signed-rank tests on the AUROC values of each method across four simulated cell-type datasets (mCAD, HSC, VSC, and GSD) under 0%, 50%, and 70% missing rates (Figure 4). The detailed AUROC values under different missing rate conditions are provided in Table S5, and the corresponding paired Wilcoxon signed-rank test results are summarized in Table S6. The results showed that, in the mCAD dataset, no significant differences were observed between missing-rate conditions for any method (all p > 0.05), indicating that the performance of all models remained generally stable as the missing rate increased in this dataset. For the HSC and VSC datasets, neither ZINB-GRAN nor DeepRIG exhibited significant differences across any pairwise comparison of missing-rate conditions (all p > 0.05), suggesting strong robustness. In contrast, PPCOR, GRNBoost2, and DeepSEM showed significant differences in some comparisons. For the GSD dataset, most methods were more sensitive to missingness. In particular, PPCOR and GRNBoost2 showed significant differences in all three pairwise comparisons, whereas both ZINB-GRAN and DeepRIG exhibited significant differences between 0% and higher missing rates (50% or 70%), but not between 50% and 70%. Overall, ZINB-GRAN maintained stable performance across varying missing rates in three of the four simulated datasets (mCAD, HSC, and VSC), with significant fluctuation observed only in the GSD dataset, demonstrating a robustness pattern similar to DeepRIG and superior to that of most baseline methods.

Display Full Size

Download

Figure 4. Model performance on computationally simulated scRNA-seq data with varying expression missingness rates. (a)-(d) Each subplot corresponds to a manually selected regulatory network model: mCAD, HSC, VSC, and GSD, respectively. Each box plot illustrates the distribution of AUROC values obtained through 10 repeated experiments for each method under a specific expression missingness rate. In the color scheme, blue, yellow, and green represent simulated scRNA-seq datasets with expression missingness rates of 0% (complete data), 50%, and 70%, respectively. The red dashed line indicates the reference baseline for the theoretical optimal AUROC value of 1.0. mCAD: mammalian cortical area development; HSC: hematopoietic stem cells development; VSC: ventral spinal cord; GSD: gonadal sex development; scRNA-seq: single-cell RNA sequencing; AUROC: Area Under the Receiver Operating Characteristic curve.

4.1.2 ZINB-GRAN enables accurate reconstruction of GRN structures from real scRNA-seq data

To further validate ZINB-GRAN in biologically realistic settings, we evaluated it on six real scRNA-seq datasets. These datasets, which capture the complexity of actual cellular systems, pose a greater challenge than simulated networks with predefined rules. For reliable benchmarking, we established cell type-specific reference GRNs based on ChIP-seq data as the ground truth. Among the compared methods, ZINB-GRAN, DeepRIG, scGREAT, CNNC, and 3DCEMA are supervised models trained using portions of the ground truth, while the others are unsupervised and operate independently of such prior information. As shown in Figure 5a, the node degree distributions across all six reference networks follow a long-tail pattern. Notably, the networks for mESC, mHSC-L, mHSC-E, and mHSC-GM also exhibit a near-normal distribution of node degrees. This pattern indicates a concentrated connectivity density, reflecting typical topological properties inherent to biological regulatory networks.

Display Full Size

Download

Figure 5. Performance evaluation of models on real scRNA-seq datasets. (a) Node degree distribution of cell type-specific reference networks; (b) Performance comparison of ZINB-GRAN and other GRN inference algorithms across different evaluation metrics. Each row presents the results of different models on specific datasets under the same evaluation metric. scRNA-seq: single-cell RNA sequencing; GRN: gene regulatory network.

As shown in Figure 5b, the ZINB-GRAN model achieved AUROC values exceeding 0.6 across all cell types, demonstrating robust discriminatory performance. Detailed quantitative results, including AUROC, AUPRC, and EPR, are provided in Table S7. It performed best in the mHSC-GM cell type, with an AUROC of 0.8211, while its performance was comparatively lower in mDC cells, attaining an AUROC of 0.6141. Notably, ZINB-GRAN outperformed other models, particularly in the mHSC-E and mHSC-GM cell types. In terms of AUPRC, the model also showed a clear advantage in these two cell types, achieving AUPRC values of 0.8277 and 0.8448, respectively, indicating a strong ability to identify positive samples. Although performance was weaker in mDC cells, with an AUPRC of 0.0968, ZINB-GRAN still exhibited better overall performance compared to other methods. In summary, the model demonstrates excellent performance in specific cell types and is recommended for prioritized application in these contexts. Further optimization targeting cell types with lower performance may enhance its overall predictive capability.

To further validate the accuracy of ZINB-GRAN in high-confidence predictions, we evaluated the alignment between its top-k predicted regulatory interactions, where k equals the number of edges in the ground truth network, and the actual regulatory network. As shown in Figure 5b, ZINB-GRAN consistently achieved EPR values greater than 1 across all six real scRNA-seq datasets, indicating that its predictive performance consistently surpasses that of a random predictor. Notably, on the four datasets whose ground truth networks exhibit near-normal degree distributions, such as mESC and mHSC-L, ZINB-GRAN outperformed all baseline methods in terms of EPR. In particular, for the hESC dataset, ZINB-GRAN attained an EPR of 1.36, significantly higher than the 1.10 achieved by the second-best method. This result suggests that the method is particularly suitable for single-cell expression data conforming to a normal distribution assumption, which is widely adopted as a fundamental premise in single-cell gene expression analysis. Furthermore, Figure 5b shows that GAE-based models, ZINB-GRAN and DeepRIG, generally perform better than the variational autoencoder-based method DeepSEM. This improvement can be attributed to the fact that GAE models not only learn regulatory patterns from expression data but also effectively integrate prior gene graph information derived from WGCNs, confirming the critical role of prior topological structure in GRN inference. Further comparison reveals that ZINB-GRAN outperforms DeepRIG on hESC and mHSC-type datasets. This advantage stems from ZINB-GRAN’s use of adversarial learning to align the latent expression distribution from the encoder with a ZINB-derived prior, thereby explicitly incorporating this biological statistical distribution as a regularization constraint. This finding indicates that introducing a ZINB-derived prior, which aligns with the characteristics of single-cell data, can provide essential statistical guidance for inferring GRNs.

Moreover, among the five supervised learning models, ZINB-GRAN achieved AUROC values comparable to those of CNNC. However, in tasks with significant class imbalance, such as GRN inference, AUROC can be a misleading metric. A model may exhibit a strong ROC curve by overfitting to the majority negative class, namely non-regulatory interactions, while performing poorly in identifying the minority positive class, namely true regulatory relationships. Therefore, AUPRC and EPR, which better reflect a model’s ability to detect positive instances, serve as more appropriate evaluation metrics in this context. When comparing the five supervised methods, ZINB-GRAN significantly outperformed CNNC, scGREAT, DeepRIG, and 3DCEMA in both AUPRC and EPR. These results confirm that the superior performance of ZINB-GRAN does not stem from overfitting to negative samples but rather from its stronger capability to identify genuine regulatory interactions. This demonstrates a substantial improvement in predictive performance for the task of GRN inference.

As shown in the ablation experiments in Figure 6, the graph adversarial autoencoder ZINB-GRAN outperforms the GAE-based DeepRIG model. This advantage stems from the synergistic effect between the continuous prior distribution derived from the ZINB-derived model and the adversarial network framework.

Display Full Size

Download

Figure 6. GAE model improves GRN prediction performance. Each color represents a different model. GAE: graph autoencoder; GRN: gene regulatory network.

Across all real scRNA-seq datasets, ZINB-GRAN and DeepRIG consistently and significantly outperform WGCNA. This indicates that GAE-based models can learn global regulatory structures from prior regulatory graphs and identify latent regulatory patterns through adversarial training, thereby substantially improving the inference performance of GRNs. Notably, WGCNA achieves results comparable to those of LEAP and PPCOR across all evaluation metrics, suggesting that correlation-based WGCNs can effectively capture gene regulatory patterns.

4.2 Using ZINB-GRAN to infer cell type-specific GRNs from peripheral blood mononuclear cells (PBMCs)

Cell type-specific GRNs determine expression patterns through transcription factor interactions, directly influencing cellular function. Recent large-scale single-cell immune atlases have revealed substantial heterogeneity among human immune cell populations and provided valuable references for immune cell-type annotation and regulatory program analysis^[23]. To validate ZINB-GRAN’s performance in inferring cell-specific GRNs, we applied it to scRNA-seq samples of over 8,000 human PBMCs from the PBMC8k dataset provided by 10x Genomics.

After preprocessing the data and performing clustering analysis using Scanpy^[14], we manually annotated the cells to determine their types and identified marker genes for each cell type based on literature research. We focused particularly on four major cell populations: CD4 T cells, CD14 monocytes, CD8 T cells, and B cells. We constructed specific GRNs for these cell types. To obtain a gold standard network for PBMCs, we acquired a comprehensive network from the hTFtarget^[24] database, which includes ChIP-seq data, transcription factor binding sites, and epigenetic modification information. Subsequently, we trained the ZINB-GRAN model on each cell type and inferred their respective specific GRNs. In the GRNs inferred by the ZINB-GRAN model, we retained regulatory relationships with edge weights exceeding 0.2. Additionally, we ranked genes based on their degree in the regulatory network and identified those with node degrees greater than 20 as core genes bound by TFs.

We reconstructed GRNs for four cell types, containing 50 for CD4 T cells, 2,000 for CD14 monocytes, 511 for CD8 T cells, and 250 for B cells, regulatory relationships, respectively. Table S8 lists the regulatory factors predicted and identified by the ZINB-GRAN model in various cell types. The model predicted the most regulatory factors in CD14 monocytes, corresponding to the number of associations predicted. In CD14 monocytes, 25 core genes bound by TFs were predicted, with four identified as specific regulatory factors including JUND, CEBPB, RPS4Y2, and BCL11A. Among these, CEBPB, encoding C/EBP, has been reported to directly bind to and activate the monocyte-specific CD14 promoter, thereby enhancing CD14 gene expression during monocytic differentiation^[25]. This supports the potential role of CEBPB-mediated regulatory programs in monocyte-associated gene expression. Further analysis revealed that CEBPB regulates CD68, TYROBP, LST1, S100A11, and other monocyte marker genes. These four genes are highly expressed only in CD14 monocytes, with low expression in other cells. Research shows that SPI1 absence or inhibition significantly affects monocyte and macrophage development, influencing their immune response roles, particularly in regulating responses to external stimuli and antigen presentation^[26]. Additionally, GATA2 was identified as a key factor associated with monocyte development. GATA2 deficiency, caused by heterozygous mutations in the GATA2 gene, is characterized by the loss of multiple immune cell populations, including monocytes, B cells, dendritic cells, and natural killer cells^[27]. In CD14 monocytes, 55 genes have node degrees greater than 20, while in other cell types, these genes have node degrees less than 20.

In B cells, one of the genes bound by TFs, MEF2C, was also identified as a B cell marker gene, as shown in Figure 7b. Research indicates that MEF2C is a calmodulin-dependent transcription factor crucial for B cell proliferation and antibody responses, with its transcriptional activity regulated by calcium signaling pathways and various kinases^[28]. Another core transcription factor, EGR1, has been implicated in B-cell immune regulation. A recent single-cell study in esophageal squamous cell carcinoma identified senescent EGR1⁺ B cells as an immunologically relevant subset associated with impaired anti-tumor immunity and poor immunotherapy response^[29]. Two EGR1 target genes, IL4R and BIRC3, are highly expressed in B cells but show low expression in other cell types, as shown in Figure 7. Additionally, IL4R has been identified as a B cell marker gene, as shown in Figure 7b. PAX5 is a critical transcription factor in early B cell differentiation and lineage commitment. Recent studies have shown that SIRT7-dependent deacetylation regulates PAX5 stability and transcriptional activity, thereby coordinating B cell development, and that PAX5 alterations are closely associated with B cell acute lymphoblastic leukemia^[30]. Furthermore, RUNX3 was identified by the ZINB-GRAN model as a cell type-specific regulatory factor. In human B cells, RUNX3 has been reported to transcriptionally cross-regulate RUNX1, suggesting a potential role of RUNX-family TFs in B-cell regulatory programs^[31]. These results demonstrate that the cell type-specific gene regulatory predictions from the ZINB-GRAN model can assist in identifying regulatory factors and help annotate cell type marker genes.

Display Full Size

Download

Figure 7. The cell type-specific GRNs predicted by the ZINB-GRAN model assist in identifying cell types. (a) tSNE visualization of the four most abundant cell types: CD4 T, CD14 monocytes, CD8 T, and B cells; (b) 30 marker genes identified for the four cell types; (c) tSNE visualization analysis of the four cell types colored by the expression levels of marker genes such as CD68, TYROBP, LST1, S100A11, LGALS2, PSAP, IL4R, and BIRC3. GRNs: gene regulatory networks; t-SNE: t-distributed stochastic neighbor embedding.

We compare the GRNs constructed by ZINB-GRAN and WGCN, focusing on the differences between the inferred GRNs of CD14 monocytes and CD8 T cells from PBMCs, as shown in Figure 8. The inferred GRNs exhibit distinct network topologies and topological property distributions, such as node degree, closeness centrality, hub score, and edge betweenness centrality, which differ significantly from WGCNs. The inferred GRNs primarily follow a power-law distribution, resembling scale-free networks, where key hub nodes have extensive connections, a pattern not seen in WGCNs. This scale-free distribution is common in real biological networks, highlighting that inferred GRNs are more effective at identifying hub genes and key regulatory factors.

Display Full Size

Download

Figure 8. The final predicted GRN in PBMCs exhibits distinct topological features compared to WGCN. This study focuses on CD14 monocytes and CD8 T cells. (a) Regulatory network inferred by ZINB-GRAN in CD14 monocytes; (b) Regulatory network inferred by WGCN in CD14 monocytes; (c) Regulatory network inferred by ZINB-GRAN in CD8 T cells; (d) Regulatory network inferred by WGCN in CD8 T cells. For each subfigure, the left and right sections respectively display the network visualization and the distribution of four topological properties. GRN: gene regulatory network; PBMCs: peripheral blood mononuclear cells; WGCN: weighted gene co-expression network.

We also compared the differences between the GRNs inferred from CD14 monocytes isolated from PBMCs using various methods such as PPCOR, GRNBOOST2, DeepSEM, and DeepRIG, and the GRN predicted by the ZINB-GRAN method. Through visualization and analysis of the network structures inferred by these methods and the ZINB-GRAN model, this section reveals some interesting findings. As shown in Figure 9, the GRN inferred by GRNBOOST2 is similar to the one predicted by the ZINB-GRAN model but exhibits a more dispersed network structure. On the other hand, the GRN predicted by the correlation-based PPCOR method is closer to the weighted gene WGCNs. Notably, there are significant differences between these methods in terms of topological properties such as node degree, closeness centrality, hub score, and edge betweenness centrality. These findings suggest that the ZINB-GRAN model is more effective in capturing regulatory relationships in real data and excels in identifying genes bound by TFs.

Display Full Size

Download

Figure 9. GRNs predicted by different methods on CD14 monocytes. (a) ZINB-GRAN; (b) DeepRIG; (c) GRNBOOST2; (d) PPCOR; (e) DeepSEM. The left panels show the visualizations of the GRNs inferred by each method, while the right panels display the distributions of four topological properties, highlighting the differences. GRNs: gene regulatory networks.

4.3 Using ZINB-GRAN model to infer new regulatory factors in TNBC cells

To validate whether ZINB-GRAN can infer new regulatory factors, we applied it to predict GRNs from scRNA-seq data of triple-negative breast cancer (TNBC), a subtype characterized by the lack of common tumor markers. TNBC accounts for 10-20% of all breast cancers and has limited treatment options. We used scRNA-seq data from 1,534 cells of six untreated TNBC tumors and normal breast tissue, obtaining expression data for 13,280 genes. The gold standard network from the hTFtarget database^[24] was used for the breast model, filtering out non-breast tissue regulatory associations. After training ZINB-GRAN on normal breast tissue data, we applied it to infer the TNBC GRN, resulting in a network with 178,038 regulatory relationships involving 12,493 genes. We evaluated the network by calculating node degree, betweenness centrality, and PageRank scores, identifying key hub genes and regulatory factors based on their centrality and importance in the network.

To further prioritize candidate key genes, we integrated network topology analysis with GO enrichment results, as shown in Figure 10. First, significantly enriched GO biological process, cellular component, and molecular function terms were identified using an FDR threshold of 0.05. Genes involved in these significantly enriched GO terms were then intersected with highly connected genes in the reconstructed GRN. These genes were further prioritized based on network centrality metrics, including node degree, betweenness centrality, and PageRank scores. Based on this network-based prioritization strategy, we identified 10 candidate key genes, including MYC, ESR1, RARA, GATA3, HIF1A, KLF4, EP300, BRCA1, HDAC2, and E2F1. Among these 10 genes, 5 have been confirmed in multiple studies to play important roles in TNBC.

Display Full Size

Download

Figure 10. ZINB-GRAN model reveals potential regulatory factors in the triple-negative breast cancer network. (a) Network visualization with node size proportional to degree; high-degree genes are colored and form a central circular pattern, indicating their hub status; (b) Top 10 genes ranked by three network metrics: node degree, betweenness centrality, and PageRank scores. MYC and KLF4 show significant topological features, highlighting their importance in the regulatory network; (c) Gene Ontology enrichment analysis of high-connectivity genes with node degree greater than 1,000 shows overrepresentation of triple-negative breast cancer-related functions.

The GRN analysis of TNBC identifies several key regulatory factors. TNBC is a highly aggressive and heterogeneous breast cancer subtype with poor prognosis, and recent advances in molecular subtyping and targeted therapy have highlighted the importance of identifying reliable biomarkers and therapeutic targets^[32]. MYC dysregulation in TNBC has been associated with tumor proliferation, disease progression, poor prognosis, and therapeutic resistance^[33]. Both MYC and KLF4 function as major hub and bottleneck genes based on connectivity, betweenness centrality, and PageRank metrics. HIF1A, a key mediator of hypoxic responses in breast cancer, contributes to tumor progression, metastasis, immune evasion, and drug resistance, and may therefore represent a potential therapeutic target^[34]. E2F1 has been reported to promote TNBC migration, invasion, and metastasis through the PSMD14/E2F1/CENPF axis, partly by inhibiting ferroptosis^[35], while BRCA1 mutation status is associated with homologous recombination repair deficiency and therapeutic sensitivity to poly(ADP-ribose) polymerase inhibitors in TNBC^[36]. KLF4 may contribute to TNBC aggressiveness through the adipocyte-induced ANGPTL4/KLF4 axis, which promotes glycolysis and metastasis, and elevated ANGPTL4 and KLF4 expression is associated with poor prognosis^[37]. Gene Ontology enrichment analysis based on the Gene Ontology knowledgebase^[38] of highly connected genes with connectivity greater than 1,000 confirmed their significance in TNBC-related biological processes, demonstrating ZINB-GRAN’s utility for identifying potential regulatory factors and treatment targets.

5. Discussion

ZINB-GRAN is a GRN prediction framework based on graph adversarial autoencoders. It transforms scRNA-seq data into WGCN matrices and treats GRN inference as a link prediction task, using GCNs as encoders and decoders to reconstruct network topology. Through adversarial training, it optimizes gene representations. Testing on real datasets from PBMCs and TNBCs demonstrates that ZINB-GRAN effectively predicts GRNs and identifies genes bound by TFs. These results have been validated in the literature, confirming ZINB-GRAN’s utility for cell-specific GRN inference.

ZINB-GRAN builds GRNs via adversarial networks. Its training process distinguishes latent representation Z from predefined distributions, forcing the model to learn global patterns over local memorization. This creates robust GRNs with better noise tolerance for sparse data. Discriminator feedback optimizes network structure, enhancing performance on complex datasets. ZINB-GRAN surpasses eleven comparison algorithms in GRN inference.

Despite its advantages, ZINB-GRAN has limitations. While scRNA-seq data provides cellular gene expression insights, it alone cannot fully elucidate specific gene regulation mechanisms. Recent single-cell multi-omics methods, such as SCENIC+, have demonstrated that integrating gene expression with chromatin accessibility can link enhancers to target genes and improve the reconstruction of context-specific GRNs^[39]. Therefore, incorporating multi-omics single-cell data, particularly single-cell assay for transposase-accessible chromatin using sequencing, data, would offer more comprehensive information for accurately inferring GRN topology and regulatory mechanisms. This integration will be our future research focus.

6. Conclusion

In this study, we proposed ZINB-GRAN, a graph adversarial autoencoder-based framework for inferring GRNs from scRNA-seq data. Unlike many existing methods that mainly focus on pairwise regulatory associations, ZINB-GRAN incorporates the global topological structure of GRNs by constructing a weighted gene co-expression network and reformulating GRN inference as a link prediction task. By combining graph convolutional encoding, topology reconstruction, and adversarial regularization with a ZINB-derived prior, ZINB-GRAN learns robust latent gene representations that better capture complex regulatory relationships in sparse and noisy single-cell transcriptomic data.

Comprehensive benchmarking on simulated and real scRNA-seq datasets demonstrated that ZINB-GRAN achieves competitive or superior performance compared with existing GRN inference methods across multiple evaluation metrics. In particular, the model showed strong capability in reconstructing regulatory networks and identifying high-confidence TF–target relationships. Case studies on human PBMCs further showed that ZINB-GRAN can infer cell type-specific GRNs and identify biologically meaningful regulatory factors associated with immune cell functions. In addition, application to triple-negative breast cancer data revealed potential key regulatory genes and pathways related to tumor progression, metastasis, therapeutic response, and disease-associated regulatory programs.

Overall, ZINB-GRAN provides an effective and flexible framework for single-cell GRN inference by integrating graph-based representation learning with biologically informed distributional regularization. The results suggest that modeling global network topology is important for improving the accuracy and interpretability of GRN reconstruction. In future work, we will further extend ZINB-GRAN by incorporating single-cell multi-omics data, such as chromatin accessibility and epigenetic information, to more comprehensively characterize regulatory mechanisms and improve the inference of context-specific GRNs.

Supplementary material

The supplementary material for this article is available at: Supplementary materials.

Acknowledgements

The authors sincerely thank Zile Wang for her valuable suggestions on the construction of the ZINB-GRAN model and the development of the code.

The authors declare that ChatGPT based on GPT-4 was used solely for language polishing during the manuscript revision process. All research content, including study design, model development, data analysis, interpretation of results, experimental results, figures, tables, and conclusions, is original and was produced, checked, and approved by the authors. No AI tools were used to generate core scholarly content. The authors take full responsibility for the integrity, originality, and accuracy of the work.

Authors contribution

Zhang H: Conceptualization, supervision, methodology, software, formal analysis, visualization, writing-original draft, writing-review & editing.

Wang H, Ma H: Software, formal analysis, visualization, writing-review & editing.

Zhao J: Conceptualization, supervision, funding acquisition, writing-review & editing.

Xia J, Zheng C: Conceptualization, supervision, writing-review & editing.

Conflicts of interest

Not applicable.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

All data underlying the findings are fully available without restriction. The gene expression matrices of the scRNA-seq datasets analyzed in this study were downloaded from the Gene Expression Omnibus under the following accession numbers: GSE75748 for the hESC dataset, GSE98664 for the mESC dataset, GSE48968 for the mDC dataset, and GSE81682 for the mHSC dataset. The ground truth data for the six real scRNA-seq datasets, derived from ChIP-seq experiments, are available through ChIP-Atlas and BEELINE. The PBMC8k scRNA-seq dataset was obtained from 10x Genomics, and the datasets for normal breast cells and TNBC were downloaded from GEO under accession number GSE118390. The ground truth regulatory networks for PBMCs and normal breast cells were acquired from the hTFtarget database. The ZINB-GRAN package is publicly available at https://github.com/Hyz65/ZINB-GRAN.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 62362062), the Talent Program of Xinjiang Autonomous Region–Youth Outstanding Talent and Youth Innovative Talent (2023TSYCCX0104), and the Xinjiang Autonomous Region Tianshan Innovation Team Program (2025D14015).

Copyright

© The Author(s) 2026.