Inferring gene regulatory network from single-cell transcriptomes with graph autoencoder model

The gene regulatory structure of cells involves not only the regulatory relationship between two genes, but also the cooperative associations of multiple genes. However, most gene regulatory network inference methods for single cell only focus on and infer the regulatory relationships of pairs of genes, ignoring the global regulatory structure which is crucial to identify the regulations in the complex biological systems. Here, we proposed a graph-based Deep learning model for Regulatory networks Inference among Genes (DeepRIG) from single-cell RNA-seq data. To learn the global regulatory structure, DeepRIG builds a prior regulatory graph by transforming the gene expression of data into the co-expression mode. Then it utilizes a graph autoencoder model to embed the global regulatory information contained in the graph into gene latent embeddings and to reconstruct the gene regulatory network. Extensive benchmarking results demonstrate that DeepRIG can accurately reconstruct the gene regulatory networks and outperform existing methods on multiple simulated networks and real-cell regulatory networks. Additionally, we applied DeepRIG to the samples of human peripheral blood mononuclear cells and triple-negative breast cancer, and presented that DeepRIG can provide accurate cell-type-specific gene regulatory networks inference and identify novel regulators of progression and inhibition.


Authors' response:
We thank the reviewer for raising this important point.The selection of a threshold is crucial in WGCNA.However, there is no consensus on the optimal threshold for all tasks, as it depends on the specific situation and analysis objective.In this study, we chose to retain all the correlation coefficients instead of adopting other predefined strategies.This approach was motivated by the expectation that it would capture a broader range of potential regulatory information than other strategies.We further compared the model performance under different thresholds of WGCN to verify our expectations (below and in Supplementary Fig 4).The model has the highest median EPR value when preserving all the edges of WGCN compared to other soft thresholds.This suggests that a WGCN that retains all weights does provide more potential regulatory information in our model.

Manuscript changes:
-Addition of Supplementary Fig 4 : -Addition of the following text: "The selection of a threshold is crucial in WGCNA.However, there is no consensus on the optimal threshold for all tasks, as it depends on the specific situation and analysis objective.In this study, we chose to retain all the correlation coefficients instead of adopting other predefined strategies.This approach was motivated by the expectation that it would capture a broader range of potential regulatory information than other strategies (Fig S4)." 2. While the authors have summarized the six real-cell scRNA-seq datasets, their details and the corresponding ground truth of TF-target interactions still need to be well described, such as the number of TF and target genes in each dataset and the number of links in the corresponding ground truth.

Authors' response:
We thank the reviewer for raising this important point.We have added the detailed clarifications for each real-cell scRNA-seq dataset, including the number of TFs and target genes, the number of links, average degrees per gene and TF, in the supplementary Table 1.

Manuscript changes:
-Addition of Supplementary Table 1 lower compared to the RMSE between the transferred results from B cells and the results of CD14 monocytes trained from scratch, as well as the results of the other two cell types trained from scratch.This analysis suggests that the transferred results from B cells align more closely with the predicted results obtained from training B cells from scratch, rather than being heavily influenced by the training on CD14 monocytes.
Furthermore, we have provided a comprehensive list of the identified regulators from both the transferred model and the model trained from scratch (below and Supplementary Table 8).In the top-5 regulators predicted from the transferring, three regulators were also predicted from scratch, among which PAX5 and EGR1 have already been validated as cell-type-specific regulators for B cells (see Supplementary Table 7).
Additionally, one of the regulators identified through the transferring process, EBF1, has been previously reported as an essential TF for maintaining B cell identity and preventing alternative fates in committed cells [1].

Manuscript changes:
-Addition of Supplementary Fig. 5 (below): -Addition of Supplementary -Addition of the following text: "We conducted further investigations to assess the transferability of our model across different cell types.We trained our model on CD14 Monocytes and subsequently transferred it to three other cell types.And we found that our model presents good transferability across different cell types when comparing the predicted results from transferring to those obtained from each cell type from scratch.
For instance, in the case of B cells, we calculated the Root Mean Square Error (RMSE) between the predicted results obtained from transferring and those derived from training B cells from scratch (Fig. S5).Notably, the RMSE between the two models within B cells was significantly lower compared to the RMSE between the transferred results from B cells and the results of CD14 monocytes trained from scratch, as well as the results of the other two cell types trained from scratch.This analysis suggests that the transferred results from B cells align more closely with the predicted results obtained from training B cells from scratch, rather than being heavily influenced by the training on CD14 monocytes.Furthermore, we have provided a comprehensive list of the identified regulators from both the transferred model and the model trained from scratch (Table S8)." 3. In the TNBC dataset it seems to me that 14,915 interactions are used to train the model, while it predicted 173,399 interactions, which is 10+ times higher.Are the trained interactions all included in the predicted interactions?How confident are the newly predicted interactions?Is there any example of predicted interactions that are not in the training data validated through other experiments?
Authors' response: We thank the reviewer for raising this important point.We generated a Venn diagram (Supplementary Fig 8A ) to visualize the overlap between the inferred GRNs and the interactions used for training the model.About 37.5% of the ground truth interactions (5598 interactions) used for model training were present in the final inferred GRNs.DeepRIG ranked the inferred regulations based on their regulatory scores.The interactions with higher scores are considered more reliable and confident.Furthermore, we displayed the distribution of their regulatory scores using a histogram, providing an overview of the confidence levels of the predicted interactions (Supplementary Fig S8B).We found that most of the newly inferred regulations are assigned the same or higher predicted scores compared to the overlapped interactions, indicating greater confidence in the newly predicted associations.Among the newly predicted interactions, an interesting finding was the regulation of the HIF1α pathway by XBP1.This prediction was validated by Chen et al. [2] through a series of experiments, including invasion assays, subcutaneous xenograft experiments, motif analysis on in-house ChIP-seq data, survival analysis, RT PCR analysis and immunoblot analysis.Furthermore, DeepRIG inferred that ZNF217 regulates the expression of KLF4, FTO, and ALKBH5, and is, in turn, regulated by HIF1α.Remarkably, Zhang et al. [3] performed RT-qPCR analysis, immunoblot assay and immunohistochemical analysis, and demonstrated that ZNF217 cointeracts with ALKBH5, functioning as a mediator to control the expression of KLF4 and NANOG involved in the post-transcriptional mechanism of HIF1α.

Manuscript changes:
-Addition of Supplementary .Furthermore, we found that most of the newly inferred regulations are assigned the same or higher predicted scores compared to the overlapped interactions, indicating greater confidence in the newly predicted associations (Fig S8B).Among the newly predicted interactions, an interesting finding was the regulation of the HIF1α pathway by XBP1.This prediction was validated by Chen et al. [2] through a series of experiments, including invasion assays, subcutaneous xenograft experiments, motif analysis on in-house ChIP-seq data, survival analysis, RT PCR analysis and immunoblot analysis.Furthermore, DeepRIG inferred that ZNF217 regulates the expression of KLF4, FTO, and ALKBH5, and is, in turn, regulated by HIF1α.Remarkably, Zhang et al. [3] performed RT-qPCR analysis, immunoblot assay and immunohistochemical analysis, and demonstrated that ZNF217 co-interacts with ALKBH5, functioning as a mediator to control the expression of KLF4 and NANOG involved in the post-transcriptional mechanism of HIF1α.The findings emphasize the critical role of these altered regulatory interactions in the development and progression of breast cancer." 4. When using simulated and real data for the performance comparison, my understanding is that the underlying networks and Chip-seq validated interactions are used to train deepRIG.This makes the comparison not completely fair since other methods may not be trained on network structures but only on expression data (e.g.deepSEM).This should be addressed or clearly explained in the result section.In addition, the authors may consider using the simulated data to demonstrate what portion of all the interactions should be used to train the model in order for it to make accurate predictions.

Authors' response:
We thank the reviewer for raising this important point.We acknowledge that comparing supervised and unsupervised methods directly may raise concerns about fairness, as they have different access to prior knowledge.It is indeed true that supervised learning methods, such as DeepRIG, CNNC, and 3DCEMA, benefit from the utilization of ground truth ChIP-seq data, which can significantly enhance their performance in GRN inference.We humbly recognize the advantage that prior knowledge provides to supervised methods in terms of performance improvement.
At the same time, it is important to highlight that this advantage does not diminish the value of unsupervised methods.Unsupervised approaches, despite lacking the benefit of ChIP-seq data, still play a crucial role in inferring GRNs and have their own merits in handling diverse biological scenarios where prior knowledge may be limited or unavailable.By presenting the comparison in a balanced and objective manner, we aim to underscore the complementary strengths of both supervised and unsupervised methods in GRN inference.
We have clearly stated and explained it in the result section.

Manuscript changes:
-Addition of following text: "We acknowledge that comparing supervised and unsupervised methods directly may raise concerns about fairness, as they have different access to prior knowledge.It is indeed true that supervised learning methods, such as DeepRIG, CNNC, and 3DCEMA, benefit from the utilization of ground truth ChIP-seq data, which can significantly enhance their performance in GRN inference.We humbly recognize the advantage that prior knowledge provides to supervised methods in terms of performance improvement.At the same time, it is important to highlight that this advantage does not diminish the value of unsupervised methods.Unsupervised approaches, despite lacking the benefit of ChIP-seq data, still play a crucial role in inferring GRNs and have their own merits in handling diverse biological scenarios where prior knowledge may be limited or unavailable." In addition, the authors may consider using the simulated data to demonstrate what portion of all the interactions should be used to train the model in order for it to make accurate predictions.

Authors' response:
We appreciate the valuable suggestion from the reviewer.To ensure what portion of all the interactions should be used to train the model for making more accurate predictions, we further investigated on mESC scRNA-seq dataset.We used different portions of all the ground truth interactions to train our model and utilized the rest to test the model (below and Supplementary Fig 3).We found that DeepRIG reached highest EPR median value when using 67% ground truth interactions for training, compared to other portions.Notably, the reason we used mESC scRNA-seq data instead of simulated data is that the number of ground truth interactions from mESC is largest among all benchmark datasets, which can ensure there is enough much interactions to test the model no matter what portion of ground truth is used for training.

Manuscript changes:
-Addition of Supplementary Fig 3: -Addition of following text: "To ensure what portion of all the interactions should be used to train the model for making more accurate predictions, we further investigated on mESC scRNA-seq dataset.We used different portions of all the ground truth interactions to train our model and utilized the rest to test the model (Fig S3).We found that DeepRIG reached highest EPR median value when using 67% ground truth interactions for training, compared to other portions." 5. In figure 8B and related paragraphs, the authors should briefly explain the meanings of betweenness and pagerank scores and what biological properties do they describe about the genes.

Authors' response:
We thank the reviewer for raising this point.Bwtweenness and pagerank scores are two important network analytics metrics that assess the importance or influence of nodes in a network.Specifically, betweenness centrality measures the intermediacy of a node in the network, while PageRank measures the global importance of a node based on the network's link structure and importance propagation.Genes with high betweenness and PageRank scores are more likely to be considered as important hub genes or regulators.We have added above descriptions into figure 8 and related paragraphs.

Manuscript changes:
-Addition of following text: "Bwtweenness and pagerank scores are two important network analytics metrics that assess the importance or influence of nodes in a network.Specifically, betweenness centrality measures the intermediacy of a node in the network, while PageRank measures the global importance of a node based on the network's link structure and importance propagation.Genes with high betweenness and PageRank scores are more likely to be considered as important hub genes or regulators." Reviewer #3: The While the manuscript introduces a novel method to improve GRN inference, it lacks clarity on how the WGCNA network and the final GRN differ (and whether those differences are biologically meaningful at all).There are also concerns in the way the comparison to other methods is presented, as discussed below.
Major comments: 1.There is no information on how the WGCNA network and the final GRN differ, especially in the real data presented in figures 7 and 8, where most of the networks have a single central node.The authors should presented clear data on how the GRN differs between the "prior regulatory graph" defined by the WGCNA and the final inferred GRN.Please present what nodes and edges change in figure 7 and 8 if one were to reconstruct the network just using the weights of WGCNA.

Authors' response:
We appreciate the valuable suggestion from the reviewer.We explored how the WGCN and the final inferred GRNs differ in the application of PBMC and TNBC in terms of both predictive performance and network topologies, as followed: (1) The difference between WGCN and the final inferred GRNs in inferring performance on six real scRNAseq datasets.We compared the final GRNs inferred by DeepRIG with the straightforward WGCNA model that simply calculated the Spearman's correlation coefficients.As the weights of WGCN range from 0 to 1, we consider those weights bigger than 0.5 as potential regulations, and reconstruct the GRNs from WGCN.We evaluated two models on six real scRNA-seq datasets.The results demonstrated that the final GRNs inferred by DeepRIG present better performance than reconstructed GRNs from WGCN in terms of AUROC, AUPRC, and EPR values (Fig 6).
(2) We also demonstrated the contrasting network topologies between the WGCN and the final inferred GRNs from PBMCs and TNBC.
Specifically, we visualized the WGCN and the final inferred GRNs for CD14+ Monocytes and CD8 T cells in PBMCs, and analyzed the distributions of various topological attributes (below and Supplementary Fig 6).Our observations revealed distinct characteristics of the final inferred GRNs, including the presence of sub-regulatory networks with central regulators, which were not observed in the WGCN.Furthermore, significant differences in the distributions of several topological attributes such as degree, closeness, hub score, and edge betweenness were observed between the WGCN and the final inferred GRNs.The visualization and analysis of these network structures and topological attributes suggest that the final inferred GRNs are more effective in identifying hub genes or key regulators.Notably, in CD14+ Monocytes, we identified several regulators, such as CEBPB and STAT1, that exhibited sharp degree changes between the WGCN and the final predicted GRNs.CEBPB was observed to have 101 degrees in the final inferred GRNs but 0 in the WGCN.By comparing its regulatory structures in the GRNs of four cell types of PBMC, we found that CEBPB exhibits a unique regulatory pattern and complex regulatory structure specific to CD14+ monocytes (Fig 7B).This celltype-specific regulatory structure suggests that CEBPB gene may play a crucial functional role in CD14+ monocytes, as it has been previously demonstrated as a critical TF for monocytic cell proliferation and differentiation [4].Similarly, STAT1, sharp degree changes between the WGCN (0 degree) and the final predicted GRN (178 degrees), has been reported as a critical factor driving peak pSTAT1 levels in CD14+ monocytes [5].In CD8 T cells, we identified JUND as a key regulator, with a degree of 50+ in the final predicted GRN compared to 0 in the WGCN.In fact, JUND is known for its ubiquitous expression during CD8 T cell differentiation [6,7] and it plays an important role in the response of CD8 T cells to acute infection and cancer [8].In B cells, we observed a significant degree change between the WGCN and the final inferred GRN for MEF2C, indicating its central role as a transcription factor in the regulatory network.MEF2C has been reported to regulate the calcineurin signaling pathway and is essential for B cell proliferation [9].The study revealed that MEF2C expression is specific to receptor stimulation in B cells, underscoring its specific regulatory role in receptor-mediated signal transduction.

Manuscript changes:
-Addition of Supplementary Fig 6: -Addition of following text: "We then explored the differences between the final inferred GRNs and the WGCN from CD14 Monocytes and CD8T cells in PBMCs (Fig S6).As expected, they showed contrasting network topologies in visualizations and the distributions of various topological attributes.Our observations revealed distinct characteristics of the final inferred GRNs, including the presence of subregulatory networks with central regulators, which were not observed in the WGCN.Furthermore, significant differences in the distributions of several topological attributes such as degree, closeness, hub score, and edge betweenness were observed between the WGCN and the final inferred GRNs.The visualization and analysis of these network structures and topological attributes suggest that the final inferred GRNs are more effective in identifying hub genes or key regulators." -Addition of Supplementary Fig 9 (and below): -Addition of following text: "We further compared the distributions of network topologies between the final predicted GRNs and WGCN of the TNBC scRNA-seq dataset.The results showed significant differences in node degree, closeness, hub score, and edge betweenness (Fig S9).As expected, the distribution of network topologies from the final inferred GRNs exhibited long-tailed characteristics.Ten potential regulators identified in the final inferred GRNs (Fig 8 ) showed sharp degree changes compared to those in WGCN.In addition, genes including ZNF217, CREBBP, XBP1, RARG, RXRA, E2F4, RAC3, PBX1, and ELF1, also exhibited significant degree changes." 2. All the real data comparisons in Figure 5B are unfair because DeepRIG is a supervised model trained on ground truth ChIP-seq data and then evaluated using ChIP-seq data while most of the other methods are unsupervised and do not use ChIP-seq information.This should be clearly stated in the figure and in the main text for this section.Indeed, CNNC, the other supervised method has really similar performance, bringing the question of whether DeepRIG poses a real practical improvement other than being a supervised method.

Authors' response:
We thank the reviewer for this important comment.We thank the reviewer for raising this important point.We acknowledge that comparing supervised and unsupervised methods directly may raise concerns about fairness, as they have different access to prior knowledge.It is indeed true that supervised learning methods, such as DeepRIG, CNNC, and 3DCEMA, benefit from the utilization of ground truth ChIP-seq data, which can significantly enhance their performance in GRN inference.We humbly recognize the advantage that prior knowledge provides to supervised methods in terms of performance improvement.
At the same time, it is important to highlight that this advantage does not diminish the value of unsupervised methods.Unsupervised approaches, despite lacking the benefit of ChIP-seq data, still play a crucial role in inferring GRNs and have their own merits in handling diverse biological scenarios where prior knowledge may be limited or unavailable.By presenting the comparison in a balanced and objective manner, we aim to underscore the complementary strengths of both supervised and unsupervised methods in GRN inference.
We have clearly stated and explained it in the figure and in the main text for this section.

Manuscript changes:
Addition of following text: "We acknowledge that comparing supervised and unsupervised methods directly may raise concerns about fairness, as they have different access to prior knowledge.It is indeed true that supervised learning methods, such as DeepRIG, CNNC, and 3DCEMA, benefit from the utilization of ground truth ChIP-seq data, which can significantly enhance their performance in GRN inference.We humbly recognize the advantage that prior knowledge provides to supervised methods in terms of performance improvement.At the same time, it is important to highlight that this advantage does not diminish the value of unsupervised methods.Unsupervised approaches, despite lacking the benefit of ChIP-seq data, still play a crucial role in inferring GRNs and have their own merits in handling diverse biological scenarios where prior knowledge may be limited or unavailable." Indeed, CNNC, the other supervised method has really similar performance, bringing the question of whether DeepRIG poses a real practical improvement other than being a supervised method.

Authors' response:
We thank the reviewer for this comment.While CNNC showed comparable performance in terms of AUROC, our model significantly outperformed it in terms of AUPRC and EPR, which are more suitable metrics for evaluating imbalanced sample classification tasks like GRNs inference.Unlike AUROC, which can lead to misleading evaluations in such scenarios, AUPRC and EPR provide a more accurate assessment of the model's performance on the minority class (i.e.positive samples).Through a meticulous comparison of the three supervised learning methods, DeepRIG exhibited superior performance, proving its practical superiority over CNNC and 3DCEMA, particularly in capturing the regulatory relationships with high precision and recall.

Manuscript changes:
Addition of following text: "While CNNC showed comparable performance in terms of AUROC, our model significantly outperformed it in terms of AUPRC and EPR, which are more suitable metrics for evaluating imbalanced sample classification tasks like GRNs inference.Unlike AUROC, which can lead to misleading evaluations in such scenarios, AUPRC and EPR provide a more accurate assessment of the model's performance on the minority class (i.e.positive samples).Through a meticulous comparison of the three supervised learning methods, DeepRIG exhibited superior performance, proving its practical superiority over CNNC and 3DCEMA, particularly in capturing the regulatory relationships with high precision and recall." 3. Although the authors explain in lines 565-569 that they only considered undirected regulatory networks because the ChIP-seq data they use does not specify activation or inhibition, single-cell gene expression networks do indirectly contain that information.No directionality nor inhibitory relationships are considered in the current model.In fact, many of the nodes in the presented networks from Figure 3 have been shown to have mutually inhibitory relationships, such as Gata1 and Pu1 in the HSC network.The authors could improve the potential of their manuscript by commenting/expanding their conclusions to positive and negative regulation, which is widespread in real biological networks.

Authors' response:
We appreciate the valuable suggestion from the reviewer.Among nine baselines, several GRNs inference methods, including PIDC, GRNBOOST2, and DeepSEM, output unsigned regulatory networks, making it challenging to differentiate between activating and inhibitory regulations.Moreover, the ground truth ChIP-seq data also provides unsigned regulatory associations, further complicating the evaluation.To address this limitation, we focused on evaluating unsigned GRNs of three types of ground truth for all baseline methods.
In contrast, our model provides a significant advantage by outputting signed regulatory networks, representing both activating and inhibitory regulations, which closely mirror the complexity of real biological networks.A positive score in the predicted interaction indicates an activating regulation, while a negative score signifies an inhibitory one.To thoroughly assess the performance of DeepRIG, PPCOR, LEAP, SCODE, SCRIBE, and CNNC in predicting activating and inhibitory regulations, we conducted an indepth investigation on four curated networks (Supplementary Fig 2).We designed an approach where activating regulations were assigned label 1, inhibitory regulations were labeled -1, while 0 denotes no regulation.We then calculated the root mean square error (RMSE) between the predicted signed GRNs and the corresponding ground truth labels.The results indicated that DeepRIG achieved the smallest RMSE on HSC, VSC, and GSD networks, and the second-smallest RMSE on the mCAD network.These findings suggest that our model outperformed the other baseline methods in predicting both activating and inhibitory regulations.

Manuscript changes:
-Addition of Supplementary Fig 2 and below: -Addition of following text: "Several GRNs inference methods, including PIDC, GRNBOOST2, and DeepSEM, output unsigned regulatory networks, making it challenging to differentiate between activating and inhibitory regulations.Moreover, the ground truth ChIP-seq data also provides unsigned regulatory associations, further complicating the evaluation.To address this limitation, we focused on evaluating unsigned GRNs of three types of ground truth for all baseline methods.
In contrast, our model provides a significant advantage by outputting signed regulatory networks, representing both activating and inhibitory regulations, which closely mirror the complexity of real biological networks.A positive score in the predicted interaction indicates an activating regulation, while a negative score signifies an inhibitory one.To thoroughly assess the performance of DeepRIG, PPCOR, LEAP, SCODE, SCRIBE, and CNNC in predicting activating and inhibitory regulations, we conducted an indepth investigation on four curated networks (Fig S2).We designed an approach where activating regulations were assigned label 1, inhibitory regulations were labeled -1, while 0 denotes no regulation.
We then calculated the root mean square error (RMSE) between the predicted signed GRNs and the corresponding ground truth labels.The results indicated that DeepRIG achieved the smallest RMSE on HSC, VSC, and GSD networks, and the second-smallest RMSE on the mCAD network.These findings suggest that our model outperformed the other baseline methods in predicting both activating and inhibitory regulations." 4. Line 154: "to better mimic the zero-flame".Are the authors referring to "Zero-inflation"?It would be most relevant to model the number of UMIs detected per single-cell too.The number of UMIs (or number of genes, which is positively correlated with the former) is, together with the number of single-cells, an important driver of the sparsity in single-cell datasets.

Authors' response:
We thank the reviewer for bringing the spelling mistake to our attention.We apologize for the error and appreciate the reviewer's clarification regarding the term "Zero-inflation".As mentioned by the reviewer, sparsity is a prominent feature of single-cell datasets.Due to the inherent heterogeneity and sparsity among cells, many genes may not be detected in certain cells or may have very few unique molecular identifiers (UMIs) detected.This leads to the presence of zero values in the single-cell data, indicating the zero-inflation characteristic of the dataset.Therefore, in order to accurately model the number of UMIs detected per cell, it is essential to account for the zero-inflation characteristics observed in single-cell data.
5. Figure 4. Please change the legend from "no. of cells" to "Dropout rate".

Authors' response:
We appreciate the reviewer for bringing to our attention the mistake in Fig 4 .We sincerely apologize for the oversight and have promptly made the necessary correction, revising the legend to "Dropout rate".

Manuscript changes:
Modification of Fig 4.
6.The numbers of cells in figure 5C is very low for current datasets.The authors should benchmark or show how their method scales up this with tens of thousands or hundreds of thousands of cells, which are now routine in single-cell experiments.Similarly, the number of genes detected in an experiment can reach 8000 or more.

Authors' response:
We thank the reviewer for raising this important point.We applied our model to the whole PBMC8k dataset with 8339 cells to demonstrate its scalability in handling larger cell numbers.While we acknowledge that the PBMC8k dataset is not in the tens of thousands or millions, this analysis allowed us to evaluate the time complexity and scalability of our approach and other baseline methods on relatively larger datasets.Additionally, generating a suitable ground truth for data sets of tens of thousands or millions of cells is somewhat challenging.We also compared our model to PIDC, PPCOR, GRNBOOST2, DeepSEM, CNNC which do not require peseudotime or stamped information.As expected, PPCOR only took 1 minute to predict GRNs, and DeepRIG and DeepSEM completed the GRN inference in just 4 minutes and 10 minutes, respectively.On the other hand, PIDC and CNNC required 60 to over 100 minutes to generate the outputs, while GRNBOOST2 took thousands of minutes.Regarding gene scalability, we evaluated our model and five baselines on the TNBC scRNA-seq dataset containing 7185 genes.Notably, PPCOR completed predictions in 22 minutes, while DeepRIG and DeepSEM took 61 and 140 minutes, respectively, for GRN inference.GRNBOOST2 required 1383 minutes, and PIDC failed to provide predictions.This comparison highlights the scalability of our model in handling datasets with a large number of cells and genes.Additionally, we observed that the time complexity of most methods is primarily influenced by the number of genes.

Manuscript changes:
-Addition of following text: "We further applied our model to the whole PBMC8k dataset with 8339 cells to demonstrate its scalability in handling larger cell numbers.We also compared our model to PIDC, PPCOR, GRNBOOST2, DeepSEM, CNNC which do not require peseudotime or stamped information.As expected, PPCOR only took 1 minute to predict GRNs, and DeepRIG and DeepSEM completed the GRN inference in just 4 minutes and 10 minutes, respectively.On the other hand, PIDC and CNNC required 60 to over 100 minutes to generate the outputs, while GRNBOOST2 took thousands of minutes.This comparison clearly demonstrated the scalability of our model when dealing with datasets containing a large number of cells.
Regarding gene scalability, we evaluated our model and five baselines on the TNBC scRNA-seq dataset containing 7185 genes.Notably, PPCOR completed predictions in 22 minutes, while DeepRIG and DeepSEM took 61 and 140 minutes, respectively, for GRN inference.GRNBOOST2 required 1383 minutes, and PIDC failed to provide predictions.This comparison highlights the scalability of our model in handling datasets with a large number of cells and genes.Additionally, we observed that the time complexity of most methods is primarily influenced by the number of genes." 7. The authors should also include the results from the other network inference methods in figure 7, and comment on how they differ.

Authors' response:
We thank the reviewer for raising this important point.We compared the performance of several other GRN inference methods, including PIDC, PPCOR, GRNBOOST2, DeepSEM, and CNNC, with our proposed DeepRIG method on CD14 Monocytes of PBMCs.Visualizing and analyzing the network structures of the inferred GRNs from these methods, as well as DeepRIG, revealed interesting findings (below and Supplementary Fig. 7).We observed that the GRNs inferred by CNNC and GRNBOOST2 exhibited similarities to the predictions of our model, albeit with a more dispersed network structure.On the other hand, the GRNs predicted by the correlation-based methods, PIDC and PPCOR, were found to be more akin to the WGCN.Notably, significant differences were observed among these methods in terms of topological attributes, such as degrees, closeness, hub scores, and edge betweenness.A particularly intriguing observation was that the inferred GRNs from our model and CNNC exhibited a power-law distribution, resembling scale-free networks.In scale-free networks, a few key hub nodes have rich connections to other nodes, while most nodes have limited connections.This distribution pattern is commonly observed in various real biological networks, including metabolic networks [14,15], proteinprotein interaction networks [16], and transcriptional regulation networks [17].Our analysis suggests that the GRNs inferred by our model better capture the regulatory relationships present in the real data and are more effective in identifying hub genes or regulators.

Manuscript changes:
-Addition of supplementary Fig 7: -Addition of the following text: "We compared the performance of several other GRN inference methods, including PIDC, PPCOR, GRNBOOST2, DeepSEM, and CNNC, with our proposed DeepRIG method on CD14 Monocytes of PBMCs.Visualizing and analyzing the network structures of the inferred GRNs from these methods, as well as DeepRIG, revealed interesting findings (Fig. S7).We observed that the GRNs inferred by CNNC and GRNBOOST2 exhibited similarities to the predictions of our model, albeit with a more dispersed network structure.On the other hand, the GRNs predicted by the correlation-based methods, PIDC and PPCOR, were found to be more akin to the WGCN.Notably, significant differences were observed among these methods in terms of topological attributes, such as degrees, closeness, hub scores, and edge betweenness.A particularly intriguing observation was that the inferred GRNs from our model and CNNC exhibited a power-law distribution, resembling scale-free networks.In scale-free networks, a few key hub nodes have rich connections to other nodes, while most nodes have limited connections.This distribution pattern is commonly observed in various real biological networks, including metabolic networks, proteinprotein interaction networks, and transcriptional regulation networks.Our analysis suggests that the GRNs inferred by our model better capture the regulatory relationships present in the real data and are more effective in identifying hub genes or regulators." 8. Why are only two TFs discussed in figure 7? How about the other TFs expressed in this dataset?
Authors' response: We appreciate the reviewer's valuable comment.Our approach focuses on uncovering the distinct gene regulation mechanisms exhibited by cells with varying transcription levels.In Figure 7, our goal was to provide cell-type-specific GRN inference and identify regulators that are specific to each cell type.In the manuscript, we highlighted two regulators, CEBPB and SPI1, which were found to exhibit cell-type-specific regulatory structures in CD14 Monocytes, as well as two regulators, MEF2C and EGR1, in B cells.
We acknowledge the reviewer's suggestion that additional discussion is needed to cover more cell-typespecific regulators identified by our model, as they present significant differences in regulatory structures across cell types.Therefore, we have included introductions of several additional cell-type-specific regulators identified by our model, such as GATA2 We have also provided a comprehensive list of all identified cell-type-specific regulators in the corresponding predicted GRNs (Supplementary Table 7).
We would like to express our gratitude for the reviewer's insightful suggestion, which has enhanced the clarity and completeness of our manuscript.

Manuscript changes:
-Addition of the following text: "We have also identified GATA2 and STAT1 as two specific regulators that exhibit cell-type-specific regulatory structures in CD14 Monocytes.GATA2 demonstrated 201 degrees in CD14 Monocytes, while only 3 degrees in B cells.Furthermore, GATA2 exhibited the highest Betweenness Centrality value (0.59), indicating its significant influence and control over other genes.GATA2 is recognized as an important TF involved in regulating various biological processes, including ontogeny in monocytes.Similarly, STAT1, which showed sharp degree changes across cell types, has been reported as a critical factor driving peak pSTAT1 levels in CD14+ monocytes.In the inferred GRNs of B cells, two TFs, PAX5 and RUNX3 were also identified as a cell-type-specific regulator by our model.Indeed, the wellestablished transcription factor PAX5 plays a pivotal role in governing the identity and function of B cells throughout B lymphopoiesis and has been associated with human B cell malignancies.RUNX3 is known for its cross-regulation mechanism with RUNX1 in B cells, and its expression is required for efficient B cell proliferation.The comprehensive list of identified regulators specific to each cell type can be found in Table S7." -Addition of Supplementary Table 7.

Fig 8 (
below): -Addition of the following text: "We then investigated the overlap between the inferred GRNs and the interactions used for training the model (Fig S8).About 37.5% of the ground truth interactions (5598 interactions) used for model training were present in the final inferred GRNs (Fig S8A) :

Table 8
study introduces DeepRIG, a graph-based deep learning model for inferring gene regulatory networks (GRNs) from single-cell RNA sequencing data.Unlike traditional methods that only consider pairwise gene regulatory relationships, DeepRIG captures the global regulatory structure by transforming gene expression data into a co-expression mode to build a prior regulatory graph.It then employs a graph autoencoder model to embed the global regulatory information into gene latent embeddings and reconstruct the GRN.Benchmarking results show that DeepRIG outperforms existing methods in accurately reconstructing GRNs on completely synthetic and simulated regulatory networks.The model was applied to human peripheral blood mononuclear cells and triple-negative breast cancer samples.