We evaluated a diverse set of state-of-the-art deep learning models for their ability to predict the effects of genetic variants on enhancer activity in the human genome. These models differ markedly in architecture, parameter count, and the types of data used during training or pre-training. Notably, they were originally developed with distinct objectives in mind. For instance, models like Nucleotide Transformer [22] and DNABERT-2 [21] were designed for broad genomic representation learning, whereas others—such as Geneformer [26] and TREDNet [18]—were tailored for more specific applications.

Geneformer, in particular, is a Transformer-based model initially built to analyze single-cell gene expression profiles. However, it can be fine-tuned for a range of tasks, including chromatin state prediction and in silico perturbation modeling [26]. In our study, we adapted and fine-tuned Geneformer and other models to predict the effects of regulatory variants, leveraging their capacity to capture biological context for accurate assessment of enhancer activity changes (see Section 4, Model Adaptation and Fine-tuning). Importantly, this work compares not only model architectures but also their training foundations—highlighting how differences in pre-training data and objectives can impact downstream performance. To better disentangle these factors, we additionally assessed how the same model performs when applied across different experimental datasets. Despite their original design differences, all selected models were repurposed to predict enhancer variant effects (see Section 4).

Our benchmark spans 54,859 single-nucleotide polymorphisms located in enhancer regions, compiled from nine datasets across four human cell lines: K562, HepG2, NPC, and HeLa. This evaluation provides insights into how both model design and training data influence predictive accuracy for regulatory variant interpretation.

We performed regression analyses between model predictions and experimental log2-fold-changes for variant effects in enhancers across the human genome sequence. For each dataset, we conducted a separate analysis and averaged the results to provide a robust comparison across experimental settings (Pearson Correlation Figure 2, Spearman Correlation Figure S1, Section 4). We acknowledge that our classification groups together all models that are not purely CNN- or Transformer-based, despite significant architectural differences, particularly in the case of HyenaDNA, which diverges substantially from attention-based models. Nonetheless, our primary objective is to establish a broad guideline for enhancer variant prediction models, rather than to define rigid model categories.

CNN-based models, such as TREDNet, ChromBPNet, and SEI, demonstrate the best alignment with experimentally recorded variant effects, achieving Pearson correlation coefficients of 0.297, 0.289 and 0.276, respectively, while the most accurate Transformer-based model, Nucleotide Transformer v2, scores only at 0.105, well below the level of statistical uncertainty. Note that the reported values represent averages across nine Pearson correlation coefficients (one per dataset). Correlation strength varies by dataset, with some reaching values above 0.8 (Tables S1 and S2).

These findings highlight the effectiveness of CNNs in recognizing spatial and structural genomic patterns essential for predicting enhancer variant effects. Convolutional layers identify genomic features, such as local motifs and epigenetic markers, that influence genetic variant effects [29]. This capability is crucial for capturing the complex genomic relationships underlying gene regulation and enhancer activity. In contrast, Transformer-based models fine-tuned for specific cell lines show lower correlation values (0.019–0.105). Despite their ability to model long-range relationships and complex sequences, they struggle to detect the subtle genomic patterns linked to regulatory variant effects.

This challenge likely stems from the complexity of genomic data, lower granularity in feature detection compared to CNNs, and the extensive data required to fully exploit Transformer capabilities. To bridge the gap between local pattern recognition and long-range dependency modeling, hybrid models combine components such as CNNs, Transformers, and LSTMs. Models like the HyenaDNA series and Borzoi show moderate correlations, with HyenaDNA medium-160k-seqlen-hf achieving the highest among hybrids at 0.124. The integration of diverse components enhances performance compared to Transformer-based models while optimizing computational costs relative to CNN-based models. However, despite their architectural versatility, these hybrid models do not surpass the performance of CNN-based models, highlighting a trade-off between flexibility and focused efficiency (Figure 2 and Figure S1).

The classification task provides a deeper understanding of model performance by evaluating true positive and true negative rates in predicting variants that upregulate or downregulate gene expression (Figure 3). Among 54,859 experimental SNPs, approximately half are linked to upregulation and half to downregulation.

CNN-based models, like TREDNet, ChromBPNet, and SEI, excel in predicting variant effects. TREDNet achieves balanced performance with correctly predicted upregulation (0.60) and downregulation (0.60) rates, while SEI leads in upregulation detection (0.62) but struggles with downregulation (0.56). ChromBPNet provides an upregulation of 0.59 and a downregulation of 0.58. Transformer-based models show mixed results. Nucleotide Transformer v2-500m-multi-species matches SEI’s upregulation prediction rate (0.62) but underperforms in downregulation (0.52). DNABERT-2 performs steadily (0.54, 0.55), while Geneformer shows the weakest outcomes. Borzoi achieves a strong upregulation prediction rate (0.61) but falters in downregulation (0.47). Enformer shows moderate, balanced results (0.58, 0.55), while HyenaDNA remains consistent but unremarkable.

A crucial aspect of each model is its ability to identify causal SNPs within their respective Linkage Disequilibrium (LD) blocks. This is arguably the most important function of variant classification models, as it enables resolution of causal variants from the large sets of associated SNPs in LD blocks linked to GWAS tag SNPs. To investigate this, we curated a dataset of 14,183 causal SNPs specific to the HepG2 cell line (Figure 4, Dataset 2, see Section 4).

The choice of HepG2 was motivated by three factors: (i) the high quality and comprehensiveness of available enhancer and variant annotations in this cell line; (ii) its extensive use in prior studies, which provides a solid basis for benchmarking and comparison; and (iii) the availability of experimentally supported causal SNPs, which allows us to extract and validate true causative variants for model evaluation. For each causal SNP, we identified associated LD-block SNPs with an r²-value above 0.8 within a 500 kbp window, resulting in 263,286 SNPs. This design enables a direct and rigorous assessment of each method’s ability to recover the causal SNP from within an average ~26:1 pool of associated variants.

For each identified SNP, we generated two sequences of 1 kbp: one containing the altered variant and the other containing the reference nucleotide. In both cases, the SNP or reference nucleotides were positioned at the center of the sequence. These sequences were then input into several models to assess their ability to predict causal SNPs by calculating the log2 fold-change between the scores of the alternative sequence and the reference sequence (Section 4). This procedure was repeated for all LD-block associated SNPs and the log2 fold-change score of the known causal SNP was compared to its associated variants. Finally, we computed the percentage of causal SNPs within an LD-block correctly predicted as causal based on their score being either the highest (“top-1” test) or within the top two and three highest scores (“top-2” and “top-3” tests).

The Borzoi model achieved the highest accuracy, correctly identifying 42.5% of causal SNPs, followed by the SEI, Enformer, ChromBPNet and TREDNet models, with 40.7%, 39.4%, 35.2%, and 33.8% in the top1 test, respectively, (Figure 4). The causal variant detection rate increased substantially in the top-2 and top-3 test with Borzoi (44.9%, 48.9), SEI (43.6%, 47.0%), Enformer (42.4%, 46.2%), ChromBPNet (43.4%, 45.1%), and TREDNet (36.1%, 40.1%). Among the Transformer-based models, the Gena LM bert-based-t2t-multi demonstrated the best performance (top-1 30.4%, top-2 32.7%, and top-3 37.5%), marginally surpassing other Nucleotide Transformer and Gena LM variants, which achieved accuracies between 28.0% and 30.0% in top-1 test. The HyenaDNA models, categorized as hybrid architectures, displayed moderate performance, with accuracies ranging from 24.3% to 26.5% in top-1 test, and larger sequence lengths, such as those used in the HyenaDNA large model, led to improved results. The Geneformer model, however, had the lowest accuracy, identifying 24.3% of causal SNPs in top-1 test. Overall, hybrid and CNN architectures outperformed Transformer-based models in this application, underscoring their suitability for tasks involving causal SNP detection.

Overall, models like TREDNet and SEI consistently outperform others in both regression analysis and classification tasks, highlighting the strength of CNN-based architectures in capturing local genomic patterns critical for regulatory variant effect prediction. Hybrid models, such as Borzoi and Enformer, provide balanced performance by combining local pattern recognition with long-range dependency modeling. While they do not surpass CNNs in overall predictive accuracy, they are particularly well-suited for tasks requiring the detection of causal variants within LD-blocks, where modeling long-range context is essential. Transformer-based models, though powerful for learning complex dependencies across extended sequences, currently lack the resolution of CNN and hybrid models for tasks centered on fine-scale motif recognition. However, when optimized, they are best applied to tasks that demand integration of distal regulatory interactions, such as enhancer–promoter contact or tissue-specific regulatory architecture. Based on these observations, we suggest a two-step strategy: first, apply CNN-based models to efficiently identify high-risk variants across the genome, and then leverage hybrid models such as Borzoi to resolve the causal SNPs within the corresponding LD blocks, where long-range context plays a decisive role.

The results of experimental assays of regulatory variants depend on the selected degree of certainty in separating significant and insignificant variant effects. To investigate the impact of experimental data significance on correlation with modeling results, we analyzed the top models for each architecture type: TREDNet (CNN), HyenaDNA (hybrid), and Nucleotide Transformer (Transformer). The experimental data was binned using the following significance thresholds, p < 0.5, p < 0.1, p < 5 × 10⁻⁵, p < 10⁻⁵.

The results demonstrate that models align better with experimental findings when applied to variants with higher statistical significance. For instance, data points with greater experimental confidence show increased clustering in the first and third quadrants, representing cases where both experimental and predicted results are either positive (upregulation) or negative (downregulation) (Figure 5A–C and Figures S2 and S3a–d).

Among the evaluated architectures, HyenaDNA exhibited the greatest improvement in predictive accuracy with the increase in the statistical significance of experimental results (from all data 0.124 to 0.226 with p < 10⁻⁵, Figure 5C). However, TREDNet consistently displayed superior overall performance, with a higher density of points aligning with experimental outcomes (from all data 0.298 to 0.421 with p < 10⁻⁵, Figure 5B). This trend is further illustrated in the density plots, where TREDNet exhibits increased clustering in the first and third quadrants, demonstrating stronger correlations between model predictions and experimental values compared to both hybrid and Transformer models (Figure 5D). Moreover, the correlation coefficients across varying significance thresholds highlight TREDNet’s robustness, achieving the highest values under all conditions. True positives and true negatives increase with the significance of the experimental results (positive/positive 0.63, negative/negative 0.66, with p < 0.05; positive/positive 0.72, negative/negative 0.78, with p < 10⁻⁵), with the hybrid model and the Transformer model following in performance (Figure 5E). The rest of models’ performance relative to experimental data significance can be found in Figures S2 and S3a–d. The correlation between predictions and experimental data improves notably as values around zero are removed, confirming that high-uncertainty regions reduce the interpretability and strength of correlation (Figure S4).

Given the strong cell-line specificity of enhancers [31], determining whether model performance depends on the type of cell line used is essential for understanding their predictive capabilities and limitations. For this, we evaluated whether a model fine-tuned for a specific cell line demonstrates superior performance compared to models from other architectures, such as CNNs, which often rely on generalization rather than cell-line-specific adaptation. Indeed, Transformer-based models, pre-trained on extensive datasets and fine-tuned for targeted tasks, can offer both flexibility and accuracy, making them well-suited for capturing the unique regulatory landscape of cell lines. In contrast, CNN-based models, such as TREDNet and SEI, though highly effective for specific tasks, could require additional steps, such as retraining or the addition of specialized layers, to adapt their outputs for cell-line-specific predictions. However, the results indicate that TREDNet and SEI remain the most effective models for predicting variant effects across different cell lines, even after model fine-tuning (Pearson Correlation Table 1, Spearman Correlation Table S1, Dataset-wise Sup Table S2 and Sup Table S3). For example, the Nucleotide Transformer v2-250m-multi-species model aligns with experimental results in K562 and HepG2 (0.1998 and 0.11) but performs less effectively in NPC and HeLa cells. This variability could stem from differences in pre-training datasets—such as multi-species versus human-specific data—and the fine-tuning approach.

Notably, the hybrid model Enformer, which integrates features from Transformers and CNNs, performs particularly well in HeLa cells (0.245), ranking second only to SEI (0.279). As observed in the previous test, CNN-based models demonstrate superior performance, showcasing robustness across cell-line variants. We also assess the impact of the certainty of experimental results, SEI, which does not require training or fine-tuning for specific cell lines, outperforms fine-tuned models. TREDNet, likely benefiting from its second-phase training tailored to the specific cell line, achieves the best overall performance, further highlighting its ability to leverage cell-line-specific information effectively.

Models designed to predict variant effects rely primarily on the change in confidence when detecting the enhancer region after introducing the variant.

Next, we investigated whether improving the accuracy of enhancer sequence detection enhances predictions of variant effects. While SEI, Borzoi, and Enformer do not require fine-tuning, TREDNet incorporates a two-phase training despite being a CNN-based model. The second phase can be considered fine-tuning [18].

Transformer-based architectures, particularly Nucleotide Transformer models pre-trained on multi-species datasets, emerge as strong performers in enhancer classification (Figure 6). Models such as Nucleotide Transformer v2-250m-multi-species and Nucleotide Transformer 2.5b-multi-species consistently demonstrate high accuracy, as indicated by their auROC and auPRC metrics (auROC: 0.967, 0.963; auPRC: 0.958, 0.951). However, the correlation between the models’ outputs and the MPRA experimental results, which measure the fold-change effect between reference and alternative sequences, remains relatively low (ranging from 0.112 to 0.15). SEI stands out as the top performer in detecting enhancer activity, achieving the highest auROC (0.974) and auPRC (0.971) scores among all models. It is also the second best for the correlation between model outputs and experimental results (0.295). TREDNet demonstrates strong performance with a higher correlation (0.318), while balancing high auROC (0.872) and auPRC (0.814) scores. The HyenaDNA models show competitive classification metrics, with the Hyenadna medium-160k-seqlen-hf and Hyenadna large-1m-seqlen-hf versions attaining notable auROC (0.933, 0.904) and auPRC (0.921, 0.893) scores. However, their Pearson correlation coefficients are moderate (0.152, 0.151), indicating that while these models classify enhancer activity effectively, further optimization is needed to improve their predictive accuracy for variant-specific effects. Borzoi performs the weakest in detecting enhancer regions (auROC: 0.672, auPRC: 0.65), while Geneformer fares poorly in detecting variants (auROC: 0.027). Yet, this weak enhancer detection performance doesn’t translate in inaccurate detection of causal variants, as Borzoi outperforms all other methods in this critical test (Figure 6).

Fine-tuned models consistently outperform their one-shot counterparts, achieving F1 scores between 0.85 and 0.92, compared to 0.70 to 0.80 for one-shot models (Figure 7, top). This improvement in performance is reflected in the correlation between model predictions and experimental results, particularly in detecting enhancer variants (Figure 7, bottom). The performance enhancement is consistent across most architectures, with CNN-based models showing remarkable results even without fine-tuning. However, SEI, which is available only in its ‘one-shot’ implementation, achieves a high F1 score and correlation, surpassing many fine-tuned Transformer models also when the models were trained on the difference cell-line than the tested dataset (Figure S6). Additionally, some Transformer-based models, such as Nucleotide Transformer 2.5 b–1000 g, show no improvement in detecting enhancer variants.

TREDNet and SEI achieve superior performance in regression (Figure 2) and classification tasks (Figure 3). TREDNet and ChromBPNeT demonstrate balanced predictive capability for both upregulated and downregulated variants, while SEI excels in detecting upregulated variants. CNNs effectively capture local genomic patterns like motifs and epigenetic markers, making them well-suited for enhancer activity analysis. Their simpler architecture yields high Pearson correlation coefficients (0.297 for TREDNet, 0.289 for ChromBPNet, and 0.276 for SEI). However, their inability to model long-range dependencies limits performance in tasks requiring broader sequence context.

Transformers, including the Nucleotide Transformer and DNABERT-2, capture long-range dependencies but struggle with fine-grained regulatory signals. They underperform in enhancer variant effect prediction, with accuracies between 28.0% and 30.0% in causal SNP identification (Figure 4). These models require substantial training data, and their performance could improve with expanded datasets or refined focus on regulatory signals.

Hybrid models, such as HyenaDNA and Borzoi, integrate CNNs, Transformers, and LSTMs to balance local and long-range dependencies. While HyenaDNA shows moderate correlation with experimental results, Borzoi outperforms other models in causal SNP detection (Figure 4), making it valuable for GWAS analysis. However, hybrid models do not consistently surpass CNNs in overall predictive accuracy and should be assessed based on task-specific requirements.

The reliability of model predictions is closely tied to the quality of experimental data. Model alignment with experimental results improves as statistical confidence in variant effects increases (Figure 5 and Figure 6). TREDNet, for instance, maintains high predictive accuracy even under varying significance thresholds, emphasizing the robustness of CNN-based models. When working with high-quality data—for example, datasets with clean sequencing reads, low background noise, strong signal-to-noise ratios, and highly significant peaks or variants (low p-values, e.g., <0.001)—CNN models are the preferred choice, as they excel at capturing local regulatory patterns with high accuracy. In contrast, when datasets consist of a mixture of quality levels—such as variable read depth, moderate noise, less stringent peak calls, or broader distributions of statistical confidence—hybrid models like HyenaDNA are recommended, as they maintain robust performance under more heterogeneous conditions. Finally, in cases where very large datasets are available (e.g., >100,000 sequences), fine-tuning Transformer-based models such as Nucleotide Transformer v2 becomes advantageous, as these architectures are able to exploit scale to learn complex dependencies and integrate long-range regulatory interactions more effectively. These findings suggest that enhancing experimental rigor and statistical refinement can reduce uncertainty in variant effect assessments, leading to more reliable predictions.

In addition to data quality, model generalization across different cell lines is a critical factor. Transformer models fine-tuned for specific cell lines can capture regulatory landscapes but often underperform compared to CNNs. For instance, the Nucleotide Transformer v2-250m-multi-species achieves reasonable results in K562 and HepG2 but struggles in NPC and HeLa, likely due to pre-training limitations. In contrast, CNN models generalize well across cell lines without requiring fine-tuning, further demonstrating their robustness in variant effect prediction. Enhancer detection models, such as the Nucleotide Transformer v2-250m-multi-species, achieve high classification accuracy (auROC: 0.967, auPRC: 0.958), yet their predictions show weak correlation with experimental results (Figure 6). CNN-based models, particularly SEI, perform well in both enhancer detection and variant effect prediction. While fine-tuning improves Transformer performance, SEI, despite its one-shot implementation, still outperforms many fine-tuned Transformers (Figure 7).

A notable advantage of Transformer architectures lies in their potential for zero-shot applications (Figure 7), where models can be applied directly without additional training or fine-tuning. In real-case scenarios, this property makes Transformers highly attractive for tasks where training data are scarce or unavailable. At present, the zero-shot models tested here do not yet match the performance of CNN-based approaches, but we consider this a temporary limitation. As Transformer models continue to grow in size and scope—incorporating more comprehensive pre-training across diverse cell types, tissues, and species—they are expected to increasingly capture the regulatory information needed for accurate variant effect prediction in zero-shot settings.

In conclusion, our findings suggest that SEI and TREDNet are the most effective models for predicting the effects of genetic variants on enhancer activity. Hybrid models like Borzoi demonstrate strong potential in extracting causal SNPs from GWAS data, showcasing their value in applications that require the integration of both local and global sequence contexts. Transformer-based models, on the other hand, are particularly effective at identifying enhancer regions due to their ability to capture long-range dependencies, although they still face challenges with fine-grained genomic features and often require substantial data for optimal performance. These results highlight the strengths of different model architectures for various aspects of regulatory variant prediction and underscore the need for further optimization to enhance their performance. Future work should aim to improve feature detection in Transformer models and refine the flexibility of hybrid architectures to better capture both local and long-range genomic dependencies.

The datasets from HepG2, K562, NPC, and HeLa cell lines underwent a standardized preprocessing pipeline to ensure consistency and feature extraction for downstream analysis. Two main pipelines were applied: one for training and testing enhancer classifier models, such as Nucleotide Transformers and DNABERT-2, and another for mutagenesis experiments. For the enhancer classification pipeline, 1kb sequences were extracted for both enhancer (positive) and control regions. To maintain consistency across models, all nucleotides were capitalized, as some models are case-sensitive. Unknown nucleotides marked as “N” were removed from the sequences. The datasets were then split into training, validation, and test sets, and tokenized using the corresponding model’s tokenizer with a maximum length of 512 tokens. Labels were generated to distinguish between positive and control sequences. In the mutagenesis experiment pipeline, the reference and alternative nucleotide locations were extracted for each row of the dataset. For each nucleotide of interest, a 1kb sequence was generated with the target nucleotide positioned in the center. Chromosome and position information, along with reference and alternative alleles, were also extracted. While longer input sequences generally lead to improved correlation with experimental measurements—indicating that additional genomic context enhances predictive power—the signal-to-noise ratio tends to decrease at higher lengths, suggesting that the additional input may also introduce irrelevant or noisy features that affect model robustness (Figure S5). Log2 activity metrics, such as fold-change values or effect sizes, were calculated to assess regulatory activity. Additionally, signed p-values were derived from raw or adjusted p-values to capture both the statistical significance and directionality of the regulatory effects.

The preprocessing step ensured that variations in dataset structure did not interfere with downstream analysis. To further enhance consistency, hg19 genomic coordinates were used for variant mapping, and enhancer activity metrics were uniformly applied across datasets. Differences in expression or activity levels between experimental conditions, such as mutant versus control or modern versus archaic sequences, were computed when applicable. Ultimately, all datasets were converted into a standardized format that included chromosome information, positions, reference and alternative alleles, log2 activity ratios, signed p-values, raw p-values, and other relevant regulatory features.

The objective of fine-tuning is to adapt models originally designed as enhancer classifiers to assess the impact of SNPs on enhancer activity. By introducing variant sequences, we evaluate how nucleotide changes affect model scores. To achieve this, we fine-tuned pre-trained DNA language models—including Geneformer (initially developed for scRNA-seq data)—as enhancer classifiers (Table 2).

Sequences were preprocessed as described earlier, with a length of 1 kbp, balancing genomic context and model compatibility. This length also reflects the average enhancer size of ~400 bp. When adapting models like Enformer (196 kb input) and Borzoi (512 kb input), which do not accept arbitrary sequence lengths, we retained their full input size. The 1 kb variant-centered sequence was inserted into representative genomic contexts of the appropriate size, and predictions were averaged over multiple such contexts. Yet, for evaluation settings like Borzoi’s sliding window analysis where the method natively supports shorter sequences, we used 1 kb inputs directly.

Most models were loaded from the Hugging Face Hub using AutoTokenizer and AutoModelForSequenceClassification. Exceptions such as TREDNet, Enformer, SEI, and Borzoi were sourced from GitHub and required custom pipelines for preprocessing and inference.

Fine-tuning was performed using the Hugging Face Transformers library for flexible training and evaluation. Models such as HyenaDNA, Geneformer, Gena-LM, DNABERT-2, and Nucleotide Transformers were trained via the Trainer API on four biosample-specific datasets: HeLa (BioS2), neural progenitor cells (NPC, BioS45), HepG2 (BioS73), and K562 (BioS74). Sequences were cleaned to remove ambiguous bases (Ns) and split into training, validation, and test sets. Model-specific tokenizers processed sequences with a maximum input length of 512 tokens. For non-fine-tuned models like SEI, ChromBPNet, and Enformer, the same 1 kbp sequence length was used with padding when appropriate.

Key hyperparameters included a batch size of 8, a learning rate of 1 × 10⁻⁵, and a maximum of 200 epochs. Early stopping with a patience of 5 epochs was used to avoid overfitting, with progress logged every 500 steps. Fine-tuned models were saved in output directories named by model checkpoint and biosample ID and uploaded to the Hugging Face Hub for reproducibility. GPU memory was cleared after each session for efficiency. This workflow enabled the adaptation of general-purpose DNA language models for cell-type-specific regulatory prediction while ensuring consistency with non-fine-tuned models.

Training and fine-tuning were performed on the NIH BioWulf high-performance computing clusters equipped with A100 GPUs, which enabled efficient processing of large-scale datasets (e.g., ~130,000 sequences of positive and control enhancer regions). Full training runs typically completed within a few hours using 8 A100 GPUs, while fine-tuning on cell line–specific subsets required proportionally less time and computational resources. CNN-based models (SEI, ChromBPNet): These models were the most computationally efficient, with training converging within 2–4 h on ~130 k sequences. Fine-tuning on smaller, cell line–specific subsets (10–30 k sequences) often required less than 2 h on a single A100 GPU, highlighting their suitability for rapid retraining and deployment. Hybrid models (Borzoi, HyenaDNA, Caduceus): Hybrid architectures combining convolutional and attention/state-space components required moderately more resources than CNNs but remained efficient. Full training typically required 4–6 h with 8 GPUs, while fine-tuning was feasible within 1–2 h using 2–4 GPUs. Transformer-based models (Nucleotide Transformer, DNABERT-2, GENA-LM, Geneformer) required substantially greater computational resources due to their parameter size and sequence length handling. Fine-tuning the larger variants (e.g., Nucleotide Transformer multi-species, GENA-LM large) could take 6–8 h across 4 GPUs. Long-range attention models, Enformer and Borzoi, designed to model >100 kb genomic contexts, were the most computationally demanding, 4 GPUs and ~8–10 h for large datasets.

During the fine-tuning process, the models were optimized to distinguish enhancer sequences from control sequences. In contrast, the validation process focused on evaluating the difference between the reference enhancer sequence and the alternative sequence, where the middle nucleotide was mutated. The difference in output scores between these two sequences, obtained from the fine-tuned model, allowed for the estimation of the variant’s effect on the enhancer. If the output scores for both sequences were identical, the variant had no effect. Conversely, a change in output scores indicated either up- or down-regulation.

The prediction process involved generating regulatory effect predictions for each dataset by calculating the log2 ratio of alternative to reference allele probabilities for each variant. This metric provided a standardized approach to assess the relative regulatory effects of alternative alleles compared to reference alleles. Since each fine-tuned model produced logits for each class (enhancer or not), predictions for both reference and alternative alleles were processed by first calculating margin logits, which represent the difference in predicted class probabilities. These logits were then transformed into probabilities using the sigmoid function. The log2 ratio of these probabilities was computed to quantify the regulatory differences between the two alleles. Additionally, precomputed predictions stored in pickle files were loaded for specific models and datasets. The predictions for reference and alternative alleles were used to compute a ratio, representing the alternative-to-reference prediction ratio. This ratio was subsequently log2-transformed to ensure consistency in the output format. Additional preprocessing steps were applied to specific datasets as needed, such as loading specialized prediction files.

The prediction pipeline was applied across all datasets by iterating through experiments defined in a metadata table. For each experiment, data were preprocessed, and predictions were either computed or loaded accordingly. The results were compiled into a unified dictionary containing log2 prediction ratios for all variants across experiments. As in the fine-tuning process, the validation of variant/reference sequences was performed using 1 kb sequences, ensuring consistency in genomic context and input length across models. Several metrics were employed to evaluate model performance, including the Area Under the Curve (AUC), Precision-Recall Curve (PRC), F1-score, and correlation coefficients such as Pearson and Spearman. These metrics provided a comprehensive assessment of model performance, encompassing both classification accuracy and the correlation between predicted and observed values.

This study utilizes several deep learning models designed for genomic sequence analysis, with a particular focus on applications to the Homo sapiens genome. The models represent the forefront of innovation in deep learning for genomics, leveraging diverse architectures and pre-training paradigms to handle the complexity of human genomic data (Table 3). To maintain focus on models explicitly designed for the human genome, this study excluded models like EVO and GPN, which are primarily aimed at other species or biological domains.

We constructed a dataset of positive and control enhancer sequences from four human cell lines: HepG2, Hela, K562, and Neural Progenitor Cells (NPCs). These cell lines were selected for their relevance in genomic research and the availability of comprehensive epigenomic data, such as DNase-Seq and histone modification profiles, sourced from the ENCODE project using genome annotations for hg19 and hg38.

Positive enhancer sequences were defined through a two-step process. First, open chromatin regions were identified using DNase-Seq data. These regions were then intersected with H3K27ac histone modification signals, a hallmark of active enhancers. To refine the dataset, we excluded non-enhancer regions such as exonic regions (coding sequences), promoter regions (near transcription start sites), and low-confidence regions listed in the ENCODE Blacklist (e.g., repetitive sequences). Specific blacklisted peaks from UCSC browser tables were also removed.

Control sequences were selected using a similar filtering process but specifically excluded regions overlapping both DNase-Seq peaks and H3K27ac signals to ensure they represented non-active genomic regions. Both positive and control sequences underwent identical exclusion criteria to eliminate confounding genomic features such as exons and promoters.

All sequences were standardized to 1 kb in length to balance sufficient coverage of regulatory regions with computational efficiency. This size is widely used in enhancer studies as it captures relevant regulatory elements while avoiding extraneous genomic features. The dataset includes an equal number of positive and control sequences for each cell line: HepG2 (14,062 each), NPC (12,466 each), K562 (19,968 each), and Hela (31,247 each). Metadata from the ENCODE project was used to retrieve biosample IDs and associated DNase-Seq and histone modification data. Reference files provided coordinates for exons, blacklist regions, and promoters, ensuring only valid enhancer or control regions were included in the final dataset.

To systematically characterize the functional effects of genetic variants on gene regulation, we compiled a diverse set of high-throughput reporter assay datasets spanning multiple cell types, genomic contexts, and experimental designs (Table 4).

Datasets 1 and 2 were sourced from the study by [35], which utilized the Survey of Regulatory Elements (SuRE) reporter assay to identify reporter assay quantitative trait loci (raQTLs)—SNPs influencing regulatory element activity. For K562 cells, 19,237 raQTLs were identified with an average allelic fold-change of 4.0-fold, while HepG2 cells yielded 14,183 raQTLs with a higher average fold-change of 7.8-fold. Most raQTLs were cell-type-specific, showing limited overlap between the two cell lines.

Dataset 3 was derived from [36] and focuses on a ~600 bp enhancer region at the SORT1 locus (1p13). This region includes SNP rs12740374, which creates a C/EBP binding site that influences SORT1 expression. The enhancer was tested in HepG2 cells using a luciferase reporter assay. Data for flipped sequences were excluded due to the orientation-independent nature of enhancers, and deletion datasets were omitted due to computational model constraints, which focus on single-nucleotide variants.

Dataset 4, from [37] used MPRA to investigate 2756 variants linked to red blood cell traits in K562 cells. The study identified 84 highly significant SNPs with functional effects, enriched in open chromatin regions shared between K562 cells and primary human erythroid progenitors (HEPs). This dataset provides insights into erythroid-specific regulatory activity.

Dataset 5, generated by [38], includes MPRA results for 284 SNPs tested in HepG2 cells. These SNPs, located within regions of high linkage disequilibrium identified through eQTL analyses, offer a detailed map of allele-specific regulatory effects relevant to liver-specific gene regulation and genotype-phenotype associations.

Dataset 6, from [39], includes lentiMPRA results for 14,042 single-nucleotide variants fixed or nearly fixed in modern humans but absent in archaic humans (e.g., Neanderthals). Tested across three cell types—embryonic stem cells, neural progenitor cells, and fetal osteoblasts—the study identified 1791 sequences with regulatory activity, including 407 with differential expression between modern and archaic humans. This dataset focuses on NPC data.

Datasets 7, 8, and 9, from [40], include MPRA results for three coding exons—SORL1 exon 17, TRAF3IP2 exon 2, and PPARG exon 6—tested in HeLa cells. Respectively, these datasets cover 1962, 1614, and 1665 SNPs. The study explored enhancer activity within these coding exons (eExons), finding that ~12% of mutations altered enhancer activity by at least 1.2-fold, with ~1% causing changes of less than 2-fold.

Code, data, and analysis repository: https://github.com/tanoManzo/AI4Genomic.git. Fine-tuned models: https://huggingface.co/tanoManzo. (Access data: 5 September 2025).