Introduction

Deciphering gene regulatory networks (GRNs) is crucial for understanding the complex biological functions of cells, revealing how precise spatiotemporal regulation of genes maintains functions in normal cells and dysfunctions in disease states [1][2][3]. scRNA-seq is a powerful technique for studying cellular dynamics, enabling measurement of subtle differences in gene expression between cells. This is especially useful for study of those participating in cellular processes such as cell differentiation or cell cycle progression [4]. By analyzing gene expression profiles, we can effectively capture the dynamic dependencies between genes, facilitating the reverse engineering of GRNs. Consequently, these inferred networks provide crucial biological information for the relationships between transcription factors (TFs) and their target genes [5].

Dozens of computational methods to infer GRNs from scRNAseq data have been proposed, uncovering the potential regulatory dependencies that underlie cellular phenotypes from a genomic perspective [5][6][7]. Most early methods focused on static scRNA-seq data, including GENIE3 [8], GRNBoost2 [9], SCENIC [10], NLNET [11], PIDC [12], PPCOR [13], and SINCERA [14]. Another type of GRN inference methods utilize time-ordered or pseudotime-ordered gene relationships across different stages of cell development, including LEAP [15], AR1MA1-VBEM [16], SINCERITIES [17], SCIM-ITAR [18], NME [19], Normi [20], Velorama [21], SINGE, [22] and scTIGER [23]. Most of these approaches begin by using probabilistic or regression models to fit the time or pseudotime series gene expression data and primarily rely on correlation metrics to predict gene regulatory relationships. For example, LEAP [15] uses Pearson's correlation and permutation testing to estimate pairwise co-expression between genes. MMFGRN leverages advanced matrix factorization techniques to integrate multi-modal data to infer regulatory relationships [24]. Additionally, some methods use differential equations to simulate changes of gene expression values over pseudotimed cell orders, including dynGENIE3 [25], Inference Snapshot [26], SCODE [27], SCOUP [28] and GRISLI [29]. Although those methods attempt to incorporate time-series features, they frequently rely on heuristic optimization or linear assumptions, limiting their ability to model complex, non-linear regulatory relationships.

Recently, the increasing availability of verified TF-gene pairs for specific cell types has facilitated the development of supervised deep learning methods for reconstructing GRNs, including CNNC [30], DeepDRIM [31], dynDeepDRIM [32], GENELink [33], GNNLink [34], TDL [35], SCGRNs [36], STGRNS [37], DeepIMAGER [38], and scGeneRAI [39]. Furthermore, several GRN inference methods, such as scMEGA [40], Dicts [41], LINGER [42], and SCENIC+ [43], utilize not only TF-binding motif enrichments but also other single-cell multi-omics data (e.g., scRNA-seq and scATAC-seq). Although the training accuracies of these supervised deep learning methods are notably high [37,38], the models trained on partially known TF-gene pairs remain biased due to the absence of a comprehensive global ground truth for the biological processes under investigation. Moreover, the TF-gene pairs used by these supervised methods are primarily derived from specific cell types, limiting the effectiveness of transfer learning when applied to different cell types [31,38,44]. As a result, unsupervised methods for inferring GRNs from transcriptomics data remain the dominant approach [45], serving as a crucial strategy for cells with limited availability of omics data.

The BEELINE framework [46] was designed for evaluating dozens of GRN inference tools on scRNA-seq data and found that tools incorporating pseudotime orders of cells showed better accuracy on synthetic datasets, but less accuracy on real datasets. This is because the scRNA-seq has considerable noise and sparsity, making it challenging to accurately infer GRNs from time-series scRNA-seq data [39,46]. Furthermore, the performance of GRN inference algorithms is often compromised by high rates of false positive regulations [47], especially in many correlation-based methods [10,15,27,48]. To address these limitations, we propose an unsupervised GRN prediction method by employing GRANGER recurrent autoencoders to network gene expression regulations (referred to as GRANGER). GRANGER mainly integrates GRANGER causality [49] into a recurrent variational autoencoder (VAE) [50] to infer GRNs. In addition, a negative binomial (NB) distribution is also incorporated into the loss function to correct for the distribution of gene expression profiles across all cells, as incorporating the NB distribution into deep learning frameworks has been demonstrated to improve both data imputation and the precision of cell clustering [51]. We evaluated GRANGER's inference accuracy using both simulation and real scRNA-seq datasets, and found it demonstrated superior performance compared to other unsupervised GRN inference methods.

Materials and methods

Overview of GRANGER

To infer GRNs from a scRNA-seq dataset of m genes, GRANGER translates the expression profiles of each gene g in all cells into a time-series sequence

where T represents the number of time points, adhering to the underlying principles of GRANGER causality during time series generation process [49] (Fig. 1A). GRANGER applies a causal recurrent VAE to all sequences X 1 , . . . , X g , . . . , X m to predict regulatory relationships among genes. The causal recurrent VAE in GRANGER consists of an encoder and a multi-head decoder. Each decoder head is responsible for generating X g , the time-series expression of a gene. The encoder reduces the multivariate time-series of each gene to a lower-dimensional feature space. In the inferred GRN, the edge between two genes is directed, representing a regulatory relationship between a gene (typically a TF) and its target genes (Fig. 1B).

Here the inferred GRN is represented as a GRANGER causal graph [50] G = (V, E), where the set of nodes V consists of the set of time series X 1 , . . . , X m , and each edge h, g ∈ E connecting X h to X g , which represents X h causes X g in time instant t with a lag, if the past profiles X h t-of X h provide statistically significant information about the prediction of X g t at time t (note that self-loops are allowed in graph G). The graph G = (V, E) can be characterized by an adjacency matrix A = A gh m×m ,

Therefore, the main objective in GRANGER is to infer the matrix A, each of whose elements A gh represents whether gene h regulates gene g.

T be the feature matrix of all sequences from time a to time b. As shown in Fig. 1B, a time segment X (t-2l-1):t is extracted from X 1:T of the m time-series sequences to set as a batch sample of training data. The first half X (t-2l-1):(t-l-1) of the time segment is taken as the input to encoder, and the second half X (t-l):(t-1) of the time segment and the latent space output by encoder are taken as the input of decoder, whose aim is to predict or reconstruct the segment X (t-l):t with time lag l. The decoder contains m heads and the g-th head predicts the g-th sequence X g , g = 1, 2, . . . , m.

A recurrent neural network (RNN) with gated recurrent units (GRUs) [52] is separately employed in the encoder and each head of decoder to form a recurrent VAE.

Data preprocessing and time series generation

The Python package SCANPY (version1.9.3) [53] was used for preprocessing and filtering the raw scRNA-seq data by: (i) removing genes expressed in fewer than 10% of cells, (ii) eliminating low-quality cells that express fewer than three genes or exhibit high mitochondrial gene expression, (iii) normalizing and logtransforming the filtered matrix, and (iv) selecting a subset of highly variable genes for further analysis. If the scRNA-seq data did not contain inherent time-series or pseudotime information, the PAGA method [54] was used through the 'scanpy.tl.paga' function to predict the pseudotime order of cells.

Variational autoencoder based on GRANGER causality

Suppose gene regulatory relationships can be modeled as a GRANGER causal graph G for a set of time series X 1 , . . . , X m . Let C (X g ) be the set of causes (or parents) of X g in graph G, i.e., C (X g ) = X h : A gh = 1 . At time point t, X g t can be predicted by its past profiles X g t-and the past profiles C(X g ) t-of its cause variables C (X g ) . So X g t can be represented by a function of X g t-and C(X g ) t- with additive independent noise ε g t . That is,

To learn the function f g for each X g and infer the matrix A, we used a VAE whose encoder and decoder are both composed of RNNs. If we input the first half segment X (t-2l-1):(t-l-1) to the encoder and use the decoder to reconstruct the second half segment X (t-l):t , the VAE with time lag l is:

where D ϕ and E ψ represent decoder and encoder with parameter sets ϕ and ψ, respectively, and ε t is an additive term without any specific assumption. Furthermore, the decoder has multiple heads, in which the g-th head is used to learn Xg t = f g X g t-, C(X g ) t- in Equation (2). Finally, the full vector Xt is formed by stacking the output of all m heads. Thus, the VAE based on GRANGER causality can learn the functions f g : g = 1, 2, . . . , m .

For X τ ∈ R m , τ = t -2l -1 , . . . , tl -1 , let h τ be the corresponding hidden state in the encoder. The encoder is formulated by the following equations from (4) to (9), the first four of which are the GRU equations:

∼

where r τ is the reset gate vector of GRU, z τ is the update gate vector of GRU, represents Hadamard product, and each b denotes the corresponding bias. W r , W hr , W z , W hz , W in , W h , W μ , W σ ⊂ ψ are weight parameters, where W in and W h are the weight matrices for input and hidden state in the encoder, respectively, and W μ and W σ are the weights to calculate the mean μ and the standard deviation σ of the learned Gaussian distribution.

For X τ ∈ R m , τ = tl , . . . , t, let H τ denote the corresponding hidden state in the decoder and H g τ the hidden state of the gth head. The initial state of the decoder is sampled from the Gaussian distribution which is parameterized with μ and σ . That is,

where z ∼ N (0, I), V re is weight for reparameterization, and b re is the bias. The decoder is formulated by Equations (11) (18) for the g-th head, the first four of which are the GRU equations:

Downloaded from https://academic.oup.com/bib/article/26/2/bbaf089/8068119 by guest on 25 January 2026

where R τ is the reset gate vector of GRU, Z τ is the update gate vector of GRU and each b denotes the corresponding bias.

⊂ ϕ are weight parameters, where V g in , V g h and V g out are weight matrices for input, hidden state, and output layer of the g-th head in the encoder, respectively. To capture the characteristics of scRNA-seq data, X g t at time install t is assumed to be an NB distribution parameterized with mean μ g t and dispersion θ g t . Two NB output layers with Equations ( 17) and ( 18) were added to the last hidden state of the g-th head H g t to estimate the mean μ g t and the dispersion θ g t of the learned NB distribution, respectively. V g μ and V g θ are the learned weights from the last hidden layer of the two NB output layers, respectively. Âg is the g-th row of the estimated adjacency matrix  of the GRANGER causal graph.

Loss function and network optimization

It is assume that each time-series variable X g t fits an NB distribution, i.e., X g t ∼ NB μ g t , θ g t , where μ g t is the mean and θ g t is the dispersion. That is,

Then the NB-based loss function is formulated as the negative log-likelihood of NB distributions:

As the adjacency matrix Âof the predicted GRNs is assumed to be sparse, a sparsity-inducing penalty R  to  is appended to the loss function. The final loss function is defined as

where the first term is the expected negative reconstruction error quantifying how well the model can reconstruct X g (t-l):t for each g and the second term is the relative Kullback-Leibler (KL) divergence between the learned Gaussian distribution in the encoder and the prior distribution z ∼ N (0, 1). The first term in Equation ( 21) is estimated by computing the cross-entropy between the training sample distribution and the reported probabilities. We used ℓ 1 norm to compute R (• ), because ℓ 1 norm is easier to obtain sparse solutions than ℓ 2 norm and helps reduce the risk of overfitting. GRANGER was trained by minimizing the loss function in Equation (21) with the proximal gradient descent (PGD) on the weights of input layers V g in and stochastic gradient descent (SGD) on all other parameters [55]. A two-stage training strategy is proposed: (i) we first trained the model by computing the gradients of the convex part of L (ψ, ϕ) that is the first and second term in Equation ( 21), updating ϕ and ψ except V in using SGD and updating V in using PGD; (ii) we fixed all zero elements in V in based on obtained causal matrix Â, and continue training the model with SGD on L (ψ, ϕ). In practice, we used the iterative shrinkage thresholding algorithm (ISTA) with constant step size to do PGD [56]. To further improve the accuracy of predicted GRNs, two-layer GRUs are used in the recurrent VAE, where the input Y t of the second layer is the hidden state h (1) t of the first layer.

Benchmark datasets and case study for GRANGER evaluation

To fairly evaluate the performance of GRANGER, we utilized the six synthetic and four curated datasets provided by the GRN evaluation framework BEELINE [46]. The six synthetic datasets simulate a linear trajectory (dyn-LI), a cyclical trajectory (dyn-CY), a long linear trajectory (dyn-LL), a bifurcating trajectory (dyn-BF), a trifurcating trajectory (dyn-TF), and a bifurcating converging trajectory (dyn-BFC). The four curated datasets focus on hematopoietic stem cell differentiation (HSC), mammalian cortical area development (mCAD), ventral spinal cord development (VSC) and gonadal sex determination (GSD). Each of the six synthetic datasets contains five sub-datasets with different numbers of cells (100, 200, 500, 2000, and 5000 cells), and each subdataset consists of 10 samples. Each of the four curated datasets contains 2000 cells and 10 samples. For each of curated datasets, BEELINE also provides samples containing 50% and 70% dropout zeros. All these benchmark datasets contain pseudotime-ordered information that were predicted by Slingshot [57] and preprocessed by the BEELINE framework. In addition, we conducted performance validation on four real scRNA-seq datasets, each separately sequenced from human embryonic stem cells (hESC), mouse dendritic cells (mDC), human mature hepatocytes (hHEP), and erythroid-lineage mouse hematopoietic stem cells (mHSC-E). We predicted the GRNs for the TFs with corresponding GRNs provided by BEELINE, including hESC (eight TFs and 310 targeted genes), hHEP (two TFs and 378 targeted genes), mDC (two TFs and 31 targeted genes), and mHSC-E (six TFs and 459 targeted genes). As a real application, GRANGER was applied to predict GRNs from whole mouse brain data obtained from the Allen Brain Cell Atlas [58]. We downloaded the hippocampal formation (HPF) scRNA-seq data from the batch that was generated by 10Xv2 chemistry. We focused on Glutamatergic neurons (Glut) in the CA1 and prosubiculum (ProS) in HPF, labeled as CA1-ProS-Glut_1, which has 1055 sequenced neurons and predicted GRNs for five TFs, Sox10, Onecut2, Gbx1, E2f7, and Prox1. The data filtering and preprocessing of 1055 neurons was performed by using SCANPY (version1.9.3) [53]. PAGA was used to generate pseudotime series for the data [54]. For each TF, the top 20 highly correlated genes and top 100 highly variable genes were selected for GRN prediction by GRANGER. The functional enrichment of all genes in the merged GRN network were performed by using the PANTHER database [59][60][61].

Methods comparison and evaluation metrics

We compared GRANGER with popular unsupervised algorithms: GRNBoost2, SINCERITIES, PIDC, PPCOR, SCODE, GENIE3, and SINGE, which have shown high performance in validations within the comprehensive BEELINE framework [46] and other benchmarking studies [6,7,62]. We also included Normi, a recently published unsupervised algorithm that claims outstanding performance. The technical details of these eight methods are listed in Table 1 and Table S1. To ensure the consistency of The weights of the recurrent VAE in GRANGER were updated using the Adam optimizer in PyTorch [63], with a learning rate of 10 -3 . The default batch size and epoch number were set to 8 and 500, respectively. The parameters for GRNBoost2, PIDC, PPCOR, SCODE, SINCERITIES, GENIE3, and SINGE were same as those used in the BEELINE framework (https://github.com/Murali- group/Beeline/tree/master/Algorithms). The parameters of Normi were also set to be default. All experiments were conducted on a Linux workstation equipped with an Nvidia GTX4090 GPU. We evaluated GRANGER and the other methods by calculating the area under the receiver operating characteristic curve (AUROC) and the area under the precision-recall curve (AUPRC) of each sample. As GRN inference is performed on datasets with imbalanced positive and negative samples, we primarily rely on the AUPRC metric, which is recommended for evaluating such tasks [64]. For real datasets (hESC, mDC, hHEP, and mHSC-E), we calculated the AUPRC ratio as the AUPRC divided by the expected AUPRC of a random predictor, providing a normalized measure of performance relative to randomness. The early precision ratio (EPR) was calculated as the fraction of true positives among the top-k predicted edges compared to a random predictor.

Results

The NB-based loss function improves the performance of GRANGER causal recurrent VAE

We systematically designed the loss function of GRANGER, by including not only the similar reconstruction loss and KL divergence loss as traditional VAE, but also an NB-based loss L NB (ψ, ϕ) and a spare-inducing penalty R Â . The NB distribution has been employed to address the dropout events in scRNA-seq data [51,[65][66][67], and recent studies suggest that the NB model is more appropriate than other models for scRNA-seq data sequenced by droplet-based platforms [68,69].

We evaluated the impact of the NB-based loss hyperparameter λ NB in four datasets (dyn-CY, dyn-LI, HSC, and mCAD), each containing 2000 cells. We tested different values of λ NB , ranging from 0 to 5, and assessed the performance of GRANGERs using the AUROC and AUPRC metrics. Across all datasets, we observed that when λ NB was set to 1, GRANGER consistently achieved the highest AUROC and AUPRC scores compared to other values of λ NB (Fig. S1). Specifically, in the dyn-CY dataset, the AUPRC reached around 0.90 when λ NB was 1, demonstrating its effectiveness in improving predictive accuracy. Meanwhile, the AUPRC of the dyn-LI dataset was even higher at 0.98. In contrast, when λ NB was set to 0 (i.e., without the NB-based loss), both AUROC and AUPRC scores were significantly lower, indicating the importance of including the NB-based loss to handle missing data and improve GRN prediction accuracy. Similarly, the trend was observed across the parameter testing of the HSC and mCAD datasets. The results confirm that including the NB-based loss consistently improves GRANGER's ability to infer GRNs.

Proper sparsity-inducing penalty improves the performance of GRANGER

We tested the performance of GRANGER with different values of the scale factor λ A on four samples of 2000 cells from different datasets, fixing the other two hyperparameters with l = 200 cells and λ NB = 1. We ran GRANGER on the four datasets with different λ A settings ranging from 0.1 to 1. We found that the AUROC and AUPRC scores reached their peak when λ A was set to values between 0.2 and 0.4 (Fig. S2). For example, in the dyn-CY dataset, λ A = 0.4 resulted in the best performance, with an AUROC score of ∼ 0.84 and an AUPRC score around 0.9. Similarly, for the dyn-LI dataset, the optimal performance occurred with λ A = 0.3, where both AUROC and AUPRC scores were consistently high (AUROC ∼0.83 and AUPRC ∼0.98). In contrast, when λ A was increased to 0.8 or higher, the model showed a significant drop in performance, with some experiments failing to complete within the allocated 500 epochs. This degradation is likely due to the model becoming overly sparse, leading to the learned weights of the model tending to be smaller or even close to zero. Based on these results, we set the default value of λ A to 0.3 for subsequent experiments, as it provided superior results across multiple datasets.

Optimizing time lag and GRU structure of GRANGER

Besides the hyperparameters λ NB and λ A , the time lag l of each training sample also has a significant impact on causality inference. We tested how an appropriate division of time lag l could help the model to improve GRANGER causality inference. We first ran GRANGER on four samples from four datasets (dyn-CY, dyn-LI, HSC, and mCAD), each containing 2000 cells with settings of λ NB = 1, λ A = 0.3,and different time lag values, l = 20, 50, 100, 200, 300, 400, 500, and 1000. Among all time lags, the AUPRC and AUROC scores where l = 200 rank in the top two, consistently outperforming most other time lags in the four samples (Fig. 2A). For example, in the dyn-CY-2000 dataset, the highest AUROC and AUPRC scores (around 0.9 and 0.86, respectively) were observed when the time lag was set to 200. Similar trends were seen in dyn-LI-2000 and mCAD-2000, with the best performance observed at a time lag of 200 or 300. In addition, we conducted experiments on samples with different cell numbers for dyn-CY and dyn-LI with varying time lags. As shown in Fig. S3, we found that different optimal settings of l values depend on the sample size. For samples containing 100 cells in dyn-CY and dyn-LI datasets, l = 15 is a relatively optimal setting; for samples with 500 cells,l = 200 is a relatively optimal setting; and for larger samples of 5000 cells, l = 200 ∼ 2, 000 are relatively optimal settings. To reduce computational complexity, we suggest that the setting l = 200 for large samples containing 500 or more cells.

In the encoder and each head of the decoder, two-layer GRUs are employed instead of single-lay GRUs. We conducted experiments on both single-layer GRU and two-layer GRUs on four datasets (dyn-CY, dyn-LI, HSC, and mCAD) each containing 2000 cells. For each of the four datasets, the median AUPRC and AUROC scores of GRANGER with two-layer GRUs are significantly higher than those of the model with single-layer GRUs (Fig. 2B). In dyn-CY-2000, For example, the AUPRC score for the two-layer GRU structure reaches around 0.91, whereas the single-layer model achieves only ∼ 0.65. Similarly, in dyn-LI-2000, the two-layer GRU achieves an AUPRC score of around 0.95, compared to the singlelayer's lower score 0.89. These results support that, compared with single-layer GRUs, two-layer GRUs preserve more contextual information, enabling the model to capture more complex patterns from the input data.

Impact of pseudo-temporal information on GRN prediction accuracy

In GRANGER, we integrated the PAGA method [54] to infer pseudo-time orders for cells. To evaluate the impact of different pseudo-time trajectory inference methods on GRN prediction, we compared the performance of GRANGER using three widely recognized methods, Slingshot [57], PAGA [54], SCORPIUS, [70] and randomly shuff led pseudo-time series (as control) across six datasets (dynCY-2000, dynLI-2000, dynLL-2000, GSD, HSC, and mCAD). As shown in Fig. S4, the results demonstrate that pseudotime trajectory information generated by Slingshot and PAGA significantly enhances GRN prediction performance, achieving comparable and high AUPRC values across all datasets. SCORPIUS also yielded consistent results, though slightly lower than those achieved by Slingshot and PAGA. In contrast, GRANGER's performance on datasets with randomly shuff led pseudo-time series was slightly reduced, underscoring the importance of accurately inferred pseudo-temporal trajectories for robust GRN inference. Notably, the synthetic datasets (dynCY-2000, dynLI-2000, dynLL-2000) showed a greater dependency on precise pseudo-time information compared to the curated datasets (GSD, HSC, and mCAD), ref lecting the complexity of GRNs in real biological systems. Additionally, the minimal differences between Slingshot and PAGA suggest their equivalence in predicting both tree and linear topological trajectories, consistent with findings from previous comparative analysis [71]. Overall, this analysis emphasizes the critical role of pseudo-temporal information in improving GRN inference and provides guidance on selecting suitable trajectory inference methods for scRNA-seq studies.

Comparing performance of GRANGER with other methods

We compared GRANGER with eight unsupervised GRN inference methods (GRNBoost2, PIDC, PPCOR, SCODE, SINCERITIES, GENIE, Normi, and SINGE). We ran GRANGER with the default hyperparameters trained in the previous sections (i.e., λ NB = 1, λ A = 0.3 and time lag l = 200 for the samples of 2000 cells). For six synthetic datasets (dyn-LI, dyn-CY, dyn-LL, dyn-BF, dyn-BFC, and dyn-TF), the results show that GRANGER consistently outperforms other methods across different datasets and cell numbers (Fig. 3). For example, in the dyn-LI and dyn-CY synthetic datasets, GRANGER achieves the highest median AUROC and AUPRC scores for cell numbers ranging from 100 to 5000 cells. For 100 and 200 cells in dyn-LI, GRANGER's AUROC scores range from 0.85 to 0.9, surpassing other methods, where the AUROC scores for GRNBoost2, PIDC, PPCOR, and SCODE are around 0.6-0.7 and the AUROC scores for SCODE f luctuate around 0.5-0.6. In dyn-CY, GRANGER achieves AUROC scores of 0.85-0.9 across all cell numbers, significantly outperforming methods like SCODE, which yields AUROC scores below 0.6 at smaller cell numbers. GRANGER's AUPRC scores show even greater improvement for all six datasets compared to the other methods. For complex trajectories (e.g., dyn-BF, dyn-TF, and dyn-BFC), GRANGER also maintains superior performance with minimal decline, whereas other methods exhibit a noticeable drop in AUPRC scores.

The performance improvement is most pronounced at smaller cell numbers (100 and 200 cells), where GRANGER demonstrates a significant advantage over competitors. For example, for dyn-LI with 100 cells, GRANGER achieves an AUPRC of ∼ 0.95, whereas SCODE shows much lower AUROC values around 0.55. Similarly, in terms of AUPRC, GRANGER outperforms other methods with values approaching 0.98 for small cell sizes, while some methods (i.e., GRNboots2, PIDC, and PPCOR) struggle to reach 0.4. These results highlight GRANGER's reliability in handling smaller datasets, which present greater challenges for accurate GRN inference.

For four curated datasets (HSC, VSC, mCAD, and GSD), GRANGER continues to demonstrate superior performance (Fig. 3, last row). On the HSC and mCAD datasets, GRANGER consistently achieves higher median AUROC and AUPRC scores compared to other methods. While its performance on the GSD dataset is more similar to that of SCODE and SINGE, GRANGER still manages to maintain competitive results. We also evaluated the performance of GRANGER on four real scRNA-seq datasets (hESC, hHEP, mDC, and mHSC-E) using ground truth GRNs derived from ChIP-seq datasets curated by BEELINE. The results are summarized in Tables S2 and S3, highlighting GRANGER's consistent superiority in GRN inference. Table S2 presents the AUPRC ratio values, where GRANGER achieved the highest scores across all datasets, including 1.57 for hESC and 1.91 for hHEP, outperforming other methods. Similarly, Table S3 showcases GRANGER's outstanding EPR values, such as 1.69 for hESC and 1.78 for hHEP. Taken together, these findings demonstrate that GRANGER maintains high predictive accuracy across diverse scRNA-seq datasets and varying cell numbers, highlighting its versatility and reliability in biological data analysis.

Robustness of GRANGER against dropout events

We conducted benchmarking using simulation datasets with different dropout levels to evaluate GRANGER's performance in handling dropout noise. For each of the four curated datasets, we compared GRANGER with eight methods on samples containing 50% and 70% dropout zeros. Figure 4 demonstrates that GRANGER's performance shows minimal degradation at both dropout rates compared to the original data. GRANGER consistently achieves the highest AUPRC scores for different dropout rates compared to the other eight methods. For example, GRANGER achieves AUPRC scores around 0.8 in the HSC dataset and 0.7 in the GSD dataset, whereas the other methods' AUPRC scores are all below 0.5 in the HSC dataset and below 0.3 in the GSD dataset. Although GRANGER performs similarly AUROC scores to other methods in several datasets, it achieves higher AUROC in the HSC and mCAD datasets across different dropout rates. These results show the strong resilience of GRANGER in predicting GRNs from sparse scRNA-seq data, reinforcing its suitability for handling real-world biological data with dropout challenges.

GRANGER detected potential GRNs in mouse brain cells

To access the effectiveness of GRANGER in predicting GRNs in real-world scenarios, we applied it to scRNA-seq datasets of the whole mouse brain [58]. As a case study, we focused on 1055 CA1-ProS-Glut_1 neurons. We specifically predicted GRNs for five TFs, Sox10, Onecut2, Gbx1, E2f7, and Prox1, which play critical regulatory roles in neurons [72,73]. GRANGER successfully predicted regulatory interactions for each TF: 46 regulations for E2f7, 19 for Sox10, 20 for Gbx1, 24 for Prox1, and 20 for Onecut2. These regulated genes are involved in several functional categories, including metabolism, transcription, cell adhesion/structure, the cell cycle, and other biological processes (Fig. 5A, see Table S4 for detailed annotations).

We identified that the GRNs of Sox10, Onecut2, Gbx1, E2f7, and Prox1 form a highly interconnected network comprising 95 genes, including 24 co-regulated genes (marked in pink in Fig. 5A). Several of these co-regulated genes, such as Fos [74], Bdnf [75,76], Ncam2 [77], and Penk [78], stand out for their significant regulatory roles in critical cellular functions. These co-regulated genes highlight the interplay of TFs in regulating diverse yet interconnected mechanisms in hippocampal neurons. To further investigate the biological mechanisms, we then performed functional enrichment analysis for these 95 genes and revealed multiple biological processes enriched in critical regulatory functions (Table S5). Notably, nervous system development (GO:0007399) and cellcell signaling (GO:0007267) were the most significantly enriched processes. The process of negative regulation of neuroblast proliferation (GO:0007406) had an exceptionally high log-fold change of 57.67. Additionally, hormone-related processes such as hormone secretion (GO:0046879) and hormone transport (GO:0009914) were notably enriched.

To validate the predicted regulatory interactions in the five GRNs, we conducted comprehensive literature reviews and analyzed ChIP-seq data to assess binding activities. First, a substantial proportion of the predicted interactions were supported by literature evidence indicating correlations or interactions between the TFs and their targets: 16 out of 24 interactions in the Prox1 GRN (Table S6), 14 out of 19 in the Sox10 GRN (Table S7), 14 out of 20 in the Onecut2 GRN (Table S8), 29 out of 46 in the E2f7 GRN (Table S9), and 13 out of 20 in the Gbx1 GRN (Table S10). Second, ChIP-seq data analysis provided additional confirmation of the predictions. ChIP-seq datasets were identified for Prox1, Sox10, and Onecut2 in mice; however, no suitable datasets were currently available for E2f7 and Gbx1. Notably, strong binding peaks were detected at critical genomic regions, including promoter regions, 3' UTRs, and introns, for all 24 target genes in the Prox1 GRN, 19 target genes in the Sox10 GRN, and 20 target genes in the Onecut2 GRN. Moreover, these binding peaks frequently co-localized with active regulatory markers such as K4me2, K4me3, K27ac, and ATAC-seq signals. Specifically, co-occupancy with these markers was observed for 20 of 24 targets in the Prox1 GRN, 13 of 19 in the Sox10 GRN, and 11 of 20 in the Onecut2 GRN. An example is shown in Fig. 5B, illustrating Prox1's regulatory activity over the Limd1 gene, where Prox1 ChIP-seq data in mouse neurons exhibit strong binding peaks at the Limd1 promoter region. Furthermore, these peaks overlap with active regulatory markers, including K4me3, K27ac, and ATAC-seq signals, highlighting Prox1's role in directly modulating Limd1 expression. These results demonstrate the reliability of GRANGER in accurately identifying novel regulatory interactions across diverse TFs with strong support by both literature evidence and ChIP-seq signals.

One of the key insights of GRANGER is that it can effectively detect gene regulatory relationships, even in cases where gene coexpression and correlation are not strong. For example, among the 46 genes in the E2f7 GRN, there are four subgroups showing high expression levels and strong correlations within each group, while many other genes exhibit low correlations (Fig. 5C). Similar patterns were observed in the GRNs of Sox10, Onecut2, Gbx1, and Prox1 (Fig. S5). The ability of GRANGER to infer regulatory interactions without always depending on high co-expression among genes allows it to identify meaningful biological relationships that other correlation-based methods might miss.

Discussion and conclusion

In this study, we introduced GRANGER, an unsupervised deep learning method to infer GRNs from scRNA-seq data. The combination of recurrent VAEs, GRANGER causality and NB-based loss function offers several advantages over traditional GRN inference methods, particularly in terms of capturing temporal dependencies, addressing data sparsity and improving prediction accuracy. One of the notable advantages of using recurrent VAEs is their ability to capture sequential patterns and temporal dynamics in time-series data. Traditional autoencoders compress data into a lower-dimensional space, but recurrent VAEs extend this by using memory cells (such as GRUs) to model time-dependent variations. By leveraging this capability, GRANGER captures gene expression changes over time, which is essential for accurately reconstructing the underlying GRNs in dynamic biological processes. Furthermore, by combining recurrent VAEs with GRANGER causality, GRANGER offers a causal interpretation of gene regulation that is often missing in many correlation-based approaches. Considering that scRNA-seq data are extremely sparse and noisy due to experimental limitations [51,[65][66][67], GRANGER incorporates an NB-based loss function to well leverage statistical characteristics of the data. This feature, combined with a sparsity-inducing penalty, enables GRANGER to achieve superior performance and robustness against dropout noise over eight competing methods.

We applied GRANGER to scRNA-seq data from the BRAIN Initiative project, specifically targeting the identification of GRNs for five TFs (E2f7, Gbx1, Sox10, Prox1, and Onecut2) in neurons. GRANGER recalled known regulations but it also identified novel insights that are functional important and could be experimentally validated in future studies. For example, additional ChIP-seq experiments on mouse cell lines or tissues can confirm TF binding signals and provide direct evidence of these TF regulating their targeted genes under different conditions. Such experimental validations would validate GRANGER's predictions and uncover novel insights for the functional roles of TFs in complex biological systems.

Despite its advantages, GRANGER has certain limitations that need further exploration and improvement. First, the reliance on linear GRANGER causality may not fully capture non-linear dependencies between genes [79]. Integrating non-linear causality measures or hybrid models could enhance GRANGER's ability to detect more complex regulatory interactions. Furthermore, GRANGER can be expanded by incorporating additional omics data, such as epigenomics or proteomics, to build multi-modal models capable of providing a more comprehensive view of gene regulation [45]. We plan to enhance GRANGER's interpretability and accuracy by incorporating attention mechanisms [80], which can dynamically adjust the weighting of gene interactions based on the biological context. In summary, GRANGER provides an unsupervised deep learning framework capable of inferring robust and interpretable GRNs from noisy, high-dimensional scRNA-seq data, offering a powerful tool to discover novel regulatory interactions.

Key Points

• GRANGER is an unsupervised deep learning method that has superior performance in cell-specific GRN inference from scRNA-seq dataset.

• GRANGER uses recurrent VAEs and GRANGER causality to infer directional regulatory relationships among genes, offering a more comprehensive view of gene interactions compared to correlation-based methods. • GRANGER integrates NB-based loss functions and sparsity-inducing penalties to address the sparsity and dropout noise in scRNA-seq data.