Deciphering Cell Lineage Gene Regulatory Network via MTGRN

We appreciate the reviewer’s high-quality comments. Below, we will address each of the questions raised.

1. In response to the reviewer's question that NicheNet and the ground truth network have a significant overlap. It is a very important point, we will discuss below. First, we proposed a method named Random in Table 1 for baseline comparison. In Random, we randomly selected the top-k edges (where k is equal to the number of edges in the ground truth) from the prior knowledge as the inferred GRN. As for results, we can find that the AUROC for Random across five datasets is around 0.5, indicating there is no significant overlap between the prior knowledge and the ground truth. Second, we analyzed the overlap between the edges provided by NicheNet and the edges in ground truth. The result is presented below:

   • Prior Knowledge: 5318181 edges
   • mESC ground truth network: 24557 edges (overlap: 16695)
   • mHSC-E ground truth network: 24726 edges (overlap: 14732)
   • mHSC-GM ground truth network: 16198 edges (overlap: 9818)
   • mHSC-L ground truth network: 4705 edges (overlap: 2494)

   we can find that NicheNet includes 5318181 edges, while the ground truth edge counts in our datasets are 24557, 24726, 16198, and 4705, with intersections of 16695, 14732, 9818, and 2494 edges, respectively. In each dataset, the effective regulatory edges provided by NicheNet account for only 0.31%, 0.27%, 0.18%, and 0.04% of the entire prior knowledge. This indicates that NicheNet provides a coarse prior network containing considerable noise, including potential false-positive regulatory relationships from databases or experiments. This is also why the AUROC for random is around 0.5 in Table 1. From above, there is no significant overlap between NicheNet and the ground truth network, which is why CEFCON [1] also uses it as prior knowledge. MTGRN should identify approximately 10,000 true regulatory edges from a noisy network of over 5 million edges, making this task challenging. For the second question about GENIE3 and GRNBoost2, we chose to compare our model with them for two key reasons. First, these methods are considered classic in GRN inference, so comparing against them allowed us to benchmark our model's performance against established standards. Second, this comparison emphasize incorporating prior knowledge can enhance GRN inference. Except for GENIE3 and GRNBoost2, the other three baseline methods (NetREX, CEFCON, and Celloracle) all rely on prior networks, in comparison with these three methods, our model consistently demonstrated significant performance gains, showing the effectiveness of our method. If reviewer has any recommendations for models suitable for comparison, we would be glad to compare our model with them as well

2. In selecting the top k edges, we did not incorporate any gene enrichment information. Using the mESC dataset as an example, MTGRN needs to select the top 24557 edges out of a total of 5318181 edges available in the NicheNet network, the ground truth edges is a tiny fraction of the entire search space, which includes numerous false-positive regulatory interactions that could interfere with network inference. Due to this noise, some unrelated genes may appear more enriched in the prior network than actual regulatory factors. Because MTGRN will predict the gene expression for future M time points based on the former N time points, if an unrelated gene highly enriched in the prior network, it will not effectively fit the expression in the M future cells, so during attention map computation, our model will assigns lower scores to such genes, meaning that even if a gene appears highly enriched, it will still receive a low score if it does not contribute to accurate predictions (beacause they are not the true regulatory factors for the target genes). We can refer to the method Random in Table 1 for a clearer explanation. Since Random randomly selects edges from the prior network, it is more likely to pick enriched genes. However, it achieved an AUROC close to 0.5 across all five datasets. If our model’s edge selection were also based on gene enrichment, our scores would be similar to those of Random. Instead, MTGRN’s accuracy across all five datasets is significantly higher than Random, demonstrating that our selection of the top-k edges is not based on gene enrichment.

[1] Wang P, Wen X, Li H, et al. Deciphering driver regulators of cell fate decisions from single-cell transcriptomics data with CEFCON[J]. Nature Communications, 2023, 14(1): 8459.

Due to the character limit, we will continue addressing the remaining questions in the next comment.

3. In response to reviewer’s comment that several variables are not defined in the paper. We sincerely apologize for the oversight in explaining some of the symbols, we will provide a detailed explanation below. $X_{\text{input}} \in \mathbb{R}^{G \times W \times d}$ will go through three independent linear transformation to obtain Q, K, V, they are all have shape of $G \times W \times d$. To calculate the attention matrix, we use the formula:$\frac{Q K^\top}{\sqrt{d}}$, the result is an attention matrix with dimensions $G \times W \times W$. As for the second quesiton, in Section 3.3 of the paper, we mention that

the output of temporal block $X_{output} \in \mathbb{R}^{G \times W \times d_{model}}$ will be transposed to $X_{output} \in \mathbb{R}^{W \times G \times d_{model}}$, which will be input into Spatial attention module alongside prior knowledge

The difference in dimensions between the spatial attention module and the temporal attention module, as mentioned by the reviewer, is due to a transposition we applied to the output of the temporal module. We switched the positions of G and W, allowing us to compute the attention between genes in the spatial attention module. Consequently, the attention map in the spatial attention module has dimensions of $G \times G$. 4. In response to reviewer’s comment that already have work using transformers for GRN inference, we want to emphasize that MTGRN’s contribution is not applying transformer to infer GRN (we did not mention this as a contribution in our introduction). Our key contribution is reframing GRN inference as a multivariate time series forecasting problem. MTGRN transforms single-cell data into time sequence data using trajectory inference and then infers GRNs through a time-series prediction approach. The proposed time and spatial modules are designed to learn Granger causality [1] in the cell development process. In the time module, each cell can only attend to cells from prior time points, allowing each gene to observe all genes in preceding cells. However, only a small subset of genes influences the expression of a given target gene. To address this, we use the spatial module to filter out gene pairs without regulatory relationships in the prior knowledge, ensuring that each gene’s expression is inferred based only on historical expression data of genes with known interactions in prior network. Transformer is merely a tool we use for GRN inference, our main contribution lies in restructuring the GRN inference task as a new paradigm of multivariate time series prediction problem. 5. In response to question that why not use time-series scRNA-seq datasets where the time points are given rather than learned? The datasets we use for benchmaring come from BEELINE [2]. For example, in the mESC dataset, each cell is named in a format like “RamDA_mESC_00h_C01,” which indicates its collection time, here at 0 hours. The mESC data were collected at five distinct time points (0, 12, 24, 48, and 72 hours). As for why we don’t use these specific time points directly, there are three main reasons. First, in each collection time point, multiple cells are sequenced to obtain their gene expression data, so these cells will have the same time label (e.g., 00h), which prevents us from establishing a precise temporal order among them. Second, we think cells collected at the same time still have an inherent time order. By assigning a pseudotime to cells collected at same time point, we can establish a more accurate order within each time point. Third, previous methods such as CellOracle [3] inferred GRN by clustering cells of the same type and then constructed a GRN for each cluster. However, we believe that even within the same cell type, there is an inherent developmental progression. Therefore, we use trajectory inference to assign pseudotimes to cells, allowing us to capture this developmental order more effectively.

We are very grateful for your valuable suggestions and hope our responses clarify any concerns. We greatly value this opportunity to discuss our work. If there are any further questions, please feel free to ask, and we will be happy to provide further details. We hope this discussion will help improve the score of our paper, and thank you again!

[1] Shojaie A, Fox E B. Granger causality: A review and recent advances[J]. Annual Review of Statistics and Its Application, 2022, 9(1): 289-319.

[2] Pratapa A, Jalihal A P, Law J N, et al. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data[J]. Nature methods, 2020, 17(2): 147-154.

[3] Kamimoto K, Stringa B, Hoffmann C M, et al. Dissecting cell identity via network inference and in silico gene perturbation[J]. Nature, 2023, 614(7949): 742-751.

We thank the reviewer for their recognition of our paper and for the valuable feedback provided. We sincerely appreciate these insights! Below, we will address each of the questions raised.

1. In response to reviewer's comments that how did the predicted gene expression metrics such as spearmanR or pearsonR look like? We conducted the relevant experiments and included the results in the supplementary material. The results demonstrate that MTGRN accurately models gene expression, with Spearman R and Pearson R scores both exceeding 0.98.
2. In response to reviewer's comments that does the model understand genes that are co-regulated by multiple transcription factors? We sincerely thank the reviewer for raising such an excellent question, which will guide the next steps in optimizing MTGRN. Currently, we are sorry that MTGRN has not delved deeply into this aspect, but we wiil try this in the future.
3. In response to reviewer's comments that Does the model show any new gene-gene interactions? Actually, MTGRN is capable of handling this task. Once trained, during the inference stage, we can obtain the embedding for each gene before calculating the attention map. By simply computing the similarity between gene embeddings, we can predict new gene interactions. After obtaining the gene embeddings, of course more complex methods can be employed to predict new gene interactions, we have proposed a straightforward approach to demonstrate that MTGRN is capable of handling this task.

We want to emphasize that we have also included the results of MTGRN identifying dynamic GRNs in the supplementary material. In fact, MTGRN is able to infer dynamic GRN, as we assign each cell a pseudotime and use time-series prediction to infer the GRN. During inference, we can segment cells along the differentiation trajectory into different time segments (e.g., grouping every $n$ consecutive cells into one segment). By inputting cells from each segment into MTGRN, we can obtain the GRN for that specific time segment. For the entire differentiation trajectory, assuming we generate $m$ time segment, we can apply Gaussian smoothing to the regulatory edge scores across these $m$ GRNs to construct a dynamic GRN. This approach is similar to Dictys [1].

We are very grateful for your valuable suggestions and hope our responses clarify any concerns. We greatly value this opportunity to discuss our work. If there are any further questions, please feel free to ask, and we will be happy to provide further details. We hope this discussion will help improve the score of our paper, and indeed we will cite the extraordinary paper reviewer mentioned in comments, thank you again!

[1] Wang L, Trasanidis N, Wu T, et al. Dictys: dynamic gene regulatory network dissects developmental continuum with single-cell multiomics[J]. Nature Methods, 2023, 20(9): 1368-1378.

We carefully read the article “Constructing the dynamic transcriptional regulatory networks to identify phenotype-specific transcription regulators” which focuses on gene temporal dynamics. This approach highlights how dynamic characteristics of transcription regulators can reveal phenotype-specific transcription factors (TFs) and pathways, addressing limitations in static TRN models. The framework’s use of graph autoencoders and statistical methods for identifying dynamic interactions let me learn a lot and we cited it in the Related Work section in our latest submitted version.

the reference links are listed below:

[1] Pratapa A, Jalihal A P, Law J N, et al. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data[J]. Nature methods, 2020, 17(2): 147-154.

[2] Nestorowa S, Hamey F K, Pijuan Sala B, et al. A single-cell resolution map of mouse hematopoietic stem and progenitor cell differentiation[J]. Blood, The Journal of the American Society of Hematology, 2016, 128(8): e20-e31.

[3] Hayashi T, Ozaki H, Sasagawa Y, et al. Single-cell full-length total RNA sequencing uncovers dynamics of recursive splicing and enhancer RNAs[J]. Nature communications, 2018, 9(1): 619.

[4] Camp J G, Sekine K, Gerber T, et al. Multilineage communication regulates human liver bud development from pluripotency[J]. Nature, 2017, 546(7659): 533-538.

[5] Wang P, Wen X, Li H, et al. Deciphering driver regulators of cell fate decisions from single-cell transcriptomics data with CEFCON[J]. Nature Communications, 2023, 14(1): 8459.

[6] Kamimoto K, Stringa B, Hoffmann C M, et al. Dissecting cell identity via network inference and in silico gene perturbation[J]. Nature, 2023, 614(7949): 742-751.

[7] Wang L, Trasanidis N, Wu T, et al. Dictys: dynamic gene regulatory network dissects developmental continuum with single-cell multiomics[J]. Nature Methods, 2023, 20(9): 1368-1378.

We are very grateful for your valuable suggestions and hope our responses clarify any concerns. We greatly value this opportunity to discuss our work. If there are any further questions, please feel free to ask, and we will be happy to provide further details. We hope this discussion will help improve the score of our paper, and thank you again!

We sincerely appreciate the reviewer’s thoughtful comments. we will address each of the questions below.

1. In response to reviewer's comments that the experimental analysis is based on data from different cell lines rather than dynamic or developmental data. We want to emphasize that all the data we use is dynamic developmental data, reviewer could refer to BEELINE [1] for more context. It provide detailed descriptions of the five datasets and we summarize the main details in BEELINE below :

• mHSC dataset is sourced from [2] and includes 1,656 hematopoietic stem and progenitor cells (HSPCs), which can be divided into three lineages: mHSC-E (erythroid lineage), mHSC-L (lymphoid lineage), and mHSC-GM (myeloid lineage). Each lineage includes all cells in the progression from the starting cell to the endpoint cell.
• mESC dataset is sourced from [3], contains scRNA-seq results from 421 cells that track the development of mouse embryonic stem cells into primitive endoderm cells. These cells were collected at five distinct time points: 0, 12, 24, 48, and 72 hours.
• hHep dataset is from [4] and includes scRNA-seq results from an experiment where induced pluripotent stem cells (iPSCs) were differentiated into hepatocyte-like cells. This dataset includes 425 scRNA-seq measurements taken at various time points: day 0, day 6, day 8, day 14 and day 21.

In summary, all five datasets we use involve dynamic data related to cell development.

2. In response to reviewer's comments that we did not compare with latest state-of-the-art methods. In fact, the comparison methods CEFCON [5] and CellOracle [6] are GRN inference models published in Nature Communications and Nature in 2023, which, in our view, represent the best-performing models in the GRN inference filed. If reviewer has any recommendations for models suitable for comparison, we would be glad to compare our model with them as well.
3. In response to reviewer's comments that a more thorough explanation of the model's architecture and hyperparameter settings would help. We will discuss the model and hyperparameter settings below. The main contribution of MTGRN is to infer GRNs through multivariate time-series prediction, where single-cell sequencing data is transformed into time-series data, and GRNs are inferred by predicting gene expression in the next M time points based on the gene expression in the previous N time points. MTGRN is composed of temporal attention module and spatial attention module. In the time module, each cell can only attend to cells from prior time points, allowing each gene to observe all genes in preceding cells. However, only a small subset of genes influences the expression of a given target gene. To address this, we use the spatial module to filter out gene pairs without regulatory relationships in the prior knowledge, ensuring that each gene’s expression is inferred based only on historical expression data of genes with known interactions in prior network. For hyperparameter settings, we applied grid search to optimize parameters, and the final hyperparameter values are as follows:

• input_length (number of cells in the previous N time points): 16
• predict_length (number of cells in the following M time points): 16
• d_model (dimension of Transformer): 128
• d_ff (dimension in feedforward): 512
• heads (number of attention heads in Transformer): 4

Additional configurations, such as GPU setting, learning rate, and the number of training epochs, are described in detail in Section 4 under Reproducibility paragraph, reviewer could refer to this section for more information. We hope our response solve your questions.

4. In response to reviewer's comments that we should incorporate dynamic gene expression data to infer dynamic networks. In fact, MTGRN is able to infer dynamic GRN, as we assign each cell a pseudotime and use time-series prediction to infer the GRN. During inference, we can segment cells along the differentiation trajectory into different time segments (e.g., grouping every $n$ consecutive cells into one segment). By inputting cells from each segment into MTGRN, we can obtain the GRN for that specific time segment. For the entire differentiation trajectory, assuming we generate $m$ time segment, we can apply Gaussian smoothing to the regulatory edge scores across these $m$ GRNs to construct a dynamic GRN. This approach is similar to Dictys [7]. We did not include this experiment in the main text of our submission, but we have now uploaded the results as supplementary material. Reviewer could examine these results in supplementary material, hoping this clarifies your questions!

Due to the character limit, we will give the reference links in the next comment.

We greatly value the feedback provided by reviewers and sincerely apologize for any lack of clarity in presenting our work during the initial submission, which may have caused confusion. In the rebuttal phase, we have addressed the issues you raised in detail and the additional experiments you mentioned have been included in the supplementary material for your reference. We hope these clarifications and supplementary analyses will help you evaluate the significance of our study comprehensively.

We are pleased to note that during the rebuttal phase, we successfully addressed the concerns of other reviewers, which led to a positive reevaluation and corresponding score improvements. We genuinely hope our detailed and thoughtful response will also resolve your concerns and help improve the score of our paper. We sincerely look forward to your feedback and thank you again!

We greatly value the feedback provided by reviewers and sincerely apologize for any lack of clarity in presenting our work during the initial submission, which may have caused confusion. In the rebuttal phase, we have addressed the issues you raised in detail and the additional experiments you mentioned have been included in the supplementary material for your reference. We hope these clarifications and supplementary analyses will help you evaluate the significance of our study comprehensively.

We are pleased to note that during the rebuttal phase, we successfully addressed the concerns of other reviewers, which led to a positive reevaluation and corresponding score improvements. We genuinely hope our detailed and thoughtful response will also resolve your concerns and help improve the score of our paper. We sincerely look forward to your feedback and thank you again!

We are looking forward to receive feedback from the reviewer G2id.

We sincerely appreciate the reviewer’s thoughtful comments. we will address each of the questions below.

1. Regarding the question of whether prior knowledge plays a significant role in our model's performance, actually we proposed a method named Random in Table 1 for baseline comparison. In Random, we randomly selected the top-k edges (where k is equal to the number of edges in the ground truth) from the prior knowledge as the inferred GRN. As for results, we can find that the AUROC for Random across five datasets is around 0.5, significantly lower than that of our model (MTGRN), indicating that prior knowledge alone is not the reason for our model’s high performance.

2. Methods we compared against in Table 1 also rely on prior knowledge to infer GRNs. For example, CellOracle [1] uses ATAC-seq data to construct a base GRN, providing prior knowledge of TF-target relationships. CEFCON [2] also relies on NicheNet as prior information, applying a graph neural network to infer the GRN. For fairness, the baselines we chose also incorporate prior knowledge, and our model achieved higher performance in these comparisons, further supporting that our model’s performance is not due to the use of prior knowledge.

3. We analyzed the overlap between the edges provided by NicheNet and the edges in ground truth. The result is presented below:

   • Prior Knowledge: 5318181 edges
   • mESC ground truth network: 24557 edges (overlap: 16695)
   • mHSC-E ground truth network: 24726 edges (overlap: 14732)
   • mHSC-GM ground truth network: 16198 edges (overlap: 9818)
   • mHSC-L ground truth network: 4705 edges (overlap: 2494)

   we can find that NicheNet includes 5318181 edges, while the ground truth edge counts in our datasets are 24557, 24726, 16198, and 4705, with intersections of 16695, 14732, 9818, and 2494 edges, respectively. In each dataset, the effective regulatory edges provided by NicheNet account for only 0.31%, 0.27%, 0.18%, and 0.04% of the entire prior knowledge. This indicates that NicheNet provides a coarse prior network containing considerable noise, including potential false-positive regulatory relationships from databases or experiments. This is also why the AUROC for random is around 0.5 in Table 1. MTGRN should identify approximately 10,000 true regulatory edges from a noisy network of over 5 million edges, making this task challenging and countering the concern that our high performance is due to the strength of prior knowledge (actually the information provided by prior is little).

4. In response to the reviewer's question that what are the technical aspects of the model architecture that contribute most to the performance? we emphasize that MTGRN’s contribution is reframing GRN inference as a multivariate time series forecasting problem, learning Granger causality [3] in cell development. Our model’s ability stems from effectively utilizing the temporal information inherent in cell development. As shown in Figure 2, we first apply a Temporal Attention Module that allows each cell to focus on gene expression information from earlier developmental stages. At this module, each gene of the cell will attend the expression of all genes. However, each gene is typically influenced by only a subset of regulatory factors. We therefore apply the Spatial Attention Module, which leverages prior knowledge to filter out gene pairs without regulatory relationships (although this filtering is somewhat coarse due to noise in the prior knowledge). Combining these two modules allows us using the historical expression information of candidate TFs to predict a gene’s future expression. This is the core factor improving MTGRN’s performance.

5. In response to the reviewer's question that what aspects should others try to build on? we believe building GRNs using methods related to RNA velocity or single-cell trajectory inference could yield additional insights. Incorporating other modalities, such as ATAC-seq data like CellOracle, could also strengthen the accuracy of prior knowledge, which avoids the noise present in NicheNet.

We are very grateful for your valuable suggestions and hope our responses clarify any concerns. We greatly value this opportunity to discuss our work. If there are any further questions, please feel free to ask, and we will be happy to provide further details. We hope this discussion will help improve the score of our paper, and thank you again!

[1] Wang P, Wen X, Li H, et al. Deciphering driver regulators of cell fate decisions from single-cell transcriptomics data with CEFCON[J]. Nature Communications, 2023, 14(1): 8459.

[2] Kamimoto K, Stringa B, Hoffmann C M, et al. Dissecting cell identity via network inference and in silico gene perturbation[J]. Nature, 2023, 614(7949): 742-751.

[3] Shojaie A, Fox E B. Granger causality: A review and recent advances[J]. Annual Review of Statistics and Its Application, 2022, 9(1): 289-319.

We would like to clarify the perturbation experiments mentioned by reviewer in weakness section. Actually, it serve as an important technical validation of the inferred gene regulatory network. Similar approaches can be found in studies like Celloracle [1]. In fact, without perturbation analysis, it would be difficult to validate the inferred network’s accuracy. It is through the perturbation of key factors that the robustness of the network can be definitively proven.

In Figure 7 (a), our TF perturbation analysis is entirely based on the GRN inferred by our model (MTGRN). If the inferred network were incorrect, the calculated offset vectors for each cell would deviate in the wrong direction. However, by performing perturbations on the network constructed by MTGRN, we observed that knocking out Gata1 (a well-known transcription factor promoting erythroid development) causes the offset vectors of developing cells to shift in the opposite direction of the developmental trajectory. This strongly supports the accuracy of our model’s GRN inference and serves as a robust technical validation. It is an important technical discussion for analyzing model performance, and we do not fully agree with the reviewer’s perspective on this point. We hope to have further discussion.

We greatly appreciate the opportunity for this rebuttal and hope our responses address the reviewer’s concerns. Thanks once again!

[1] Kamimoto K, Stringa B, Hoffmann C M, et al. Dissecting cell identity via network inference and in silico gene perturbation[J]. Nature, 2023, 614(7949): 742-751.

We greatly appreciate the reviewers’ feedback and hope they respect others’ work with careful consideration. During the rebuttal phase, we value constructive and evidence-based discussions rather than dismissive or careless remarks. We have thoroughly addressed all the concerns you raised, and if you have further questions, we kindly request your response. Thank you!