Benchmarking DNA foundation models for genomic and genetic tasks

Domain: physics
Published: 2025-11-28
Editorial outcome: minor_revision

Abstract
The rapid evolution of DNA foundation models promises to revolutionize genomics, yet comprehensive evaluations are lacking. Here, we present a comprehensive, unbiased benchmark of five models (DNABERT-2, Nucleotide Transformer V2, HyenaDNA, Caduceus-Ph, and GROVER) across diverse genomic and genetic tasks including sequence classification, gene expression prediction, variant effect quantification, and topologically associating domain (TAD) region recognition, using zero-shot embeddings. Our analysis reveals that mean token embedding consistently and significantly improves sequence classification performance, outperforming other pooling strategies. Model performance varies among tasks and datasets; while general purpose DNA foundation models showed competitive performance in pathogenic variant identification, they were less effective in predicting gene expression and identifying putative causal QTLs compared to specialized models. Our findings offer a framework for model selection, highlighting the impact of architecture, pre-training data, and embedding strategies on performance in genomic and genetic tasks.

Introduction
Led by the advances in Natural Language Processing (NLP) in recent years, foundation language models through self-supervised pre-training have been the paradigm of decoding information in sequences. By representing sequences as numerical embeddings, foundation language models can outperform previous methods in many downstream tasks such as sequence classification and sequence generation. As natural language-based foundation models like GPT-4, Llama 3, and Qwen 3 have been proven successful, similar ideas have been extended to other domains by interpreting domain-specific languages with unique semantic rules, and examples include foundation models on programming codes, protein sequences and single-cell sequencing. With the long-lasting interests in decoding DNA sequences to understand the epigenetic patterns, transcriptional regulations, and disease associations, DNA foundation language models have also emerged recently including DNABERT-2, Nucleotide Transformer, HyenaDNA, Caduceus-Ph, and GROVER. These models are pre-trained on large genomic datasets such as the human reference genome, human whole-genome sequencing datasets like 1000 Genomes project datasets, and multi-species genome datasets. After fine-tuning, they have shown promising results in DNA sequence classification tasks.

A critical aspect of DNA foundation models is the method used to generate sequence embeddings, with sentence-level summary token, mean token embeddings, and maximum pooling being three primary approaches. The comparative efficacy of these embedding methods in DNA sequence analysis remains understudied, despite their potential impact on model performance across varying sequence lengths and biological contexts.

With the rapid evolution of DNA foundation models of various architectures and the wide range of genomic analysis tasks to be solved, there is a pressing need for effectively evaluating these models. However, most of the current evaluations on DNA foundation models are biased, as they are conducted after fine-tuning, which may introduce biases in model performance comparison. For instance, different models may have various levels of overfitting depending on which layers are selected to update during fine-tuning. The use of advanced parameter-efficient fine-tuning methods further complicates this issue by introducing additional hyperparameters that could impact model fitting. A recent work leveraged this problem by directly comparing DNA foundation models based on their output embeddings, where the weights in all layers were frozen and a trainable convolutional neural network (CNN) was appended to the last layer. While this approach mitigates fine-tuning biases, the study scope was limited to several human genome analysis tasks. There is yet to be a comprehensive benchmark which addresses crucial challenges including model performance across diverse genomic tasks and species, scalability with sequence length, and ability to capture biologically significant features.

In this study, we first provide a comprehensive and unbiased evaluation of five state-of-the-art DNA foundation language models across 57 diverse datasets spanning four major categories: human genome region classification, multi-species genome region classification, human epigenetic trait classification, and multi-species epigenetic trait classification. Our benchmarking assesses models’ performance across multiple species, sequence lengths, and genomic tasks including promoter identification, enhancer classification, transcription factor binding site prediction, and epigenetic modification detection. We systematically evaluate three different pooling methods (summary token, mean token, and maximum pooling) for each model, revealing substantial and consistent performance differences. Beyond classification tasks, we extend our evaluation to assess models’ capabilities in gene expression prediction, where we evaluate how well zero-shot embeddings can predict tissue-specific expression levels from genomic sequences using GTEx v8 data. We further assess variant effect quantification, examining models’ ability to distinguish pathogenic from common variants and identify functional quantitative trait loci (QTLs). Additionally, we investigate topologically associating domain (TAD) recognition to determine whether models inherently learn higher-order chromatin structures. We also conduct novel experiments to understand factors influencing model performance, including a controlled pre-training experiment comparing multi-species versus human-only training data. Our evaluation framework provides actionable insights for model selection, the strengths and limitations of different architectural choices, embedding strategies, and pre-training approaches across diverse genomic applications. All code and datasets are available at https://github.com/ChongWuLab/dna_foundation_benchmark.

Results
Sequence classification: pooling methods benchmark
We evaluated the zero-shot embeddings from DNA foundation models on 57 sequence classification datasets, including 52 binary classification datasets and 5 multi-class datasets (Section Methods: Benchmarking Datasets). For each dataset, we first generated zero-shot embeddings for all sequences, then split the samples into training and testing sets, and finally trained a classifier and reported its performance predicting the labels of each sequence in the test set from their zero-shot embeddings (Section Methods: Benchmarking Methods). The detailed workflow of sequence classification benchmarking can be found in Fig. 1.

First, we evaluated three output pooling methods: sentence-level summary-token ([CLS] or [SEP]) embedding, mean token embedding, and maximum pooling. We sought to determine the optimal method for generating sequence-level representations, as the choice of pooling method can significantly alter how genomic information is captured and passed to downstream analyses. The results revealed significant differences in the performance of pooling methods across all DNA foundation models. Using DeLong’s test for statistical significance (p < 0.01), we systematically assessed which pooling method performed better for each dataset. Our analysis showed that mean token embedding consistently delivered superior performance against others (Supplementary Table 1). Specifically, mean token embedding delivered higher AUROC (AUC) with statistical significance than both summary-token embedding and maximum pooling, in 41 out of the 52 binary sequence classification datasets for DNABERT-2, 42 for NT-v2, 35 for HyenaDNA, 37 for Caduceus-Ph, and 41 for GROVER. On the contrary, maximum pooling or summary-token pooling rarely outperformed the other methods, being optimal in only a few specific datasets.

Figure 2 illustrates the distribution of AUC scores for all models when using different pooling methods. The average increase in AUC when switching from summary token to mean token embedding was 4.0% (interquartile range: 2.0%−5.5%) for DNABERT-2, 6.8% (interquartile range: 3.7%−9.6%) for NT-v2, 8.7% (interquartile range: 4.6%−12.9%) for HyenaDNA, 5.9% (interquartile range: 2.8%−9.1%) for Caduceus-Ph, and 1.4% (interquartile range: 0.7%−1.9%) for GROVER across all binary classification tasks. This consistent enhancement with statistical significance underscores the superiority of mean token embedding.

The improved performance with mean token embedding suggests that, by averaging the embeddings of all non-padding tokens, mean token embedding may provide a more comprehensive representation of the entire DNA sequence as opposed to relying on a single summary token. This finding is particularly relevant for some DNA sequence classification tasks, such as promoter and enhancer identification, where the discriminative features may be distributed throughout the sequence rather than concentrated in a specific region. For example, in the promoter identification task for the GM12878 cell line, mean token embedding improved the AUC from 0.964 to 0.986 for DNABERT-2, a 2.3% increase of statistical significance (p < 0.01 by DeLong’s test). More strikingly, for the B. amyloliquefaciens genome, the improvement was from 0.689 to 0.864 for HyenaDNA, representing a 25.4% increase of statistical significance (p < 0.01 by DeLong’s test). These examples illustrate how mean token embedding can capture distributed features more effectively across the entire sequence.

Additionally, we observed that the performance differences among the models were reduced when using mean token embedding, implying that this pooling method helps to mitigate the architectural variations across the models. Specifically, the range of average AUC scores across models decreased from 0.708 to 0.799 with summary-token pooling, to 0.795 to 0.822 with mean pooling. This observation further reinforces the value of mean pooling as a more robust approach that can help standardize model evaluation.

With the optimal pooling method established, we also evaluated our choice of the downstream classifier. We selected random forest as our primary model, as it requires minimal hyperparameter tuning and can inherently handle high-dimensional inputs without dimension reduction, in both ways helps avoid the introduction of additional evaluation bias. It also has the capacity to capture the complex, non-linear relationships potentially encoded in genomic sequences. To evaluate our approach, we run random forest, naïve Bayes, and elastic-net logistic regression on all binary sequence classification datasets (Supplementary Fig. 1). The results reconfirmed that mean token embedding was the superior pooling strategy across all foundation models and all classifiers. Within mean pooling configuration, both random forest and elastic-net proved to be the more competitive classifiers. Given its strong performance and well-documented generalization capabilities, we selected random forest as the standard classifier for all analyses in our sequence classification benchmark, providing a reliable estimate of each foundation model’s performance on similar tasks.

Sequence classification: human genome regions
We evaluated the performance of the five DNA foundation models, DNABERT-2, NT-v2, HyenaDNA, Caduceus-Ph, and GROVER, on a diverse set of human genome sequence region classification tasks. For these tasks, as shown in Table 1, when using mean token pooling, all five models achieved AUC scores above 0.8 on the majority of tasks, indicating their ability to capture meaningful semantic information from human DNA sequences. These results show that zero-shot embeddings generated by these models are sufficiently informative for supervised learning models, even without fine-tuning. Among the five models, Caduceus-Ph exhibited superior overall performance across multiple human genome classification tasks. For promoter identification in cell lines GM12878, HUVEC, Hela-S3 and NHEK, both DNABERT-2 and Caduceus-Ph achieved statistically significant superior performance. DNABERT-2 showed particular strength in splice site prediction, significantly outperforming other models in both donor and acceptor identification tasks with AUCs of 0.906 and 0.897, respectively. For transcription factor binding site (TFBS) prediction, Caduceus-Ph consistently outperformed all other models, demonstrating its ability to capture complex regulatory patterns in the human genome.

When compared to our baseline CNN model trained directly on DNA sequences, the DNA foundation models demonstrated clear advantages in several key tasks (Supplementary Data 1, 2). For example, for GM12878 and HUVEC cell lines promoter identification, DNABERT-2, GROVER, Caduceus-Ph, and NT-v2 all showed significant improvements over the baseline CNN. In enhancer identification tasks. Overall, these four DNA foundation models demonstrated statistically significant performance gains compared to the baseline CNN.

Sequence classification: multispecies genome regions
In multispecies genome classification tasks using mean token pooling (Table 2), HyenaDNA demonstrated unexpected advantage in Arabidopsis promoter identification, achieving statistically significant superior performance for both TATA and Non TATA promoters with AUC scores of 0.961 and 0.955 respectively, despite being pre-trained exclusively on human genomes. This suggests HyenaDNA’s architecture enables effective transfer of learned representations to plant genomes. For the human versus worm classification task, Caduceus-Ph and GROVER demonstrated the strongest performance (AUCs: 0.992 and 0.984), with both models achieving statistically significant advantages over other models. In bacterial promoter identification, GROVER showed exceptional performance for R. Capsulatus (AUC: 0.715) while no model achieved clear statistical superiority for B. amyloliquefaciens, suggesting challenges in generalizing to more distantly related prokaryotic genomes.

Compared to the baseline CNN model, DNA foundation models generally underperformed for these multispecies tasks (Supplementary Data 1, 2). This pattern was particularly evident for the almost all datasets in Table 2, as most of the times all foundation models were outperformed by the baseline CNN. Such performance gap could be expected, as the baseline CNN is specifically optimized for each dataset with all parameters updated during training, while foundation models rely solely on their pre-training to generate embeddings. In tasks involving multispecies genomes, such discrepancy may have been amplified.

Sequence classification: human & multispecies epigenetic modification
We further investigated the performance of the foundation models on several epigenetic modification prediction tasks. For human epigenetic modifications, specifically the detection of 5-methylcytosine (5mC) and N6-methyladenosine (6mA), using mean token pooling, NT-v2, Caduceus-Ph, and GROVER consistently demonstrated strong performance (Table 3). All three models achieved statistically significant superiority for 5mC detection, with AUCs of 0.738 (NT-v2), 0.783 (Caduceus-Ph), and 0.744 (GROVER) respectively. A similar pattern was observed for 6mA detection, where NT-v2, Caduceus-Ph, and GROVER again outperformed the other models.

On yeast epigenetic mark prediction, among the foundation models, DNABERT-2 showed especially robust performance across the full range of yeast epigenetic marks and performance on the 4mC detection tasks across multiple species is detailed in Table 3. As noted in our Methods, the interpretation of results for the eukaryotic 4mC datasets (A. thaliana, C. elegans, D. melanogaster) warrant particular consideration due to the nature of their annotations, and these tasks are thus viewed as exploratory for model capabilities on such data.

With these considerations, the 4mC detection tasks (Table 4) showed that NT-v2 achieved notable performance for A. thaliana (AUC: 0.633), C. elegans (AUC: 0.649) and D. melanogaster (AUC: 0.652). Caduceus-Ph demonstrated strong results for the E. coli dataset (AUC: 0.628), where 4mC annotations are more directly supported. For C. elegans, both NT-v2 and GROVER showed significantly better performance than other models.

Compared to the baseline CNN, DNA foundation models generally underperformed on the human and multi-species epigenetic tasks (Supplementary Data 1, 2). However, this trend was reversed for yeast epigenetic mark prediction, Multiple DNA foundation models demonstrated a clear and significant advantage over the baseline CNN. In contrast to other epigenetic tasks, for most yeast datasets, several foundation models outperformed the task-specific CNN. For instance, in predicting Yeast H3, H3K4me1, and H3K79me3 marks, five and four foundation models respectively showed superior performance. This suggests that the features learned during pre-training are particularly effective for capturing the signals related to yeast histone modifications. Overall, epigenetic prediction proved more challenging for all models, with AUC scores typically 10–15% lower than for tasks distinguishing functional genomic regions.

Lastly, in the rest five classification tasks with more than two classes (Supplementary Table 2), we observed notable performance variations among the foundation models. HyenaDNA demonstrated exceptional performance in the Regulatory Region Type classification task, achieving an accuracy of 0.83, significantly outperforming the rest. Similarly, for Splice Site Type classification using NT’s dataset, HyenaDNA achieved the highest accuracy of 0.563. For Covid Variants classification, DNABERT-2 and GROVER take the lead, while for Enhancer Strength prediction, GROVER, DNABERT-2, and Caduceus-Ph all performed comparably with accuracies over 0.7. These results suggest that different model architectures have varying strengths in multi-class discrimination tasks, with HyenaDNA’s architecture appearing particularly well-suited for capturing complex patterns necessary for discriminating between multiple regulatory region types.

Gene expression prediction
We evaluated DNA foundation models and revealed several key findings about their effectiveness in predicting gene expression from genomic sequences. Across all foundation models tested, random forest regression significantly outperformed XGBoost, achieving better correlation and lower mean squared error (MSE) between predicted and actual gene expression values (all p < 0.001) (Supplementary Table 3), regardless of whether they used shorter or longer input sequences. Based on these results, we use random forest regression to predict gene expression from DNA foundation model zero-shot embeddings, and our subsequent conclusions are all based on random forest regression outcomes.

On average, the Pearson correlation between predicted and actual gene expression was modest, ranging from 0.114 to 0.123 across models using 6k bp sequences inputs (or corresponding center-cropped sequences for GROVER) (Table 4). This suggests that zero-shot embeddings from DNA foundation models, without fine-tuning, provide limited information for predicting subject-specific tissue-level gene expression. Within the short sequence models, Caduceus-Ph and HyenaDNA exhibited the highest average correlations, significantly outperforming DNABERT-2 and GROVER. GROVER, which used a shorter 2048 bp input, achieved a slightly lower average correlation than the other short-sequence models. Among long sequence models, Enformer achieved similar performance to Caduceus-Ph. HyenaDNA-450K demonstrated significantly better performance than Enformer (p < 0.01). When comparing HyenaDNA and Caduceus-Ph with their longer sequence input counterparts (HyenaDNA-450K and Caduceus-Ph with longer sequence inputs) on the same set of genes, we observed that extending the sequence length significantly improved performance for HyenaDNA (p < 0.001), while the improvement for Caduceus-Ph was not statistically significant (Supplementary Table 4).

The distribution of gene prediction correlations showed similar patterns across all models (Fig. 3). While average correlations were modest, certain genes were consistently predicted with substantially higher accuracy, with correlations exceeding 0.4 and reaching up to 0.9. This indicates that while zero-shot embeddings may not capture the regulatory complexity for all genes, they successfully identify and model genes whose expression is tightly controlled by their local sequence.

Analysis of the best-predicted genes revealed a remarkable consistency in model performance, identifying a core set of highly predictable genes regardless of model architecture (Supplementary Data 3). With 6k bp inputs, the gene CUTALP was consistently the top-predicted gene across all five short-sequence models, with correlations exceeding 0.89. Other highly predictable genes that consistently ranked among the top 10 across these models included CPNE1, NT5C3B, DDX11, and PEX6. For the long sequence models, which analyzed an extended genomic context of up to ~196 Kbp, DDX11 and HLA-DRB5 consistently achieved the highest correlations across all three models. Other genes showing strong predictability with longer sequences included DHFR and CD151. Of particular relevance to blood tissue, CD151 is a tetraspanin family member that plays important roles in platelet activation and aggregation, making it highly relevant to blood function. DHFR (dihydrofolate reductase) is essential for DNA synthesis and is targeted by several drugs used to treat blood disorders, further highlighting the biological relevance of these predictions.

While our results identify specific genes whose expression patterns appear more predictable from sequence data, the overall performance of these zero-shot embeddings remains modest. Therefore, these findings should be considered preliminary indicators of model capabilities.

Variant effect quantification
Our analysis of zero-shot embeddings for variant effect prediction reveals that model performance is highly task-dependent, creating a clear distinction between performance on pathogenic variant identification, and putative causal QTL classification. For the task of distinguishing pathogenic from common SNPs, the transformer-based foundation models NT-v2 and Caduceus-Ph surprisingly emerged as the top performers. NT-v2 was particularly dominant, achieving the highest average test AUC of 0.73 and the largest effect size (Cohen’s d: 0.88), substantially outperforming all other models, including those trained on functional tracks like Enformer (AUC: 0.69, Cohen’s d: 0.73) and Sei (Table 5). This suggests that the pre-training objectives of these specific foundation models are exceptionally well-suited for capturing the subtle, sequence-level patterns that differentiate pathogenic variants.

In contrast, for the putative causal QTL variant effect benchmark to distinguish putative causal QTLs from non-causal variants, models explicitly trained to predict genomic tracks (AlphaGenome, Enformer, and Sei) held a clear and consistent advantage (Table 6). AlphaGenome was the standout model, achieving the highest AUC and Cohen’s d scores across all four QTL types (e.g., AUC = 0.80 for sQTLs, AUC = 0.86 for ipaQTLs). Enformer and Sei also significantly outperformed the general DNA foundation models; for example, in eQTL prediction, Enformer achieved an AUC of 0.77, compared to the best-performing general foundation model, Caduceus-Ph, at 0.65. A notable finding for both Enformer and Sei was that their internal hidden states were often as predictive as their final processed output tracks; for instance, in sQTL prediction, Sei’s hidden state representation achieved an AUC of 0.65 compared to 0.63 for its output tracks.

Among the general-purpose foundation models, Caduceus-Ph again demonstrated the most consistent and robust performance, serving as a reliable baseline. The performance of the other foundation models was more erratic. This instability was particularly pronounced on the QTL tasks with the smaller sample sizes (ipaQTL, paQTL and sQTL) where several models, including DNABERT-2, HyenaDNA, Caduceus-Ph with long sequence input, and GROVER, may yield an average AUC slightly below 0.5. This result highlights the challenges of benchmarking on small datasets with strict chromosome-based splits. Our nested cross-validation design makes this variance explicit: the detailed results for each test chromosome group (Supplementary Data 4) reveal that performance on these small datasets is highly sensitive to the choice of holdout chromosomes.

Finally, a direct comparison between the short- and long-sequence versions of HyenaDNA and Caduceus-Ph didn’t suggest a clear advantage for using longer inputs in this context. For the pathogenic SNP task, the short-sequence version of Caduceus-Ph (AUC: 0.70) clearly outperformed its long-sequence counterpart (AUC: 0.62). This suggests that for a single nucleotide change, a larger genomic context may dilute the local signal when using this embedding subtraction methodology.

Pre-training experiment
We investigated the impact of pre-training dataset diversity by re-pretraining HyenaDNA on DNABERT-2’s multi-species dataset comprising 135 species across 6 taxonomic categories, and examined their performance on the 57 sequence classification benchmark datasets. We maintained the hyperparameters comparable to the original Hyena-1K checkpoint referring to their descriptions, though our implementation may slightly differ from the original model’s training scale in number of training steps and batch size as these were not clearly stated.

The newly pre-trained model, which leveraged the diverse multi-species dataset, demonstrated significant performance gains over the original human-genome-trained HyenaDNA-1K on most datasets. It achieved statistically significant improvements in 14 of the 49 evaluated datasets, particularly in areas requiring cross-species generalization and diverse epigenetic pattern recognition (Supplementary Data 5). For instance, in human epigenetic modification tasks, the AUC for 5mC detection increased from 0.707 to 0.749, and in the cross-species human versus worm genome classification, the AUC improved from 0.968 to 0.984. This advantage also extended to multi-species epigenetic datasets such as C. elegans 4mC detection and various promoter identification tasks. On the other hand, the original HyenaDNA-1K achieved superior AUC on 3 out of the 49 datasets, related to human enhancer identification, human open chromatin region identification and yeast epigenetic mark prediction. This suggests that specialized human genome pre-training still retain a few benefits from certain tasks.

These findings provide evidence that the architecture of HyenaDNA is fundamentally robust on training datasets, and multi-species pre-training may enhance its generalizability for diverse tasks and cross-species generalization. Our results suggest that multi-species pre-training should be considered a critical design choice. Future DNA foundation models can benefit from exploring optimal species sampling strategies and taxonomic diversity in pre-training.

TAD region recognition
A key question for foundation models is whether their internal representations, such as self-attention patterns, learn biologically meaningful, higher-order genomic structures without explicit training. To investigate this, we assessed NT-v2’s capability to recognize topologically associating domains (TADs) by analyzing attention patterns in the self-attention mechanism. Unfortunately, such analysis is currently constrained by model architecture; other foundation models evaluated in our study either lack accessible attention mechanisms (HyenaDNA and Caduceus-Ph), accessible attention matrices (DNABERT-2), or have input length limitations that preclude meaningful analysis of large genomic structures like TADs. We processed 1500 TAD-centered sequences and 1500 randomly selected background sequences through NT-v2, extracting and averaging attention matrices across all layers and attention heads. The attention weight at any coordinate (x, y) indicates how strongly the query token at position y attends to the key token at position x during self-attention computation. Therefore, if NT-v2 can recognize TAD structures inherently, we would expect the attention weights to increase in the central region (approximately token positions 300–700), resulting in a vertical band of positive values in the average attention matrix when key tokens fall within the TAD region. But for random background sequences, the average attention matrix would not show any clear pattern. Hence the differences between these two attention matrices will also show a vertical band.

Supplementary Fig. 2 displays the heatmap of the difference between attention matrices for TAD-centered versus background sequences. We observed no distinct patterns in the attention difference matrix. The heatmap shows predominantly uniform values near zero, with only minimal variation along the diagonal where tokens attend to themselves. This suggests that NT-v2’s self-attention mechanism does not recognize TAD boundaries without specific fine-tuning.

Our analysis represents the current extent of attention mechanism interpretability for DNA foundation models associated with chromatin structures, and reveals current DNA foundation models’ limitations on interpretability. This finding highlights a critical challenge and a vital area for future research, and extending such interpretability analyses is crucial for moving these models from “black box” predictors to tools that can generate essential biological insights.

Runtime analysis
Our runtime analysis revealed computational efficiency profiles across the DNA foundation models on a single A100 GPU (Fig. 4). Despite the logarithmic scale of sequence lengths, none of the models demonstrated clear quadratic scaling, suggesting that attention computation was not the dominant factor within our experiment of batch size equal to 1. However, testing with larger batch sizes to potentially observe quadratic scaling would require multiple and more powerful GPUs, which we leave for future work. HyenaDNA models exhibited the most favorable scaling characteristics, particularly for longer sequences. The HyenaDNA-160K model maintained relatively stable performance for shorter sequences (below 2K nucleotides), while HyenaDNA-1M showed excellent scaling even at sequence lengths approaching 500K nucleotides. This efficiency can be attributed to HyenaDNA’s architecture that leverages long convolutions with implicit parameterization and data-controlled gating. DNABERT-2 demonstrated competitive runtime for sequences approaching its recommended length limit (~2K nucleotides), beyond which we observe a dramatic runtime increase, indicated by the spike in the orange line. NT-v2 consistently required the highest runtime across all sequence lengths, ~2.5s for 100 replications, which reflects its substantially larger model size of 500M parameters, yet the runtime is highly stable. GROVER demonstrated the fastest runtime among all tested models for sequences within its supported length range (2K nucleotides), while Caduceus-Ph-131K showed moderate performance with gradually increasing runtime as sequence length increased.

Methods
To evaluate DNA foundation language models comprehensively, we identified the three most recent state-of-the-art DNA foundation language models, including DNABERT-2, Nucleotide Transformer version-2, HyenaDNA, Caduceus, and GROVER. These foundation models take DNA sequence as input, tokenize into sequence of tokens, and generate embeddings of fixed dimension for each token after passing multiple layers. In the following, we will briefly describe these three models.

DNABERT-2 has the network architecture similar to Bidirectional Encoder Representations from Transformers (BERT), which usually contains a positional embedding layer added to input embeddings, and a series of encoders each consisting of a multi-head self-attention layer and a feedforward network. It is pre-trained using the masked language modeling approach on genomes from 135 species, including the human reference genome. DNABERT-2 tokenizes DNA sequences by the Byte Pair Encoding (BPE) method, which is an iterative algorithm that searches for nucleotides combinations and builds the vocabulary at the same time; it makes no assumption on fixed words and grammars, so each input sequence is independently tokenized merely based on its pattern. It is worth noting that the number of tokens in the tokenized sequence is not fixed in DNABERT-2. DNABERT-2 modifies the architecture of BERT by using Attention with Linear Biases (ALiBi) instead of positional embedding layer. DNABERT-2 has about 117 million trainable parameters, the output embedding dimension is 768. There is no hard limit on the input sequence length, although the runtime is still quadratically increasing with sequence length.

Nucleotide Transformer Version 2 (NT-v2) is also based on the BERT architecture, and it is pre-trained using the masked language modeling approach on genomes from 850 species, including the human reference genome. To tokenize DNA sequence, NT-v2 employs the 6mers tokenization method that uses a sliding window of size 6 and reads every 6 nucleotides; if there are leftover elements at the end of sequence, nucleotides will be tokenized individually into {A, T, C, G, N}. Therefore, the number of tokens produced by the tokenizer will be ~1/6 of DNA sequence length. NT-v2 modifies BERT by replacing the learned positional embeddings with the rotary embeddings, which rotates the embeddings output by each attention layer based on the token’s position, and the Swish activation without bias. These modifications reduce the number of model parameters in Nucleotide Transformer Version 1, and thus reduce the computation cost. The largest NT-v2 model has around 500 million trainable parameters, the output embedding dimension is 1024, and the input sequence length limit is 12,000 nucleotides.

HyenaDNA differs from the architectures of DNABERT-2 and NT-v2 by eschewing the attention mechanism in favor of a decoder-based architecture. HyenaDNA is pre-trained exclusively on the human reference genome using an next nucleotide prediction approach. The key component of this model is the Hyena operators, which integrate long convolutions with implicit parameterization and data-controlled gating. Benefiting from this architecture, HyenaDNA can process very long DNA sequences with fewer model parameters than attention-based transformer. This enables a straightforward tokenization approach in HyenaDNA, where each nucleotide is treated as an individual token. HyenaDNA can also perform in-context learning such as soft-prompting, and details can be found in its original article. The largest HyenaDNA model has around 30 million trainable parameters, the output embedding dimension is 256, and the input sequence length limit is one million nucleotides.

Caduceus-Ph is a long-range DNA foundation language model that leverages the Mamba DNA architecture—an extension of selective state space models (SSMs)—to capture the unique structural and contextual properties of genomic sequences. It introduces two key innovations: bi-directional sequence modeling via the BiMamba block, which enables the incorporation of both upstream and downstream genomic context, and reverse complement (RC) equivariance achieved through a dedicated module that respects the intrinsic reverse complementary nature of DNA. In our implementation, RC handling is performed via a post-hoc conjoining approach, simplifying embedding extraction and facilitating the ensembling of predictions across strands at inference time. Unlike traditional attention-based architectures, Caduceus-Ph operates at nucleotide resolution and scales efficiently to long sequences without incurring quadratic computational overhead. The largest Caduceus-Ph configuration comprises ~35 million trainable parameters, features an output embedding dimension of 256, and maximum input sequence of 131 thousand nucleotides.

GROVER is a DNA foundation language model designed to learn the intrinsic “grammar” of the human genome by capturing both local token properties and long-range contextual dependencies. Built on a transformer encoder architecture akin to BERT, GROVER comprises 12 transformer layers that include multi-head self-attention, feedforward networks, and normalization layers. Its novelty lies in the use of byte pair encoding (BPE) to tokenize genomic sequences—the vocabulary is iteratively optimized over 600 cycles to balance token frequency and capture the heterogeneous sequence composition of the genome. Trained with a masked token prediction objective, GROVER not only reconstructs masked regions but also develops token embeddings that encapsulate key genomic features such as nucleotide frequency, GC content, and sequence length. The GROVER configuration comprises ~117 million trainable parameters, features an output embedding dimension of 768. It supports an input sequence length limit of 512 tokens. Due to the variable-length nature of BPE tokens which can range from single nucleotides to longer 6mers, the corresponding string nucleotide sequence length can vary, typically averaging around two thousand.

Benchmarking datasets
Sequence classification datasets. To unbiasedly evaluate the foundation models, we collected and curated datasets from five distinct sequence classification sources. These datasets were selected to reflect a wider range of potential downstream tasks, ensuring they are both challenging and achievable. While these datasets are not genome-wide in terms of exhaustive coverage, they have been carefully curated by focusing on sub-sampled genomic regions that capture local regulatory signals within controlled windows, each with specific sampling strategies to reduce noise and redundancy.

DNase-I hypersensitive sites detection: The datasets used in the study consist of positive DNA sequences for the 280 Dnase I hypersensitive sites (DHS), and negative sequences for the 737 non-Dnase I hypersensitive sites, originally collected by ref. and has been widely used for benchmarking the detection of DNase I hypersensitive sites. The positive sequences are experimentally identified DNase I hypersensitive sites (DHSs), which are markers for open chromatin and active regulatory elements. The 737 negative sequences are genomic regions verified as non-DHS. Identification of the DNA sequences containing DHS is crucial for detecting DNA regulatory regions, as DHS is indicative of genomic regulatory regions. The datasets consist of a curated set of human genomic regions representing a specific challenge set designed for direct comparison against prior computational methods. sequence length ranges from 225−275 base pairs.

5mC and 6mA modifications detection in human: We use the two benchmark datasets to predict 5-methylcytosine (5mC) and N6-methyladenosine (6mA) modifications in human (Homo sapiens) DNA in this study. The datasets were originally collected by ref. which have been subsequently adopted by relevant tasks. These datasets consist of short, 41bp sequences centered on specific cytosine or adenine sites in the human genome. Both the positive (methylated) and negative (non-methylated) samples are derived from experimental profiling studies. Sequences with high similarity were computationally filtered out to reduce redundancy and potential bias.

Promoter identification in multiple species: We use 8 datasets from a study that established a benchmark for identifying promoters across multiple species, covering 4 distinct organisms including human (4 different cell lines of GM12878, NHEK, HeLa-S3, HUVEC), B. amyloliquefaciens, R. capsulatus, and Arabidopsis (TATA and non-TATA). These datasets focus on genomic regions proximal to transcription start sites (TSS). The positive samples are promoter sequences derived from experimentally-based ENCODE annotations for human cell lines, as well as established promoter sequences for bacterial and plant species. The negative samples for this task were computationally generated by identifying non-promoter genomic regions that have the highest sequence similarity to the corresponding positive promoter sequences, creating a challenging classification problem. Sequence lengths vary significantly, from under 100bp in bacteria to over 2000bp in human cell lines.

N4-methylcytosine (4mC) site detection in multiple species: For this task, we utilized six distinct datasets that the authors collected to benchmark, corresponding to E. coli, C. elegans, G. pickeringii, G. subterraneus, D. melanogaster, and A. thaliana, consist of 41bp sequences centered on potential 4mC sites. The positive samples are experimentally identified 4mC sites collected from public databases. The negative samples (non-4mC sites) were generated by randomly sampling sequences from the respective genomes. It is important to note that while 4mC is a well-established modification in prokaryotes like E. coli and the Geobacter species, for eukaryotic species (A. thaliana, C. elegans, D. melanogaster), 4mC site annotations are often reliant on computational predictions rather than direct, genome-wide high-stringency chemical mapping due to ongoing research into its prevalence and detection []. This distinction is important for the interpretation of model performance on predicting sequences containing 4mC sites.

Genomic Benchmarks dataset collection: We also downloaded a suite of datasets from the Genomic Benchmarks collection, which provides benchmarks for classifying key regulatory elements. The tasks from this collection include identifying promoters, enhancers, and open chromatin regions in humans, as well as distinguishing between coding and intergenic DNA and differentiating sequences between species (human and C. elegans). For these datasets, positive samples are derived from experimentally-validated sources like the FANTOM5 project and the Ensembl Regulatory Build. The negative samples consist of randomly sampled, length-matched sequences from non-overlapping genomic regions, and sequence lengths typically range from 200 to over 2000 base pairs.

We also adopted a large collection of datasets from the downstream task benchmarks in DNABERT-2 and NT-v2. These datasets cover a diverse range of tasks including promoter (proximal and core), transcription factor binding site, splice site, and epigenetic mark prediction across species like human, mouse, and yeast. The majority of these tasks are region-specific, focusing on sequences of a few hundred base pairs centered on the functional element. A common paradigm in their construction is the use of experimentally-derived positive samples (e.g., from ENCODE ChIP-seq), while negative samples are sometimes generated through methods like dinucleotide shuffling, random sampling with matched GC-content, or even adversarial example generation to increase task difficulty. Since these datasets were originally pre-processed into specific splits, we re-shuffled them to ensure a fair and rigorous comparison in our framework.

With all the datasets we collected, there are in total 57 datasets included in this study, including 52 binary classification datasets and 5 multi-class datasets. A key characteristic of these 57 classification tasks is their focus on short sequences (41 to ~2000 bp), a design intended to test local signal recognition. Detailed descriptions and our naming of all 57 datasets can be found in Supplementary Note 1. In Supplementary Note 2, we also annotated each dataset on whether the sequences are relied on experimentally validated or golden standard sources. The specific training size, testing size, class label distribution, and details of sequence lengths for all datasets used in our study can be found in Supplementary Data 6.

To facilitate a systematic analysis, we categorized these 57 datasets into four distinct categories based on the nature of their respective classification tasks: (1) Human Genome Sequence Region Classification: This category encompasses tasks such as identification of transcription factor binding sites, promoter regions, and other functional elements within the human genome. (2) Multi-Species Genome Sequence Region Classification: These tasks involve distinguishing genomic regions across different species, for example, differentiating between human and Caenorhabditis elegans (worm) genome sequences. (3) Human Genome Epigenetic Trait Classification: This group includes tasks related to identifying epigenetic modifications specific to the human genome, such as detection of N4-methylcytosine (4mC) sites. (4) Multi-Species Genome Epigenetic Trait Classification: These tasks focus on identifying and classifying epigenetic traits across multiple species’ genomes. We present and analyze our findings separately for each of these four dataset categories, allowing for a comprehensive assessment of the models’ capabilities and limitations in various genomic classification scenarios. It is important to note that the input sequences for all 57 classification datasets, which range from 41bp to ~2000bp, fall well within the maximum input length of all foundation models evaluated in these tasks. Therefore, a direct comparison was possible without sequence truncation, ensuring that each model was evaluated on identical, biologically relevant inputs for these tasks.

Gene expression prediction datasets. To evaluate DNA foundation models in predicting gene expression, we used the GTEx v8 dataset, which provides gene expression measurements across multiple human tissues. We focused on whole blood (610 subjects, 21,004 genes) tissue. We generated subject-specific DNA sequences by incorporating individual genetic variants from whole-genome sequencing VCF data into the human reference genome (GRCh38/hg38). For each gene, we extracted sequences spanning 3000 base pairs upstream and downstream of the transcription start site (TSS), resulting in ~6000-nucleotide sequences that includes subject-specific SNPs and insertion/deletion variants. The variants have been phased so both alleles are generated. Consequently, for each gene we have generated {subject-specific sequences on both alleles, subject gene expression} pairs.

To consider models which are capable of processing longer sequences, we created an extended-context dataset with sequences spanning TSS ± 98K base pairs. We included sequences of such length for two key reasons: (1) Some versions of HyenaDNA and Caduceus can process sequences significantly longer than 6000 bps, and we wanted to assess whether they generate improved zero-shot embeddings with extended genomic context; and (2) Enformer, another transformer model that has demonstrated success in predicting gene expression, accepts sequences up to 196,608 bps long, allowing us to benchmark it against DNA foundation models.

For each gene, we regressed the gene expression on covariates provided by GTEx v8, including demographic variables such as sex, probabilistic estimation of expression residuals (PEER) factors, and top genotype principal components used to account for population structure. The residuals from this regression, representing covariate-corrected gene expression levels, served as the quantitative trait for our prediction models to learn from the sequence data. The 6000 bps dataset includes all genes. For the long sequence dataset, due to computational constraints, we randomly selected 1000 genes and excluded genes not present in whole blood tissue dataset or genes whose TSS ± 98K exceed chromosome boundaries, making the final total of 768 genes.

Variant effect quantification datasets. To directly assess the capacity of DNA foundation models to represent genetic variants within their zero-shot embeddings, we developed two distinct benchmarks: one focused on differentiating pathogenic from common variants, and another on identifying functional quantitative trait loci (QTLs).

Our first benchmark utilizes data from the Genomics Long-Range Benchmark, which provides two distinct categories: pathogenic single nucleotide polymorphisms (SNPs) associated with disease phenotypes and common SNPs frequently observed in the population. Following the methodology of generating paired sequences, we constructed a reference sequence and a corresponding alternative sequence for each SNP by replacing the central nucleotide. This process resulted in a dataset of {reference sequence, alternative sequence} pairs, where each pair corresponds to a single SNP and is labeled as either pathogenic or common. To accommodate models with different input length capabilities, we created two versions of the dataset. For models supporting longer contexts, we extracted 196,608bp sequences, resulting in 22,222 pathogenic and 17,374 common sequences. For models with shorter input limits, we used 6000bp sequences, yielding 22,239 pathogenic and 17,398 common sequences. In both cases, a small number of SNPs were omitted where the required sequence length exceeded chromosome boundaries. The purpose of this benchmark is to evaluate if the change between the reference and alternate allele embeddings is distinct enough to distinguish pathogenic from benign genetic variation.

As a second benchmark, we adopted the curated putative causal QTL dataset from the Borzoi study. The dataset contains high-confidence causal variants, was originally derived from GTEx v8 and has been statistically fine-mapped. We focus on QTLs from the whole blood tissue. In total, there are 1896 expression QTLs (eQTLs), 540 splicing QTLs (sQTLs), 116 intronic-polyadenylation QTLs (ipaQTLs), and 142 polyadenylation QTLs (paQTLs). A key strength of this dataset is its rigorously constructed negative set, where each putative causal QTL is paired with a carefully matched non-causal variant based on characteristics such as distance to the relevant functional site and the expression level of the associated gene. For each variant, this resulted in a {reference sequence, alternative sequence} pair, labeled according to its specific QTL classification (e.g., eQTL, sQTL) or as a non-causal variant. Similar to the pathogenic versus common SNP benchmark, we constructed sequence datasets of both long (196,608bp) and short (6000bp) sequences, again filtering out a small number of SNPs due to chromosome length limits.

TAD region dataset. To evaluate DNA foundation models’ ability to recognize topologically associating domains (TADs), we utilized the processed dataset of TAD boundaries measured in IMR90 cells following the Enformer study. The dataset contains TAD regions of 2400 base pairs across the human genome, each with an estimated boundary strength score reflecting the magnitude of interaction changes at the boundary. We filtered the dataset to include only TAD regions with boundary strength in the top 5%, and then selected 1500 TAD regions with the strongest boundaries, similar to the workflow in Enformer study. For each selected TAD region, we extracted 6000bp sequences from the human reference genome, positioning the 2400bp TAD region in the center with flanking context on both sides. This generated a total of 1500 TAD-centered sequences. As a control, we randomly sampled 1500 background sequences of identical length (6000bp) from the reference genome, representing genomic regions without known TAD boundaries. This balanced design allowed for direct comparison between TAD-centered sequences and background sequences.

The objective of this dataset was to compare the attention matrices generated by foundation models when processing TAD-centered sequences versus background sequences. Since NT-v2 is the only model in our evaluation that both returns attention matrices and support sufficiently long input sequences, we focused exclusively on NT-v2 for this task.

Benchmarking methods
Sequence classification benchmark. For each classification task, we first established a clear data partitioning strategy. We maintained the training and testing splits from the original works where they were available and clearly defined. For datasets lacking pre-defined partitions, we applied a random 70:30 split for training and testing. For datasets adopted from the downstream task benchmarks in DNABERT-2 and NT-v2, we performed a new random 70:30 split maintaining the original class label distributions, ensuring a standardized and fair comparison across all models in our framework. This re-shuffling is always performed unless there was clear evidence that the original datasets were not manually processed by the authors of DNA foundation models. It should also be noted that due to the pre-processed nature of many classification datasets lacking complete chromosomal information, a uniform whole-chromosome holdout strategy was not feasible.

On the choice of supervised classifiers, we focus on classifiers that require minimal hyperparameter tuning to ensure that performance reflects the quality of DNA foundation model embeddings rather than classifier optimization. Additionally, we select classifiers that inherently handle high-dimensional inputs without requiring dimensionality reduction, as such preprocessing introduces another layer of tuning and potential confounding. Therefore, we included random forest, Naïve Bayes, and elastic-net logistic regression. During training, we performed 5-fold cross-validation that divides the training set into non-overlapping train-validation pairs for hyperparameters tuning, and then reported the testing performance on the test set. The hyperparameter grids for all these classifiers are detailed in the Supplementary Table 5. We evaluated model performance on the test set using AUROC (AUC) and accuracy. AUC serves as our primary measure of performance throughout this work, with the accuracy providing supplementary information.

On the choice of baseline models for comparison against DNA foundation models, we noted that the short-sequence nature of the benchmark datasets is by definition incompatible with genomic models requiring longer and fixed-length inputs, such as Enformer, Sei, and AlphaGenome. Forcing these short sequences into a long-context model via excessive ‘N’ padding would create additional artificial signal structures to handle. Therefore, to establish a relevant, task-specific baseline, we implemented a simple CNN consisting of three convolution layers with pooling and a classification head. The CNN architecture takes 5-dimensional one-hot encoded input corresponding to nucleotides A, T, C, G, N and comprises three 1D convolutional layers with 64, 128, and 256 channels respectively, max pooling after the first two layers, adaptive global maxpooling, and a final linear layer for classification. Unlike the foundation models where embeddings remained frozen during training, the CNN was trained directly on one-hot encoded DNA sequences with all parameters updated during training. This comparison aimed to identify when pre-trained DNA foundation models could outperform task-specific neural networks.

To evaluate statistical differences in performance between models, we applied the DeLong’s test for comparing AUC values. For each dataset (task), pairwise comparisons were conducted between all model pairs using one-sided tests with a significance threshold of α = 0.01. A model’s performance was considered significantly better only when its AUC value was higher and the corresponding test yielded p < 0.01. For the five datasets involving multi-class classification tasks, where DeLong’s test become less adaptable, we used classification accuracy as the primary metric instead.

Sequence classification benchmark: pooling methods. For our model-wise comparison in sequence classification benchmark, we ensured fair evaluation by using the same pooling method across all models. However, given the lack of comprehensive studies on optimal pooling methods for DNA foundation models, we also conducted a systematic comparison of different pooling methods for each individual model. In foundation language models, output pooling methods refer to techniques for generating a single, fixed-dimensional embedding to represent an entire sequence. Despite their importance, there has not yet been a comprehensive study on which output pooling methods are most effective for DNA foundation models. All DNA foundation models in this study create additional special tokens during tokenization, aligning with common practice in natural language processing. We focused on three common output pooling strategies:

Sentence-level summary token method: In the original BERT architecture, a special “CLS” token is appended to the start of every tokenized sequence to serve as a sentence-level summary. After processing through multiple self-attention layers, this token captures the context of the entire input sequence. DNABERT-2, NT-v2, and GROVER implement this approach. Similarly, HyenaDNA and Caduceus-Ph append an end-of-sequence “SEP” token that can be used to represent the whole sequence.

Mean pooling method: This approach calculates the average of all non-padding token embeddings in the sequence to create a unified representation. This potentially captures information from all parts of the sequence equally.

Maximum pooling method: This method selects the maximum value across each dimension from all token embeddings, potentially highlighting the most salient features.

For these comparisons, we employed the same statistical testing framework described earlier, using the DeLong’s test with a significance threshold of 0.01 to determine significant differences in performance.

Gene expression prediction. To evaluate DNA foundation models in predicting gene expression, we employed a regression-based approach using data from the GTEx v8 dataset. As described in the dataset subsection, we worked with two sequence lengths: shorter 6000 nucleotides (applied to all 610 subjects and all 21,004 genes) and longer sequences of up to 196K nucleotides (applied to a subsample of 768 genes).

For the shorter sequences, we generated zero-shot embeddings using DNABERT-2, NT-v2, HyenaDNA, Caduceus-Ph, and GROVER. Due to GROVER’s architecture, its input is limited to 512 tokens generated via Byte Pair Encoding (BPE). Since the average BPE token length is approximately four nucleotides, this corresponds to an input of roughly 2048 base pairs. To accommodate this, we extracted the central 2048bp from each 6000bp sequence, ensuring the TSS remained centered. GROVER’s performance on this shorter input was then directly compared against the other models that utilized the full 6000bp sequence. For the longer sequences, we generated embeddings using HyenaDNA-450K, Caduceus-Ph, and Enformer. Caduceus-Ph was provided with the central 131K nucleotides due to its input length limit. For Enformer, we used the model’s final hidden state as its sequence embedding and performed mean pooling over all regions.

For each subject and gene across all models, we first generated zero-shot embeddings on both alleles separately, then calculated the mean of these embeddings to represent the joint genetic information from the two alleles.

Since gene expression is represented as a continuous value, we approached this as a regression task. We experimented with random forest and XGBoost as regression models. For each gene, we train a separate regression model, with the subjects split into training and testing sets at a ratio of 75:25, respectively. During training, we performed 5-fold cross-validation using mean squared error as the optimization metric for hyperparameter tuning, with the complete hyperparameter grid detailed in Supplementary Table 6. This regression framework was applied to the whole blood tissue with 610 subjects and 21,004 genes.

For evaluation, we used both mean squared error (MSE) and Pearson correlation coefficient between predicted and actual expression values for each gene, which quantifies how effectively each DNA foundation model could capture the genomic features that influence gene expression levels. We conducted multiple comparative analyses: First, we compared performance metrics between different DNA foundation models on short sequences to assess their relative effectiveness. Second, we compared metrics between Enformer and DNA foundation models on longer sequences to evaluate the benefits of specialized architectures. Third, we performed paired comparisons on the same genes between HyenaDNA (shorter sequences) versus HyenaDNA-450K (longer sequences), and Caduceus-Ph on shorter versus longer sequences, to determine whether extended genomic context improves prediction accuracy. To assess statistical significance in these comparisons, we applied the paired Wilcoxon signed-rank test, which allows for robust non-parametric comparison of paired samples. Additionally, we identified and analyzed the subset of genes that were best predicted by each model to understand potential model-specific strengths in capturing particular genomic features or regulatory mechanisms.

Variant effect quantification benchmark. To directly assess the capacity of DNA foundation models to represent genetic variants within their zero-shot embeddings, we developed a unified robust quantitative framework. Our general methodology, inspired by the Enformer study, was to first compute a unified effect vector for each variant by subtracting the zero-shot embedding of the reference sequence from that of the alternative sequence (embedding(alt) - embedding(ref)). This standardized approach was applied across all models in our comparison, including the general DNA foundation models (DNABERT-2, NT-v2, GROVER, HyenaDNA, Caduceus-Ph), as well as Sei, Enformer, and, AlphaGenome. We then trained a random forest classifier on these high-dimensional effect vectors to distinguish between the variant labels (e.g., pathogenic vs. common).

To prevent data leakage from effects such as linkage disequilibrium, we employed a strict chromosome-based holdout strategy for evaluation. However, recognizing that a single arbitrary chromosome split could lead to high variance in performance estimates, especially for smaller datasets like ipaQTL and paQTL, we implemented a nested cross-validation scheme to ensure a more robust and stable evaluation. We established three distinct test sets comprised of non-overlapping chromosome groups: Group 1 (chromosomes 3, 6, 9, 12, 16, 18, 19, 21), Group 2 (chromosomes 2, 5, 11, 14, 17, 20, 22, and X), and Group 3 (chromosomes 1, 4, 7, 8, 10, 13, 15). For each of these three main folds, the corresponding group served as the test set. The remaining chromosomes served as a training-validation set, on which we performed a 4-fold inner cross-validation, also partitioned by chromosome, to tune the random forest hyperparameters. Once the optimal parameters were selected, the classifier was retrained on the entire training-validation set and evaluated on the held-out test chromosomes. The final reported test AUC and Cohen’s d scores are the average of the performances across these three independent test sets.

For the short sequence (6000bp) datasets, we compared DNABERT-2, NT-v2, GROVER, Caduceus-Ph, HyenaDNA, and Sei. Due to model-specific constraints, GROVER received the central 2048 nucleotides and Sei received the central 4096 nucleotides. For Sei, we evaluated two distinct representations: (1) its direct output across 21,907 predictive tracks, and (2) its final hidden state (dimension 15,360) after the spline transformation layer. For the long sequence (196,608bp) datasets, we compared Caduceus-Ph, HyenaDNA (450K checkpoint), and Enformer, with Caduceus-Ph using the central 131,072 nucleotides. For Enformer, we assessed both its final hidden state and its 5313 human-specific output tracks. For the long-sequence QTL datasets, we also included AlphaGenome, providing it the central 131,072 nucleotides of the input and averaging its output tracks over the central 2048bp to derive the effect vector.

For all experiments, no dimensionality reduction was performed on the effect vectors. However, to ensure stable and robust training for models generating relatively high dimensional embeddings (e.g., the hidden states of Sei and Enformer), we adjusted the hyperparameter search space for the random forest. Specifically, we removed the “sqrt” option from the max features tuning grid (Supplementary Table 7), thereby prioritizing the more restrictive “logarithm” option. This standard practice for high-dimensional data helps decorrelate the trees, which stabilizes the training process, reduces the risk of overfitting, and ensures a more reliable performance estimate.

TAD region recognition benchmark. We assessed NT-v2’s capability to recognize topologically associating domains by analyzing attention patterns across TAD-centered and background sequences. We processed 1500 TAD-centered sequences and 1500 control sequences through NT-v2 and extracted attention matrices from all layers and heads. For each sequence type, we averaged these matrices across all layers, heads, and sequences to produce composite attention patterns.

Our analysis focused on the center 400 tokens (corresponding to the 2400bp TAD region) in each averaged attention matrix. We compared attention values within this central region between TAD-centered and background sequences to determine whether NT-v2 naturally allocates higher attention to genomic regions containing strong TAD boundaries.

Pre-training dataset experiment. To investigate the effect of pre-training dataset diversity on zero-shot embedding quality, we pre-trained HyenaDNA using DNABERT-2’s multi-species dataset. This was made possible thanks to the HyenaDNA’s reproducible pre-training pipeline, allowing us to pre-train the model on different datasets while maintaining the same training settings for fair comparison.

For this experiment, we used DNABERT-2’s multi-species dataset, which encompasses genomes from 135 species across 6 categories, totaling 32.49 billion nucleotide bases (approximately 12 times larger than the human genome dataset). Since the DNABERT-2 dataset contains processed sequences of 1000bps, we used the HyenaDNA architecture with input length 1000 (HyenaDNA-1K). We strictly maintained the pre-training settings as specified in the original work, including all model hyperparameters, optimizer configurations. We roughly maintained their training steps, and batch sizes.

We compared our multi-species pre-trained model against the official human-genome pre-trained HyenaDNA-1K on 49 of the 57 sequence classification tasks (excluding 8 datasets with sequences >1000bps). This comparison allowed us to directly assess whether species diversity in pre-training data improves zero-shot embedding quality while controlling for all other variables, including model architecture and training methodology.

Runtime benchmark. To evaluate the computational efficiency of DNA foundation models, we measured the time required to generate zero-shot embeddings across varying sequence lengths. Experiments were conducted on a single A100 GPU, measuring the average forward pass time for randomly generated DNA sequences of length 2k (k ranging from 6−19, spanning 64 to ~500,000 nucleotides). We recorded performance for each model until it reached its sequence length limit or exceeded the A100 memory capacity.

All measurements used a consistent batch size of 1 to accommodate the substantial memory requirements of longer sequences. To reduce measurement error, we performed 100 replications and recorded the total runtime. For DNABERT-2 and GROVER, which use Byte Pair Encoding in tokenization with variable token counts, we estimated the number of tokens as (sequence length / 4). This approximation may introduce slight measurement variations due to potential sequence truncation or padding requirements.

Model configuration selection
For the DNA foundation models we benchmarked, we selected specific model checkpoints based on their capabilities and our experimental requirements. For NT-v2, we used the largest NT-v2-500M model, which showed optimal performance in the original study. With HyenaDNA, we included different checkpoints for different tasks: Hyena-160K for default choice of benchmarks, Hyena-450K in specifically gene expression prediction and variant effect quantification, to leverage its ability to process longer sequences, and Hyena-1K for pre-training experiments to match our pre-processing format. For Caduceus-Ph, we consistently employed Caduceus-Ph-131K across all tasks to maximize context length while maintaining efficiency. DNABERT-2 and GROVER only have one model checkpoint. The embedding dimensions of all DNA foundation models are detailed in Supplementary Table 8. For Enformer, we used the widely-accepted official PyTorch implementation. Sei only has one model checkpoint. For AlphaGenome, we used their provided API to access the regular model checkpoint.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper presents a collection of genomic sequence classification benchmark tasks designed to measure the performance of DNA foundation language models. The study aims to assist future researchers in selecting and optimizing DNA language models. The results showcase the strengths and weaknesses of the three models tested and the comparative efficacy of using different sequence embedding methods. The benchmark is well-conceived and covers four DNA classification tasks over 57 multi-species datasets. This breadth ensures that the benchmark provides a holistic view of the capabilities of DNA language models across various classification challenges. Such a diverse and extensive benchmark is unparalleled in recent genomic benchmarks. Several detailed comments about the work are listed in the following:\nThe scope of models used in this comparison seems limited. The most recent DNA foundation models benchmark paper, Bend (Marin, F. I. et al., 2024), used 15 versions of various models. Including more models, such as GROVER, DEEPSEA, or even newer models like MAMBA, could further validate the evaluation methods and conclusions of the study. Moreover, adding a simple downstream baseline, such as a shallow two-layered CNN, could help validate the proposed benchmarking and provide a direct and obvious comparison to the complex language models.\nThe benchmark performance of models on DNA modification seems to be much lower in accuracy than other methylation prediction/imputation models (e.g., DeepCpG). So, in general, how do foundation models compare to other task-specific models in performance?\nLack of comparison with other classifiers: The study argues that random forest should be utilized due to its minimal tuning requirement (compared with CNN used for zero-shot embeddings in Bend (Marin, F. I. et al, 2024), which requires much more hyperparameters) and its superior performance over XGBoost. However, other simple machine learning models, such as SVM, logistic regression, KNN or Naïve Bayes also generally require minimal parameter tuning. Providing a detailed comparison of more classifier models could help to fully capture the representation in the last hidden embedding, as well as solidify the argument in the superiority of choosing Random Forest over other models, including the CNN used in past papers.\nModels are inconsistently pre-trained on different datasets. To better understand why one model excels another on certain tasks, it could be useful to train all language models from scratch using the same dataset. So we can gain a better insight into the origins of advantages and disadvantages. DNABERT2 and NT-V2 are originally pre-trained on different multi-species sequence datasets. HyenaDNA is pre-trained purely on the human reference genome. Training all models on the same dataset under the same conditions would allow for a more controlled comparison and better assessment of the generalizability of the evaluation methods.\nAlternative pooling methods should be tested and used in a context-dependent manner. The study thoroughly compared the effects of two pooling methods over different models. It demonstrates that mean pooling is generally more effective than summary-level token pooling in generating zero-shot embeddings. However, some models might be intrinsically designed to favor certain pooling methods (for example, NT-v2 with its final averaging pooling layer in pre-training). Since the pooling method can be model-specific, task-specific, and sensitive to sequence length, testing more pooling methods (such as max-pooling, attention-based pooling, or hybrid approaches) could facilitate more intriguing analyses.\nThe study could be more impactful and cutting-edge by extending the test datasets to include longer sequences. This would allow in-depth exploration of the relationship between model performance and sequence length. Transformer-based models are known to be inefficient with long sequences. Test-time training models and MAMBA that scale up to million-bps-long sequences may be compared and discussed in this context.\nAlthough this study focuses on zero-shot embedding to avoid potential bias introduced by fine-tuning, fine-tuning is still widely used and could significantly alter or improve the model’s performance. Future research on how fine-tuning affects performance across different tasks could provide a more complete picture of each model’s capability.\nWhile the study is comprehensive in its evaluation of classification tasks, real-life research questions are not limited to this one type of task. Incorporating more types of tasks, such as the more complex regression-based tasks, could contribute to the generalizability of the proposed method and make this work more impactful.", "reviewer_2": "The manuscript benchmarks three DNA foundation models—DNABERT-2, Nucleotide Transformer version-2 (NT-v2), and HyenaDNA—for genomic sequence classification using zero-shot embeddings. While the authors' approach to comparing the performance of different sequence embeddings is generally appropriate and well-designed, its effectiveness is undermined by the poor quality of datasets chosen for model evaluation. Specifically, the human genome datasets are very limited in size, and some are not based on experimental data. I strongly advise the authors to carefully review the datasets selected for evaluation and select genome-wide, high-quality experimental data. Additionally, to make the manuscript's contribution more impactful, the model comparison should include other models that focus on regulatory signal predictions from genomic sequences.\nMajor points:\nAlthough the datasets cover a diverse range of tasks, the size of each individual human dataset is small and inadequate for a comprehensive comparison of DNA foundation models. A more appropriate approach would be to use whole-genome datasets, as this would provide a more accurate assessment of the models' ability to generalize across diverse genomic sequences.\nThe authors should clarify why their performance results differ from those in previous studies. For example, the original NT-v2 paper (The Nucleotide Transformer: Building and Evaluating Robust Foundation Models for Human Genomics) showed that NT-v2 predicts promoters (All, non-TATA, and TATA) better than DNABERT-2. But in this paper, the authors found that DNABERT-2 performed better for some of the promoter classification tasks. The authors should explain what might have caused this difference.\nOne of the datasets the authors used is 4-methylcytosine (4mC), which is not widely accepted as an epigenetic modification. The paper the authors cited relies on predicted 4mC data rather than experimentally validated annotations, making it a poor choice for model evaluation.\nTo make the paper more impactful, it would be important to extend the comparison of these three models to include models like Enformer and Sei. These would provide a more complete picture of the models' capabilities.\nMinor points:\nThe HyenaDNA model processes longer sequences than the other two models, so the appropriate way to compare them would be to generate encodings from sequences of equal length. This could be achieved by concatenating the encodings generated by the models with shorter input sizes to match the input sequence length of HyenaDNA.\nThe authors mention that in some cases, they maintained the original training and testing split of datasets, but in other instances, this was not feasible. In such cases, a more appropriate approach for comparing models would be to always randomly split the data, otherwise it could unfairly advantage certain models over others.", "reviewer_3": "In this manuscript, the authors presented an unbiased evaluation of existing DNA foundation models. The evaluation focuses on the performance of zero-shot embeddings in genomic classification tasks such as promoter/enhancer/TFBS identification, epigenetic modification prediction, as well as application to model organisms. The authors also briefly discussed the effect of pooling methods on model performance. While the topic proves critical in the context of applying LLM to genomics, several key aspects of model evaluation are missing from the current manuscript. These limitations needs to be properly addressed before the paper can be considered for publication.\nMajor comments:\nWhile the authors presented a comprehensive evaluation of the predictive performance of DNA foundation models, the evaluation on model interpretability is largely missing. Since these models were trained with millions of data points from reference genome, they should be able to capture the important sequence patterns or motifs associated with transcriptional regulation. While the foundation models are generally hard to interpret due to complex architecture of deep neural networks, the authors could focus on the convolutional kernel or attention weights (position weight matrix analysis) and see if they can recapitulate known sequence motifs (for reference, check this review doi: 10.1146/annurev-genom-021623-024727).\nWhile it is straightforward to evaluate model performance on reference genome-based prediction of epigenetic tracks (TF binding, methylation, etc), these tasks are not informative to most biologists who are interested in using DNA foundation models. At the end of the day, these epigenetic tracks can be measured with high-throughput sequencing techniques such as CHIP-seq, ATAC-seq, Bisulfite-seq, etc. Biologists are more likely to use foundation models for predicting variant effect on transcriptional regulation (comparing predictions of reference allele and alternative allele). Therefore, I think the authors should design experiments to compare the performance of DNA foundation models on variant effect prediction. For instance, whether models can accurately predict the expression-change direction of eQTL variants (the authors can check the Enformer paper for datasets and experimental design. doi:10.1038/s41592-021-01252-x)\nAnother question the audience will be interested in is the comparison between DNA foundation models in general versus other widely-used sequence-based models such as Enformer, Sei, etc. Do they generally perform better or worse in terms of the tasks evaluated in this manuscript? While the deep learning in genomics community would like to see a direct comparison between these models in terms of the specific tasks in Result section, I understand adding these analyses can be time-consuming. I encourage the authors to (at minimal) address the question in the Discussion section, discussing the pros and cons of using foundation models versus traditional sequence-based models.\nMinor comments:\nInstead of random splitting of training and testing sets, I suggest the authors adopt whole chromosome holdout (e.g. chromosome 8 and 9 as testing set) to avoid overlap of regulatory region between training and testing set.\nFigure 2, please add asterisks on top of each boxplot denoting the significance level of comparing mean pooling versus CLS pooling.\nThe README file is clear. I think users should be able to reproduce the results from the GitHub repo.\nMinor comment:\nSome file paths are still based on authors' institutional computing cluster (such as \"/rsrch4/home/biostatistics/hfeng3/review_datasets/results_final/ntv2/\" in \"classify_ntv2.py\"). Please revise those paths in the context of GitHub directory"}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors addressed my comments very well. I don't have any more comments and commend the authors for the good work!", "reviewer_2": "I have reviewed the revised version of the manuscript “Benchmarking DNA Foundation Models for Genomic Sequence Classification” as well as the authors’ point-by-point response. While I appreciate that the authors have implemented some important improvements, I find that most of the major concerns outlined in my initial review remain unaddressed.\nThe addition of the GTEx v8 gene expression prediction task, which measures model prediction accuracy of gene expression across individuals, is a well-considered choice for the dataset, however, it improves the evaluation of only one task. For a benchmarking paper, it is critical that the appropriate evaluation of DNA foundation models be conducted using genome-wide, high-quality experimental data for the majority, if not all, of the tasks.\nI appreciate that the authors included Enformer and a simple CNN model as references for model performance comparison. However, Enformer was only added to the gene expression task, whereas it should also be included in other human genome tasks. Additionally, the authors state that Sei's architecture and training are \"highly optimized for variant effect prediction,\" which is not correct. If there are technical difficulties in applying the Sei model as a reference, the authors should instead use other appropriate models available for each specific task.\nMy final concern is that the authors chose to preserve the original training and testing splits of certain datasets, but not others. The authors explain that this choice was made to avoid potential data leakage. While I understand this reasoning, I believe it is the authors’ responsibility, particularly in the context of a benchmarking study, to design datasets in a way that prevents data leakage and allows for fair model comparison. It cannot be overstated that finding new comprehensive datasets, as well as the careful design of previously used datasets, is what makes a benchmarking study valuable and impactful.", "reviewer_3": "The authors presented a revision to the original manuscript. I applaud the authors for their efforts addressing all three reviewers’ comments. Most of my comments were properly addressed except for one: comment 2 regarding the performance of variant effect prediction among DNA foundation models. I think some major issues still need to be resolved before final publication. I encourage the authors to make the following revisions.\nComments about the “variant effect quantification” section:\nGROVER was excluded from this evolution task for no apparent reasons. I noticed that the authors conducted evaluations on all five foundation models except for the task of “variant effect quantification” and “Gene expression prediction”. Please explain in the manuscript why GROVER was excluded from these two tasks.\nThe authors adopted p-value as a metric to evaluate the performance, which is problematic for the following two reasons. (1): The authors did not specify which test was used to derive the p-value: parametric or non-parametric. If parametric test was used, the test p-value is not a good metric for performance evaluation given the scale of L2 distance varies by foundation model (see Figure 5). (2): The test p-value also depends on the sample size. Often times, a relatively small test statistics can result in significant p-value if the sample is large enough. In the case of variant effect prediction, a good model should be able to clearly separate the pathogenic variants from neutral variants. Therefore, I suggest the authors adopt the following two metrics for this specific task: (1) AUROC/AUPRC of the classification task pathogenic vs neutral; (2) Effect size (different in distance between the pathogenic and neutral variants). Note that for effect size to be comparable across models, the pathogenic variant scores need to be first normalized against the neutral variant scores.\nIt is unclear how many pathogenic/neutral variants are included in the benchmark dataset. Please add the count in the method section. If the set is relatively small, the authors should consider all other benchmark variant sets for the purpose of a comprehensive evaluation. For example, the authors can follow the Figure 4d of Nucleotide Transformer paper, adding the eQTL, meQTL, and Clinvar variant sets."}, "author_response": "Response: Thank you for your nice summary and positive evaluation of our manuscript. Your encouraging comments motivate us to further refine and enhance our study. Please find our point-by-point response to your comments below.\n(1) The scope of models used in this comparison seems limited. The most recent DNA foundation models benchmark paper, Bend (Marin, F. I. et al., 2024), used 15 versions of various models. Including more models, such as GROVER, DEEPSEA, or even newer models like MAMBA, could further validate the evaluation methods and conclusions of the study. Moreover, adding a simple downstream baseline, such as a shallow two-layered CNN, could help validate the proposed benchmarking and provide a direct and obvious comparison to the complex language models.\nResponse: We thank the reviewer for helpful and insightful suggestions. In response to the reviewer's suggestion to broaden our model selection, we have expanded our benchmark to include two additional state-of-the-art DNA foundation models, bringing our total to five models evaluated across all 57 datasets:\nWe have added Caduceus-Ph (detailed in the Methods section), a long-range DNA foundation model leveraging the MambaDNA architecture which is an extension of selective state space models. Caduceus-Ph introduces two key innovations: bi-directional sequence modeling via the BiMamba block, and reverse complement equivariance through a dedicated module. To generate zero-shot embedding, we take the average of the embeddings of a DNA sequence and its reverse complement sequence. Caduceus-Ph tokenizes in single-nucleotide resolution, and can handle very long sequences like HyenaDNA. It comprises approximately 35 million parameters with an embedding dimension of 256 and supports sequences up to 131 thousand nucleotides.\nWe also included GROVER (described in the Methods section), a DNA foundation model built on a transformer encoder architecture like BERT, and comprises 12 transformer layers with multi-head self-attention mechanisms. It employs byte-pair encoding (BPE) to tokenize genomic sequences, with a vocabulary optimized over 600 cycles to capture the heterogeneous composition of the genome. GROVER includes approximately 117 million trainable parameters with an output embedding dimension of 768 and supports sequences up to 512 tokens, corresponding to approximately 2 thousand nucleotides due to its variable-length tokenization.\nAs suggested, we implemented a simple downstream baseline CNN to provide a direct comparison with the complex foundation models. While the reviewer suggested a shallow two-layered CNN, we designed our baseline to be as simple as possible while maintaining sufficient capacity for meaningful comparison. Our CNN takes 5-dimensional one-hot encoded input corresponding to nucleotides A, T, C, G, N and consists of three 1D convolutional layers with 64, 128, and 256 channels, max pooling after each layer, adaptive global max pooling, and a final linear layer for classification. This architecture allows for a fair comparison between DNA foundation models' zero-shot embeddings paired with simple classifiers versus task-specific models with all parameters updated during training, thereby evaluating the representational power of pre-trained DNA foundation models.\nThe expanded model comparison is reflected throughout our Results section (Tables 1-6), providing comprehensive evaluation across all 57 datasets and strengthening our findings regarding embedding strategies and architectural trade-offs in DNA sequence analysis. ... [Tables omitted for brevity] ...\n(2) The benchmark performance of models on DNA modification seems to be much lower in accuracy than other methylation prediction/imputation models (e.g., DeepCpG). So, in general, how do foundation models compare to other task-specific models in performance?\nResponse: Thank you for your helpful comments. First, for epigenetic modification prediction tasks, our results show that foundation models indeed underperform compared to specialized models. As shown in Table 3, for 5-methylcytosine (5mC) prediction, the best performing model Caduceus-Ph achieved an AUC of 0.7752, while for N6-methyladenosine (6mA) prediction, GROVER reached an AUC of 0.7744. These values are lower than those reported for specialized models like iDNA-ABF (0.95 for 5mC and 0.91 for 6mA), which is specifically designed for methylation prediction.\nFurthermore, when comparing foundation models to our simple CNN baseline (described in the Methods section), we found that \"all DNA foundation models underperformed in most epigenetic modification tasks compared to baseline CNN, likely reflecting the subtle nature of these modifications that may require full parameter updates rather than frozen embeddings\" (from Results: Human & Multispecies Epigenetic Modification section).\nHowever, for certain genomic classification tasks based on the human genome, zero-shot embeddings from foundation models demonstrated strong performance. As shown in Table 1, for promoter identification across different cell lines, models like DNABERT-2 and Caduceus-Ph achieved AUCs above 0.98, and for transcription factor binding site prediction, Caduceus-Ph demonstrated robust performance with AUCs ranging from 0.8217 to 0.933 across different datasets. These observations align with the expectation that epigenetic patterns are more subtle than those distinguishing functional genomic regions. Specifically, functional genomic regions like promoters typically contain recognizable sequence motifs and characteristic nucleotide compositions that create strong DNA sequence signals. In contrast, epigenetic modifications often occur in diverse sequence contexts and are influenced by complex factors beyond DNA sequence, including chromatin structure, cellular differentiation state, and environmental conditions, thus may be benefited more from task-specific models and careful training.\nLastly, we anticipate that with reasonable fine-tuning, DNA foundation models could potentially achieve performance levels comparable to task-specific models, but this remains to be studied in future work. We have not conducted extensive fine-tuning of these foundation models for specific tasks in this study, as our study specifically evaluates DNA foundation models using zero-shot embeddings with frozen weights, deliberately assessing their inherent representational capacity rather than optimized performance after task-specific fine-tuning. We have added related discussions on Page 24-26.\n(3) Lack of comparison with other classifiers: The study argues that random forest should be utilized due to its minimal tuning requirement (compared with CNN used for zero-shot embeddings in Bend (Marin, F. I. et al, 2024), which requires much more hyperparameters) and its superior performance over XGBoost. However, other simple machine learning models, such as SVM, logistic regression, KNN or Naïve Bayes also generally require minimal parameter tuning. Providing a detailed comparison of more classifier models could help to fully capture the representation in the last hidden embedding, as well as solidify the argument in the superiority of choosing Random Forest over other models, including the CNN used in past papers.\nResponse: We appreciate the reviewer's suggestion to broaden our comparison of classifiers.\nIn our evaluation methodology (Methods section \"Sequence Classification Benchmark\"), we aimed to use a classifier that requires minimal hyperparameter tuning to ensure that performance reflects the quality of DNA foundation model embeddings rather than classifier optimization. Additionally, we focused on models that inherently handle high-dimensional inputs without requiring dimensionality reduction. Following this idea, in our revision, we have explored multiple classifiers such as SVM, Naïve Bayes, XGboost, and elastic-net regression. We applied elastic-net instead of logistic regression due to the high dimensionality of DNA sequence embeddings. We found that random forest outperformed these alternatives across all evaluation metrics (Supplementary Figure 1). This finding aligns with previous research demonstrating that tree-based models often excel on tabular data with complex feature interactions (as cited in reference [22] Grinsztajn et al., NeurIPS 2022). The high-dimensional nature of DNA sequence embeddings, with potential intricate interactions between features, makes random forest particularly well-suited for this application. We have updated the Methods descriptions to justify our choice.\nFor comparison with more complex approaches, we implemented a simple CNN baseline consisting of three convolution layers with pooling and a classification head that was trained directly on one-hot encoded DNA sequences with all parameters updated during training (Methods section). This baseline provides an estimate for evaluating when pre-trained embeddings outperform task-specific baseline models. The results comparing foundation models against the CNN baseline are presented in Supplementary Tables 6 and 7 and discussed throughout the Results section. ... [Tables omitted for brevity] ...\n(4) Models are inconsistently pre-trained on different datasets. To better understand why one model excels another on certain tasks, it could be useful to train all language models from scratch using the same dataset. So we can gain a better insight into the origins of advantages and disadvantages. DNABERT2 and NT-V2 are originally pre-trained on different multi-species sequence datasets. HyenaDNA is pre-trained purely on the human reference genome. Training all models on the same dataset under the same conditions would allow for a more controlled comparison and better assessment of the generalizability of the evaluation methods.\nResponse: We acknowledge this important point about pre-training dataset variability. As noted in our Methods section, the models in our benchmark indeed have different pre-training backgrounds: DNABERT-2 was pre-trained on genomes from 135 species, NT-v2 on 850 species, and HyenaDNA exclusively on the human reference genome.\nWhile ideally, we would pre-train all models on identical datasets for the most controlled comparison, this proved technically unfeasible despite our considerable efforts. These include the large computational resources required for pre-training each model, the extensive hyperparameter optimization, and the significant engineering effort to adapt to model-specific pipeline and to reproduce training environments.\nTo address this concern within our constraints, we conducted a controlled pre-training experiment specifically described in the \"Pre-Training Dataset Experiment\" section of our Methods. We chose HyenaDNA for this controlled experiment because of its reproducible pre-training pipeline, allowing us to pre-train the model on different datasets while maintaining the same training settings for fair comparison. Even with HyenaDNA, we encountered technical challenges requiring careful code modifications to accommodate our hardware infrastructure and newer CUDA versions, but these were ultimately resolvable through extensive efforts. While we also investigated the feasibility of re-training other foundation models, these efforts were precluded by the lack of accessible, documented pre-training scripts and clearly defined training environments in their respective official repositories. Using DNABERT-2's multi-species dataset (135 species across 6 categories, totaling 32.49B nucleotides - approximately 12 times larger than the human genome dataset), we pre-trained a new version of HyenaDNA-1K while \"strictly maintaining the pre-training settings as specified in the original work, including all model hyperparameters and optimizer configurations. Our results, detailed in the \"Pre-Training Experiment\" subsection in Results, show that the newly pre-trained model improved the AUC of human-genome pre-trained HyenaDNA-1K with statistical significance across diverse genomic tasks, particularly in epigenetic modification detection and cross-species classification. For example, in 5mC detection, AUC improved from 0.745 to 0.771, and in human versus worm classification from 0.960 to 0.975.\nInterestingly, the original HyenaDNA-1K maintained superior AUC on specific human regulatory element tasks, particularly those related to enhancer regions and open chromatin structures, suggesting that human-specific pre-training provides advantages for certain human genomic tasks.\nThese findings demonstrate that the architecture of HyenaDNA is fundamentally robust on training datasets, and multi-species pre-training may enhance its generalizability, particularly for cross-species generalization and epigenetic modification detection. These results suggest that multi-species pre-training should be considered a critical design choice for future DNA foundation models.\nWe believe that our controlled experiment provides valuable insights into the impact of pre-training data composition on model performance and have summarized our findings and experiences on Page 22.\n(5) Alternative pooling methods should be tested and used in a context-dependent manner. The study thoroughly compared the effects of two pooling methods over different models. It demonstrates that mean pooling is generally more effective than summary-level token pooling in generating zero-shot embeddings. However, some models might be intrinsically designed to favor certain pooling methods (for example, NT-v2 with its final averaging pooling layer in pre-training). Since the pooling method can be model-specific, task-specific, and sensitive to sequence length, testing more pooling methods (such as max-pooling, attention-based pooling, or hybrid approaches) could facilitate more intriguing analyses.\nResponse: We acknowledge the suggestion to investigate additional pooling approaches. In our revised study, we have expanded our analysis to include maximum pooling alongside the previously evaluated mean pooling and summary-token methods. As described in the \"Sequence Classification Benchmark: Pooling Methods\" section of our Methods, we systematically compared three common output pooling strategies:\n1. Sentence-level summary token method (using [CLS] or [SEP] tokens)\n2. Mean pooling method (averaging all non-padding token embeddings)\n3. Maximum pooling method (selecting maximum values across dimensions)\nOur results, presented in Figure 2 and detailed in the \"Sequence Classification: Pooling Methods Benchmark\" section, demonstrate that mean token embedding consistently delivered statistically superior performance across all models. Using DeLong's test for statistical significance (p < 0.01), we found that mean pooling was the statistically optimal choice for a majority of the binary sequence classification datasets: 37 out of 52 (71.2%) for DNABERT-2, 39 for NT-v2, 32 for HyenaDNA, 34 for Caduceus-Ph, and 35 for GROVER.\n... [Figure 2 description omitted] ...\nThe average increase in AUC when switching from summary token to mean token embedding ranged from 4.3% (DNABERT-2) to 9.7% (HyenaDNA) across all classification tasks. This consistent enhancement underscores the superiority of mean token embedding.\nInterestingly, our results show that maximum pooling rarely outperformed the other methods, being optimal in only a handful of specific datasets, primarily for HyenaDNA and Caduceus-Ph in certain epigenetic modification detection tasks. This suggests that while maximum pooling can occasionally be beneficial for specific task-model combinations, it generally does not match the performance of mean pooling across our comprehensive evaluation.\nWe also observed that the performance differences among the models were reduced when using mean token embedding, with the range of AUC score differences decreasing from 0.063 with summary-token pooling to 0.032 with mean pooling. This finding suggests that mean pooling helps mitigate the architectural variations across the models and provides a more robust approach for standardizing model evaluation.\nBased on these comprehensive results, we recommend mean token embedding as the default pooling method for zero-shot embeddings in DNA sequence classification tasks. We have added relevant results on Page 17-18.\n(6) The study could be more impactful and cutting-edge by extending the test datasets to include longer sequences. This would allow in-depth exploration of the relationship between model performance and sequence length. Transformer-based models are known to be inefficient with long sequences. Test-time training models and MAMBA that scale up to million-bps-long sequences may be compared and discussed in this context.\nResponse: We thank the reviewer for this valuable suggestion. It is indeed essential to explore the capabilities of DNA foundation models on longer sequences, as this reflects real genomic contexts where regulatory elements can act at considerable distances from their targets\nTo address this comprehensively, we have substantially expanded our evaluation framework to include longer sequences through our newly added gene expression prediction benchmark, as detailed in the \"Gene Expression Prediction Datasets\" section of our Methods. This new benchmark enables systematic assessment of how sequence length impacts predictive performance across different model architectures. Systematically applying our long-sequence analysis (e.g., TSS ± 98kbp contexts) across all benchmark tasks was not feasible. This was due to the intrinsic biological focus of many tasks being on more localized genomic features, and the datasets consequently providing sequences of corresponding, often shorter, pre-defined lengths. For example, epigenetic modification datasets (5mC, 6mA, 4mC) are inherently centered on predicting events at single nucleotide positions based on their immediate, fixed-length contexts (typically 41bp. Similarly, many functional element classification tasks, such as promoter identification, rely on experimentally defined sequence windows whose lengths are chosen to capture the specific biological elements under study. To address the impact of sequence length on gene expression prediction, our benchmark utilized the GTEx v8 dataset to create two distinct analytical setups:\n1. Short-Length Analysis: Sequences spanning 6,000 base pairs (TSS ± 3,000 bp) were generated for all genes (e.g., 21,004 genes in whole blood tissue). This allowed evaluation of DNABERT-2, NT-v2, HyenaDNA (6K version), and Caduceus-Ph (6K version).\n2. Long-Length Analysis: For the long sequence dataset, due to computational constraints encountered in Enformer, we randomly selected 1000 genes and excluded genes whose TSS ± 98K exceeds chromosome start or end limit, making the final total of 768 genes. For these 768 genes, we generated sequences spanning up to 196,000 base pairs (TSS ± 98k bp). This enabled us to assess models capable of processing longer inputs, specifically HyenaDNA-450K, Caduceus-Ph (using its 131K bp maximum input), and Enformer (up to 196,608 bp), which served as a key reference given its design for long-range genomic predictions. This setup was crucial for investigating whether extended genomic context improves zero-shot embedding quality for gene expression.\nOur results, detailed in the \"Gene Expression Prediction\" section of our manuscript (Page 19, Tables 6 & 7, Figure 3), indicate that while overall prediction correlations for most genes were modest (average ~0.121-0.137). Importantly, when directly comparing short versus long sequence inputs for the same 768 genes (Table 7), extending the context from 6K to ~196K bp yielded mixed results. HyenaDNA-450K (0.137 correlation) showed a statistically significant improvement over its 6K counterpart (0.122 correlation, p=0.0005). However, Caduceus-Ph with 131K input (0.127 correlation) did not significantly outperform its 6K version (0.124 correlation, p=0.975). Enformer achieved a correlation of 0.129 on this subset. These findings suggest that while architectures like Caduceus-Ph and HyenaDNA efficiently handle longer inputs, the benefit of extended context for gene expression prediction using zero-shot embeddings is still unclear. Further architectural innovations may be necessary to fully leverage very long sequences.\n(7) Although this study focuses on zero-shot embedding to avoid potential bias introduced by fine-tuning, fine-tuning is still widely used and could significantly alter or improve the model’s performance. Future research on how fine-tuning affects performance across different tasks could provide a more complete picture of each model’s capability.\nResponse: We appreciate the reviewer's insightful comment regarding fine-tuning. We fully agree that fine-tuning represents a critical approach for maximizing model performance in specialized genomic tasks.\nAs noted in our Results section on gene expression prediction, we found that zero-shot embeddings alone cannot perfectly capture the relationship between genetic variants and gene expression, only for certain genes they can capture significant sequence-based regulatory information, so still reasonable fine-tuning is needed. This observation applies broadly across the more complex tasks in our benchmark. Furthermore, in our Sequence Classification benchmark, when comparing foundation models to our baseline CNN (which is fully updated during training), we observed that DNA foundation models generally underperformed for multispecies tasks and underperformed in most epigenetic modification tasks. These results highlight specific areas where fine-tuning would likely provide substantial benefits.\nOn the other hand, our zero-shot embedding based evaluation provides several valuable insights that may directly inform fine-tuning strategies. For example, our finding that mean token embedding consistently outperforms summary token embedding across all models suggests that fine-tuning approaches might benefit from focusing on mean token representations.\nSimilarly, our newly added pre-training experiment demonstrates that multi-species pre-training may enhance model generalizability, particularly for cross-species generalization and epigenetic modification detection, providing valuable guidance for selecting appropriate base models prior to fine-tuning.\nWe have expanded our Discussion section (Page 25) to address how our findings inform future fine-tuning research, emphasizing that comprehensive evaluation of fine-tuned models across diverse genomic tasks represents an important next step. Such research would complement our zero-shot benchmark by revealing how architectural differences influence adaptation efficiency and downstream performance.\n(8) While the study is comprehensive in its evaluation of classification tasks, real-life research questions are not limited to this one type of task. Incorporating more types of tasks, such as the more complex regression-based tasks, could contribute to the generalizability of the proposed method and make this work more impactful.\nResponse: Thank you for your helpful suggestion. To address scientific questions beyond classification tasks, we have enhanced our benchmark to include a more complex, regression-based task (i.e., predicting gene expression using genomic sequence) and real world research applications such as variant effect quantification and TAD region recognition.\nOne extension beyond classification tasks is a comprehensive gene expression prediction benchmark using the GTEx v8 dataset. This regression task evaluates the models' ability to predict continuous gene expression values from subject-specific genomic sequences. We focused on whole blood tissue. This benchmark fundamentally differs from classification and allows us to investigate how well foundation models capture regulatory signals influencing gene expression.\nOur findings from this regression task (Results, Pages 18-19; Table 6; Figure 3; Supplementary Table 8) revealed that while the average prediction performance across all genes was modest (mean Pearson correlations typically ranging from 0.121 to 0.137 depending on the model and sequence length), some genes demonstrated notably higher predictability, such as CUTALP, CPNE1 and PEX6. This indicates that zero-shot embeddings from these models can, for certain genes, capture significant sequence-based regulatory information.\nWe also introduced a \"Variant Effect Quantification Benchmark\" (Methods) to assess models' ability to distinguish between clinically significant and benign genetic variants. Using 6000 bp sequences containing either pathogenic or common SNPs, we quantified how dramatically each variant altered the sequence's representation in embedding space. The results, shown in Figure 5, demonstrate that all models successfully differentiated between variant types with statistical significance, with NT-v2 showing the strongest differentiation capability (p=1.25e-36).\nAdditionally, we included a \"TAD Region Recognition Benchmark\" (Methods) to evaluate models' ability to recognize higher-order chromatin structures from genomic sequence alone. By analyzing attention patterns in NT-v2, we found that NT-v2's self-attention mechanism does not recognize TAD boundaries without specific fine-tuning (Figure 6), highlighting an important limitation of current foundation models.\nThese diverse task types provide a more comprehensive evaluation of DNA foundation models' capabilities beyond simple classification, enhancing the generalizability and impact of our benchmark. ... [Figure 5 description omitted] ...\nResponse to Reviewer #2:\nThe manuscript benchmarks three DNA foundation models—DNABERT-2, Nucleotide Transformer version-2 (NT-v2), and HyenaDNA—for genomic sequence classification using zero-shot embeddings. While the authors' approach to comparing the performance of different sequence embeddings is generally appropriate and well-designed, its effectiveness is undermined by the poor quality of datasets chosen for model evaluation. Specifically, the human genome datasets are very limited in size, and some are not based on experimental data. I strongly advise the authors to carefully review the datasets selected for evaluation and select genome-wide, high-quality experimental data. Additionally, to make the manuscript's contribution more impactful, the model comparison should include other models that focus on regulatory signal predictions from genomic sequences.\nResponse: We thank the reviewer for the insightful and constructive feedback on our manuscript. In response to these critical points, we have made every effort to strengthen our study design. The key enhancements, which are detailed in our point-by-point responses below, include:\n1. We have explicitly addressed the concerns regarding dataset quality by contextualizing the results from the 4-methylcytosine (4mC) datasets, which are based on computational predictions rather than direct experimental validation. We now clearly describe the limitation of these benchmarking datasets and present them separately to avoid misinterpretation.\n2. We have incorporated a new, large-scale gene expression prediction task using the experimentally-derived GTEx v8 dataset. This new benchmark assesses model performance on a genome-wide scale (21,004 genes) for predicting gene expression using DNA sequence, directly addressing the concerns about dataset size and quality.\n3. We have expanded our model comparison to include Enformer, a state-of-the-art model specifically designed for predicting regulatory signals from genomic sequence. This provides a vital reference point for evaluating the DNA foundation models.\n4. We have further broadened the scope of our evaluation to include a \"Variant Effect Quantification\" benchmark and a \"TAD Region Recognition\" benchmark, adding depth and practical relevance to our findings.\nWe believe these substantial additions directly address the reviewer's concerns and have significantly strengthened the manuscript, providing a more comprehensive and impactful evaluation of DNA foundation models. We have detailed the specific implementation and results of these new analyses in our responses below.\nMajor points:\n1. Although the datasets cover a diverse range of tasks, the size of each individual human dataset is small and inadequate for a comprehensive comparison of DNA foundation models. A more appropriate approach would be to use whole-genome datasets, as this would provide a more accurate assessment of the models' ability to generalize across diverse genomic sequences.\nResponse: We appreciate the reviewer's concern regarding dataset size limitations and the need for genome-wide evaluations. To address this, we have substantially expanded our benchmark to include a comprehensive gene expression prediction experiment, as detailed in the \"Gene Expression Prediction Datasets\" section of our Methods (Pages 6-7). This new benchmark leverages the GTEx v8 dataset and involves two key analyses:\n1. Standard-Length Genome-Wide Analysis: We generated subject-specific DNA sequences (6,000 bp, centered around the TSS) for all available genes in the whole blood tissue, which involved 21,004 genes across 610 subjects. This represents a genome-scale regression task, assessing the models' ability to predict gene expression from local genomic context across thousands of genes.\n2. Extended-Length Context Analysis: To evaluate performance on longer sequences, we created a dataset spanning up to 196,000 base pairs (TSS ± 98K bp). For this, we initially selected 1,000 genes, which after filtering genes whose TSS ± 98K exceeds chromosome start or end limit resulted in a final set of 768 genes for analysis. This allowed us to specifically assess models like HyenaDNA-450K, Caduceus-Ph, and Enformer on their capacity to utilize extended genomic context.\nThis comprehensive, genome-wide regression benchmark using GTEx data significantly broadens the scope of our evaluation beyond classification tasks. It provides critical insights into how well current DNA foundation models can capture the complex regulatory signals influencing gene expression from sequence data, both from local and more distal genomic regions.\nBesides gene expression prediction tasks, we have added a few additional tasks including “Variant Effect Quantification Benchmark” and “TAD Region Recognition Benchmark\" to breath the width of the benchmarks. The “variant effect quantification benchmark” evaluates models' ability to distinguish between pathogenic and common SNPs by analyzing latent space distances between reference and variant alleles across the genome, providing direct assessment of clinical relevance; and “TAD Region Recognition Benchmark” assesses whether models inherently capture higher-order chromatin organization by analyzing attention patterns across 1,500 topologically associating domains and 1,500 control regions.\nWe believe these benchmark tasks provide a comprehensive evaluation framework that balances specificity and scale, allowing us to thoroughly assess the strengths and limitations of current DNA foundation models. We have added relevant descriptions and results on page 10-18.\n2. The authors should clarify why their performance results differ from those in previous studies. For example, the original NT-v2 paper (The Nucleotide Transformer: Building and Evaluating Robust Foundation Models for Human Genomics) showed that NT-v2 predicts promoters (All, non-TATA, and TATA) better than DNABERT-2. But in this paper, the authors found that DNABERT-2 performed better for some of the promoter classification tasks. The authors should explain what might have caused this difference.\nResponse: Thank you for highlighting this discrepancy between our findings and those reported in the original NT-v2 paper. The different performance results can be attributed to several key methodological differences in our evaluation approach:\nZero-shot embeddings vs. fine-tuning: Our study deliberately evaluates models using zero-shot embeddings with frozen weights, whereas the original NT-v2 paper employed fine-tuning approaches. As we note in our Methods section, most current evaluations are conducted after fine-tuning, which may introduce biases in model performance comparison. For instance, different models may have various levels of overfitting depending on which layers are selected to update during fine-tuning. Our approach mitigates these biases by directly comparing the quality of the embeddings themselves.\nPooling method differences: As demonstrated in our \"Sequence Classification: Pooling Methods Benchmark\" section and Figure 2, the choice of pooling method substantially impacts model performance. We found that mean token embedding consistently outperforms summary token embedding across all models, with performance improvements of 4.3-9.7%. Therefore, for a fair comparison, we used mean token embedding for all the DNA foundation models.\n... [Figure 2 description omitted] ...\nStandardized evaluation pipeline: To ensure fair comparison, we implemented a consistent evaluation framework with 1) identical data preprocessing procedure, 2) random forest classifier with identical hyperparameter optimization for all embeddings from DNA foundation models, 3) statistical significance testing using DeLong's test (p<0.01). The NT-v2 paper's evaluation framework differed in classifier choice, hyperparameter optimization strategy, and significance testing approach.\nIn summary, our benchmark prioritizes a fair, side-by-side comparison of the pre-trained representations themselves, by eliminating variations in fine-tuning and output pooling. While this leads to different results compared to task-specific, fine-tuned results reported in original papers, it provides a clearer understanding of the models' foundational capabilities and generalizability. We have elaborated on the importance of standardized evaluation in our manuscript (e.g., Introduction, page 3; Methods, page 8).\n3. One of the datasets the authors used is 4-methylcytosine (4mC), which is not widely accepted as an epigenetic modification. The paper the authors cited relies on predicted 4mC data rather than experimentally validated annotations, making it a poor choice for model evaluation.\nResponse: Thank you for your critical observation regarding the 4-methylcytosine (4mC) datasets. The biological significance and prevalence of 4mC in higher eukaryotes are indeed less established compared to 5mC or 6mA, and a substantial portion of the available annotation for 4mC sites, particularly in multicellular organisms, is derived from computational predictions rather than direct, high-stringency experimental validation like chemical mapping.\nRecognizing these limitations, we have taken several steps in our manuscript to ensure these results are contextualized appropriately. We have incorporated a discussion in our manuscript Datasets and Results section explicitly stating: “It is important to note that while 4mC is a well-established modification in prokaryotes like E. coli and the Geobacter species, for eukaryotic species (A. thaliana, C. elegans, D. melanogaster), 4mC site annotations are often substantially reliant on computational predictions rather than direct, genome-wide high-stringency chemical mapping due to ongoing research into its prevalence and detection. This distinction in annotation methodology is important for the interpretation of model performance on these 4mC datasets.”\n“As noted in our Methods, the interpretation of results for the eukaryotic 4mC datasets (A. thaliana, C. elegans, D. melanogaster) warrants particular consideration due to the predictive nature of their annotations, and these tasks are thus viewed as exploratory for model capabilities on such data.” Also, in our Results section, the findings for the 4mC datasets are presented separately in Table 5. This allows readers to consider them separately from epigenetic marks.\nWe believe these changes ensure that the performance on 4mC datasets is understood as exploratory and interpreted with the necessary caution due to the nature of the underlying annotations. ... [Tables omitted for brevity] ...\n4. To make the paper more impactful, it would be important to extend the comparison of these three models to include models like Enformer and Sei. These would provide a more complete picture of the models' capabilities.\nResponse: Thank you for your helpful and insightful suggestions. In our revised manuscript, we have expanded our benchmark beyond classification tasks to include a comprehensive gene expression prediction benchmark that allows us to evaluate models on longer sequences and more complex regression tasks and compare DNA foundation models with Enformer.\nA significant addition is our comprehensive gene expression prediction benchmark, which includes a comparison with Enformer. As detailed in our Methods section (\"Gene Expression Prediction Datasets”) and in response to your major comment 1, we developed an extensive regression framework using GTEx v8 data, focusing on whole blood tissue. Our findings from this regression task (Results, page 18-19) revealed that while the average prediction performance across all genes was modest (mean Pearson correlations typically ranging from 0.121 to 0.137), some genes demonstrated notably higher predictability. For example, using 6,000 bp sequences, genes such as CUTALP, CPNE1, and PEX were consistently among the highly predictable genes across various foundation models (Supplementary Table 8), indicating that zero-shot embeddings can capture significant sequence-based regulatory information for certain loci. When comparing long-sequence models on the 768-gene subset (Table 7), extending context length from 6K to ~196K bp yielded mixed results for average correlation improvement; HyenaDNA showed a modest significant increase, while Caduceus-Ph did not.\nFor the highly predictable MYB gene on this long-sequence subset, Caduceus-Ph achieved the highest correlation (0.5822), followed by Enformer (0.5626) and HyenaDNA (0.5400), though these specific values represented only modest gains over 6K sequence predictions for MYB. This suggests that while current architectures can process longer inputs, effectively leveraging very distal information for gene expression prediction via zero-shot embeddings remains a challenge.\nBesides adding Enformer to our new gene expression prediction benchmark, we also implemented a simple Convolutional Neural Network (CNN) baseline consisting of three convolutional layers with pooling and a classification head that was trained directly on one-hot encoded DNA sequences with all parameters updated during training. This comparison (Supplementary Tables 6 & 7) revealed that foundation models significantly outperformed the CNN baseline on several human genome tasks, such as promoter identification (e.g., all five foundation models showed statistically significant improvements for HeLa-S3 cell line promoters) and enhancer classification (DNABERT-2, GROVER, and NT-v2 showed improvements). Conversely, the CNN baseline often outperformed foundation models on multispecies tasks (e.g., B. amyloliquefaciens promoter, mouse TFBS datasets), highlighting current limitations in the cross-species generalization of the foundation models' zero-shot embeddings.\nWhile we considered including Sei in our evaluation, several factors constrained its integration: (1) Sei's architecture and training are highly optimized for variant effect prediction by predicting chromatin profiles, rather than generating general-purpose genomic sequence representations like the foundation models we benchmarked. (2) Its embedding dimensionality, tokenization approach, and output types differ substantially, complicating a direct, fair comparison of zero-shot embeddings. The inclusion of Enformer and our CNN baseline provides valuable reference points for contextualizing foundation model performance, which also enhance our benchmark's comprehensiveness. We have added these results throughout the Results section.\nMinor points:\n5. The HyenaDNA model processes longer sequences than the other two models, so the appropriate way to compare them would be to generate encodings from sequences of equal length. This could be achieved by concatenating the encodings generated by the models with shorter input sizes to match the input sequence length of HyenaDNA.\nResponse: Thank you for your helpful comments regarding the comparison of models with differing maximum input sequence lengths. For the 57 sequence classification datasets, the input sequence lengths, as provided by their original sources, were handled by all five evaluated foundation models. These lengths typically range from 41 bp for epigenetic modification datasets to approximately 2,000 bp for some promoter identification tasks, all falling within the processing capabilities of the models used. Therefore, for these classification tasks, no truncation or manual extension like concatenation of embeddings was necessary, allowing for a direct comparison. We have clarified this on Page 6-7 to avoid confusions.\nIn our newly added gene expression prediction (GTEx) and variant effect quantification benchmarks, we selected a sequence length of 6K bp for the standard analysis. To specifically assess the impact of longer contexts and evaluate models designed for such inputs, we designed separate long-sequence experiments for the GTEx benchmark. These experiments utilized inputs up to approximately 196K bp for HyenaDNA-450K and Enformer, and 131K bp for Caduceus-Ph, as detailed in our Methods (\"Gene Expression Prediction Datasets\") and Results (\"Gene Expression Prediction\"). This approach allowed us to evaluate the benefits of extended genomic context where applicable, rather than compare short-context models to longer-context models by concatenating embeddings, which could introduce their own interpretational challenges. The original HyenaDNA paper (Nguyen et al., NeurIPS 2023) noted that longer sequences do not necessarily yield better performance, with optimal length varying by task. Our experimental design (using common fixed lengths for some benchmarks and long-sequence setups for others) reflects this understanding.\n6. The authors mention that in some cases, they maintained the original training and testing split of datasets, but in other instances, this was not feasible. In such cases, a more appropriate approach for comparing models would be to always randomly split the data, otherwise it could unfairly advantage certain models over others.\nResponse: Thank you for your suggestions on dataset partitioning. We considered applying random splits for all datasets. However, for datasets already partitioned in their original publications (31 out of 57), we chose to preserve these existing splits. Our primary concern with re-splitting these was the potential for data leakage; if an original split was carefully designed to avoid overlaps (e.g., due to sequence similarity or potential overlap), a new random split could inadvertently compromise this separation and lead to inflated performance metrics. Preserving established splits prevents this potential leakage, ensures a fairer model comparison, and maintains consistency with prior evaluations in the field, as original splits often account for inherent biological dependencies. For the remaining datasets without pre-existing partitions, we applied a consistent random 70:30 split. Importantly, our overall results (Tables 1-5) do not suggest any systematic bias favoring models on datasets where their original splits were used, with performance leadership varying across tasks irrespective of dataset origin. The original source of each of the 31 already-partitioned datasets can be found in the Supplementary Information: Dataset Naming and Supplementary Table 2.\nResponse to Reviewer #3:\nIn this manuscript, the authors presented an unbiased evaluation of existing DNA foundation models. The evaluation focuses on the performance of zero-shot embeddings in genomic classification tasks such as promoter/enhancer/TFBS identification, epigenetic modification prediction, as well as application to model organisms. The authors also briefly discussed the effect of pooling methods on model performance. While the topic proves critical in the context of applying LLM to genomics, several key aspects of model evaluation are missing from the current manuscript. These limitations needs to be properly addressed before the paper can be considered for publication.\nResponse: We sincerely thank the reviewer for your nice summary and insightful evaluation of our manuscript. We have undertaken a substantial revisions to strengthen our manuscript. Please find our point-by-point response to your comments below.\nMajor comments:\n1. While the authors presented a comprehensive evaluation of the predictive performance of DNA foundation models, the evaluation on model interpretability is largely missing. Since these models were trained with millions of data points from reference genome, they should be able to capture the important sequence patterns or motifs associated with transcriptional regulation. While the foundation models are generally hard to interpret due to complex architecture of deep neural networks, the authors could focus on the convolutional kernel or attention weights (position weight matrix analysis) and see if they can recapitulate known sequence motifs (for reference, check this review doi: 10.1146/annurev-genom-021623-024727).\nResponse: We appreciate the reviewer's suggestion regarding model interpretability. We fully agree that interpretability is a critical dimension of foundation model evaluation, particularly for genomic models where biological insights are paramount.\nWe have conducted an interpretability experiment examining whether foundation models could capture higher-order genomic structures, specifically topologically associating domains (TADs). As described in our TAD Region Recognition experiment (see Results section), we assessed NT-v2's capability to recognize TAD boundaries through analysis of its attention patterns. We compared attention matrices between 1,500 TAD-centered sequences and 1,500 background sequences, examining whether the model would show differential attention in regions containing TAD boundaries. Figure 6 in our manuscript illustrates this analysis, showing the difference between attention matrices for these two sequence types.\nOur findings revealed that NT-v2's self-attention mechanism does not inherently recognize TAD structures in its zero-shot embeddings. The heatmap in Figure 6 shows predominantly uniform values near zero, with only minimal variation along the diagonal, suggesting the model does not encode these higher-order chromatin structures without specific fine-tuning.\n... [Figure 6 description omitted] ...\nWe acknowledge that examining convolutional kernels or position weight matrices for known transcriptional motifs would provide valuable insights. This represents an important direction for future work, particularly for models like DNABERT-2 and NT-v2 which may capture sequence motifs related to transcription factor binding sites. We have added relevant discussions on Page 26.\n2. While it is straightforward to evaluate model performance on reference genome-based prediction of epigenetic tracks (TF binding, methylation, etc), these tasks are not informative to most biologists who are interested in using DNA foundation models. At the end of the day, these epigenetic tracks can be measured with high-throughput sequencing techniques such as CHIP-seq, ATAC-seq, Bisulfite-seq, etc. Biologists are more likely to use foundation models for predicting variant effect on transcriptional regulation (comparing predictions of reference allele and alternative allele). Therefore, I think the authors should design experiments to compare the performance of DNA foundation models on variant effect prediction. For instance, whether models can accurately predict the expression-change direction of eQTL variants (the authors can check the Enformer paper for datasets and experimental design. doi:10.1038/s41592-021-01252-x)\nResponse: We appreciate the reviewer's insight regarding variant effect prediction as a critical application area for biologists. We fully agree that predicting the functional impact of genetic variants represents one of the most valuable applications of DNA foundation models. To address this important dimension, we have incorporated two new benchmarks specifically focused on variant effects:\nThe first benchmark, detailed in the \"Variant Effect Quantification\" section of our Methods and Results (Page 7, 11 and 21), was specifically designed to assess how DNA foundation models distinguish between pathogenic and common genetic variants. Understanding this capability is essential for determining whether these models can capture meaningful biological signals relevant to disease mechanisms without explicit training on variant classification tasks.\nFor this benchmark, we utilized variant effect quantification data from the Genomics Long-Range Benchmark hosted on HuggingFace, which provides two distinct categories: pathogenic SNPs (associated with disease phenotypes) and common SNPs (frequently observed in the population and typically considered benign). Following methodology introduced in the Nucleotide Transformer paper, we constructed paired sequences for each SNP: a 6000 bp reference sequence based on the reference genome, and a corresponding variant sequence with its center nucleotide replaced by the alternative genotype. This design allows for direct comparison between nearly identical sequences differing by only a single nucleotide.\nFor each sequence pair, we calculated four distance metrics between their zero-shot embeddings to quantify how the models represent genetic variation: L1 distance, L2 distance, cosine similarity, and dot product (cosine similarity without normalization). We then applied the Mann-Whitney U test to compare distributions of distances caused by pathogenic versus common SNPs, testing whether pathogenic SNPs induce larger changes in embedding space.\nAs shown in Figure 5, all models successfully differentiated between pathogenic and common variants with statistical significance. NT-v2 demonstrated superior performance with the most significant differentiation (𝑝 = 1.25×10!\"#), followed by DNABERT-2 (𝑝 = 2.25×10!$#), while HyenaDNA (𝑝 = 5.29×10!%%) and Caduceus-Ph (𝑝 = 9.29×10!&) showed progressively less pronounced, though still significant, differentiation. These patterns remained consistent across all four metrics, with NT-v2 consistently producing the lowest or second lowest p-values. The consistency across multiple distance metrics provides strong evidence that these models inherently capture biologically meaningful representations of variant effects.\n... [Figure 5 description omitted] ...\nOur second benchmark, described in the \"Gene Expression Prediction\" section (Methods, pages 6-7; Results, pages 18-19). This evaluates the models' ability to predict gene expression levels from subject-specific DNA sequences incorporating individual genetic variants. This benchmark directly tests the models’ capacity to interpret the functional consequences of genetic variants on gene expression. Success in this task is a critical prerequisite for downstream applications, such as using DNA foundation model based predicted expression to conduct causal inference and identify significantly associated genes for complex diseases by testing association between predicted gene expression and trait of interest. Using the GTEx v8 dataset, we focused primarily on whole blood tissue (610 subjects, 21,004 genes). For the standard analysis, we generated 6,000 bp sequences (TSS ± 3,000 bp) for each gene, incorporating subject-specific variants. Zero-shot embeddings were generated, and Random Forest regression models (which outperformed XGBoost in our tests, as noted on page 18) were trained to predict gene expression. While the average Pearson correlation between predicted and actual expression values across all genes was modest (typically ranging from 0.121 to 0.137, Table 6), a consistent subset of genes, including CUTALP, CPNE1 and PEX6, showed significantly higher predictability.\nWe further extended this gene expression analysis to much longer genomic sequences (up to ~196K bp, using TSS ± 98K bp) for a filtered subset of 768 genes, evaluating HyenaDNA-450K, Caduceus-Ph (131K input), and Enformer. Interestingly, while HyenaDNA showed a statistically significant, albeit modest, average improvement with longer context, Caduceus-Ph did not (Table 7). These findings suggest that while zero-shot embeddings capture some regulatory signals, fully leveraging distal genomic context for precise gene expression prediction without fine-tuning remains a challenge for current models. The consistent performance on specific genes, however, indicates that these models do capture certain strong, sequence-encoded regulatory features.\nWe appreciate the reviewer's valuable suggestion to specifically evaluate eQTL variant effect prediction. We sincerely think predicting the expression-change direction of eQTL variants and identifying causal eQTL are important to the community and we plan to thoroughly investigate this in a separate project. We sincerely believe that DNA foundation models hold great potential to complement traditional statistical methods for identifying causal eQTLs. We have added relevant discussion on Page 26: The DNA foundation models hold great potential in identifying causal variants for omics and complex diseases. The sequence-based representations learned by foundation models may capture functional genomic context that complements statistical signals from fine-mapping and colocalization, potentially enhancing our ability to distinguish causal regulatory variants from those in linkage disequilibrium. We leave such exciting direction to future research.\n3. Another question the audience will be interested in is the comparison between DNA foundation models in general versus other widely-used sequence-based models such as Enformer, Sei, etc. Do they generally perform better or worse in terms of the tasks evaluated in this manuscript? While the deep learning in genomics community would like to see a direct comparison between these models in terms of the specific tasks in Result section, I understand adding these analyses can be time-consuming. I encourage the authors to (at minimal) address the question in the Discussion section, discussing the pros and cons of using foundation models versus traditional sequence-based models.\nResponse: Thank you for your helpful and insightful suggestions. Following your suggestions, we have expanded our benchmark to include a comparison with Enformer in our gene expression prediction benchmark (Figure 3 and Table 7). For these long-sequence gene expression prediction tasks, we found Enformer showed comparable performance to Caduceus-Ph and HyenaDNA, with no model clearly outperforming the others across all genes. Specifically, Enformer achieved the second-highest correlation (0.5626) for the MYB gene, behind Caduceus-Ph (0.5822).\nTo systematically evaluate the tradeoffs between foundation models and specialized architectures, we implemented a simple Convolutional Neural Network (CNN) baseline with three convolutional layers, pooling operations, and a classification head trained directly on one-hot encoded DNA sequences. This task-specific model revealed important patterns (Supplementary Tables 6 & 7): Foundation models significantly outperformed the CNN baseline on several human genome tasks, such as promoter identification (e.g., all five foundation models showed statistically significant improvements for HeLa-S3 cell line promoters) and enhancer classification (DNABERT-2, GROVER, and NT-v2 showed improvements). Conversely, the CNN baseline often outperformed foundation models on multispecies tasks (e.g., B. amyloliquefaciens promoter, mouse TFBS datasets), highlighting current limitations in the cross-species generalization of the foundation models' zero-shot embeddings.\nA key trade-off exists between the versatility of foundation models and the task-specific optimization of specialized models. Foundation models excel in zero-shot settings and show strong potential for transfer learning across tasks, while specialized models may provide superior performance on their target applications. Moving forward, hybrid approaches that combine foundation model pre-training with specialized architecture components could potentially offer the best of both paradigms. We have also added relevant discussions on Page 24-26: Our comparison between DNA foundation models and specialized architecture reveals a fundamental trade-off between versatility and task-specific optimization in genomic sequence analysis. While foundation models demonstrate remarkable adaptability across diverse tasks through extensive pre-training, specialized models incorporate inductive biases that can yield superior performance within their targeted domains. This complementarity suggests promising directions for hybrid architectures that combine pre-trained genomic representations with task-specific components to simultaneously achieve broad applicability and domain-specific excellence.\nMinor comments:\n4. Instead of random splitting of training and testing sets, I suggest the authors adopt whole chromosome holdout (e.g. chromosome 8 and 9 as testing set) to avoid overlap of regulatory region between training and testing set.\nResponse: Thank you for your helpful suggestions. We fully agree that chromosome-based holdout evaluation represents the gold standard for genomic models, as it provides a more stringent assessment of generalization by preventing regulatory element overlap between training and testing sets. Unfortunately, we faced practical constraints with this approach. Many datasets in our study were obtained in post-processed format without chromosome information. Therefore, we maintained train-test splits from original works where available and used random 70:30 splitting elsewhere as described in our Methods section. This approach prevents potential data leakage, maintains consistency with prior evaluations, and ensures fair model comparison. We have cautioned this limitation on Page 8, “Sequence Classification Benchmark” under Methods section.\nFor our newly added gene expression prediction task, we eliminated the concern of regulatory region overlapping by defining a separate prediction question for each gene, which effectively avoids overlapping regulatory regions between training and testing samples since all subjects roughly capture the same region.\n5. Figure 2, please add asterisks on top of each boxplot denoting the significance level of comparing mean pooling versus CLS pooling.\nResponse: We have now included three pooling methods (mean, summary token, and maximum pooling) in our comparison and performed systematic pairwise testing using DeLong's test with a p-value threshold of 0.01 to determine significant differences in performance across all 52 binary classification datasets. Briefly, mean token pooling demonstrated statistically superior performance for the majority of datasets across all models (71.2% for DNABERT-2, 75.0% for NT-v2, 61.5% for HyenaDNA, 65.4% for Caduceus-Ph, and 67.3% for GROVER). The full pairwise statistical comparisons are presented in Supplementary Table 5. We have also updated the manuscript on Page 12-13 to describe these statistical significance results.\nReviewer #3 (Remarks on code availability):\nThe README file is clear. I think users should be able to reproduce the results from the GitHub repo.\nResponse: Thank you for your supportive comments.\nMinor comment:\n6. Some file paths are still based on authors' institutional computing cluster (such as \"/rsrch4/home/biostatistics/hfeng3/review_datasets/results_final/ntv2/\" in \"classify_ntv2.py\"). Please revise those paths in the context of GitHub directory.\nResponse: Thank you for your careful code review. We have thoroughly addressed the path dependencies in our codebase. Additionally, we have significantly improved the quality of our codes, allowing researchers to easily validate our findings and extend our benchmark to new models.\n\nA Benchmarking DNA Foundation Models for Genomic Sequence Classification\nResponse to Reviewer #1:\nThe authors addressed my comments very well. I don't have any more comments and commend the authors for the good work!\nResponse: Thank you for your support and recognition of our work. We really appreciate it.\nReviewer #1 (Remarks on code availability):\nThe README should be clear for readers to reproduce results.\nResponse: Thank you for checking our codes. We appreciate it.\nResponse to Reviewer #2:\nI have reviewed the revised version of the manuscript “Benchmarking DNA Foundation Models for Genomic Sequence Classification” as well as the authors’ point-by-point response. While I appreciate that the authors have implemented some important improvements, I find that most of the major concerns outlined in my initial review remain unaddressed.\nResponse: Thank you for your critical comments, which have pushed us to think harder and further improved our work. We have made every effort to address your remaining comments. Specifically, we have added further description for each dataset we used (whether experimental evidence or computationally generated, whether genome-wide) and added several new tasks regarding variant effect quantification. Furthermore, we have added Sei (Chen et al. Nature Genetics, 2022, PMID: 35817977) to our benchmark tasks in relevant tasks and also added recent AlphaGenome (Avsec et al., Preprint from DeepMind; https://www.biorxiv.org/content/10.1101/2025.06.25.661532v2). We have also benchmarked the tracks and embeddings for different models. In addition, we have added results for random data splitting, which remains similar and does not change our main conclusions. For your convenience, we first state your comments and then provide our responses. Changes in the manuscript are highlighted in blue. We have provided our detailed responses below and sincerely hope you will find our revision satisfactory.\n1. The addition of the GTEx v8 gene expression prediction task, which measures model prediction accuracy of gene expression across individuals, is a well-considered choice for the dataset, however, it improves the evaluation of only one task. For a benchmarking paper, it is critical that the appropriate evaluation of DNA foundation models be conducted using genome-wide, high-quality experimental data for the majority, if not all, of the tasks.\nResponse: Thank you for this insightful comment regarding benchmarking tasks. To address this critical point, we have made several revisions to the manuscript.\nFirst, to enhance the clarity of our dataset origins and scope, we have thoroughly revised the Methods: Benchmarking Datasets section (pages 20-22). For all five major dataset collections, we now explicitly detail their origins, indicating whether the positive samples are derived from direct experimental evidence or gold standard sources (e.g., ENCODE ChIP-seq, GTEx) or are computationally generated in Supplementary Text 2 (also attached below). We also clarify the scope of each dataset, noting whether it represents specific genomic regions or a genome-wide collection. We believe this important information will help readers better understand the contexts and implications of our results. ... [Supplementary Text 2 omitted] ...\nSecond, following your suggestions and further enhancing the contribution of our work, we introduced an entirely new benchmark focused on Variant Effect Quantification, which directly evaluates model performance on genome-wide experimental data. This new section is detailed in Methods: Variant Effect Quantification Datasets (page 23). This benchmark includes two genome-wide tasks: our expanded pathogenic vs. common SNP task using long-range sequences (196,608 bp) and the incorporation of a new high-quality dataset of fine-mapped causal quantitative trait loci (QTLs), derived from GTEx v8 whole-blood tissue. This QTL task requires models to distinguish putative causal variants from carefully matched non-causal negatives, providing a direct test of their ability to capture subtle regulatory signals. To ensure this new benchmark is rigorous, we have also included and compared several models considered state-of-the-art for genome-wide prediction tasks, namely Enformer (Avsec et al. Nature Methods, 2021, PMID: 34608324), Sei (Chen et al. Nature Genetics, 2022, PMID: 35817977), and AlphaGenome (Avsec et al., Preprint from DeepMind; https://www.biorxiv.org/content/10.1101/2025.06.25.661532v2). We also implemented a robust data partitioning strategy for our benchmarks. Since complete chromosomal information was available for these datasets, we were able to move beyond a simple random split and employ a nested cross-validation scheme with strict chromosome-based holdouts. This approach ensures that our performance estimates are robust and stable by averaging results across multiple, independent chromosome-based test sets.\nOur methodology for applying these models is detailed in Methods: Variant Effect Quantification Benchmark (pages 27-28), and the results are presented in the new Results: Variant Effect Quantification section (pages 12-13).\nThe results from this new benchmark reveal a key distinction: while specialized models like AlphaGenome and Enformer excel at the fine-grained QTL prediction tasks, the foundation models proved effective at distinguishing pathogenic from common SNPs. Notably, NT-v2 emerged as the top performer on the pathogenic SNP task, suggesting its diverse pre-training captures fundamental disease-relevant signals, while Caduceus-Ph served as the most consistent foundation model across these new variant effect benchmarks. The results are also briefly described in our response to your Comment 2.\nThe results reflect their distinct training objectives. Specialized foundation models like AlphaGenome and Enformer were explicitly designed to predict functional genomic annotations and regulatory tracks such as transcription factors and chromatin accessibility, which capture the regulatory mechanisms that underlie quantitative trait variation. On the other hand, general foundation models such as NT-v2 were pre-trained on phased DNA sequence to understand general genomic “grammar” itself and learn fundamental patterns of sequence constraint and evolutionary pressure that correlate with pathogenicity. We have added relevant discussion on pages 17-18.\n2. I appreciate that the authors included Enformer and a simple CNN model as references for model performance comparison. However, Enformer was only added to the gene expression task, whereas it should also be included in other human genome tasks.\nResponse: We thank the reviewer for this excellent suggestion. We agree that expanding the evaluation of state-of-the-art models like Enformer across multiple human genome tasks is crucial for a comprehensive benchmark. Following your suggestions, we have now incorporated Enformer, into two new, large-scale variant effect quantification benchmarks, which we believe are highly appropriate for evaluating its capabilities and provide performance reference for DNA foundation models.\nFirst, we have carefully and very seriously considered applying Enformer to the human genome datasets among the 57 sequence classification datasets. However, we determined this was not ideal due to a fundamental architectural constraint: Enformer requires a fixed input length of 196,608 bp. The sequences in our classification benchmark are significantly shorter, mostly ranging around several hundred nucleotides. The standard approach for handling such a discrepancy would be to pad each sequence with hundreds of thousands of 'N' characters (https://github.com/google-deepmind/deepmind-research/issues/287). This practice is not only computationally extensive but more importantly would also unfairly penalize Enformer's performance, particularly given its use of large 128 bp windows for processing DNA sequence.\nTherefore, we designed these new variant effect benchmarks with different sequence lengths, specifically to create an appropriate and rigorous setting where the full capabilities of Enformer and other long-range models could be leveraged. Our design, which includes both short and long sequence versions of the variant datasets, was specifically chosen to probe the impact of genomic context. We believe this approach effectively addresses the concern of expanding Enformer's evaluation.\nSpecifically, we now included Enformer in our expanded pathogenic vs. common SNP benchmark, which now uses long-range input sequences of 196,608 bp. Enformer demonstrated strong performance (AUROC ≈ 0.69 using its hidden states), confirming its ability to leverage large genomic context to identify disease-relevant variants. However, it was notably outperformed by the foundation model NT-v2 (AUROC = 0.73), a surprising and significant finding that underscores the power of pre-training on diverse species for this task.\nMore significantly, we benchmarked Enformer on our newly incorporated QTL task. This dataset consists of high-confidence, fine-mapped putative causal variants from GTEx v8 whole blood tissue. Enformer showcased its strength as an architecture explicitly trained to predict genomic tracks. It delivered exceptional performance on distinguishing putative causal expression QTLs (eQTLs) from non-causal variants, achieving an AUROC (AUC) of 0.77, and also performed well on splicing QTLs (sQTLs). These results position it as a top-tier model for predicting the impact of variants on specific regulatory functions, though the newly introduced AlphaGenome consistently achieved the highest performance across all QTL categories. Notably, while our uniform evaluation framework produces different absolute result values than reported in the original AlphaGenome manuscript as we handle tracks differently to ensure unbiased comparison across all models, the conclusions are highly consistent with AlphaGenome manuscript. ... [Tables 9 and 10 omitted] ...\nIn summary, we want to keep our benchmarking fair and unbiased. We have added relevant discussion on pages 17-18.\n3. Additionally, the authors state that Sei's architecture and training are \"highly optimized for variant effect prediction,\" which is not correct. If there are technical difficulties in applying the Sei model as a reference, the authors should instead use other appropriate models available for each specific task.\nResponse: We sincerely thank the reviewer for this important correction and apologize for the inaccurate characterization of the Sei model. Sei takes input sequences of 4096 nucleotides and outputs predictions for 21,907 chromatin profiles, which serve as comprehensive feature representations for each sequence. These predictions can capture diverse regulatory activities in genomic sequence.\nFurthermore, following your guidance, we have added open-source Sei model and API-based AlphaGenome into our variant effect quantification benchmark. To provide a thorough evaluation of Sei, we extracted two distinct high-dimensional representations for each variant sequence: the final output tracks (21,907 features) and the final hidden state (15,360 features) after the spline transformation layer. For AlphaGenome, which we included in the QTL tasks, we utilized its 131k bp architecture and derived an effect vector by averaging its output tracks over the central 2,048 bp. For both benchmarks, these representations were treated as zero-shot embeddings, applying the same methodology used for the other foundation models to ensure fair comparisons.\nThe inclusion of these powerful models yielded several important insights, now detailed in the Results: Variant Effect Quantification section (pages 12-13, Tables 9 and 10).\nIn the pathogenic SNP task, Sei performed robustly with an AUC of approximately 0.67 using its output tracks, placing it ahead of most general-purpose foundation models. In the more fine-grained causal QTL tasks, AlphaGenome was the clear standout performer, achieving the highest AUC across all four QTL types, including 0.80 for sQTLs and 0.86 for ipaQTLs. Sei also demonstrated strong, competitive performance, often securing the second-best results. We have added relevant discussion on pages 17-18.\nA notable finding with Sei was that its internal hidden states were often as predictive as, and in some cases (such as for sQTLs) even slightly more predictive than, its final processed output tracks, suggesting rich information is encoded before the final projection.\nWe clarify that the application of these specialized models was restricted to the variant effect benchmark due to fundamental architectural and practical constraints. Models such as Sei and Enformer are designed for very long, fixed-length inputs (>100,000 bp), making them unsuitable for our much shorter classification sequences without resorting to massive, methodologically unsound padding that would be computationally expensive and unfairly penalize their performance. Importantly, Sei does not recognize the padding tokens, making it even more challenging to directly apply to short sequence. Similarly, the use of the AlphaGenome model was also limited by its long-sequence design and further constrained by practical API quota issues. We therefore designed the new variant effect benchmarks specifically to provide a fair and rigorous evaluation ground for these long-range architectures.\nFinally, we wish to emphasize the core goal of this study. For the sequence classification benchmark, our primary aim is to accurately benchmark the zero-shot performance among DNA foundation models to understand their intrinsic differences, rather than to exhaustively investigate whether they can outperform every existing task-specific model. To this end, we provide a simple CNN to establish a reasonable performance baseline, a strategy consistent with recent large-scale benchmarking efforts in related fields that also compare foundation models against simple baselines such as simple linear models (Ahlmann-Eltze et al., Nature Methods, 2025, PMID: 40759747). We hope these benchmarking results will provide useful information to the community to decide whether and which DNA foundation models to be used for downstream tasks.\n4. My final concern is that the authors chose to preserve the original training and testing splits of certain datasets, but not others. The authors explain that this choice was made to avoid potential data leakage. While I understand this reasoning, I believe it is the authors’ responsibility, particularly in the context of a benchmarking study, to design datasets in a way that prevents data leakage and allows for fair model comparison. It cannot be overstated that finding new comprehensive datasets, as well as the careful design of previously used datasets, is what makes a benchmarking study valuable and impactful.\nResponse: We are grateful to the reviewer for raising this critical point concerning data partitioning. We completely agree that the careful and consistent handling of data splits is fundamental to the integrity of a benchmarking study, ensuring fair comparisons and preventing data leakage. It is indeed our responsibility to guarantee this, and we have taken significant steps in the revised manuscript to address this concern thoroughly.\nIn response, we undertook a comprehensive re-evaluation of our data partitioning strategy across all of our benchmarks. For the 57 sequence classification tasks, we have now re-partitioned all datasets that previously relied on prior splits (which varies differently across different datasets such as 90:10, 95:5, and 70:30 train-test split), particularly those adopted from the DNABERT-2 and NT-v2 benchmarks, by combining the original train, test, and/or validation sets and applying a uniform, random 70:30 train-test split.\nThe only exceptions to our re-splitting are the datasets from the \"Genomic Benchmarks\" datasets collected by us, and the \"Enhancer” and “Enhancer Strength\" datasets adopted from NT-v2 benchmarks. We retained their original splits because: (1) The Genomic Benchmarks collection is a publicly established resource specifically designed for standardized benchmarking, and has not been processed by any target-oriented models. (2) The enhancer datasets training set is a mixture of experimentally verified sequences and synthetic data from a generative model, which is intractable to separate. The test set, however, consists solely of experimental data to ensure an unbiased evaluation of real-world performance (Geng et al., Biophysical Chemistry, 2022, PMID: 35605495). Therefore, any re-split would risk contaminating the test set with synthetic data, compromising the integrity of the benchmark. Crucially, we double checked that these datasets were published independently and were processed by any foundation model (e.g. NT-v2 did not process “Enhancers dataset”) we evaluate, thereby posing no risk of inherent model bias.\nUpon re-running our benchmark with these newly standardized data splits, while some AUCs changed in decimal places, we found that our main conclusions remain largely robust, especially for the comparisons among DNA foundation models. One minor change in conclusion is the relative comparisons between DNA foundation models and baseline CNN as detailed in Results: Sequence Classification: Human & Multispecies Epigenetic Modification section (pages 8-9). On the newly partitioned yeast epigenetic mark prediction data, multiple DNA foundation models now demonstrate a performance advantage over the baseline CNN. We have updated the entire manuscript to reflect these new results.\nResponse to Reviewer #3:\nAll my comments were properly addressed by the author. I have no further comments. I commend the authors for their good work.\nReviewer #3 (Remarks on code availability):\nThe GitHub repo is well documented. Results should be reproducible.\nResponse: Thank you so much for supporting our work and careful review. We truly appreciate it."}]

Decision letter