Learning to Discover Regulatory Elements for Gene Expression Prediction

Thanks for your valuable comments, and we understand your concern about reproduction. Therefore, we provide our implementation in the following link for the reproduction: https://anonymous.4open.science/r/Seq2Exp-1E7E/README.md The datasets are currently too large to be hosted on an anonymous GitHub repository or provided as a zip file. Therefore, we only release our code at rebuttal. And we promise that the datasets and codes will be made publicly available after the paper is open. In addition, we provide our replies as follows.

W1 Quality & Q1:

• W1 Quality: We provide our implementation in the link above.
• W1 Quality & Q1: For the hyperparameters, we summarize them in the Appendix A.4 in our revised paper. Specifically, we set α3​ and β3​ to be 10 and 90, respectively. In practice, rather than directly assigning these values, we first set the mean of the beta distribution to 0.1, for example, α=1 and β=9. We then scale these parameters by a constant factor of 10, resulting in α=10 and β=90.
• W1 Quality: For the running time, we report the running time for one iteration of all the methods in the table below. For EPInformer, it uses part of the sequence and dilated CNN which can conduct an efficient training, and they only focus on predefined regions. And token-level sequence models like Enformer and Caduceus need longer time to encode the full sequence. Seq2Exp needs longer training time, since our framework needs to learn the positions of regulatory elements, i.e., a generator and predictor. Therefore, the training time is doubled to discover the relevant sub-sequences for the prediction.

W2 Clarity & Q2 & Q3. Thanks for pointing it out.

• W2 Clarity & Q3: We have added a complete and precise description about these mRNA Half-life features and promoter activity in the Appendix A.3 highlighted as red to make it clear. These features, previously used in studies (Lin et al., 2024), are included in our implementation and baselines to ensure a fair comparison. Specifically, these features are incorporated into the final output linear layer of the gene expression predictor $g_{\phi}$​, concatenated with the hidden features as input of the final layer.
• W2 Clarity & Q2: Hi-C data is inherently represented as a matrix with a resolution of 5000 base pairs. Here, we specifically focus on the TSS, meaning we are only interested in the links between the bin containing the TSS and other bins. Although the resolution is 5000 bp, we still obtain the contact frequencies between each individual base pair and the TSS base pair. However, the contact values for a specific base pair will be identical to those of the surrounding base pairs within the same 5000 bp bin. Here we follow the implementations of ABC model [1] and EPInformer [2].

W3 Significance: Thanks for providing your concerns.

Regarding the independence of mask selection between base pairs, we assume independence to facilitate the derivation of the information bottleneck and KL loss in the current framework. However, exploring "lumpy masks," which consider dependencies across base pairs, could be an interesting direction for future work. This approach has the potential to enhance the discovery of regulatory elements by capturing more complex patterns and interactions.

If you ask about the Assumption 1: DNA sequences and signals are conditional independent to each other given the ground truth regulatory elements, this assumption could also benefit from further discussion. If the ground truth positions of CpG islands are known, the observed methylation peaks can be understood as a result of these positions. Moreover, DNA sequences can be viewed as a combination of regulatory and non-regulatory elements, leading to a scenario where signal measurements are independent of the DNA sequences themselves.

Q4. Learning α2​ and β2​ could be an interesting direction for future work. For simplicity in our current approach, we set α2​ equal to the signal value and β2​ as a constant to incorporate the signal information. We also experimented with making β2​ learnable parameters but did not observe much differences. Future work could explore more effective ways to leverage the signal values.

[1] Activity-by-contact model of enhancer–promoter regulation from thousands of crispr perturbations
[2] EPInformer: a scalable deep learning framework for gene expression prediction by integrating promoter-enhancer sequences with multimodal epigenomic data

Thank you for your response and valuable comments. We understand your concern about the training time, which is primarily due to our framework using the long DNA sequence. DNA language models, such as Caduceus and HyenaDNA, have demonstrated their effectiveness across various downstream tasks while maintaining reasonable training time with long DNA sequence. However, the integration of these models with epigenomic signals remains underexplored. Incorporating such signals could provide valuable insights and pave the way for powerful models capable of handling multiple cell types in future research. Therefore, while our method requires higher computational costs due to its use of long DNA sequence, we believe it is crucial to develop a learnable way to consider these epigenomic signals modeling to discover the regulatory elements. Moreover, developing more efficient and trainable methods to achieve this remains an important and promising direction for further exploration. Once again, thank you for your response and your willingness to support us with a score of 7. We truly appreciate it.

Thank you for your valuable comments. We have provided the replies as follows.

W1. Thanks for your suggestion. We have added this citation to support this sentence.

W2. We have added the detailed implementation about how to select the window around the target gene in Section 2.1 task description as follows.

Specifically, in our implementation, each example contains one target gene. We first identify the transcription start site (TSS) of the target gene, then select input sequences $X_{seq}$ and $X_{sig}$ consist of $L = 200,000$ base pairs, centered on the TSS. Then, the entire sequences provide sufficient contextual information for accurate prediction of the target gene expression value $Y$.

W3. Thanks for pointing it out. We have added a complete and precise description about these mRNA Half-life features in the Appendix A.3 highlighted as red to make the list complete.

W4. Thanks for pointing it out. We removed the $R^2$ metric and retained the Pearson correlation in our experiment in our revised paper shown in Table 1 and Table 2.

W5. Description about the size of the field for each method is summarized as follows. The size of the field of view of Enformer, HyenaDNA, Mamba, Caduceus are the full sequence. And the size of the field of view for EPInformer is the extracted potential enhancer candidates based on DNase-seq measurement following ABC model [5]. Our proposed Seq2Exp also has the view of full sequence, since the generator will take the full sequence to extract the relevant sub-sequences for the prediction.

Q1. Thanks for posting this question. Currently, our model uses single base pair masks instead of continuous lumpy masks. Exploring lumpy masks could be an interesting direction for future work, as it can provide a valuable prior and enhance the discovery of regulatory elements.

Q2. Thanks for posting this question. We would like to clarify that the proposed Seq2Exp and other baselines are all under the same train/val/test dataset to make a fair comparison. Specifically, chromosomes 3 and 21 are used as the validation set, and chromosomes 22 and X are reserved for the test set. And all other chromosomes are used in training. You can find this description in Section 5.1.4.

Q3. We would like to clarify that we train a predictor model based on the mask obtained from MACS3 and compare its performance to that of a predictor model trained on regions discovered by our framework. The motivation of this experiment is to determine which method—MACS3 or our framework—can discover regions more closely related to target gene expression. Statistical peak-calling methods, such as MACS3, can be considered a straightforward approach for measuring these elements. This methodology has been utilized by models like the ABC model [5] and EPInformer [4]. Our final prediction results demonstrate that the regions found by our framework enable the predictor model to achieve better performance shown in Table 2, indicating that our framework discovers better regions than those identified by statistical methods.

[1] Effective gene expression prediction from sequence by integrating long-range interactions

[2] Mamba: Linear-Time Sequence Modeling with Selective State Spaces

[3] HyenaDNA: Long-Range Genomic Sequence Modeling at Single Nucleotide Resolution

[4] EPInformer: a scalable deep learning framework for gene expression prediction by integrating promoter-enhancer sequences with multimodal epigenomic data

[5] Activity-by-contact model of enhancer–promoter regulation from thousands of crispr perturbations

Thanks for your suggestions. We have included the information about the response to W5 in Section 5.1.2 of the new revised paper. Once again, we sincerely appreciate your support.

Thank you for your valuable comments. We have revised our paper to make it more clear for the readers.

W1. We have changed the description from ‘target Y’ or ‘$\mathbf{Y \in R}$’ to target gene expression $Y$ or add some additional sentences to further explain these symbols. We highlight these modifications in red in the revised paper. Moreover, we add the description about the model input and target gene expression in Section 2.1 task description to provide these details at the beginning of our paper to make it more clear.

W2. Thanks for pointing it out. We have reorganized the content to section 5.1.4, providing additional details as follows: Both generator $p_\theta$ and predictor $q_\phi$ are based on Caduceus architecture [1], an advanced long-sequence model considering the bi-direction and RC-equivariance for DNA. Specifically, we train for 50,000 steps on a 4-layer Caduceus architecture from scratch with a hidden dimension of 128, and more hyperparameters can be found in the Appendix A.4.

Q1. Thanks for pointing it out. We have added these hyperparameters in Appendix A.4.

[1] Caduceus: Bi-Directional Equivariant Long-Range DNA Sequence Modeling

Thanks for your comments, and we provide more information about the weaknesses and questions.

W1. Thank you for sharing your concerns regarding the contributions of our proposed method compared to EPInformer. EPInformer builds upon the annotated enhancer positions identified by the Activity-By-Contact (ABC) model [1] and then leverages deep learning to predict gene expression. Its primary contribution lies in effectively utilizing these pre-annotated enhancer regions. In the genomics field, signal peaks derived from measurements like DNase-seq or H3K27ac often serve as indicators of enhancer positions. Essentially, EPInformer’s methodology relies heavily on these signal measurements, denoted as Rm​ in Figure 1 of our paper.
However, EPInformer does not account for Rg​, the possibility that measured enhancer peaks may not interact with the gene of interest. This limitation motivates our approach, which integrates DNA sequence information and utilizes deep learning rather than solely relying on statistical measurements to discover regulatory elements.
One of the core contributions of Seq2Exp lies in its ability to discover regulatory elements (including enhancers) using the principle of information bottleneck. To the best of our knowledge, we are the first to employ deep learning for this purpose, whereas prior studies focus on the results of signal measurements. As demonstrated by our updated experimental results, Seq2Exp outperforms EPInformer in performance.
Furthermore, we implement a new version: Seq2Exp-soft, which learns the soft importance values instead of the hard mask within our framework, indicating an improved performance.

Moreover, we include an additional cell type, H1, and apply both our method and the baseline EPInformer for evaluation below. Additionally, we ran experiments with five different random seeds for K562 and GM12878 as part of the rebuttal. Due to time constraints, we provide the results of a single run on H1 using the soft version in the rebuttal. Results for five random seeds, the hard version, and all other baselines on H1 will be reported when the training is finished. You can find more details about experiments in our responses to Q1, Q2, and Q3. We wish to convey that the performance improvement observed across most experimental settings cannot be ignored.

Q1 & Q3. Thanks for posting your concerns, we wish to clarify that we use the codes of EPInformer to re-run the experiments, with the totally same setting as the original paper. In this way, we get the MSE, MAE, Pearson coefficient. To make the results more trustable, we additionally run the experiments with five different random seeds on K562 and GM12878 for the rebuttal and provide the results as below. Note that we are running multiple runs results on other baselines and hard versions, and we will report these in our paper when they are finished. In addition, we also provide our implementation in this link for reproduction: https://anonymous.4open.science/r/Seq2Exp-1E7E/README.md

Q2. We understand your concerns. We wish to convey that in most metrics and cell types, our proposed methods have surpassed the baselines methods greater than the std indicating a valid improvement. Specifically, on GM12878, the error reduction of Sep2Exp-hard on MSE and MAE is a little larger than or equal to the std of EPInformer experiments, indicating a valid improvement. Furthermore, to address your concern, we improve the performances of our framework with a soft mask version. Specifically, on GM12878, Seq2Exp-soft achieves an error reduction of more than 0.01 in both MAE and MSE compared to EPInformer. Note that 0.01 is about 5% and 3% improvement on MAE and MSE for GM12878, respectively. Furthermore, when the top 10% is selected from a soft mask for retraining, MAE can also achieve an error reduction of about 0.01. Moreover, in the experiment of K562, the improvement of both soft version and hard version is clear. We also provide the updated results as below and in our revised paper.

[1] Activity-by-contact model of enhancer–promoter regulation from thousands of crispr perturbations

H1

K562

GM12878

Dear reviewer hbLE,

As the discussion period is approaching its deadline. We would appreciate knowing if our responses have addressed your concerns. If you have any remaining questions or concerns, please don't hesitate to share them with us, we will be happy to respond! If we have sufficiently alleviated your concerns, we hope you'll consider raising your score. Thank you for your time and consideration!

This is definitely good. Thanks for replying to my questions. I've reviewed the updated draft and author responses, and as a result have updated my rating to this paper. As a reader, I now understand "what new" approach they bring to the table for prediction gene regulation. Please open source your code and likely provide a library for ease of use.