
==== Front
HGG Adv
HGG Adv
Human Genetics and Genomics Advances
2666-2477
Elsevier

S2666-2477(24)00087-3
10.1016/j.xhgg.2024.100347
100347
Article
Enhancing personalized gene expression prediction from DNA sequences using genomic foundation models
Ramprasad Pratik 1
Pai Nidhi 1
Pan Wei panxx014@umn.edu
12∗
1 Division of Biostatistics and Health Data Science, University of Minnesota, Minneapolis, Minneapolis, MN, USA
∗ Corresponding author panxx014@umn.edu
2 Lead contact

27 8 2024
10 10 2024
27 8 2024
5 4 1003477 3 2024
23 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Artificial intelligence (AI)/deep learning (DL) models that predict molecular phenotypes like gene expression directly from DNA sequences have recently emerged. While these models have proven effective at capturing the variation across genes, their ability to explain inter-individual differences has been limited. We hypothesize that the performance gap can be narrowed through the use of pre-trained embeddings from the Nucleotide Transformer, a large foundation model trained on 3,000+ genomes. We train a transformer model using the pre-trained embeddings and compare its predictive performance to Enformer, the current state-of-the-art model, using genotype and expression data from 290 individuals. Our model significantly outperforms Enformer in terms of correlation across individuals, and narrows the performance gap with an elastic net regression approach that uses just the genetic variants as predictors. Although simple regression models have their advantages in personalized prediction tasks, DL approaches based on foundation models pre-trained on diverse genomes have unique strengths in flexibility and interpretability. With further methodological and computational improvements with more training data, these models may eventually predict molecular phenotypes from DNA sequences with an accuracy surpassing that of regression-based approaches. Our work demonstrates the potential for large pre-trained AI/DL models to advance functional genomics.

AI models can predict gene expression from DNA sequences but struggle to account for variations across individuals. We show how AI foundation models can be used to improve performance over existing approaches. Our results highlight the potential of foundation models to enhance gene expression prediction and advance functional genomics.

Keywords

AI
DL
Enformer
Nucleotide Transformer
elastic net regression
SNP
gene expression
transformers
foundation models
==== Body
pmcIntroduction

Many artificial intelligence (AI)/deep learning (DL) models have recently been introduced for molecular genomic modeling tasks such as predicting gene expression1,2 and epigenetic features.3,4 In particular, several of these state-of-the-art AI/DL models are trained using DNA sequences as input1,2,3,5 rather than using just the observed genetic variants (usually single-nucleotide polymorphisms [SNPs]) as predictors, as in regression-based approaches.6 These in silico models shed light on biological mechanisms based on the patterns they learn from sequences,7 but more importantly, they can also predict outcomes for arbitrary sequences, including those with rare variants. In this way, they can supplement in vitro or in vivo experiments, which may be constrained by time, equipment, and labor. These sequence-to-expression models “learn” by recognizing short patterns of nucleotides, known as motifs, that are repeated throughout the genome, and outcomes are predicted based on combinations of such motifs.8

Enformer is one such DL model, which was trained in a multitask framework to predict several functional genomic outputs directly from DNA sequences.1 Recently, several groups have benchmarked the performance of Enformer and other DL models, including ExPecto,2 Xpresso,5 and Basenji2,9 in predicting gene expression directly from a sequence.8,10,11 Using paired whole-genome sequencing and RNA sequencing data from the Geuvadis consortium, Huang et al.10 computed cross-gene and cross-individual correlations between model predictions and true RNA abundance values. The average cross-gene correlation was calculated by fixing an individual, calculating rank correlation between the predicted and the observed for all genes in the study, and then averaging across individuals, whereas the average cross-individual correlation was calculated by fixing a gene, calculating rank correlation between all individuals, and averaging over genes. In terms of cross-gene correlation, Enformer outperformed all other DL models included in the study, but its cross-individual correlation, as well as the other models’, was near 0. Sasse et al.8 reached similar conclusions in that these models can be effective at predicting population average gene expression but struggle to predict individual expression for a fixed gene. Both studies also found that these models tend to mispredict the sign of the effect of an SNP. In fact, the models tended to agree with each other more often on the magnitude of the effect than the direction, indicating that they recognize the presence of important variants but struggle to determine the direction of their effect.10 Furthermore, these models are outperformed by traditional statistical models, such as regularized linear regression as implemented in the well-known PrediXcan,6 in predicting individualized impact of genetic variations. Such models use only selected SNP allele counts as inputs, ignoring the overwhelming majority of a DNA sequence. Therefore, although simpler, such statistical models based on allele counts may be less effective due to the limited information used to train them.

One explanation for the DL models’ poor performances in predicting individualized gene expression is that they are trained on one reference genome rather than diverse genomes from different individuals, and the mechanisms that differentiate expression between genes may differ from those that differentiate expression between individuals.8 Recently, Dalla-Torre et al.3 built the Nucleotide Transformer (NT), a family of transformer-based foundation models pre-trained on DNA sequences from diverse genomes. The largest of these models has 2.5 billion parameters and was trained on genomes from 850 diverse phyla, including 11 model organisms. To solve downstream tasks such as gene expression or epigenetic mark prediction, the model can either be fine-tuned, or specific layers can be probed as embeddings to be fed into a separate, smaller downstream model. The authors show that using the NT outperforms existing methods in predicting regulatory elements, splice sites, and histone modifications. Furthermore, using a zero-shot learning-based approach, the authors show that the model embeddings capture which variants are functionally important. By computing the distance in embedding space between a sequence with a variant and the reference sequence, the model embeddings can classify ClinVar “likely pathogenic” variants versus variants not in the ClinVar dataset with an AUC score of 0.8. Thus, the NT seems promising in identifying important sources of genetic variation to connect DNA sequences with molecular phenotype outcomes. While the authors have rigorously benchmarked the NT’s performance in predicting regulatory elements, splice sites, and histone modifications, its application to predicting personalized gene expression levels remains underexplored. Moreover, previous studies have generally not leveraged paired DNA sequence-gene expression data for training DL models, as is common with traditional statistical models like PrediXcan. This gap prompts an investigation into whether integrating such data into DL frameworks could enhance their predictive accuracy and robustness for personalized gene expression prediction.

Here, we use paired DNA sequence-gene expression data for personalized gene expression prediction. In particular, we leverage the pre-trained NT to obtain embeddings of DNA sequences as input to a specialized downstream transformer model to predict gene expression. We find that using embeddings from the foundation model significantly outperforms the predictions from Enformer and is competitive with the elastic net regression approach across several genes.

Material and methods

Data and modeling pipeline

Figure 1 shows the general framework we used to train our model. Since our goal is to model the effect of personal genome variation on gene expression, we use normalized RNA sequencing data from the Genotype-Tissue Expression (GTEx) project,12 measured across five different tissues (muscle skeletal, whole blood, skin sun exposed lower leg, artery tibial, and adipose subcutaneous) and paired with Whole-Genome Sequencing (WGS) data from 290 individuals. The tissues were selected based on the maximal cohort of individuals for whom comprehensive data was available across all chosen tissues. We construct training and test sets for the donors according to a random 80/20 split, and perform another random 80/20 split on the training set to derive the validation set. To ensure consistency, these partitioning schemes were implemented uniformly across all genes, tissues, and models.Figure 1 Modeling framework

Input DNA sequences for each donor-gene combination are constructed by inserting the corresponding variants into the human reference genome. The sequences are then tokenized and input into the pre-trained model to generate embeddings. Embeddings from the last layer of the model are fed into the downstream transformer model, which is trained to extract relevant features from the embeddings and predict gene expression.

Gene expression data are available across all selected tissues for 17,289 gene IDs in total, but for the purpose of this analysis, we only consider a subset of 600 gene IDs, based on the following criteria. First, from the subset of locally heritable genes (based on the models from http://gusevlab.org/projects/fusion/),13 we select the 200 shortest and 200 longest genes, with length measured as the distance between the transcription start site (TSS) and end site. Additionally, for the purpose of out-of-sample testing, we randomly select 200 genes from those that lie outside the set of locally heritable genes.

Our choice of genes is based on the following reasoning. First, as discussed below, we use an input sequence of length 4 kb for our model, and we reason that, for locally heritable genes, shorter genes are more likely to have their expression quantitative trait loci located within a closer distance from the TSS. Thus, we selected 200 short genes from the set of locally heritable genes. Second, we include the 200 longest locally heritable genes in order to measure the extent to which the shorter input sequence length of our model affects its performance on longer genes relative to Enformer, which uses much longer sequence inputs. Finally, since the set of locally heritable genes was derived using linear models, we sought to determine whether our model could capture locally heritable features that might be undetectable to linear models trained on the variant allele counts alone, such as elastic net regression.

Personal input sequences for each individual are constructed by inserting their single nucleotide variants (SNVs) into the Homo Sapiens Assembly 38 reference sequence14 around each gene’s TSS. Two sequences are constructed for each individual, corresponding to each haplotype. For computational reasons and the fact that the NT model was trained using sequences of length 6 kb, we construct our input sequences using a window length of 2 kb on either side of the TSS. The raw sequences are then tokenized and input into the (frozen) pre-trained foundation model in order to obtain the hidden-state embeddings. These pre-trained embeddings are then used as inputs for our downstream transformer blocks.

A key aspect of our approach is that, instead of training separate models for each gene-tissue combination, as is typically done in variant-to-expression prediction models proposed in the literature, we sought to train a single model by stacking the embeddings from a large set of genes. This was based on the hypothesis that the embeddings from the NT model might capture regulatory information present within the sequences that has the potential to explain variation in gene expression across individuals and across genes.

NT

The NT is an encoder-only transformer that is trained in a masked language modeling framework. The inputs for training the model are 6-kb-long DNA sequence chunks tokenized into 6-mers. Some of these tokens are randomly masked, and the model is trained to predict the most likely token at masked positions. The architecture of the model consists of an embedding layer, followed by a stack of transformer layers, each with a multihead attention block, followed by a language model head. The language model head outputs a probability distribution over tokens at each masked position in the sequence, and the model is trained using the cross-entropy loss between the predicted probabilities and ground truth token label.

While there are several different versions of the NT, each of varying size and pre-trained on one of three different datasets, in this work, we exclusively use the NTv2-Multispecies 100M model. This version was trained on 850 genomes selected across diverse phyla (including the human genome, as part of the “mammalian vertebrates” group). The model itself has 22 transformer layers with 16 attention heads and, in total, about 100 million parameters. Our choice of pre-trained model was determined by the performance results reported in Dalla-Torre et al.3 and the size of the embedding dimension (512). While larger models, such as the 250m parameter version, encode richer information and can potentially further enhance predictive performance, they are computationally more demanding and use larger embedding dimensions (up to 2,560).

In order to determine the optimal layer of the foundation model to probe for embeddings, we ran a study (Figure 2) where we trained a family of models using embeddings from each of layers 12–22 on a set of 100 locally heritable genes and evaluated their predictions on the validation set of donors. Based on this study, we determined that the most appropriate embeddings to use for our task of gene expression prediction are those from the last layer.Figure 2 Correlations between predicted and true gene expression levels for downstream models trained on different layers of the foundation model

The model trained using embeddings from the last layer (22) significantly outperformed those trained on the previous layers.

Downstream model architecture and training

Our downstream transformer model architecture is depicted in Figure 3. It consists of the following components: a sequence of 2 cross-attention blocks, 2 self-attention blocks, and a multilayer perceptron (MLP) regression head. The attention blocks use the same hidden dimension as the pre-trained embeddings; i.e., 512 with 8 attention heads. In total, the model contains about 12.9 million parameters.Figure 3 Downstream model architecture

A single input consists of two separate sets of embeddings, one for each haplotype. The haplotype embeddings are processed through a series of shared cross-attention blocks, where each set attends to the other in order to integrate complementary information. They are then passed through a series of shared self-attention blocks to enhance the representation of intra-sequence dependencies. Finally, the two sets of embeddings are fused via mean pooling and fed into an MLP regression head, which generates the final prediction.

Our architecture is loosely based on the multimodal transformer model proposed by Tsai et al.,15 with the key difference being that our transformer blocks are shared across the two sets of input haplotype embeddings, given that they both belong to the same modality. The cross-attention blocks facilitate the integration and alignment of complementary information from the two sets of embeddings, enabling the model to effectively capture interactions between them. Meanwhile, the self-attention blocks enable the model to capture dependencies within each set of haplotype embeddings, enhancing the representation of intra-sequence relationships.

The architectural configuration of our model was determined based on a study (Figure 4) where we constructed a family of models by varying the depth of the two sets of transformer blocks (the cross-encoder and the encoder blocks), trained these models on a set of 100 locally heritable genes, and measured the resulting validation loss after 100 epochs to determine the best-performing architecture. In general, we saw marginal improvements in performance as more attention layers were added to both transformer blocks, and the best-performing model was the largest one considered, where we used two attention layers in both sets of transformer blocks. This was the final configuration selected for this work.Figure 4 Validation loss vs. epoch for a family of models constructed by varying the number of attention layers (depth) in the two transformer blocks of our model

These two parameters are referred to here as the cross-encoder depth (ced) and the encoder depth (ed). The worst performance was the case where both ced and ed are 0, where we use no transformer blocks whatsoever, and the model only contains the MLP regression head. We saw a drastic improvement in performance after adding just a single attention block, followed by marginal improvements as more layers were added. The best-performing model was the largest one considered, where we used 2 layers in both the cross-attention and self-attention modules. This was the final model configuration selected for this work.

We use FlashAttention16 to minimize the computational cost of processing long sequence lengths (a 4-kb DNA sequence chunk amounts to an embedding matrix of dimension 667×512). The two sets of output embeddings were fused via mean pooling prior to the regression head, which comprises of a fully connected layer, a Tanh activation function, and a linear layer to output the predicted gene expression levels. The structure of the regression head is adapted from the official implementation of the BERT classifier.17

The model was trained on a single A100 GPU on the Minnesota Supercomputing Institute Agate cluster, using data from 232 donors on the set of 400 locally heritable genes. Training was performed over 100 epochs using the Adam optimizer, with a batch size of 64, an initial learning rate of 2×10−4, and a cosine annealing schedule.18 Regularization was implemented using a 10% dropout rate19 at the initial embedding layer, the transformer attention and MLP layers, and the output linear layer.

Enformer

The first benchmark is Enformer, the top performing model in two reviews of several sequence-to-expression DL models.8,10 As above, personalized nucleotide sequences are created for each individual by inserting their SNVs into the Homo Sapiens Assembly 38 reference sequence around each gene’s TSS, which is done separately for each of the individual’s two haplotypes. Enformer is then used to obtain separate outputs for each haplotype. Unlike the NT-based model above, the Enformer input sequence length is 393,216 bases.

Given a genetic sequence, Enformer produces predictions for 5,313 human output tracks in 986 bins. After obtaining output for all tracks and bins, we summarize output across bins by averaging the values from the center three bins (447, 448, and 449). With respect to combining values from different tracks, Sasse et al.8 performed a sensitivity test comparing using the single track best matching their data source versus weighting all tracks via an elastic net model trained to predict population average gene expression from the reference genome. Since using the elastic net model did not outperform using one track in their work, here, we use the best-performing track for each tissue to allow some additional flexibility beyond using a manually pre-specified track. More specifically, to find the best-performing track for each tissue, for every individual in the training set, for each of the 5,313 tracks, we compute the correlation across genes between the Enformer prediction in that track vs. observed gene expression. For this Enformer analysis, we use unnormalized gene transcripts per million as the measure of observed gene expression. The track with the highest cross-gene correlation averaged over individuals is chosen as the best-performing track for that tissue (Table 1). Although we are ultimately interested in correlation across individuals rather than across genes, we found that the average cross-individual correlations were too close to 0 to be able to meaningfully pick out a best-performing track. After using the best-performing track to obtain a prediction for each haplotype, these predictions are then averaged over the two haplotypes to produce one expression prediction for each individual for a particular gene and tissue.Table 1 Best-performing track for each tissue based on average cross-gene Pearson correlation

Tissue	Best track	Track name	Avg. r	
Muscle skeletal	5,231	CAGE: skeletal muscle – soleus muscle	0.198	
Whole blood	4,950	CAGE: blood, adult, pool 1	0.464	
Skin sun exposed lower leg)	4,970	CAGE: skin, adult	0.584	
Artery tibial	5,137	CAGE: CD133+ stem cells – cord blood derived, pool 1	0.591	
Adipose subcutaneous	5,137	CAGE: CD133+ stem cells – cord blood derived, pool1	0.663	
CAGE, cap analysis gene expression.

Track names are listed from https://www.github.com/calico/basenji/blob/master/manuscripts/cross2020/targets_human.txt.

To validate the Enformer predictions, we looked at the distribution of cross-individual and cross-gene correlations across the 58 test set individuals and the 400 locally heritable genes for each of the five tissues considered (Figure 5). For each gene, the cross-individual correlation is calculated as the Spearman rank correlation between Enformer predicted and observed expression values across the 58 individuals. Similarly, for each individual, cross-gene correlation is calculated as the rank correlation between predicted and observed expression across the 400 genes. Depending on the tissue, the average cross-gene correlation ranges between 0.35 and 0.45, while the cross-individual correlations are roughly symmetric about 0 for all tissues. This is analogous to the analysis described in Huang et al.,10 and these results show similar patterns as in Figure 1C in their paper; that is, Enformer is able to capture variation across genes to some extent but struggles to capture variation across individuals. One difference between the results is that the average cross-gene correlation in Huang et al.10 is about 0.55, but this is likely due to using different tissues and a different dataset.Figure 5 Distribution of cross-individual (top-row) and cross-gene Spearman rank correlations between Enformer predicted and observed gene expression for each tissue

For cross-individual correlations, each data point in the histogram is the correlation across individuals between predicted and observed gene expression for a fixed gene, and for cross-gene correlations, each data point is the correlation across genes for a fixed individual. Cross-individual correlations are roughly symmetric about 0 for all five tissues, while median cross-individual correlations are between 0.35 and 0.45, depending on the tissue. This figure is analogous to Figure 1C in Huang et al.

Elastic net linear regression

For our second benchmark, we used the elastic net linear regression model,20 a widely used method for gene expression prediction from single-nucleotide variation, especially in the context of transcriptome-wide association studies.6 In this case, the predictors were the allele counts of the SNPs lying within a 2-kb window either side of each gene’s TSS, extracted using Plink software.21 Variants with a minimum allele frequency of less than 0.01 were excluded. We used a fixed elastic net mixing parameter of 0.5, and, for a fair comparison, we trained individual models for each gene using the same training/validation split that we used for our transformer model. Models were trained using scikit-learn’s ElasticNetCV function on Python.22

Results

In this section, we compare the performance of our model to those of our two benchmarks, Enformer and elastic net regression, across the set of 400 locally heritable genes and 5 tissues considered (2,000 predictions in total). In each case, gene expression was predicted directly from the DNA sequence (or, in the case of elastic net, from the SNV allele counts) of the 58 donors in the test set, and the Pearson correlation was computed against the observed expressions for each donor. For Enformer, we did not normalize the gene expression values, which is consistent with the methodology of Huang et al.10 and does not impact correlation for a fixed gene. Since Sasse et al.8 and Huang et al.10 found that Enformer predictions often failed to correctly attribute the variants’ direction of effect but were more consistent with respect to the magnitude of the effect, we also compared the performance of our model to Enformer both in terms of the signed and absolute Pearson correlations. These performance comparisons are presented in Figures 6, 7, and 8. Some key performance comparison metrics are presented in Table 2.Figure 6 Comparison of model performance to elastic net on the set of 400 locally heritable genes across 5 tissues (muscle skeletal, whole blood, skin sun exposed lower leg, artery tibial, and adipose subcutaneous)

The metric used for evaluation is the Pearson correlation. Our model outperformed elastic net on around 34.5%, with a median underperformance of 0.02 across the 2,000 evaluations.

Figure 7 Comparison of model performance to Enformer on the set of 400 locally heritable genes across 5 tissues (muscle skeletal, whole blood, skin sun exposed lower leg, artery tibial, and adipose subcutaneous)

The metric used for evaluation is the Pearson correlation. Our model outperformed Enformer on around 87% of the 2,000 evaluations.

Figure 8 Comparison of model performance to Enformer on the set of 400 locally heritable genes across 5 tissues (muscle skeletal, whole blood, skin sun exposed lower leg, artery tibial, and adipose subcutaneous)

The metric used for evaluation is the absolute value of the Pearson correlation. Our model outperformed Enformer on around 83% of the 2,000 evaluations.

Table 2 Outperformance of our model vs. the various benchmarks across the set of the shortest 200 locally heritable genes

Benchmark	Median outperformance	Outperformance rate	
Elastic net	−0.02	0.35	
Enformer	0.29	0.87	
Enformer (absolute correlations)	0.19	0.83	
Median outperformance refers to the median pairwise difference in Pearson correlations of our model vs. each benchmark. “Outperformance rate” refers to the fraction of genes where our model outperformed the corresponding benchmark.

Additionally, we considered a few different aspects of relative performance to determine which factors were driving the performance differences between our model and elastic net regression. First, we considered the impact of tissues on relative performance (Figure 9). Here, a one-way ANOVA confirmed that we did not see much difference in relative performance across the five tissues considered (p = 0.995).Figure 9 Boxplot (ex-outliers) of pairwise correlation differences vs. elastic net regression across the five tissues considered

ANOVA F-statistic p-value of 0.995 indicates no significant difference in relative performance across tissues.

Second (Figure 10), we looked at the heritability measure (H squared) from the list of heritable genes,13 expecting to see better performance on our model for the genes with lower H-squared values, since the heritability measures were derived using linear models, and it is possible that our model could capture non-linear features in the data that are missed by the former. We saw a statistically significant negative slope (p = 1.46×10−10) on the pairwise performance difference (Pearson correlation measure for our model minus elastic net for each gene) and the heritability measure (square root of H squared) across the 400 locally heritable genes and 5 tissues considered for this study.Figure 10 Linear regression of the pairwise difference in correlations between our model and elastic net against the square root of the heritability measure H squared (sourced from the list we used to filter for heritable genes), measured across 400 locally heritable genes and 5 tissues (2,000 evaluations in total)

The negative slope indicates that our model’s performance relative to elastic net is better on genes with lower heritability.

Third (Figure 11), we considered the relation between pairwise performance difference and gene length. Here, we did not expect to see much difference, as, for the sake of consistency, we used the same window length around each gene’s TSS for both our model and elastic net. However, we did see a mild but statistically significant negative slope (p = 1.91×10−6) on the linear regression of correlation difference against the log of gene length, indicating that our model’s performance deteriorates slightly relative to elastic net as gene length increases.Figure 11 Linear regression of the pairwise correlation difference against the log of gene length

The negative slope indicates that, on average, our model’s performance relative to elastic net deteriorates slightly for longer genes.

We ran similar studies to analyze our model’s performance relative to Enformer. To summarize, the one-way ANOVA on relative performance across tissues (Figure 12) yielded a statistically significant F-statistic (p = 3.27×10−9), but this is largely due to Enformer, since our model’s performance was roughly the same across the tissues considered. We found that our model’s performance relative to Enformer increases with higher gene local heritability (Figure 13). Conversely, we found that our model’s relative performance deteriorated with increasing gene length (Figure 14). This outcome was anticipated, since Enformer utilizes much longer input sequences for prediction, which is more relevant for longer genes.Figure 12 Boxplot (ex-outliers) of pairwise correlation differences vs. Enformer across the five tissues considered

The ANOVA F-statistic p value of 3.27×10−9 indicates a significant difference in relative performance across tissues. This is largely due to Enformer, as our model’s performance is roughly the same across tissues.

Figure 13 Linear regression of the pairwise difference in correlations between our model and Enformer against the square root of the heritability measure H squared (sourced from the list we used to filter for heritable genes), measured across 400 locally heritable genes and 5 tissues (2,000 evaluations in total)

The positive slope indicates that our model’s performance relative to Enformer is better on genes with higher estimated heritability.

Figure 14 Linear regression of the pairwise correlation difference against the log of gene length

The negative slope indicates that, on average, our model’s performance relative to Enformer deteriorates for longer genes. This is expected behavior, since Enformer uses a much larger input sequence window for prediction and is therefore able to utilize sequence information that is not accessible to our model.

Finally, in order to assess our model’s ability to predict the expression levels for out-of-sample genes, we tested its performance on 200 genes randomly sampled from outside the set of locally heritable genes (Figures 15 and 16). We found that the model outperformed elastic net on around 44% of the 200 genes, with significant outperformance on several individual genes, while it only outperformed Enformer on around 48% of genes in this setting (in contrast to the 87% outperformance vs. Enformer on the set of locally heritable genes). We note that these genes were completely out of sample for our model, whereas for elastic net, we trained separate models for each gene using the training & validation sets of individuals and evaluated its test performance on the remaining individuals.Figure 15 Out-of-sample performance comparison of model to elastic net on 200 randomly selected genes lying outside the local heritability set across 5 tissues (1,000 points in total)

The model was trained on a set of 400 locally heritable genes, meaning that the genes used for testing were completely out of sample for our model (while individual models were trained for each of the genes in the case of elastic net). The metric used for evaluation is the Pearson correlation, and the tissues used were muscle skeletal, whole blood, skin sun exposed lower leg, artery tibial, adipose and subcutaneous. Our model outperformed elastic net on around 42% of the genes considered.

Figure 16 Out-of-sample performance comparison of model to Enformer on 200 randomly selected genes lying outside the local heritability set across 5 tissues (1,000 points in total)

The model was trained on a set of 400 locally heritable genes, meaning that the genes used for testing were completely out of sample for our model (while individual models were trained for each of the genes in the case of elastic net). The metric used for evaluation is the Pearson correlation, and the tissues used were muscle skeletal, whole blood, skin sun exposed lower leg, artery tibial, adipose subcutaneous. Our model outperformed Enformer on around 42% of the genes considered.

Discussion

In this work, we used pre-trained embeddings from a foundation model to improve personalized sequence-to-gene expression predictions for various tissues from the GTEx database. The embeddings were used as input to a downstream transformer model that predicted RNA expression values.

We compared the results from our model to Enformer, the top-performing model in two similar studies of sequence-to-expression DL models and an elastic net regression approach. Our model utilizing NT embeddings outperformed Enformer in correlation across individuals between observed and predicted expression values; in addition, it did not suffer as much from Enformer’s documented problem of mispredicting the sign of an effect. However, our approach was still outperformed, on average, by a much simpler elastic net model, albeit with a much smaller performance gap than in previous evaluations of such DL models. Elastic net has an advantage in prediction in that it takes only small sets of pre-determined SNP allele counts as inputs and does not need to learn patterns in DNA sequences that generalize across genes and tissues. Elastic net’s simplicity also perhaps limits its future utility in that each regression model is restricted to a single gene and can only use the set of SNPs that was selected during training. On the other hand, sequence-to-expression DL models can theoretically predict the effect of any genetic variant and their combinations, even rare and unseen ones, based on the transferable patterns they learn from genome sequences. Thus, although the performance of sequence-to-expression models for individualized expression prediction has been shown to be somewhat limited at the moment, these models have tremendous potential as more genome and paired DNA sequence-gene expression data are incorporated.

Several avenues exist to further enhance the performance of these models. We anticipate significant improvements by increasing the number of genes and tissues used for training, as well as extending the DNA sequence lengths used as inputs. In this study, we selected genes based on their length and local heritability and chose tissues based on data availability. However, the optimal selection of genes and tissues for training remains an open question. We hypothesize that a more systematic and principled approach to gene and tissue selection could lead to better predictive performance.

Future research could explore the impact of varying the size, complexity, and type of model used to generate embeddings. Although we focused on the NTv2-100m-multispecies model in this work due to its strong performance, computational efficiency, and lower embedding dimension, larger versions within the same model family or other foundation models could also be considered. For instance, HyenaDNA23 allows for much longer sequence lengths and tokens at single-nucleotide resolutions, while DNABERT-224 uses byte pair encoding to address some limitations of k-mer tokenization.

Data and code availability

Individual-level GTEx WGS data are accessible through dbGaP Project 26511, while individual-level GTEx expression data are available at https://gtexportal.org/home/downloads/adult-gtex/bulk_tissue_expression. The code used for training the models and the lists of donor and gene IDs used for this study are available at https://github.com/pratik-r/dlfse/.

Acknowledgments

This research was supported by 10.13039/100000002 NIH grants U01 AG073079 (to P.R. and W.P.), R01 AG065636 (to W.P.), R01 AG069895 (to P.R. and W.P.), RF1 AG067924 (to W.P.) and T32 GM132063 (to N.P. and W.P.). The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health and by 10.13039/100000054 NCI , 10.13039/100000051 NHGRI , 10.13039/100000050 NHLBI , 10.13039/100000026 NIDA , 10.13039/100000025 NIMH , and 10.13039/100000065 NINDS . The authors would also like to acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing resources that contributed to the research results reported within this paper.

Declaration of interests

The authors declare no competing interests.

Web resources

Code repository, https://github.com/pratik-r/dlfse/

Enformer track names, https://www.github.com/calico/basenji/blob/master/manuscripts/cross2020/targets_human.txt

GTEx bulk expression data download, https://gtexportal.org/home/downloads/adult-gtex/bulk_tissue_expression

TWAS/FUSION (for list of locally heritable genes), http://gusevlab.org/projects/fusion/
==== Refs
References

1 Avsec Ž. Agarwal V. Visentin D. Ledsam J.R. Grabska-Barwinska A. Taylor K.R. Assael Y. Jumper J. Kohli P. Kelley D.R. Effective gene expression prediction from sequence by integrating long-range interactions Nat. Methods 18 2021 1196 1203 34608324
2 Zhou J. Theesfeld C.L. Yao K. Chen K.M. Wong A.K. Troyanskaya O.G. Deep learning sequence-based ab initio prediction of variant effects on expression and disease risk Nat. Genet. 50 2018 1171 1179 30013180
3 Dalla-Torre H. Gonzalez L. Mendoza-Revilla J. Carranza N.L. Grzywaczewski A.H. Oteri F. Dallago C. Trop E. de Almeida B.P. Sirelkhatim H. The nucleotide transformer: Building and evaluating robust foundation models for human genomics Preprint at bioRxiv 50 2023 1171 10.1101/2023.01.11.523679
4 Kelley D.R. Snoek J. Rinn J.L. Basset: learning the regulatory code of the accessible genome with deep convolutional neural networks Genome Res. 26 2016 990 999 27197224
5 Agarwal V. Shendure J. Predicting mrna abundance directly from genomic sequence using deep convolutional neural networks Cell Rep. 31 2020 107663
6 Gamazon E.R. Wheeler H.E. Shah K.P. Mozaffari S.V. Aquino-Michaels K. Carroll R.J. Eyler A.E. Denny J.C. GTEx ConsortiumNicolae D.L. A gene-based association method for mapping traits using reference transcriptome data Nat. Genet. 47 2015 1091 1098 26258848
7 Novakovsky G. Dexter N. Libbrecht M.W. Wasserman W.W. Mostafavi S. Obtaining genetics insights from deep learning via explainable artificial intelligence Nat. Rev. Genet. 24 2023 125 137 36192604
8 Sasse A. Ng B. Spiro A.E. Tasaki S. Bennett D.A. Gaiteri C. De Jager P.L. Chikina M. Mostafavi S. Benchmarking of deep neural networks for predicting personal gene expression from dna sequence highlights shortcomings Nat. Genet. 55 2023 2060 2064 38036778
9 Kelley D.R. Cross-species regulatory sequence activity prediction PLoS Comput. Biol. 16 2020 e1008050
10 Huang C. Shuai R.W. Baokar P. Chung R. Rastogi R. Kathail P. Ioannidis N.M. Personal transcriptome variation is poorly explained by current genomic deep learning models Nat. Genet. 55 2023 2056 2059 38036790
11 Tang Z. Toneyan S. Koo P.K. Current approaches to genomic deep learning struggle to fully capture human genetic variation Nat. Genet. 55 2023 2021 2022 38036789
12 Lonsdale J. Thomas J. Salvatore M. Phillips R. Lo E. Shad S. Hasz R. Walters G. Garcia F. Young N. The genotype-tissue expression (gtex) project Nat. Genet. 45 2013 580 585 23715323
13 Gusev A. Ko A. Shi H. Bhatia G. Chung W. Penninx B.W.J.H. Jansen R. de Geus E.J.C. Boomsma D.I. Wright F.A. Integrative approaches for large-scale transcriptome-wide association studies Nat. Genet. 48 2016 245 252 26854917
14 McKenna A. Hanna M. Banks E. Sivachenko A. Cibulskis K. Kernytsky A. Garimella K. Altshuler D. Gabriel S. Daly M. DePristo M.A. The genome analysis toolkit: a mapreduce framework for analyzing next-generation dna sequencing data Genome Res. 20 2010 1297 1303 20644199
15 Tsai Y.-H.H. Bai S. Liang P.P. Zico Kolter J. Morency L.-P. Salakhutdinov R. Multimodal transformer for unaligned multimodal language sequences Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), Florence, Italy 2019 Association for Computational Linguistics
16 Dao T. Fu D. Ermon S. Rudra A. Flashattention C.R. Fast and memory-efficient exact attention with io-awareness Adv. Neural Inf. Process. Syst. 35 2022 16344 16359
17 Jacob D. Chang M.-W. Lee K. Toutanova K. BERT: Pre-training of deep bidirectional transformers for language understanding Burstein J. Doran C. Solorio T. Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers) 2019 Association for Computational Linguistics Minneapolis, Minnesota 4171 4186 10.18653/v1/N19-1423 https://aclanthology.org/N19-1423
18 Loshchilov I. Hutter F. Sgdr: Stochastic gradient descent with warm restarts Preprint at arXiv 50 2016 1171 10.48550/arXiv:1608.03983
19 Srivastava N. Hinton G. Krizhevsky A. Sutskever I. Salakhutdinov R. Dropout: a simple way to prevent neural networks from overfitting J. Mach. Learn. Res. 15 2014 1929 1958
20 Zou H. Hastie T. Regularization and variable selection via the elastic net J. Roy. Stat. Soc. B Stat. Methodol. 67 2005 301 320
21 Purcell S. Neale B. Todd-Brown K. Thomas L. Ferreira M.A.R. Bender D. Maller J. Sklar P. de Bakker P.I.W. Daly M.J. Sham P.C. Plink: a tool set for whole-genome association and population-based linkage analyses Am. J. Hum. Genet. 81 2007 559 575 17701901
22 Pedregosa F. Varoquaux G. Gramfort A. Vincent M. Bertrand T. Grisel O. Blondel M. Prettenhofer P. Weiss R. Vincent D. Scikit-learn: Machine learning in python J. Mach. Learn. Res. 12 2011 2825 2830
23 Nguyen E. Poli M. Faizi M. Thomas A. Wornow M. Birch-Sykes C. Massaroli S. Patel A. Clayton R. Bengio Y. Hyenadna: Long-range genomic sequence modeling at single nucleotide resolution Adv. Neural Inf. Process. Syst. 36 2024
24 Zhou Z. Ji Y. Li W. Dutta P. Davuluri R. Liu H. Dnabert-2: Efficient foundation model and benchmark for multi-species genome Preprint at arXiv 50 2023 1171 10.48550/arXiv:2306.15006
