Transfer learning using generative pretrained genomic DNA models for predicting perturbation-induced changes in gene expression

Transfer learning using generative pretrained genomic DNA models for predicting perturbation-induced changes in gene expression | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Transfer learning using generative pretrained genomic DNA models for predicting perturbation-induced changes in gene expression Takuya Shihashi, Itoshi Nikaido This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6461531/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Background Transfer learning applied to genomic DNA models has the potential to improve predictive capabilities, especially when target-domain datasets and computational resources are limited. Despite its promise, the practical effectiveness of transfer learning in genomic DNA models, particularly for predicting gene expression changes due to perturbations, has not been thoroughly investigated. This study aimed to systematically evaluate the performance and utility of transfer learning approaches using genomic DNA models to accurately predict perturbation-induced gene expression. Results We benchmarked three genomic DNA models across 12 distinct datasets containing perturbation-induced gene expression data to identify optimal conditions for effective transfer learning. Notably, perturbation-induced gene expression data were not included in the pre-training of these genomic DNA models. Among these, the Enformer model consistently generated accurate embeddings, demonstrating superior clustering performance and gene signature scoring aligned closely with observed experimental data. Additionally, we identified a phenomenon termed "genomic neighbouring gene interference," wherein partially overlapping DNA sequences of adjacent genes cause correlated predictions, resulting in both beneficial and detrimental effects on predictive accuracy. Conclusions Our findings highlight the efficacy of transfer learning in genomic DNA models for predicting perturbation-induced gene expression, particularly emphasizing the Enformer model's robust performance. Understanding genomic neighbouring gene interference offers critical insights for refining predictive accuracy in genomic applications. This study provides practical guidance for researchers developing transfer learning strategies and genomic DNA models, paving the way for more accurate and resource-efficient genomic predictions. Transfer learning Genomic DNA models Genome Language Model Gene expression prediction Perturbation Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Background In natural language processing (NLP), fine-tuning pretrained models such as large language models (LLMs) enables high performance on specific tasks, even with small datasets and limited computer resources 1 , 2 . There has been some success in applying deep learning architectures used in NLP to genomic DNA models 3 – 6 , expecting that fine-tuning with small datasets can be also effective in genomic DNA models 7 , 8 . This approach is beneficial when advancements in measurement technology yield new types of data, where it is challenging to prepare measurement data with sufficient quantity and diversity. However, research on transfer learning in this area is currently insufficient, and many uncertainties remain regarding its practical utility. Recently, advancements in technologies such as Perturb-seq and SciPlex have made it possible to comprehensively profile gene expression following external perturbations, such as gene overexpression or drug treatment 9 , 10 . The development of models capable of predicting perturbation-induced gene expression changes from DNA sequences can help elucidate the relationships between perturbations and cis-regulatory elements (CREs), deepening our understanding of transcriptional regulation and contributing to the construction of gene regulatory networks (GRNs) 11 . Moreover, combining this approach with DNA mutation simulations could reveal the relationship between noncoding DNA mutations and drug responses, with potential applications in personalized medicine. Given the current limitations in the numbers of cell types and perturbations available in the data, the use of transfer learning is a practical approach for developing predictive models. At present, Enformer is a representative genomic DNA model. This model combines a Transformer 12 with a convolutional neural network (CNN) to predict gene expression and epigenetic marker tracks from DNA sequences 3 . Other examples include genome language models (gLMs), such as HyenaDNA and Nucleotide Transformer, which predict enhancer element activity and splicing sites from DNA sequences 4 , 5 . These models have rarely been trained with perturbation-induced gene expression data, and there has yet to be an exploration of building predictive models for perturbation data via transfer learning. In this study, we benchmarked transfer learning in genomic DNA models, demonstrating that it is possible to predict perturbation-induced gene expression changes from DNA sequences. We evaluated multiple conditions across three genomic DNA models and 12 datasets, revealing that transfer learning via Enformer can effectively model these changes. By employing deep learning model interpretation methods, we confirmed that transfer learning can focus on regulatory elements and transcription factor motifs on the basis of biological principles. Moreover, we discovered that when the input DNA sequences have regions of overlap, the predicted gene expression values often strongly correlate with each other but deviate from the actual observed values. We referred to this issue as 'genomic neighbouring gene interference'. Results Design of transfer learning and overview of the evaluation workflow The overall structure of this study is illustrated in Fig. 1 . We used Enformer, HyenaDNA, and Nucleotide Transformer as the genomic DNA models. The parameters of the additional heads were trained using the embedding values from each model as input. During this process, we considered three approaches: (1) freezing the parameters of the pretrained models (feature-based), (2) fine-tuning a subset of the parameters via low-rank adaptation (LoRA) 13 , or (3) fine-tuning all of the parameters. Enformer uses the embedding values of 512 bp around the transcription start site (TSS) for the input values to the heads. HyenaDNA also uses the embedding values of 512 bp around the TSS. Since the Nucleotide Transformer, tokenizes sequences in 6 bp units, we used embedding values of 516 bp. In NLP classification tasks, the last token is often used for next-token prediction models, and the CLS token is used for masked token prediction models 4 , 5 . Therefore, for the HyenaDNA and Nucleotide Transformer, we also conducted training using these embedding values as inputs (Additional file 1: Supplementary Fig. 1a). We prepared 12 datasets and conducted training with each dataset 10 , 14 – 17 . These datasets include genome-wide screens with perturbations such as clustered regularly interspaced short palindromic repeats (CRISPR) activation, CRISPR inhibition, compound treatments, and coding DNA sequence (CDS) overexpression. To mitigate the reduction in model accuracy caused by noise originating from measurement techniques, we converted single-cell data into pseudobulk data and used it for model training. Moreover, we excluded nonexpressed genes. The numbers of perturbations and genes in each dataset are listed in Additional file 1: Supplementary Table 1. Selection of transfer learning methods for benchmarking There are multiple methods for transfer learning 1 . We examined the differences in accuracy among three learning methods: (1) Obtaining the embedding values in advance and training only the additional head (feature-based). The pretrained model combined with the head or the head with precomputed embedding values is used for prediction. (2) LoRA is used to train parts of the transformer layers along with the additional head 13 . The number of parameters updated is set by the rank (in this study, 16, 64, 256, and 512). (3) All parameters of both the pretrained model and the additional head are trained. The learning rate for the pretrained model was set either equal to or lower than that of the additional head. In the evaluation using the CD8T dataset by Jiang et al. 17 , the feature-based approach with precomputed embedding values achieved the highest accuracy in terms of correlations across genes and perturbations (Additional file 1: Supplementary Fig. 2a and b). With respect to correlations across perturbations, since genes with higher expression levels led to greater accuracy, we divided the genes into five groups on the basis of expression levels and examined the distribution of correlation coefficients. To verify the accuracy differences among training methods in correlations across genes, we compared the observed and predicted values in nonperturbed samples. Compared with the feature-based method, fine-tuning caused the model's predictions to become more uniform across different genes, leading to a decrease in correlation coefficients (Additional file 1: Supplementary Fig. 2c). Similar results were observed with HyenaDNA and the Nucleotide Transformer, where training only the feature-based approach using precomputed embedding values yielded the highest accuracy. In some cases, catastrophic forgetting appeared to occur (Additional file 1: Supplementary Fig. 2d-g). On the basis of these results, we conducted further investigations using the feature-based approach for the transfer learning method. Performance of transfer learning models for predicting perturbation-induced gene expression changes We applied transfer learning using three genomic DNA models, Enformer, HyenaDNA, and Nucleotide Transformer, across 12 datasets to evaluate the correlations across genes and perturbations. The Enformer model demonstrated higher accuracy, with a median gene correlation of approximately 0.6, outperforming HyenaDNA and the Nucleotide Transformer (Fig. 2 a and Additional file 1: Supplementary Fig. 3a). When the original Enformer model was compared with the transfer learning model in nonperturbed K562 cells, it was suggested that transfer learning improved gene correlations and better adapted the model to the additional data (Fig. 2 b). With respect to perturbation correlations, the Enformer model exhibited a median correlation between 0.2 and 0.5 for high-expression genes, whereas near-zero correlations were observed for low-expression genes, likely due to the limitations of measurement sensitivity (Fig. 2 c). The accuracy of the HyenaDNA and Nucleotide Transformer methods was generally lower (Fig. 2 c and Additional file 1: Supplementary Fig. 3b). To evaluate the differences among perturbations, uniform manifold approximation and projection (UMAP) embedding and Leiden clustering were performed via expression values. Enformer generated embeddings that were consistent with the observed data, whereas HyenaDNA and Nucleotide Transformer generated embeddings with inconsistencies (Fig. 3 a-d and Additional file 1: Supplementary Fig. 3c and d). Moreover, Enformer outperformed the other methods in terms of clustering concordance, as indicated by its higher adjusted Rand index (ARI) and normalized mutual information (NMI) scores (Fig. 3 e and Additional file 1: Supplementary Fig. 3e). The UMAP of observed values from the dataset generated by Jiang et al. 17 revealed good separation, and the model trained with this dataset achieved high accuracy (Fig. 2 c and Fig. 3 e-g). In contrast, the UMAP of observed values from the dataset generated by Srivatsan et al. 10 indicated poor separation, and the model trained with this dataset demonstrated low accuracy (Fig. 2 c and Fig. 3 e,h-i). In summary, the Enformer transfer learning model achieved the highest accuracy. The accuracy depends on the characteristics of the dataset used for training. Transfer learning accurately predicts biological phenotypes and gene regulatory networks We evaluated whether the Enformer transfer learning model could capture biological features. Since the number of genes that can be assessed with test data alone is limited, we evaluated the gene signatures via the prediction results of all the genes. As there was no apparent difference in the embedding results when predictions of all genes were used compared with when only the test data were used, we proceeded with this approach (Fig. 3 a-e and Additional file 1: Supplementary Fig. 4a-e). The Norman et al. 14 dataset includes biological function labels for each perturbation 14 , and by examining the distribution of these labels on UMAP plots of the observed and predicted values, we confirmed that perturbations with the same label were clustered together (Fig. 4 a and Additional file 1: Supplementary Fig. 4f). In terms of cell differentiation, the gene signature scores for erythroid and granulocyte lineages were consistent with their labels (Fig. 4 a and b). Next, we calculated activation scores for each transcription factor via gene signatures from CollecTRI 18 , which summarizes transcriptional regulatory networks on the basis of ChIP-seq data, and ranked the relative scores for each perturbation. Consistently high activation scores were observed for transcription factors such as HNF4A, IRF1, and TP73 in both the observed and predicted values when they were overexpressed via CRISPRa (Fig. 4 c and d). However, the predicted rank for transcription factors such as AHR was significantly lower, suggesting that Enformer may not have adequately learned AHR-mediated regulation during pretraining. Finally, no apparent differences in correlation were detected between the housekeeping genes and the other genes, suggesting that the model's accuracy does not depend on the housekeeping genes (Additional file 1: Supplementary Fig. 4g) 19 . Overall, these findings suggest that the Enformer transfer learning model accurately predicts cell differentiation and gene regulatory networks. Transfer learning model using Enformer attention to transcriptional regulatory regions and transcription factor motifs The Enformer transfer learning model accurately predicted gene expression changes due to perturbations and outperformed HyenaDNA and Nucleotide Transformer. To explore the differences in prediction accuracy, we identified the DNA sequences that the models focused on during predictions via attribution scores derived from Input X Gradient 3 , 20 , 21 . By subtracting the attribution scores of the control group from those of the perturbation group, we evaluated the difference in attribution scores and identified the critical sequences necessary for predicting expression changes induced by perturbations. Since we confirmed that the fold change in the predicted values correlated with the subtracted attribution score differences in certain genomic regions, we proceeded with this approach. (Additional file 1: Supplementary Fig. 5). We employed the Norman et al. 14 dataset to visualize the genomic regions highlighted by attribution scores during the prediction of TYROBP expression when CEBPA was overexpressed via CRISPRa, overlaying these regions with ATAC-seq and ChIP-seq data from the same K562 cell line (Fig. 5 a). Enformer focused on regulatory elements that activate transcription, such as ATAC-seq H3K27ac and H3K4me3. In contrast, HyenaDNA, when TSS-centred embeddings are used, focuses only on the region from the start of the input to the TSS, and when the final token is used, it fails to focus on regulatory elements. The Nucleotide Transformer, which was constrained by the short input sequence, focused only on the region surrounding the TSS. To verify whether these models focused on regulatory elements, we checked the enrichment of attribution scores in ATAC-seq peak regions. Enformer, as well as HyenaDNA and Nucleotide Transformer, when using TSS-centred embeddings, focused on ATAC-seq peaks, whereas models using the final token or CLS token did not focus on ATAC-seq peaks (Fig. 5 b). We further evaluated the Enformer transfer learning model. When we assessed whether the model focused on the ChIP-seq peaks of the same transcription factors during gene expression changes caused by transcription factor overexpression, the area under the precision–recall curve (AUPRC) ratio ranged from 3 to 6. This finding suggests that the model predicts a focus on the ChIP-seq peaks (Fig. 5 c). Furthermore, via TF-MoDISco 22 , we detected transcription factor motifs from the attribution scores for all the transcription factor perturbations included in the Norman et al. 14 dataset. Between 20% and 30% of the transcription factors had motifs detected by TF-MoDISco that matched the known motifs registered in JASPAR 23 , such as HNF4A, IRF1, CEBPA, and SPI1. Additionally, 70% of the transcription factors, such as those belonging to the same gene family, matched one of the motifs within the cluster (Fig. 5 d and e). In summary, the Enformer model achieved high predictive accuracy by focusing on CREs during prediction and recognizing transcription factor motifs at the nucleotide level. On the other hand, HyenaDNA and Nucleotide Transformer struggled to capture the sequences involved in transcriptional regulation compared with Enformer. Challenges in predicting gene expression changes from DNA sequences: Genomic neighbouring gene interference problem Identifying challenges in model predictions can guide algorithm development. In Fig. 2 c, all transfer learning models included some genes showing negative correlations between the observed and predicted values, prompting an investigation into the causes. Examining the distribution of attribution scores for genes neighbouring NFKBIB (r = -0.61) and MRPS12 (r = -0.46), which had the strongest negative correlations in Norman et al.'s 14 test data, revealed that neighbouring genes focused on each other’s TSS regions for their prediction (Fig. 6 a). The observed values of the neighbouring gene RINL were negatively correlated with both NFKBIB (r = -0.44) and MRPS12 (r = -0.63) (Fig. 6 b). However, in terms of the predicted values, RINL was positively correlated with NFKBIB (r = 0.62) and MRPS12 (r = 0.55), and other neighbouring genes also learned to correlate with RINL’s observed values (Fig. 6 c). These findings suggest that interference during training involving neighbouring genes could cause negative correlations between the observed and predicted values for NFKBIB and MRPS12. Beyond the NFKBIB region, we calculated pairwise correlation coefficients for all genes, grouping them by the distance between their TSS regions. The results indicated that neighbouring genes within ten kbp showed a strongly positive correlation with the predicted values (Fig. 6 d). We also investigated whether genes with neighbouring genes that had negatively correlated observed values displayed negative correlations in predictions. As expected, most genes with negatively correlated predictions had neighbouring genes within 200 kbp that were negatively correlated with the observed values. In these cases, one gene showed a negative correlation in the predictions, whereas the other was correctly predicted to be positively correlated (Fig. 6 e). Additionally, when neighbouring genes showed a strong positive correlation with the observed values, the predicted values tended to show similarly strong positive correlations (Fig. 6 f). To test whether this issue improves when the input sequences of neighbouring genes do not overlap, we performed transfer learning via masked input sequences, excluding regions outside the 10 kbp surrounding the TSS. Masking reduced the positive correlation between genes with TSS distances of 5–10 kbp, likely due to a reduction in the overlap of input DNA sequences to 50% or less. (Additional file 1: Supplementary Fig. 6a and b). As a supplement, masking the input sequences partially improved negative correlation of NFKBIB and MRPS12, and decreased the accuracy of correlations across genes and perturbations, so it did not address the root cause of the problem (Additional file 1: Supplementary Fig. 6c-e). In summary, the findings suggest that overlapping DNA sequences between neighbouring genes cause interference, affecting prediction accuracy both positively and negatively because of the observed correlations among these genes. We termed this phenomenon ‘genomic neighbouring gene interference’. Further improvements in model structure and training methods are needed to address these issues. Discussion In this study, we evaluated three genomic DNA models, 12 datasets, and multiple transfer learning conditions and demonstrated that transfer learning via Enformer is effective for predicting gene expression changes due to perturbations. In contrast, HyenaDNA and Nucleotide Transformer did not achieve sufficient accuracy. The datasets used in this study involved external perturbations, yet the Enformer model had not been pretrained on such data. Nevertheless, the Enformer model was able to predict gene expression changes, suggesting that it implicitly learned the principles of biological transcriptional regulation, particularly GRNs, from large-scale gene expression and epigenetic data. These GRNs drive gene expression changes in response to perturbations, allowing Enformer to make accurate predictions during transfer learning without becoming out-of-distribution. However, we identified a problem termed "genomic neighbouring gene interference," where the prediction results revealed correlations between neighbouring genes. One possible cause of this problem is that genomic features, such as CREs, are not properly associated with each individual gene 24 . To address this issue, increasing the amount of data used in pretraining, such as ChIP-seq and ATAC-seq, and the number of model parameters may naturally improve the association between genes and CREs, similar to the scaling laws observed in NLP 25 . Additionally, incorporating 3D genome structure data or explicitly including meta-information, such as cell type, in the input during pretraining could also enhance the model's performance. HyenaDNA and Nucleotide Transformer did not achieve satisfactory accuracy through transfer learning. While these models have demonstrated high accuracy in DNA sequence classification tasks that do not depend on cell type, such as predicting splicing sites and promoters. The results of this study suggest that these models are unable to learn cell type- and context-dependent information, such as transcriptional regulation. Therefore, it will be necessary to develop learning methods and datasets that capture the principles of transcriptional regulation on the basis of cell type and context 4 , 5 . Also, increasing the length of DNA sequences alone, which is the basic future task for gLM, may not improve accuracy for gene expression prediction. The feature-based method that achieved the highest accuracy in transfer learning with the Enformer model utilized less than 1 GB of GPU memory on the NVIDIA RTX 6000 Ada, completing training in a maximum of 10 minutes. In contrast, fine-tuning, including the pretrained model, requires 20–45 GB of GPU memory per batch and requires 10–24 hours of training on a single GPU for 40 epochs. In both cases, this approach significantly reduces the computational resources compared with that of the Enformer pretraining process, which requires 32 TPU v3 cores and three days for training 3 . The findings demonstrate that transfer learning is beneficial in terms of computational efficiency. Finally, this study provides an example of transfer learning in genomic DNA models and offers valuable insights for individual researchers conducting transfer learning and for those developing genomic DNA models. This study will lead to an increase in the efficiency of the development and utilization of genomic DNA models. Conclusions This study demonstrated that transfer learning with the Enformer model effectively predicts perturbation-induced gene expression changes, significantly outperforming HyenaDNA and Nucleotide Transformer. Despite lacking explicit perturbation data in pretraining, Enformer's implicit grasp of gene regulatory networks enables accurate predictions. However, genomic neighboring gene interference revealed limitations in properly associating genomic regulatory elements with their target genes. These insights will accelerate the efficient development and application of genomic DNA models. Methods Data preprocessing for perturbation data We used 12 datasets from five studies 10,14-17 . The data from Norman et al. 14 and Replogle et al. 15 were obtained from scPerturb (http://projects.sanderlab.org/scperturb/) 26 , and the data from Srivatsan et al. 10 were downloaded from the sc-pert repository (https://github.com/theislab/sc-pert), as the annotations in the scPerturb database were misaligned. The remaining data were downloaded from the original articles (GSE217460, https://data.caltech.edu/records/2cjss-wgh69). The dataset from Jiang et al. 17 was divided into CD4 T cells, CD8 T cells, B cells, and myeloid cells in the anti-CD3/CD28 treatment group. Owing to the low cell numbers, NK cells and certain other cell types or conditions in the anti-CD3/CD28 nontreatment group were not used. Because the datasets from Norman et al. 14 and Replogle et al. 15 involve CRISPR perturbations, we excluded unperturbed cells via Mixscape 27 , which is implemented in the pertpy library (v0.7.0) 28 , to clearly capture gene expression changes due to perturbations. For the Mixscape calculations, we used the highly variable genes calculated by the highly_variable_genes function in Scanpy (v1.10.2) 29 , with min_disp=0.2 for each perturbation. Each dataset was pseudobulked for each perturbation via adpbulk (v0.1.4) to remove noise and improve learning efficiency. Only perturbations with more than 100 cells were used for analysis. To prevent artificially induced gene expression profiles from affecting the learning of natural transcriptional regulatory mechanisms, we replaced the expression of target genes in the overexpression or knockdown perturbation groups with values from the nonperturbation group. The counts after pseudobulking were normalized to the CPM, converted to log2 + 1, and used for training. To exclude nonexpressed genes, we used only genes with a CPM exceeding 2 in at least one perturbation for training. Model architecture and input genome The model input consisted of genomic sequences centred on the gene's TSS. The genes were divided into training, validation, and test sets according to the classification set in the original Enformer paper, which was obtained from https://console.cloud.google.com/storage/browser/basenji_barnyard/data 3 . Only mRNAs and lncRNAs were used for training, which was based on GFF3 annotations. We used the same versions of GFF3 (v32) and FASTA (GRCh38.p13) from GENCODE 30 as in the original Enformer paper. The pretrained models used were Enformer 3 , HyenaDNA 4 , and Nucleotide Transformer 5 . These models were selected from those registered on Hugging Face 31 for their ability to handle long context lengths. Enformer We used the sequence centred on the TSS, spanning ± 98,304 bp (total 196,608 bp) as input and connected the embedding values before the final layer to an additional single-layer head. Since Enformer performs convolution operations every 128 bp and then proceeds with the transformer, we used the embedding value from the 4 bins (512 bp) surrounding the TSS as the embedding value. HyenaDNA (medium-160k) We used the sequence centred on the TSS, spanning ± 80,000 bp (total 160,000 bp) as input, using either the embedding value of the 512 bp token centred on the TSS (TSS embedding) or, following the original paper, the last token as the input for the additional head 4 . The last token was treated as a single-layer head, as in Enformer, but the TSS embedding failed to learn because of too many input elements. Therefore, we used a 2-layer head with a hidden layer of 3,072 nodes for learning. Nucleotide Transformer (2.5b-multispecies) We used the sequence centred on the TSS, spanning ± 12,288 bp (total of 24,576 bp) as input. Since Nucleotide Transformer tokenizes every six bp, we used the embedding value of 516 bp (86 tokens) as the TSS embedding value for the additional head input. We also used the embedding value of the CLS token, which is used in classification models based on bidirectional encoder representations from transformers (BERT) 32 . Similar to the approach for HyenaDNA, the CLS token was learned with a single layer, whereas the TSS embedding value was learned with a 2-layer head with a hidden layer of 3,072 nodes. Training methods for transfer learning For all the training methods, PyTorch (v2.3.1) 33 was used for training, with a learning rate of 0.0001, a gradient clip value of 0.2, and optimization with adaptive moment estimation with weight decay (AdamW), where the weight decay was set to 0.005. The loss function used was the mean square error (MSE) loss, and training was terminated when the validation loss did not improve for five epochs during training. The batch size was set to 256, and for methods that could not accommodate this batch size, the parameters were updated by accumulating gradients over 256 batches. The following five methods were implemented for training: Pretrained models were used to obtain embedding values in advance, and only the additional head was trained for 100 epochs. After training, predictions were made using only the head from the precomputed embedding values. Since there was no difference in accuracy compared with that of method (2) in terms of prediction from DNA sequences, the predictions from this method were used for evaluations, as shown in Additional file 1: Supplementary Fig. 2. Pretrained models were used to obtain embedding values in advance, and only the additional head was trained for 100 epochs. After training, predictions were made by combining the pretrained model and the head to predict from DNA sequences. LoRA 13 was used to update the parameters of the query and value in all transformer layers, as well as the parameters for the additional head. The number of epochs was 20, with the LoRA parameters set to α=2, bias=none, and the rank was varied between 16, 64, 256 and 512 for training. The entire set of parameters for both the pretrained model and the additional head were trained for 40 epochs. During this process, the learning rate for the parameters of the pretrained model was varied between 0.0001 and 1×10⁻¹⁰. Evaluation of correlation coefficients We calculated two types of Pearson correlations for model evaluation. The first type is the correlation across perturbations, which is calculated between the actual and predicted values for each gene. The second type is the correlation across genes, which is calculated between the actual and predicted values for each perturbation after computing the predicted values for all genes. Both evaluations were performed using only the genes in the test set. The predicted values for nonperturbation conditions in the original Enformer were obtained by summing the predictions of the four bins (512 bp) surrounding the TSS of the CAGE-seq (K562: CNhs12336) track of the same cell type used in the transfer learning dataset, normalizing to CPM, and then converting to log2 + 1. Embedding of gene expression We created an AnnData object to store the pseudobulked actual and predicted values as a matrix of each perturbation and gene. These values were processed via Scanpy. Then, we performed principal component analysis (PCA), k-nearest neighbours (kNN), Leiden clustering, and UMAP calculations with default parameters in Scanpy. For the evaluation of the test set, we used only the genes in the test set, and for evaluating all genes, we used all genes, including the training, validation and test sets. Clustering consistency was evaluated by performing Leiden clustering for both the observed and predicted values and then using scikit-learn (v1.5.1) to calculate the ARI and NMI. Gene signature The gene signature scores related to differentiation were calculated via the gene sets for Myeloid and Erythroid lineages provided by Norman et al. 14 , employing the score_genes function in Scanpy. TF activity was obtained by retrieving information on transcription factors and downstream genes from CollecTRI 18 and calculating a univariate linear model (ULM) with a decoupler (v1.7.0) 34 . Human housekeeping gene lists were downloaded from the HRT Atlas v1.0 19 , and correlation coefficient box plots were drawn by dividing the genes into housekeeping and nonhousekeeping categories. Attribution score calculation The attribution score for Enformer was calculated via the Input X Gradient method for the target gene and each perturbation with the InputXGradient 20 function from Captum (v0.7.0) 35 . When focusing on changes due to perturbations, we subtracted the Input X Gradient value of the nonperturbation from that of the target perturbation for the same gene. For HyenaDNA and Nucleotide Transformer, which accept tokens as input, we calculated GradientXActivation 36 via Captum's function for tokens and summed the scores for each base to calculate the attribution score for each base (https://captum.ai/tutorials/IMDB_TorchText_Interpret). Track plot The attribution score calculated by Input × Gradient was averaged every 128 bp and converted to a bigwig file via pyBigWig (v0.3.23). For the track in Fig. 6a, we calculated the Input X Gradient for all perturbations for a single gene and used the maximum absolute value of the attribution scores across all perturbations for each 128 bp bin to create the bigwig file. The ChIP-seq and ATAC-seq data were downloaded as bigwig files from ChIP-Atlas 37 . Tracks were visualized via CoolBox (v0.3.9) 38 . The gene locus information used for visualization with CoolBox was obtained from the GTF (v32) used during transfer learning. For peak enrichment, we downloaded the ATAC-seq peaks (q value < 1E-05) of the K562 cell line (SRX10184518) from ChIP-Atlas. We calculated the attribution scores for all test set genes in the group overexpressing CEBPA, converted them to absolute values, and then calculated the enrichment score via computeMatrix and plotProfile from deepTools2 39 . AUPRC ratio for transcription factor ChIP-seq peaks ChIP-seq peaks for transcription factors were downloaded from ChIP-Atlas. When ChIP-seq data for the K562 cell line were not available for specific genes, data from other cell lines were used as substitutes. For details, see Additional file 1: Supplementary Table 2. For each transcription factor, the absolute value of the attribution score was used as the score for each gene, and the ChIP-seq peak was used as the ground truth to calculate the AUPRC ratio. To minimize the influence of surrounding genes on the genome, only genes with ChIP-seq peaks located within ± 1000 bp of the TSS were used for AUPRC calculations. Motif discovery Motif detection was performed via tfmodisco-lite (v2.2.0) 22 . The parameters were set as follows: sliding_window_size = 15, flank_size = 5, target_seqlet_fdr = 0.2, trim_to_window_size = 6, initial_flank_to_add = 2, final_flank_to_add = 5, final_min_cluster_size = 20, and subcluster_perplexity = 10. For the top 150 genes by absolute fold change of the target perturbation, we calculated Input X Gradient, computed the average attribution score for each gene every 128 bp, converted these scores to z scores, and defined bins with z scores greater than 1 or less than -1 as peak bins. We collected the attribution scores for each base in these peak bins across all 150 genes to use as seqlets for TF-MoDISco input. The hypothetical scores required for TF-MoDISco calculations were obtained via the saliency function in Captum, with the attribution scores normalized to ensure a mean of zero. We used the values from the same bins as the seqlets. The detected motif patterns were compared via Tomtom via the MEME Suite (v5.3.0) 40 . Known motifs were obtained from the CORE vertebrate nonredundant PFM data of JASPAR 2024 23 in MEME format. To address cases where motifs are similar, such as those from gene families, we also used motif cluster data from JASPAR 2024 CORE vertebrates to check for membership in the same cluster. Owing to the strong dependence of q-values on motif length, we evaluated motif matches via two thresholds: when the q-value was less than 0.4 and when the motif was ranked in the top 10 according to the q-value. Input mask For the input mask, we performed transfer learning on Enformer via DNA sequences where all bases outside ± 5,000 bp (total 10,000 bp) of the TSS were replaced with N. In addition to replacing bases with N, the model structure and training method were the same as those in the full-length condition. Declarations Code availability The analysis code has been stored in the (https://github.com/rikenbit/GenPerturb) repository. Corresponding author Correspondence to: Itoshi Nikaido [email protected] Ethics declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Competing interests The authors report no competing interests. Author Contribution T.S. and I.N. designed and configured the various approaches used in this study. T.S. performed the training of transfer learning. T.S. analyzed the data. T.S.. and I.N. prepared the figures and wrote the manuscript. All authors read and approved the final manuscript. Acknowledgement We thank the members of the Bioinformatics Research Team, particularly Akihiro Matsushima, for their management of the IT infrastructure.We thank the members of the Single-cell Omics Laboratory for their discussion of single-cell RNA-sequencing technologies. Data Availability The data used in this study is based on publicly available sources, as outlined below.The perturbed scRNA-seq data were downloaded from scPerturb (http://projects.sanderlab.org/scperturb/) for the datasets by Norman et al.14 and Replogle et al.15, from the sc-pert repository (https://github.com/theislab/sc-pert) for the dataset by Srivatsan et al.10, and from the original publication (GSE217460, https://data.caltech.edu/records/2cjss-wgh69) for the remaining data. Each pretrained model was downloaded from Hugging Face (https://huggingface.co/). ATAC-seq and ChIP-seq data were downloaded from ChIP-Atlas (https://chip-atlas.org/) (Additional file 1: Supplementary Table 2). Known motif sequences were downloaded from JASPAR (https://jaspar.elixir.no/). References Kalyan KS, Rajasekharan A, Sangeetha S. Ammus: a survey of transformer-based pretrained models in natural language processing. arXiv preprint arXiv:210805542 (2021). Zhao WX et al. A survey of large language models. arXiv preprint arXiv:230318223 (2023). Avsec Ž, et al. Effective gene expression prediction from sequence by integrating long-range interactions. Nat Methods. 2021;18:1196–203. Nguyen E et al. HyenaDNA: long-range genomic sequence modeling at single nucleotide resolution. arXiv; 2023. arXiv preprint arXiv:2306.15794 (2023). Dalla-Torre H et al. Nucleotide Transformer: building and evaluating robust foundation models for human genomics. Nat Methods 1–11 (2024). Zhou Z et al. Dnabert-2: efficient foundation model and benchmark for multi-species genome. arXiv preprint arXiv:230615006 (2023). Simon E, Swanson K, Zou J. Language models for biological research: a primer. Nat Methods. 2024;21:1422–9. Li Z, et al. Applications of deep learning in understanding gene regulation. Cell Rep Methods. 2023;3:100384. Dixit A, et al. Perturb-seq: dissecting molecular circuits with scalable single-cell RNA profiling of pooled genetic screens. Cell. 2016;167:1853–e186617. Srivatsan SR, et al. Massively multiplex chemical transcriptomics at single-cell resolution. Science. 2020;367:45–51. Badia IMP, et al. Gene regulatory network inference in the era of single-cell multi-omics. Nat Rev Genet. 2023;24:739–54. Vaswani A. Attention is all you need. arXiv preprint arXiv:1706.03762 (2017). Hu EJ et al. Lora: low-rank adaptation of large language models. arXiv preprint arXiv:210609685 (2021). Norman TM, et al. Exploring genetic interaction manifolds constructed from rich single-cell phenotypes. Science. 2019;365:786–93. Replogle JM, et al. Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq. Cell. 2022;185:2559–e257528. Joung J, et al. A transcription factor atlas of directed differentiation. Cell. 2023;186:209–e22926. Jiang J et al. D-SPIN constructs gene regulatory network models from multiplexed scRNA-seq data revealing organizing principles of cellular perturbation response. bioRxiv 2023.04.19.537364 (2024). Müller-Dott S, et al. Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities. Nucleic Acids Res. 2023;51:10934–49. Hounkpe BW, Chenou F, De Lima F, De Paula EV. HRT Atlas v1.0 database: redefining human and mouse housekeeping genes and candidate reference transcripts by mining massive RNA-seq datasets. Nucleic Acids Res. 2021;49:D947–55. Sundararajan M, Taly A, Yan Q. Axiomatic attribution for deep networks. arXiv preprint arXiv:170301365 (2017). Novakovsky G, Dexter N, Libbrecht MW, Wasserman WW, Mostafavi S. Obtaining genetics insights from deep learning via explainable artificial intelligence. Nat Rev Genet. 2023;24:125–37. Shrikumar A et al. Technical note on transcription factor motif discovery from importance scores (TF-MoDISco) version 0.5. 6.5. arXiv preprint arXiv:1811.00416 (2018). Rauluseviciute I et al. JASPAR. 2024: 20th anniversary of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 52, D174–D182 (2024). Gasperini M, Tome JM, Shendure J. Towards a comprehensive catalogue of validated and target-linked human enhancers. Nat Rev Genet. 2020;21:292–310. Kaplan J et al. Scaling laws for neural language models. arXiv preprint arXiv:2001.08361 (2020). Peidli S, et al. scPerturb: harmonized single-cell perturbation data. Nat Methods. 2024;21:531–40. Papalexi E, et al. Characterizing the molecular regulation of inhibitory immune checkpoints with multimodal single-cell screens. Nat Genet. 2021;53:322–31. Heumos L et al. Pertpy: an end-to-end framework for perturbation analysis. bioRxiv 2024.08. 04.606516 (2024). Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19:15. Frankish A, et al. GENCODE: reference annotation for the human and mouse genomes in 2023. Nucleic Acids Res. 2023;51:D942–9. Wolf T. Huggingface's transformers: state-of-the-art natural language processing. arXiv preprint arXiv:1910.03771 (2019). Devlin J. Bert: pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:181004805 (2018). Paszke A et al. Pytorch: an imperative style, high-performance deep learning library. arXiv preprint arXiv:1912.01703 (2019). Badia IMP, et al. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinform Adv. 2022;2:vbac016. Kokhlikyan N et al. Captum: a unified and generic model interpretability library for pytorch. arXiv preprint arXiv:2009.07896 (2020). Dhamdhere K, Sundararajan M, Yan Q. How important is a neuron? arXiv preprint arXiv:1805.12233 (2018). Oki S, et al. ChIP-Atlas: a data-mining suite powered by full integration of public ChIP-seq data. EMBO Rep. 2018;19:e46255. Xu W, et al. CoolBox: a flexible toolkit for visual analysis of genomics data. BMC Bioinform. 2021;22:489. Ramírez F, et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 2016;44:W160–5. Bailey TL, Johnson J, Grant CE, Noble WS. The MEME suite. Nucleic Acids Res. 2015;43:W39–49. Additional Declarations No competing interests reported. Supplementary Files Additionalfile1.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6461531","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":449587554,"identity":"0c4b99ea-b288-496e-82d0-37773b77103f","order_by":0,"name":"Takuya Shihashi","email":"","orcid":"","institution":"TRIP Headquarters, RIKEN","correspondingAuthor":false,"prefix":"","firstName":"Takuya","middleName":"","lastName":"Shihashi","suffix":""},{"id":449587555,"identity":"aa752482-a891-4a30-a0ea-05cec7921e9f","order_by":1,"name":"Itoshi Nikaido","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA9ElEQVRIiWNgGAWjYBACNiBmZjAAMngYGA7wMNiAGQwMPHi0sKFqSSOshQGsBaaIh+EwftUgwCfffPhzQQFDHh/P4YcH3lSczwcyjj1gkLmDx2FsadIzDBiK2XjbDA7OOXPbso23Ld2AgecZHi08Zsw8BgyJbfwMBod5224bsPHzmEkw8BzGp8X4M0QL+weglnNALfzfCGkxkAZr4e0B2XLAgI23h42AlrQ0oBaJxDaeMwVAvyQbsPEcM5NIwOMXYIAd/szzxyZxfk/65g9vKuwM5HuSn0l87MEdYlAggcZP7DlASAsG+EG6llEwCkbBKBi2AABjJ0ZKOtCczgAAAABJRU5ErkJggg==","orcid":"","institution":"TRIP Headquarters, RIKEN","correspondingAuthor":true,"prefix":"","firstName":"Itoshi","middleName":"","lastName":"Nikaido","suffix":""}],"badges":[],"createdAt":"2025-04-16 08:53:24","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-6461531/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6461531/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":81924077,"identity":"ffb7f50c-9d72-4cd6-a1b6-b900d9de3bda","added_by":"auto","created_at":"2025-05-05 02:32:32","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":204533,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSchematic overview of the study\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003ea, Transfer learning training methods for Enformer, HyenaDNA, and Nucleotide Transformer. b, Predicted expression matrix using 12 datasets for transfer learning, including perturbation data from CRISPRa, CRISPRi, compounds, and CDS overexpression. c, Evaluation of prediction accuracy of transfer learning models via correlation coefficients and embeddings. d, Gene signature profiles of biological phenotypes and gene regulatory networks. e, Understanding the features during model predictions via attribution scores. f, Challenges in predicting the accuracy of transfer learning models.\u003c/p\u003e","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/0d6a3f74ca01ba0a320c2cc2.png"},{"id":81924417,"identity":"2745c7a7-052d-43e9-9e4f-ba80127e6688","added_by":"auto","created_at":"2025-05-05 02:40:32","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":513278,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of correlation coefficients among transfer learning models\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003ea, Correlations between predicted and observed values across genes after transfer learning. b. Scatter plots of observed values and values predicted by the Enformer transfer learning model for the nonperturbation group of the K562 cell line from Norman et al.\u003csup\u003e14\u003c/sup\u003e (left) and values predicted by the original Enformer model for CAGE-seq in the K562 cell line (right). c, Correlation between predicted and observed values across perturbations after transfer learning. For each dataset, the genes were divided into five groups on the basis of expression level: very high, high, middle, low, and very low, and the correlation coefficients for each group were plotted. The blue box plot represents the Enformer model, the orange box plot represents HyenaDNA model, and the green box plot represents the Nucleotide Transformer model.\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/de6f8dab4142c704918c5bd7.png"},{"id":81924078,"identity":"fcfa3194-9db1-4b8c-b0c6-24768292c679","added_by":"auto","created_at":"2025-05-05 02:32:32","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":863862,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of embeddings among transfer learning models\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003ea-d, Uniform manifold approximation and projection (UMAP) of observed values and predicted values for each model in Norman et al.\u003csup\u003e14\u003c/sup\u003e, showing observed (a), Enformer (b), HyenaDNA (c), and Nucleotide Transformer (d) values. The labels from the Leiden clustering of the observed values are coloured. e, Consistency evaluation\u0026nbsp; between Leiden clustering for observed and predicted values via the adjusted Rand index (ARI) and normalized mutual information (NMI). f-i, UMAPs of observed and predicted values for each dataset in Jiang et al.\u003csup\u003e17\u003c/sup\u003e CD8T (f), CD4T (g), and Srivatsan et al.\u003csup\u003e10\u003c/sup\u003e A549 (h), and K562 (i) populations. The labels from the Leiden clustering of the observed values are coloured.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/b3c163435b7d383499ac823b.png"},{"id":81924092,"identity":"673d295e-a77e-47c3-9d17-2096071cccf9","added_by":"auto","created_at":"2025-05-05 02:32:53","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":786981,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eGene signatures associated with cell differentiation and gene regulatory networks\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003ea, UMAP plot labelled with the perturbation group obtained from Norman et al.\u003csup\u003e14\u003c/sup\u003e The left plot shows the observed values, and the right plot shows the values predicted by the Enformer model. b, UMAP plot coloured according to differentiation marker gene signature scores. The upper row shows erythroid scores, and the lower row shows granulocyte scores. The left plot shows the observed values, and the right plot shows the predicted values. c, Ranking of relative TF activity score for matched transcription factor perturbation groups. d, UMAP plot coloured according to the TF activity scores of HNF4A (upper row) and IRF1 (lower row). The left plot shows the observed values, and the right plot shows the predicted values.\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/312a7c9b17d36360450aae22.png"},{"id":81924422,"identity":"6bce26fb-1a4a-4b03-b27f-a0c1d7c04a6d","added_by":"auto","created_at":"2025-05-05 02:40:32","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":583605,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDifferences in target genome regions during prediction and motif discovery\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003ea Genome track centred on the TSS of the TYROBP gene. Signals of attribution scores during gene expression changes due to CEBPA perturbation in the K562 cell line for each transfer learning model trained by Norman et al.\u003csup\u003e14\u003c/sup\u003e data are shown, along with ChIP-seq of histone modifications and ATAC-seq in the same K562 cell line. b, Enrichment of absolute attribution scores due to CEBPA perturbation from Norman et al.\u003csup\u003e14\u003c/sup\u003e data at ATAC-seq peaks of all test genes. c, Classification performance of ChIP-seq peaks for matched transcription factors on the basis of attribution scores of transcription factor perturbation predictions, calculated via the AUPRC ratio. d, Detection performance of transcription factor motifs on the basis of the attribution score by TF-MoDISco. Match rates between motifs obtained from perturbation attribution scores and transcription factor motifs registered in JASPAR (q value \u0026lt; 0.4 or within the top 10, with exact gene name matches or matches with genes in the motif cluster). e, Representative plot of motifs obtained by TF-MoDISco with matched transcription factor motifs registered in JASPAR.\u003c/p\u003e","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/fbb04b93fc9525f9dbce5d18.png"},{"id":81924419,"identity":"4825cbe3-2fe0-4e88-ba00-a52567b47ca0","added_by":"auto","created_at":"2025-05-05 02:40:32","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":1094628,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCorrelations of the predicted values due to overlap of the input DNA sequences\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003ea, Genome track centred on the TSS of the NFKBIB gene. Attribution score signals for NFKBIB and neighbouring genes from the Enformer transfer learning model, and ChIP-seq and ATAC-seq signals were plotted. b-c, Pairwise correlation coefficients of observed values (b) and correlation coefficients between observed and predicted values (c) among neighbouring genes of NFKBIB. d, Distribution of correlation coefficients divided by the distance between TSSs when calculating pairwise correlation coefficients of predicted values between all test genes. e-f, Scatter plot of the observed correlation between each gene and its neighbouring genes located within 200 kbp of the TSS in the genome (x-axis) and the prediction accuracy of the transfer learning model for each gene (y-axis). Dots represent individual genes. For the x-axis, (e) represents the correlation coefficient with the most negatively correlated neighbouring gene, and (f) represents the correlation coefficient with the most positively correlated neighbouring gene.\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/c05bd79c2bddf98d2ac07c91.png"},{"id":92249201,"identity":"5454bb45-c23a-4190-86f5-12927c3c4f6a","added_by":"auto","created_at":"2025-09-26 10:17:17","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":5150896,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/060c4c0a-a5bd-453d-bd82-293e058c7166.pdf"},{"id":81924090,"identity":"7abbef4c-5033-492c-8393-2fa7b09de97f","added_by":"auto","created_at":"2025-05-05 02:32:32","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":3196649,"visible":true,"origin":"","legend":"","description":"","filename":"Additionalfile1.docx","url":"https://assets-eu.researchsquare.com/files/rs-6461531/v1/58f2916738baf32bd3d6e5db.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Transfer learning using generative pretrained genomic DNA models for predicting perturbation-induced changes in gene expression","fulltext":[{"header":"Background","content":"\u003cp\u003eIn natural language processing (NLP), fine-tuning pretrained models such as large language models (LLMs) enables high performance on specific tasks, even with small datasets and limited computer resources\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. There has been some success in applying deep learning architectures used in NLP to genomic DNA models\u003csup\u003e\u003cspan additionalcitationids=\"CR4 CR5\" citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e, expecting that fine-tuning with small datasets can be also effective in genomic DNA models\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e,\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e. This approach is beneficial when advancements in measurement technology yield new types of data, where it is challenging to prepare measurement data with sufficient quantity and diversity. However, research on transfer learning in this area is currently insufficient, and many uncertainties remain regarding its practical utility.\u003c/p\u003e \u003cp\u003eRecently, advancements in technologies such as Perturb-seq and SciPlex have made it possible to comprehensively profile gene expression following external perturbations, such as gene overexpression or drug treatment\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e,\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e. The development of models capable of predicting perturbation-induced gene expression changes from DNA sequences can help elucidate the relationships between perturbations and cis-regulatory elements (CREs), deepening our understanding of transcriptional regulation and contributing to the construction of gene regulatory networks (GRNs)\u003csup\u003e\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. Moreover, combining this approach with DNA mutation simulations could reveal the relationship between noncoding DNA mutations and drug responses, with potential applications in personalized medicine. Given the current limitations in the numbers of cell types and perturbations available in the data, the use of transfer learning is a practical approach for developing predictive models.\u003c/p\u003e \u003cp\u003eAt present, Enformer is a representative genomic DNA model. This model combines a Transformer\u003csup\u003e\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u003c/sup\u003e with a convolutional neural network (CNN) to predict gene expression and epigenetic marker tracks from DNA sequences\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. Other examples include genome language models (gLMs), such as HyenaDNA and Nucleotide Transformer, which predict enhancer element activity and splicing sites from DNA sequences\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e. These models have rarely been trained with perturbation-induced gene expression data, and there has yet to be an exploration of building predictive models for perturbation data via transfer learning.\u003c/p\u003e \u003cp\u003eIn this study, we benchmarked transfer learning in genomic DNA models, demonstrating that it is possible to predict perturbation-induced gene expression changes from DNA sequences. We evaluated multiple conditions across three genomic DNA models and 12 datasets, revealing that transfer learning via Enformer can effectively model these changes. By employing deep learning model interpretation methods, we confirmed that transfer learning can focus on regulatory elements and transcription factor motifs on the basis of biological principles. Moreover, we discovered that when the input DNA sequences have regions of overlap, the predicted gene expression values often strongly correlate with each other but deviate from the actual observed values. We referred to this issue as 'genomic neighbouring gene interference'.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eDesign of transfer learning and overview of the evaluation workflow\u003c/h2\u003e \u003cp\u003eThe overall structure of this study is illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e. We used Enformer, HyenaDNA, and Nucleotide Transformer as the genomic DNA models. The parameters of the additional heads were trained using the embedding values from each model as input. During this process, we considered three approaches: (1) freezing the parameters of the pretrained models (feature-based), (2) fine-tuning a subset of the parameters via low-rank adaptation (LoRA)\u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e, or (3) fine-tuning all of the parameters. Enformer uses the embedding values of 512 bp around the transcription start site (TSS) for the input values to the heads. HyenaDNA also uses the embedding values of 512 bp around the TSS. Since the Nucleotide Transformer, tokenizes sequences in 6 bp units, we used embedding values of 516 bp. In NLP classification tasks, the last token is often used for next-token prediction models, and the CLS token is used for masked token prediction models\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e. Therefore, for the HyenaDNA and Nucleotide Transformer, we also conducted training using these embedding values as inputs (Additional file 1: Supplementary Fig.\u0026nbsp;1a).\u003c/p\u003e \u003cp\u003eWe prepared 12 datasets and conducted training with each dataset\u003csup\u003e\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e,\u003cspan additionalcitationids=\"CR15 CR16\" citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e. These datasets include genome-wide screens with perturbations such as clustered regularly interspaced short palindromic repeats (CRISPR) activation, CRISPR inhibition, compound treatments, and coding DNA sequence (CDS) overexpression. To mitigate the reduction in model accuracy caused by noise originating from measurement techniques, we converted single-cell data into pseudobulk data and used it for model training. Moreover, we excluded nonexpressed genes. The numbers of perturbations and genes in each dataset are listed in Additional file 1: Supplementary Table\u0026nbsp;1.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eSelection of transfer learning methods for benchmarking\u003c/h3\u003e\n\u003cp\u003eThere are multiple methods for transfer learning\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e. We examined the differences in accuracy among three learning methods: (1) Obtaining the embedding values in advance and training only the additional head (feature-based). The pretrained model combined with the head or the head with precomputed embedding values is used for prediction. (2) LoRA is used to train parts of the transformer layers along with the additional head\u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e. The number of parameters updated is set by the rank (in this study, 16, 64, 256, and 512). (3) All parameters of both the pretrained model and the additional head are trained. The learning rate for the pretrained model was set either equal to or lower than that of the additional head. In the evaluation using the CD8T dataset by Jiang et al.\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e, the feature-based approach with precomputed embedding values achieved the highest accuracy in terms of correlations across genes and perturbations (Additional file 1: Supplementary Fig.\u0026nbsp;2a and b). With respect to correlations across perturbations, since genes with higher expression levels led to greater accuracy, we divided the genes into five groups on the basis of expression levels and examined the distribution of correlation coefficients. To verify the accuracy differences among training methods in correlations across genes, we compared the observed and predicted values in nonperturbed samples. Compared with the feature-based method, fine-tuning caused the model's predictions to become more uniform across different genes, leading to a decrease in correlation coefficients (Additional file 1: Supplementary Fig.\u0026nbsp;2c). Similar results were observed with HyenaDNA and the Nucleotide Transformer, where training only the feature-based approach using precomputed embedding values yielded the highest accuracy. In some cases, catastrophic forgetting appeared to occur (Additional file 1: Supplementary Fig.\u0026nbsp;2d-g). On the basis of these results, we conducted further investigations using the feature-based approach for the transfer learning method.\u003c/p\u003e\n\u003ch3\u003ePerformance of transfer learning models for predicting perturbation-induced gene expression changes\u003c/h3\u003e\n\u003cp\u003eWe applied transfer learning using three genomic DNA models, Enformer, HyenaDNA, and Nucleotide Transformer, across 12 datasets to evaluate the correlations across genes and perturbations. The Enformer model demonstrated higher accuracy, with a median gene correlation of approximately 0.6, outperforming HyenaDNA and the Nucleotide Transformer (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea and Additional file 1: Supplementary Fig.\u0026nbsp;3a). When the original Enformer model was compared with the transfer learning model in nonperturbed K562 cells, it was suggested that transfer learning improved gene correlations and better adapted the model to the additional data (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb). With respect to perturbation correlations, the Enformer model exhibited a median correlation between 0.2 and 0.5 for high-expression genes, whereas near-zero correlations were observed for low-expression genes, likely due to the limitations of measurement sensitivity (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec). The accuracy of the HyenaDNA and Nucleotide Transformer methods was generally lower (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec and Additional file 1: Supplementary Fig.\u0026nbsp;3b).\u003c/p\u003e \u003cp\u003eTo evaluate the differences among perturbations, uniform manifold approximation and projection (UMAP) embedding and Leiden clustering were performed via expression values. Enformer generated embeddings that were consistent with the observed data, whereas HyenaDNA and Nucleotide Transformer generated embeddings with inconsistencies (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea-d and Additional file 1: Supplementary Fig.\u0026nbsp;3c and d). Moreover, Enformer outperformed the other methods in terms of clustering concordance, as indicated by its higher adjusted Rand index (ARI) and normalized mutual information (NMI) scores (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ee and Additional file 1: Supplementary Fig.\u0026nbsp;3e). The UMAP of observed values from the dataset generated by Jiang et al.\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e revealed good separation, and the model trained with this dataset achieved high accuracy (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec and Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ee-g). In contrast, the UMAP of observed values from the dataset generated by Srivatsan et al.\u003csup\u003e\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e indicated poor separation, and the model trained with this dataset demonstrated low accuracy (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec and Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ee,h-i). In summary, the Enformer transfer learning model achieved the highest accuracy. The accuracy depends on the characteristics of the dataset used for training.\u003c/p\u003e\n\u003ch3\u003eTransfer learning accurately predicts biological phenotypes and gene regulatory networks\u003c/h3\u003e\n\u003cp\u003eWe evaluated whether the Enformer transfer learning model could capture biological features. Since the number of genes that can be assessed with test data alone is limited, we evaluated the gene signatures via the prediction results of all the genes. As there was no apparent difference in the embedding results when predictions of all genes were used compared with when only the test data were used, we proceeded with this approach (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea-e and Additional file 1: Supplementary Fig.\u0026nbsp;4a-e). The Norman et al.\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e dataset includes biological function labels for each perturbation\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e, and by examining the distribution of these labels on UMAP plots of the observed and predicted values, we confirmed that perturbations with the same label were clustered together (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea and Additional file 1: Supplementary Fig.\u0026nbsp;4f). In terms of cell differentiation, the gene signature scores for erythroid and granulocyte lineages were consistent with their labels (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea and b).\u003c/p\u003e \u003cp\u003eNext, we calculated activation scores for each transcription factor via gene signatures from CollecTRI\u003csup\u003e\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e, which summarizes transcriptional regulatory networks on the basis of ChIP-seq data, and ranked the relative scores for each perturbation. Consistently high activation scores were observed for transcription factors such as HNF4A, IRF1, and TP73 in both the observed and predicted values when they were overexpressed via CRISPRa (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ec and d). However, the predicted rank for transcription factors such as AHR was significantly lower, suggesting that Enformer may not have adequately learned AHR-mediated regulation during pretraining. Finally, no apparent differences in correlation were detected between the housekeeping genes and the other genes, suggesting that the model's accuracy does not depend on the housekeeping genes (Additional file 1: Supplementary Fig.\u0026nbsp;4g)\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. Overall, these findings suggest that the Enformer transfer learning model accurately predicts cell differentiation and gene regulatory networks.\u003c/p\u003e\n\u003ch3\u003eTransfer learning model using Enformer attention to transcriptional regulatory regions and transcription factor motifs\u003c/h3\u003e\n\u003cp\u003eThe Enformer transfer learning model accurately predicted gene expression changes due to perturbations and outperformed HyenaDNA and Nucleotide Transformer. To explore the differences in prediction accuracy, we identified the DNA sequences that the models focused on during predictions via attribution scores derived from Input X Gradient\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e,\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e,\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e. By subtracting the attribution scores of the control group from those of the perturbation group, we evaluated the difference in attribution scores and identified the critical sequences necessary for predicting expression changes induced by perturbations. Since we confirmed that the fold change in the predicted values correlated with the subtracted attribution score differences in certain genomic regions, we proceeded with this approach. (Additional file 1: Supplementary Fig.\u0026nbsp;5). We employed the Norman et al.\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e dataset to visualize the genomic regions highlighted by attribution scores during the prediction of TYROBP expression when CEBPA was overexpressed via CRISPRa, overlaying these regions with ATAC-seq and ChIP-seq data from the same K562 cell line (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea). Enformer focused on regulatory elements that activate transcription, such as ATAC-seq H3K27ac and H3K4me3. In contrast, HyenaDNA, when TSS-centred embeddings are used, focuses only on the region from the start of the input to the TSS, and when the final token is used, it fails to focus on regulatory elements. The Nucleotide Transformer, which was constrained by the short input sequence, focused only on the region surrounding the TSS. To verify whether these models focused on regulatory elements, we checked the enrichment of attribution scores in ATAC-seq peak regions. Enformer, as well as HyenaDNA and Nucleotide Transformer, when using TSS-centred embeddings, focused on ATAC-seq peaks, whereas models using the final token or CLS token did not focus on ATAC-seq peaks (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eb).\u003c/p\u003e \u003cp\u003eWe further evaluated the Enformer transfer learning model. When we assessed whether the model focused on the ChIP-seq peaks of the same transcription factors during gene expression changes caused by transcription factor overexpression, the area under the precision\u0026ndash;recall curve (AUPRC) ratio ranged from 3 to 6. This finding suggests that the model predicts a focus on the ChIP-seq peaks (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ec). Furthermore, via TF-MoDISco\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e, we detected transcription factor motifs from the attribution scores for all the transcription factor perturbations included in the Norman et al.\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e dataset. Between 20% and 30% of the transcription factors had motifs detected by TF-MoDISco that matched the known motifs registered in JASPAR\u003csup\u003e\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e, such as HNF4A, IRF1, CEBPA, and SPI1. Additionally, 70% of the transcription factors, such as those belonging to the same gene family, matched one of the motifs within the cluster (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed and e). In summary, the Enformer model achieved high predictive accuracy by focusing on CREs during prediction and recognizing transcription factor motifs at the nucleotide level. On the other hand, HyenaDNA and Nucleotide Transformer struggled to capture the sequences involved in transcriptional regulation compared with Enformer.\u003c/p\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eChallenges in predicting gene expression changes from DNA sequences: Genomic neighbouring gene interference problem\u003c/h2\u003e \u003cp\u003eIdentifying challenges in model predictions can guide algorithm development. In Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec, all transfer learning models included some genes showing negative correlations between the observed and predicted values, prompting an investigation into the causes. Examining the distribution of attribution scores for genes neighbouring NFKBIB (r = -0.61) and MRPS12 (r = -0.46), which had the strongest negative correlations in Norman et al.'s\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e test data, revealed that neighbouring genes focused on each other\u0026rsquo;s TSS regions for their prediction (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea).\u003c/p\u003e \u003cp\u003eThe observed values of the neighbouring gene RINL were negatively correlated with both NFKBIB (r = -0.44) and MRPS12 (r = -0.63) (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb). However, in terms of the predicted values, RINL was positively correlated with NFKBIB (r\u0026thinsp;=\u0026thinsp;0.62) and MRPS12 (r\u0026thinsp;=\u0026thinsp;0.55), and other neighbouring genes also learned to correlate with RINL\u0026rsquo;s observed values (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec). These findings suggest that interference during training involving neighbouring genes could cause negative correlations between the observed and predicted values for NFKBIB and MRPS12.\u003c/p\u003e \u003cp\u003eBeyond the NFKBIB region, we calculated pairwise correlation coefficients for all genes, grouping them by the distance between their TSS regions. The results indicated that neighbouring genes within ten kbp showed a strongly positive correlation with the predicted values (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed). We also investigated whether genes with neighbouring genes that had negatively correlated observed values displayed negative correlations in predictions. As expected, most genes with negatively correlated predictions had neighbouring genes within 200 kbp that were negatively correlated with the observed values. In these cases, one gene showed a negative correlation in the predictions, whereas the other was correctly predicted to be positively correlated (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ee). Additionally, when neighbouring genes showed a strong positive correlation with the observed values, the predicted values tended to show similarly strong positive correlations (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ef).\u003c/p\u003e \u003cp\u003eTo test whether this issue improves when the input sequences of neighbouring genes do not overlap, we performed transfer learning via masked input sequences, excluding regions outside the 10 kbp surrounding the TSS. Masking reduced the positive correlation between genes with TSS distances of 5\u0026ndash;10 kbp, likely due to a reduction in the overlap of input DNA sequences to 50% or less. (Additional file 1: Supplementary Fig.\u0026nbsp;6a and b). As a supplement, masking the input sequences partially improved negative correlation of NFKBIB and MRPS12, and decreased the accuracy of correlations across genes and perturbations, so it did not address the root cause of the problem (Additional file 1: Supplementary Fig.\u0026nbsp;6c-e).\u003c/p\u003e \u003cp\u003eIn summary, the findings suggest that overlapping DNA sequences between neighbouring genes cause interference, affecting prediction accuracy both positively and negatively because of the observed correlations among these genes. We termed this phenomenon \u0026lsquo;genomic neighbouring gene interference\u0026rsquo;. Further improvements in model structure and training methods are needed to address these issues.\u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eIn this study, we evaluated three genomic DNA models, 12 datasets, and multiple transfer learning conditions and demonstrated that transfer learning via Enformer is effective for predicting gene expression changes due to perturbations. In contrast, HyenaDNA and Nucleotide Transformer did not achieve sufficient accuracy.\u003c/p\u003e \u003cp\u003eThe datasets used in this study involved external perturbations, yet the Enformer model had not been pretrained on such data. Nevertheless, the Enformer model was able to predict gene expression changes, suggesting that it implicitly learned the principles of biological transcriptional regulation, particularly GRNs, from large-scale gene expression and epigenetic data. These GRNs drive gene expression changes in response to perturbations, allowing Enformer to make accurate predictions during transfer learning without becoming out-of-distribution. However, we identified a problem termed \"genomic neighbouring gene interference,\" where the prediction results revealed correlations between neighbouring genes. One possible cause of this problem is that genomic features, such as CREs, are not properly associated with each individual gene\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e\u003c/sup\u003e. To address this issue, increasing the amount of data used in pretraining, such as ChIP-seq and ATAC-seq, and the number of model parameters may naturally improve the association between genes and CREs, similar to the scaling laws observed in NLP\u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e. Additionally, incorporating 3D genome structure data or explicitly including meta-information, such as cell type, in the input during pretraining could also enhance the model's performance.\u003c/p\u003e \u003cp\u003eHyenaDNA and Nucleotide Transformer did not achieve satisfactory accuracy through transfer learning. While these models have demonstrated high accuracy in DNA sequence classification tasks that do not depend on cell type, such as predicting splicing sites and promoters. The results of this study suggest that these models are unable to learn cell type- and context-dependent information, such as transcriptional regulation. Therefore, it will be necessary to develop learning methods and datasets that capture the principles of transcriptional regulation on the basis of cell type and context\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e. Also, increasing the length of DNA sequences alone, which is the basic future task for gLM, may not improve accuracy for gene expression prediction.\u003c/p\u003e \u003cp\u003eThe feature-based method that achieved the highest accuracy in transfer learning with the Enformer model utilized less than 1 GB of GPU memory on the NVIDIA RTX 6000 Ada, completing training in a maximum of 10 minutes. In contrast, fine-tuning, including the pretrained model, requires 20\u0026ndash;45 GB of GPU memory per batch and requires 10\u0026ndash;24 hours of training on a single GPU for 40 epochs. In both cases, this approach significantly reduces the computational resources compared with that of the Enformer pretraining process, which requires 32 TPU v3 cores and three days for training\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. The findings demonstrate that transfer learning is beneficial in terms of computational efficiency.\u003c/p\u003e \u003cp\u003eFinally, this study provides an example of transfer learning in genomic DNA models and offers valuable insights for individual researchers conducting transfer learning and for those developing genomic DNA models. This study will lead to an increase in the efficiency of the development and utilization of genomic DNA models.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eThis study demonstrated that transfer learning with the Enformer model effectively predicts perturbation-induced gene expression changes, significantly outperforming HyenaDNA and Nucleotide Transformer. Despite lacking explicit perturbation data in pretraining, Enformer's implicit grasp of gene regulatory networks enables accurate predictions. However, genomic neighboring gene interference revealed limitations in properly associating genomic regulatory elements with their target genes. These insights will accelerate the efficient development and application of genomic DNA models.\u003c/p\u003e "},{"header":"Methods","content":"\u003cp\u003e\u003cstrong\u003e\u003cem\u003eData preprocessing for perturbation data\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe used 12 datasets from five studies\u003csup\u003e10,14-17\u003c/sup\u003e. The data from Norman et al.\u003csup\u003e14\u003c/sup\u003e and Replogle et al.\u003csup\u003e15\u003c/sup\u003e were obtained from scPerturb (http://projects.sanderlab.org/scperturb/)\u003csup\u003e26\u003c/sup\u003e, and the data from Srivatsan et al.\u003csup\u003e10\u003c/sup\u003e were downloaded from the sc-pert repository (https://github.com/theislab/sc-pert), as the annotations in the scPerturb database were misaligned. The remaining data were downloaded from the original articles (GSE217460, https://data.caltech.edu/records/2cjss-wgh69). The dataset from Jiang et al.\u003csup\u003e17\u003c/sup\u003e was divided into CD4 T cells, CD8 T cells, B cells, and myeloid cells in the anti-CD3/CD28 treatment group. Owing to the low cell numbers, NK cells and certain other cell types or conditions in the anti-CD3/CD28 nontreatment group were not used. Because the datasets from Norman et al.\u003csup\u003e14\u003c/sup\u003e and Replogle et al.\u003csup\u003e15\u003c/sup\u003e involve CRISPR perturbations, we excluded unperturbed cells via Mixscape\u003csup\u003e27\u003c/sup\u003e, which is implemented in the pertpy library (v0.7.0)\u003csup\u003e28\u003c/sup\u003e, to clearly capture gene expression changes due to perturbations. For the Mixscape calculations, we used the highly variable genes calculated by the highly_variable_genes function in Scanpy (v1.10.2)\u003csup\u003e29\u003c/sup\u003e, with min_disp=0.2 for each perturbation. Each dataset was pseudobulked for each perturbation via adpbulk (v0.1.4) to remove noise and improve learning efficiency. Only perturbations with more than 100 cells were used for analysis. To prevent artificially induced gene expression profiles from affecting the learning of natural transcriptional regulatory mechanisms, we replaced the expression of target genes in the overexpression or knockdown perturbation groups with values from the nonperturbation group. The counts after pseudobulking were normalized to the CPM, converted to log2 + 1, and used for training. To exclude nonexpressed genes, we used only genes with a CPM exceeding 2 in at least one perturbation for training.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eModel architecture and input genome\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe model input consisted of genomic sequences centred on the gene's TSS. The genes were divided into training, validation, and test sets according to the classification set in the original Enformer paper, which was obtained from https://console.cloud.google.com/storage/browser/basenji_barnyard/data\u003csup\u003e3\u003c/sup\u003e. Only mRNAs and lncRNAs were used for training, which was based on GFF3 annotations. We used the same versions of GFF3 (v32) and FASTA (GRCh38.p13) from GENCODE\u003csup\u003e30\u003c/sup\u003e as in the original Enformer paper.\u003c/p\u003e\n\u003cp\u003eThe pretrained models used were Enformer\u003csup\u003e3\u003c/sup\u003e, HyenaDNA\u003csup\u003e4\u003c/sup\u003e, and Nucleotide Transformer\u003csup\u003e5\u003c/sup\u003e. These models were selected from those registered on Hugging Face\u003csup\u003e31\u003c/sup\u003e for their ability to handle long context lengths.\u003c/p\u003e\n\u003col\u003e\n\u003cli\u003eEnformer\u003c/li\u003e\n\u003c/ol\u003e\n\u003cp\u003eWe used the sequence centred on the TSS, spanning ± 98,304 bp (total 196,608 bp) as input and connected the embedding values before the final layer to an additional single-layer head. Since Enformer performs convolution operations every 128 bp and then proceeds with the transformer, we used the embedding value from the 4 bins (512 bp) surrounding the TSS as the embedding value.\u003c/p\u003e\n\u003col start=\"2\"\u003e\n\u003cli\u003eHyenaDNA (medium-160k)\u003c/li\u003e\n\u003c/ol\u003e\n\u003cp\u003eWe used the sequence centred on the TSS, spanning ± 80,000 bp (total 160,000 bp) as input, using either the embedding value of the 512 bp token centred on the TSS (TSS embedding) or, following the original paper, the last token as the input for the additional head\u003csup\u003e4\u003c/sup\u003e. The last token was treated as a single-layer head, as in Enformer, but the TSS embedding failed to learn because of too many input elements. Therefore, we used a 2-layer head with a hidden layer of 3,072 nodes for learning.\u003c/p\u003e\n\u003col start=\"3\"\u003e\n\u003cli\u003eNucleotide Transformer (2.5b-multispecies)\u003c/li\u003e\n\u003c/ol\u003e\n\u003cp\u003eWe used the sequence centred on the TSS, spanning ± 12,288 bp (total of 24,576 bp) as input. Since Nucleotide Transformer tokenizes every six bp, we used the embedding value of 516 bp (86 tokens) as the TSS embedding value for the additional head input. We also used the embedding value of the CLS token, which is used in classification models based on bidirectional encoder representations from transformers (BERT)\u003csup\u003e32\u003c/sup\u003e. Similar to the approach for HyenaDNA, the CLS token was learned with a single layer, whereas the TSS embedding value was learned with a 2-layer head with a hidden layer of 3,072 nodes.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eTraining methods for transfer learning\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eFor all the training methods, PyTorch (v2.3.1)\u003csup\u003e33\u003c/sup\u003e was used for training, with a learning rate of 0.0001, a gradient clip value of 0.2, and optimization with adaptive moment estimation with weight decay (AdamW), where the weight decay was set to 0.005. The loss function used was the mean square error (MSE) loss, and training was terminated when the validation loss did not improve for five epochs during training. The batch size was set to 256, and for methods that could not accommodate this batch size, the parameters were updated by accumulating gradients over 256 batches. The following five methods were implemented for training:\u003c/p\u003e\n\u003col\u003e\n\u003cli\u003ePretrained models were used to obtain embedding values in advance, and only the additional head was trained for 100 epochs. After training, predictions were made using only the head from the precomputed embedding values. Since there was no difference in accuracy compared with that of method (2) in terms of prediction from DNA sequences, the predictions from this method were used for evaluations, as shown in Additional file 1: Supplementary Fig. 2.\u003c/li\u003e\n\u003cli\u003ePretrained models were used to obtain embedding values in advance, and only the additional head was trained for 100 epochs. After training, predictions were made by combining the pretrained model and the head to predict from DNA sequences.\u003c/li\u003e\n\u003cli\u003eLoRA\u003csup\u003e13\u003c/sup\u003e was used to update the parameters of the query and value in all transformer layers, as well as the parameters for the additional head. The number of epochs was 20, with the LoRA parameters set to α=2, bias=none, and the rank was varied between 16, 64, 256 and 512 for training.\u003c/li\u003e\n\u003cli\u003eThe entire set of parameters for both the pretrained model and the additional head were trained for 40 epochs. During this process, the learning rate for the parameters of the pretrained model was varied between 0.0001 and 1×10⁻¹⁰.\u003c/li\u003e\n\u003c/ol\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eEvaluation of correlation coefficients\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe calculated two types of Pearson correlations for model evaluation. The first type is the correlation across perturbations, which is calculated between the actual and predicted values for each gene. The second type is the correlation across genes, which is calculated between the actual and predicted values for each perturbation after computing the predicted values for all genes. Both evaluations were performed using only the genes in the test set.\u003c/p\u003e\n\u003cp\u003eThe predicted values for nonperturbation conditions in the original Enformer were obtained by summing the predictions of the four bins (512 bp) surrounding the TSS of the CAGE-seq (K562: CNhs12336) track of the same cell type used in the transfer learning dataset, normalizing to CPM, and then converting to log2 + 1.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eEmbedding of gene expression\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe created an AnnData object to store the pseudobulked actual and predicted values as a matrix of each perturbation and gene. These values were processed via Scanpy. Then, we performed principal component analysis (PCA), k-nearest neighbours (kNN), Leiden clustering, and UMAP calculations with default parameters in Scanpy. For the evaluation of the test set, we used only the genes in the test set, and for evaluating all genes, we used all genes, including the training, validation and test sets. Clustering consistency was evaluated by performing Leiden clustering for both the observed and predicted values and then using scikit-learn (v1.5.1) to calculate the ARI and NMI.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eGene signature\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe gene signature scores related to differentiation were calculated via the gene sets for Myeloid and Erythroid lineages provided by Norman et al.\u003csup\u003e14\u003c/sup\u003e, employing the score_genes function in Scanpy. TF activity was obtained by retrieving information on transcription factors and downstream genes from CollecTRI\u003csup\u003e18\u003c/sup\u003e and calculating a univariate linear model (ULM) with a decoupler (v1.7.0)\u003csup\u003e34\u003c/sup\u003e. Human housekeeping gene lists were downloaded from the HRT Atlas v1.0\u003csup\u003e19\u003c/sup\u003e, and correlation coefficient box plots were drawn by dividing the genes into housekeeping and nonhousekeeping categories.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eAttribution score calculation\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe attribution score for Enformer was calculated via the Input X Gradient method for the target gene and each perturbation with the InputXGradient\u003csup\u003e20\u003c/sup\u003e function from Captum (v0.7.0)\u003csup\u003e35\u003c/sup\u003e. When focusing on changes due to perturbations, we subtracted the Input X Gradient value of the nonperturbation from that of the target perturbation for the same gene. For HyenaDNA and Nucleotide Transformer, which accept tokens as input, we calculated GradientXActivation\u003csup\u003e36\u003c/sup\u003e via Captum's function for tokens and summed the scores for each base to calculate the attribution score for each base (https://captum.ai/tutorials/IMDB_TorchText_Interpret).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eTrack plot\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe attribution score calculated by Input × Gradient was averaged every 128 bp and converted to a bigwig file via pyBigWig (v0.3.23). For the track in Fig. 6a, we calculated the Input X Gradient for all perturbations for a single gene and used the maximum absolute value of the attribution scores across all perturbations for each 128 bp bin to create the bigwig file. The ChIP-seq and ATAC-seq data were downloaded as bigwig files from ChIP-Atlas\u003csup\u003e37\u003c/sup\u003e. Tracks were visualized via CoolBox (v0.3.9)\u003csup\u003e38\u003c/sup\u003e. The gene locus information used for visualization with CoolBox was obtained from the GTF (v32) used during transfer learning. For peak enrichment, we downloaded the ATAC-seq peaks (q value \u0026lt; 1E-05) of the K562 cell line (SRX10184518) from ChIP-Atlas. We calculated the attribution scores for all test set genes in the group overexpressing CEBPA, converted them to absolute values, and then calculated the enrichment score via computeMatrix and plotProfile from deepTools2\u003csup\u003e39\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eAUPRC ratio for transcription factor ChIP-seq peaks\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eChIP-seq peaks for transcription factors were downloaded from ChIP-Atlas. When ChIP-seq data for the K562 cell line were not available for specific genes, data from other cell lines were used as substitutes. For details, see Additional file 1: Supplementary Table 2. For each transcription factor, the absolute value of the attribution score was used as the score for each gene, and the ChIP-seq peak was used as the ground truth to calculate the AUPRC ratio. To minimize the influence of surrounding genes on the genome, only genes with ChIP-seq peaks located within ± 1000 bp of the TSS were used for AUPRC calculations.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eMotif discovery\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eMotif detection was performed via tfmodisco-lite (v2.2.0)\u003csup\u003e22\u003c/sup\u003e. The parameters were set as follows: sliding_window_size = 15, flank_size = 5, target_seqlet_fdr = 0.2, trim_to_window_size = 6, initial_flank_to_add = 2, final_flank_to_add = 5, final_min_cluster_size = 20, and subcluster_perplexity = 10. For the top 150 genes by absolute fold change of the target perturbation, we calculated Input X Gradient, computed the average attribution score for each gene every 128 bp, converted these scores to z scores, and defined bins with z scores greater than 1 or less than -1 as peak bins. We collected the attribution scores for each base in these peak bins across all 150 genes to use as seqlets for TF-MoDISco input. The hypothetical scores required for TF-MoDISco calculations were obtained via the saliency function in Captum, with the attribution scores normalized to ensure a mean of zero. We used the values from the same bins as the seqlets.\u003c/p\u003e\n\u003cp\u003eThe detected motif patterns were compared via Tomtom via the MEME Suite (v5.3.0)\u003csup\u003e40\u003c/sup\u003e. Known motifs were obtained from the CORE vertebrate nonredundant PFM data of JASPAR 2024\u003csup\u003e23\u003c/sup\u003e in MEME format. To address cases where motifs are similar, such as those from gene families, we also used motif cluster data from JASPAR 2024 CORE vertebrates to check for membership in the same cluster. Owing to the strong dependence of q-values on motif length, we evaluated motif matches via two thresholds: when the q-value was less than 0.4 and when the motif was ranked in the top 10 according to the q-value.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003e\u003cem\u003eInput mask\u003c/em\u003e\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eFor the input mask, we performed transfer learning on Enformer via DNA sequences where all bases outside ± 5,000 bp (total 10,000 bp) of the TSS were replaced with N. In addition to replacing bases with N, the model structure and training method were the same as those in the full-length condition.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eCode availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe analysis code has been stored in the (https://github.com/rikenbit/GenPerturb) repository.\u003c/p\u003e\u003cp\u003e \u003cstrong\u003eCorresponding author\u003c/strong\u003e \u003cp\u003eCorrespondence to: Itoshi Nikaido [email protected]\u003c/p\u003e \u003c/p\u003e\u003ch2\u003e \u003cb\u003eEthics declarations\u003c/b\u003e \u003c/h2\u003e \u003cp\u003e \u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e \u003cp\u003eNot applicable.\u003c/p\u003e \u003c/p\u003e \u003cp\u003e \u003cstrong\u003eConsent for publication\u003c/strong\u003e \u003cp\u003eNot applicable.\u003c/p\u003e \u003c/p\u003e \u003cp\u003e \u003cstrong\u003eCompeting interests\u003c/strong\u003e \u003cp\u003eThe authors report no competing interests.\u003c/p\u003e \u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eT.S. and I.N. designed and configured the various approaches used in this study. T.S. performed the training of transfer learning. T.S. analyzed the data. T.S.. and I.N. prepared the figures and wrote the manuscript. All authors read and approved the final manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eWe thank the members of the Bioinformatics Research Team, particularly Akihiro Matsushima, for their management of the IT infrastructure.We thank the members of the Single-cell Omics Laboratory for their discussion of single-cell RNA-sequencing technologies.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eThe data used in this study is based on publicly available sources, as outlined below.The perturbed scRNA-seq data were downloaded from scPerturb (http://projects.sanderlab.org/scperturb/) for the datasets by Norman et al.14 and Replogle et al.15, from the sc-pert repository (https://github.com/theislab/sc-pert) for the dataset by Srivatsan et al.10, and from the original publication (GSE217460, https://data.caltech.edu/records/2cjss-wgh69) for the remaining data. Each pretrained model was downloaded from Hugging Face (https://huggingface.co/). ATAC-seq and ChIP-seq data were downloaded from ChIP-Atlas (https://chip-atlas.org/) (Additional file 1: Supplementary Table 2). Known motif sequences were downloaded from JASPAR (https://jaspar.elixir.no/).\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eKalyan KS, Rajasekharan A, Sangeetha S. Ammus: a survey of transformer-based pretrained models in natural language processing. arXiv preprint arXiv:210805542 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhao WX et al. A survey of large language models. arXiv preprint arXiv:230318223 (2023).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAvsec Ž, et al. Effective gene expression prediction from sequence by integrating long-range interactions. Nat Methods. 2021;18:1196\u0026ndash;203.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNguyen E et al. HyenaDNA: long-range genomic sequence modeling at single nucleotide resolution. arXiv; 2023. arXiv preprint arXiv:2306.15794 (2023).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDalla-Torre H et al. Nucleotide Transformer: building and evaluating robust foundation models for human genomics. Nat Methods 1\u0026ndash;11 (2024).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhou Z et al. Dnabert-2: efficient foundation model and benchmark for multi-species genome. arXiv preprint arXiv:230615006 (2023).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSimon E, Swanson K, Zou J. Language models for biological research: a primer. Nat Methods. 2024;21:1422\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi Z, et al. Applications of deep learning in understanding gene regulation. Cell Rep Methods. 2023;3:100384.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDixit A, et al. Perturb-seq: dissecting molecular circuits with scalable single-cell RNA profiling of pooled genetic screens. Cell. 2016;167:1853\u0026ndash;e186617.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSrivatsan SR, et al. Massively multiplex chemical transcriptomics at single-cell resolution. Science. 2020;367:45\u0026ndash;51.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBadia IMP, et al. Gene regulatory network inference in the era of single-cell multi-omics. Nat Rev Genet. 2023;24:739\u0026ndash;54.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVaswani A. Attention is all you need. arXiv preprint arXiv:1706.03762 (2017).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHu EJ et al. Lora: low-rank adaptation of large language models. arXiv preprint arXiv:210609685 (2021).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNorman TM, et al. Exploring genetic interaction manifolds constructed from rich single-cell phenotypes. Science. 2019;365:786\u0026ndash;93.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eReplogle JM, et al. Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq. Cell. 2022;185:2559\u0026ndash;e257528.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJoung J, et al. A transcription factor atlas of directed differentiation. Cell. 2023;186:209\u0026ndash;e22926.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJiang J et al. D-SPIN constructs gene regulatory network models from multiplexed scRNA-seq data revealing organizing principles of cellular perturbation response. bioRxiv 2023.04.19.537364 (2024).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eM\u0026uuml;ller-Dott S, et al. Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities. Nucleic Acids Res. 2023;51:10934\u0026ndash;49.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHounkpe BW, Chenou F, De Lima F, De Paula EV. HRT Atlas v1.0 database: redefining human and mouse housekeeping genes and candidate reference transcripts by mining massive RNA-seq datasets. Nucleic Acids Res. 2021;49:D947\u0026ndash;55.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSundararajan M, Taly A, Yan Q. Axiomatic attribution for deep networks. arXiv preprint arXiv:170301365 (2017).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNovakovsky G, Dexter N, Libbrecht MW, Wasserman WW, Mostafavi S. Obtaining genetics insights from deep learning via explainable artificial intelligence. Nat Rev Genet. 2023;24:125\u0026ndash;37.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShrikumar A et al. Technical note on transcription factor motif discovery from importance scores (TF-MoDISco) version 0.5. 6.5. arXiv preprint arXiv:1811.00416 (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRauluseviciute I et al. JASPAR. 2024: 20th anniversary of the open-access database of transcription factor binding profiles. \u003cem\u003eNucleic Acids Res.\u003c/em\u003e 52, D174\u0026ndash;D182 (2024).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGasperini M, Tome JM, Shendure J. Towards a comprehensive catalogue of validated and target-linked human enhancers. Nat Rev Genet. 2020;21:292\u0026ndash;310.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKaplan J et al. Scaling laws for neural language models. arXiv preprint arXiv:2001.08361 (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePeidli S, et al. scPerturb: harmonized single-cell perturbation data. Nat Methods. 2024;21:531\u0026ndash;40.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePapalexi E, et al. Characterizing the molecular regulation of inhibitory immune checkpoints with multimodal single-cell screens. Nat Genet. 2021;53:322\u0026ndash;31.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHeumos L et al. Pertpy: an end-to-end framework for perturbation analysis. bioRxiv 2024.08. 04.606516 (2024).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19:15.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFrankish A, et al. GENCODE: reference annotation for the human and mouse genomes in 2023. Nucleic Acids Res. 2023;51:D942\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWolf T. Huggingface's transformers: state-of-the-art natural language processing. arXiv preprint arXiv:1910.03771 (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDevlin J. Bert: pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:181004805 (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePaszke A et al. Pytorch: an imperative style, high-performance deep learning library. arXiv preprint arXiv:1912.01703 (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBadia IMP, et al. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinform Adv. 2022;2:vbac016.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKokhlikyan N et al. Captum: a unified and generic model interpretability library for pytorch. arXiv preprint arXiv:2009.07896 (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDhamdhere K, Sundararajan M, Yan Q. How important is a neuron? arXiv preprint arXiv:1805.12233 (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eOki S, et al. ChIP-Atlas: a data-mining suite powered by full integration of public ChIP-seq data. EMBO Rep. 2018;19:e46255.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXu W, et al. CoolBox: a flexible toolkit for visual analysis of genomics data. BMC Bioinform. 2021;22:489.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRam\u0026iacute;rez F, et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 2016;44:W160\u0026ndash;5.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBailey TL, Johnson J, Grant CE, Noble WS. The MEME suite. Nucleic Acids Res. 2015;43:W39\u0026ndash;49.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Transfer learning, Genomic DNA models, Genome Language Model, Gene expression prediction, Perturbation","lastPublishedDoi":"10.21203/rs.3.rs-6461531/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6461531/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003ch2\u003eBackground\u003c/h2\u003e \u003cp\u003eTransfer learning applied to genomic DNA models has the potential to improve predictive capabilities, especially when target-domain datasets and computational resources are limited. Despite its promise, the practical effectiveness of transfer learning in genomic DNA models, particularly for predicting gene expression changes due to perturbations, has not been thoroughly investigated. This study aimed to systematically evaluate the performance and utility of transfer learning approaches using genomic DNA models to accurately predict perturbation-induced gene expression.\u003c/p\u003e\u003ch2\u003eResults\u003c/h2\u003e \u003cp\u003eWe benchmarked three genomic DNA models across 12 distinct datasets containing perturbation-induced gene expression data to identify optimal conditions for effective transfer learning. Notably, perturbation-induced gene expression data were not included in the pre-training of these genomic DNA models. Among these, the Enformer model consistently generated accurate embeddings, demonstrating superior clustering performance and gene signature scoring aligned closely with observed experimental data. Additionally, we identified a phenomenon termed \"genomic neighbouring gene interference,\" wherein partially overlapping DNA sequences of adjacent genes cause correlated predictions, resulting in both beneficial and detrimental effects on predictive accuracy.\u003c/p\u003e\u003ch2\u003eConclusions\u003c/h2\u003e \u003cp\u003eOur findings highlight the efficacy of transfer learning in genomic DNA models for predicting perturbation-induced gene expression, particularly emphasizing the Enformer model's robust performance. Understanding genomic neighbouring gene interference offers critical insights for refining predictive accuracy in genomic applications. This study provides practical guidance for researchers developing transfer learning strategies and genomic DNA models, paving the way for more accurate and resource-efficient genomic predictions.\u003c/p\u003e","manuscriptTitle":"Transfer learning using generative pretrained genomic DNA models for predicting perturbation-induced changes in gene expression","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-05-05 02:32:27","doi":"10.21203/rs.3.rs-6461531/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"37ba3217-5309-4c9b-8c32-fac1ea4b0ac1","owner":[],"postedDate":"May 5th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2025-09-26T10:09:03+00:00","versionOfRecord":[],"versionCreatedAt":"2025-05-05 02:32:27","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-6461531","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6461531","identity":"rs-6461531","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}