LOL-EVE: Predicting Promoter Variant Effects from Evolutionary Sequences

We thank Reviewer 2Zmu for their thorough and constructive feedback. We provide detailed responses to the various questions and points they raised.

R_2Zmu.1: Enformer

Enformer While sequence-to-activity models like Enformer represent a valuable approach, their applicability differs fundamentally from our evolutionary model in several key aspects:

• Our model aims to capture evolutionary constraint and sequence variation patterns across mammals without requiring functional genomics data
• Recent studies by Sasse et al. [2] and Huang et al. [3] demonstrate that models like Enformer have limited performance in explaining expression variation between individuals due to cis-regulatory variants
• We focus specifically on indel effects across species, while Enformer was primarily designed and validated for SNVs in human sequences

Nonetheless, to alleviate any outstanding concern of the reviewer on that point, we are working on including Enformer on the eQTL benchmark during this rebuttal. Given the aforementioned limitations of Enformer to score indels, we will be using the reported POS of the variant as the variant effect prediction.

R_2Zmu.2: SNP extension

While indels and larger variants in noncoding regions, particularly promoters, represent a critical unmet need in variant effect prediction, we acknowledge the reviewer's request and will include SNP analyses in the supplement. This work is underway and should be completed by the end of discussion period. However, we would like to reaffirm that the paper's core contributions lie in addressing the understudied challenge of indel effects, which remains our primary focus.

R_2Zmu.3: GPN-MSA GPN-MSA, while powerful for SNVs, is not well-suited for scoring indels due to fundamental limitations. Its architecture processes whole-genome alignments with independently-treated positions, making it difficult to handle sequence shifts caused by indels. Additionally, its training objective of predicting individual nucleotides from aligned context does not naturally extend to evaluating sequence length changes.

R_2Zmu.4: Updates to Causal eQTL analysis

Thank you for the feedback! We ran additional analyses for the PIP values you suggested in Table A3 and Main Table 1. While the overall conclusions for this analysis remain unchanged, we will include these ablations in appendix to confirm the robustness of our conclusions to these parameters

R_2Zmu.5: Clarifications regarding TBFS and eQTL benchmarks

We address both points in our manuscript:

TFBS Disruption Task

• We updated the “Task” description to be more clear in the TFBS Benchmark section of the paper.

Putative Causal eQTLs

• We agree the term "putative" would add clarity and will update accordingly in the abstract and in the description of the

References:

[1] A. Sasse, B. Ng, A. E. Spiro, S. Tasaki, D. A. Bennett, C. Gaiteri, P. L. De Jager, M. Chikina, S. Mostafavi, Benchmarking of deep neural networks for predicting personal gene expression from DNA sequence highlights shortcomings. Nat. Genet. 55, 2060–2064 (2023).

[2] C. Huang, R. W. Shuai, P. Baokar, R. Chung, R. Rastogi, P. Kathail, N. M. Ioannidis, Personal transcriptome variation is poorly explained by current genomic deep learning models. Nat. Genet. 55, 2056–2059 (2023).

Thanks to the authors for their detailed response and updates to the manuscript.

While I find the “Task” description of the TFBS Benchmark section to be greatly improved in clarity, and appreciate the modification to refer to “causal eQTLs” as “putatively causal eQTLs,” I have some concerns about the updated analyses presented:

1. Updates to Causal eQTL analysis. The authors wrote in their response that “We ran additional analyses for the PIP values you suggested in Table A3 and Main Table 1.” I don’t see any analysis with different PIP cutoffs in those tables, but do see Figure A9 that has an ablation of PIP thresholds, perhaps Figure A9 is what the authors were referring to? Unless you changed the PIP cutoff for the negative set, but did not reflect this change in the methods section? In any case, the analysis presented in Figure A9 is different from what I was suggesting, which was to create a more strict negative set that is more likely to be enriched for non-causal eQTLs by using different PIP thresholds for the positive and negative sets.

2. Enformer evaluation. I appreciate the addition of Enformer to the set of models benchmarked against, but find the entropy-based scoring metrics to be non-standard. While I can appreciate that the authors may have chosen these entropy-based metrics in order to compute scores from Enformer in a similar manner to the way scores are being computed from other models, I think that Enformer scores should also be reported using similar metrics to the ones used by the Enformer authors in their paper (Avsec et al., 2021). In addition, Enformer is not currently included in the SNV expansion experiment in section A4.5, but I think it would be valuable to include Enformer as a baseline in this comparison as well. Showing that in your SNV evaluation, Enformer achieves similar performance as what is reported in their original paper (AUC ~0.7), would alleviate concerns that the reported metrics for Enformer are due to non-standard scoring methods.

   Here is the description of the scoring method from Avsec et al., 2021: “We predicted the effect of a genetic variant on various annotations by computing a forward pass through the model using the reference and alternative alleles, subtracting their difference, and summing outputs across the sequence to obtain a signed score for each training dataset. We averaged scores computed using the forward and reverse complement sequence and small sequence shifts to the left and right. […] To determine how informative different variant annotations are, we trained separate random forest classifiers for each tissue to distinguish causal from non-causal variants using eight-fold cross-validation.”

3. CADD baseline. As the authors have added CADD as a baseline in Figure A4 based on other reviewer’s suggestions, and it seems to be one of the better performing methods, it seems a bit misleading to omit it from Figure 2 and not include CADD as a baseline in the other benchmarks as well.

I have updated my score to reflect concerns about the updated Enformer analysis.

We thank Reviewer 8yAh for their thorough and constructive feedback. We provide detailed responses to the various questions and points they raised.

R_8yAh.1: Comparison with sequence-to-activity models and additional baselines

As requested CADD has been added to the background section of the paper under “Sequence-to-activity models & Meta Predictors”

CADD and prediction circularity We included CADD in our frequency-based indel prioritization task but observed limitations consistent with Grimm et al.'s [1] findings about meta predictor circularity:

• Our frequency-based results (Figure A4) show CADD performs poorly compared to evolutionary approaches
• As a meta predictor integrating many features(https://cadd.gs.washington.edu/static/ReleaseNotes_CADD_v1.7.pdf), CADD likely suffers from both type 1 and type 2 circularity as described by Grimm et al.
• Moreover, we tried to get CADD scores for our TFBS evaluation but only received 6% of the data back from their server. This suggests a potential limit in coverage for the model.

Enformer While sequence-to-activity models like Enformer represent a valuable approach, their applicability differs fundamentally from our evolutionary model in several key aspects:

• Our model aims to capture evolutionary constraint and sequence variation patterns across mammals without requiring functional genomics data
• Recent studies by Sasse et al. [2] and Huang et al. [3] demonstrate that models like Enformer have limited performance in explaining expression variation between individuals due to cis-regulatory variants
• We focus specifically on indel effects across species, while Enformer was primarily designed and validated for SNVs in human sequences

Nonetheless, to alleviate any outstanding concern of the reviewer on that point, we are working on including Enformer on the eQTL benchmark during this rebuttal. Given the aforementioned limitations of Enformer to score indels, we will be using the reported POS of the variant as the variant effect prediction.

R_8yAh.2: LOL-EVE Ablation study

• We are currently working on these ablations and will provide them by the end of the week.

R_8yAh.3: TFBS study

1. Thank you for this thoughtful suggestion. While using JASPAR's pre-computed tracks seems appealing at first glance, there are significant technical challenges in reliably distinguishing true TFBS from background noise in promoter regions, particularly within 1kb of the TSS. Position Weight Matrix (PWM) scanning typically identifies numerous potential binding sites, many of which may be biologically inactive. Recent work by Zabet & Adryan [4] and Deplancke et al. [5] has confirmed that PWM-based predictions alone have high false positive rates, especially in promoter regions which are enriched for sequence motifs.

2. Clarification of Current Metrics:

• Random baseline: Since there are 2 groups (consistent/variable expression), random accuracy would be 0.5
• "A group counts as correctly predicted when the score is lower": This is unclear phrasing. It should specify that for each TF, we expect variants in consistently expressed genes to receive lower (more deleterious) scores than variants in variably expressed genes
• The prediction is evaluated per TF, not per group. A correct prediction occurs when variants in consistently expressed genes receive lower scores than those in variably expressed genes for that TF.

R_8yAh.4: Singleton vs Common Variant Experiment This experiment has been added to the text and can be seen in Figure A4. We observe that LOL-EVE still performs the best given these cutoffs.

R_8yAh.4: Updates to Causal eQTL analysis The PIP cutoff has been updated to PIP > 0.9 vs PIP < 0.01 as well as including the non normalized AUPRC. The results for this are shown in Table A3 and Main Table 1. We observe LOL-EVE results outperform other models.

R_8yAh.5: SNP extension While indels and larger variants in noncoding regions, particularly promoters, represent a critical unmet need in variant effect prediction, we acknowledge the reviewer's request and will include SNP analyses in the supplement. This work is underway and should be completed by the end of discussion period. However, we would like to reaffirm that the paper's core contributions lie in addressing the understudied challenge of indel effects, which remains our primary focus.

We thank the reviewer for their thoughtful feedback. We address all points of feedback below, including significantly enhanced figures and text edits to address clarification requests. Please let us know if we can help clarify any other aspect of our work during the discussion period.

R_Mj1m.1: Updates to Model Architecture Figure

Workflow & Pre-training/Zero-shot Prediction: Thank you for the great suggestions. We fully revamped our Figure 1 accordingly. It now more clearly illustrates our three-stage approach as follows:

1. Data preprocessing:

• We extract 1kb promoter sequences upstream of first exons across mammals, grouping species into clades and adding control codes.

2. Pre-training:

• LOL-EVE performs next-token prediction conditioned on:
  • Previous sequence context (x<i)
  • Clade (c) , species (s) , and ESM gene embeddings (g) control codes

3. Zero-shot inference

• For each benchmark task, we:
  • Generate reference and variant sequences
  • Calculate probability scores using equation (3)
  • Compare scores to identify variants with stronger effects

R_Mj1m.2: Variant Effect Score Clarification

The score in equation (3) represents the log-likelihood ratio between variant and wildtype sequences. This captures how likely/unlikely the variant sequence is compared to wildtype, given evolutionary patterns learned during pre-training. For rare variants (4.1), lower scores indicate sequences that deviate more from evolutionarily favored patterns, which aligns with rare variants being under stronger negative selection. We edited the text directly after the introduction of equation 3 to clarify this point as suggested.

R_Mj1m.3: Fair Model Comparison

• We ensure fair comparison by:
  • Using consistent scoring methods across models (detailed in A.3)
  • Applying identical preprocessing and evaluation pipelines
  • Focusing on zero-shot performance without task-specific training
• We edited the text in the first paragraph of Section 4 to clarify this point as suggested.

We thank Reviewer SwRw for their thorough and constructive feedback. We provide detailed responses to the various questions and points they raised.

R_SwRw.1: Gene-level confounding in frequency-based indel prioritization in Frequency-based indel prioritization task and Causal eQTL prioritization.

We appreciate this comment that has the intent to improve our analysis. We agree with the reviewer about the potential risk of models exploiting gene-level selective constraints rather than variant-specific properties. The main challenge with computing metrics at the gene-level is data sparsity - there is not enough data for the majority of genes (no gene has even 100 labels per class), which would result in unreliable estimates, as shown in Figures A5 and A7. However, to address these concerns, we analyzed LOL-EVE score distributions across genes with sufficient variants (≥12) and found substantial variant-level variation within genes, with genes on average capturing 62.6% of the total score range (Figure A6), suggesting our model learns more than just gene-level features.

R_SwRw.2: Comparison with sequence-to-activity models and additional baselines

1. CADD and prediction circularity. We included CADD in our frequency-based indel prioritization task but observed limitations consistent with Grimm et al.'s [1] findings about meta predictor circularity:

• Our frequency-based results (Figure A4) show that CADD performs poorly compared to evolutionary approaches
• As a meta predictor integrating many features, CADD likely suffers from both type 1 and type 2 circularity as described by Grimm et al.
• Moreover, we tried to get CADD scores for our TFBS evaluation but only received 6% of the data back from their server. This suggests a potential limit in coverage for the model.

2. FATHMM-indel We included FATHMM-indel in the causal eQTL and TFBS disruption benchmarks, which required porting these analyses to GRCh37 coordinates [Figures A10 and Table A3]. In both cases LOL-EVE significantly outperforms this baseline. Despite coordinating with the tool authors to process variants in small batches, server limitations prevented us from including FATHMM-indel in the frequency benchmark due to the sheer size of the dataset.

3. Simple baselines As suggested by the reviewer, we also included GC content differences between the variant and wild type sequence and distance to TSS as baseline features. Please see Figure A4, Table A3, and A10. Indel length was considered in analysis for Figure A4 (variant prioritization task) and A10 (TBFS task), and included as a cutoff for Table A3 (eQTL task). Using indel cutoffs for smaller indels in the causal eQTL task was problematic due to severe class imbalance (a large majority of background variants are present in the smaller indel length bins). This is illustrated in Figure A8. Consequently, we report aggregated performance for indels with length above 5, and observe similar conclusions as reported in the original manuscript.

4. Enformer While sequence-to-activity models like Enformer represent a valuable approach, their applicability differs fundamentally from our evolutionary model in several key aspects:

• Our model aims to capture evolutionary constraint and sequence variation patterns across mammals without requiring functional genomics data

• Recent studies by Sasse et al. [2] and Huang et al. [3] demonstrate that models like Enformer have limited performance in explaining expression variation between individuals due to cis-regulatory variants

• We focus specifically on indel effects across species, while Enformer was primarily designed and validated for SNVs in human sequences

  Nonetheless, to alleviate any outstanding concern of the reviewer on that point, we are working on including Enformer on the eQTL benchmark during this rebuttal. Given the aforementioned limitations of Enformer to score indels, we will be using the reported POS of the variant as the variant effect prediction.

R_SwRw.3: SNP extensionWhile indels and larger variants in noncoding regions, particularly promoters, represent a critical unmet need in variant effect prediction, we acknowledge the reviewer's request and will include SNP analyses in the supplement. This work is underway and should be completed by the end of discussion period. However, we would like to reaffirm that the paper's core contributions lie in addressing the understudied challenge of indel effects, which remains our primary focus.

Thank you for the clarifications and comparisons to additional methods. As other reviewers have noted, there seems to be a disconnect between the new figures added to the appendix and the original figures/framing in the main paper. I assume that this is because of the lack of time and that this will be fixed before the discussion period ends. Assuming those discrepancies are addressed, I still have three primary concerns, ordered by importance.

• Comparison to CADD [high importance]: The 11/19 version of the paper showed that CADD significantly outperformed LOL-EVE on the causal eQTL classification task. This was unexplainably removed in the 11/22 version. I am not sure if this is because of the circularity issue mentioned by the authors. However, after re-reading Grimm et al., it is still unclear to me how CADD suffers from either the type 1 or type 2 circularity issue on these specific tasks (which have nothing to do with missense variant effect prediction). Correct me if I'm wrong but the regulatory predictions that CADD uses as features do not seem to come from models supervised on fine-mapped eQTL data.

• Comparison to Enformer [high importance]: As other reviewers have noted, the method to evaluate Enformer on indels is non-standard and not clearly explained. A very simple metric, if you do not want to use the random forest method used by Enformer and proposed by Reviewer 2Zmu, is to just average the |ISM| scores for all accessibility or CAGE output tracks.

• Gene-level confounding [low importance]: I appreciate the authors' willingness to try to determine if gene-specific metrics could be computed to reduce the impact of confounding by gene-level selective constraints. I disagree, however, with the authors' conclusion that it is infeasible to get rid of this confounding. Based on the number of indels per gene presented in Fig. A5, I do think you could compute frequency-based indel prioritization metrics on the always >50 genes with >= 10 variants. While the metric per gene might be noisy, you could average over genes to get a pretty stable metric. For the causal eQTL prioritization task, it seems like you could subset to putatively causal eQTLs with a negative (background) eQTL in the same gene. Then, having n_{g} positive and n_{g} negative eQTLs for each gene g, you could compute a genome-wide AUROC. Regardless, I don't think this is necessarily of great importance given Fig. A6, especially given the limited time remaining in the discussion.

Overall, I believe that the authors need to include CADD performance metrics and compute Enformer scores in a more standard way before I can recommend acceptance. I encourage the authors--regardless of the results they find--to modify their main text and figures to reflect the performance of the new tools benchmarked and to not obfuscate results that make their model look worse. For now, I will keep my score as is.