Optimization of regulatory DNA with active learning

2.1. Active learning on a synthetic genotype–phenotype landscape

We focus on optimization of regulatory DNA sequences using sequence-to-expression (STE) machine learning models as predictors. In a typical use case, an STE model is trained on $n$ measured pairs $(\mathbf{x},y)$ of genotype–phenotype associations, where $\mathbf{x}$ is a DNA sequence of length $N$ and $y$ is a readout of expression, typically quantified via fluorescence reporters or suitably designed screening assays that couple expression to fitness. In one-shot optimization, the STE model remains fixed and is iteratively queried by a sequence searching algorithm, typically using hill climbing or global optimization heuristics. We focus on an active learning approach (Fig. 1A), whereby the STE model is retrained during $M$ learning loops, with new batches of data containing $q$ new samples $(\mathbf{x},y)$ that are iteratively selected through the optimization loop. By careful design of a sampling-and-selection routine that balances exploration and exploitation of the sequence space, active learning can continuously improve model accuracy and navigate the predicted landscape in high-confidence regions of the STE model.

To first explore the ability of active learning to traverse highly non-convex expression landscapes, we focussed on synthetic data generated with a theoretical fitness model [43]. This allows sampling the landscape across the entire sequence space and computing global maxima as a ground truth baseline. We focussed on the NK model for fitness landscapes, because it has few parameters and produces landscapes with tunable ruggedness. The NK model was originally developed for gene-to-gene interactions [42], and we adapted it to model epistatic interactions between positions in a nucleotide sequence (Fig. 1B). In the adapted model, $N$ is the number of positions in a sequence, and $K$ models the number of other positions each nucleotide interacts epistatically with (i.e., the order of interaction). The parameter $K$ can be tuned between $K=0$ (no epistasis, only additive effects) and $K=N-1$ to control the ruggedness, with larger $K$ leading to increasingly complex landscapes with many local maxima and minima; more details on our implementation of the NK landscape can be found in the Methods.

We generated several NK landscapes for sequences of length $N=10$, resulting in a sequence space with a total of $4^{10}=1,048,576$ variants. To visualize the ruggedness of the landscape, we employed the t-SNE dimensionality reduction algorithm (Fig. 1C). As $K$ increases, variants with more extreme phenotypes tend to appear more dispersed across sequence space, reflecting the increased ruggedness and non-convexity of the landscape. To first assess the challenge of modelling such complex landscapes, we trained feedforward neural networks on NK fitness values (Fig. 1D). Models were trained on 2,000 sequences with Latin Hypercube Sampling (LHS) (Supplementary Figure S1), which represents a coverage of $\sim$0.2 % of the sequence space; this is substantially larger than some of the largest datasets employed in the literature. Prediction results on a held-out test set show that in the absence of epistatic interactions ($K=0$, Fig. 1D), the landscape can be regressed with reasonable accuracy. However, for $K=1$ and above, the introduction of higher-order interactions makes the landscape significantly more challenging to regress, even after landscape-specific optimization of the neural architecture.

We built an active learning loop designed to find the global optimum of the NK fitness landscape (Fig. 1A), assuming an initial set of $n=1,000$ variants for training that were randomly selected from the whole input space. We employed an ensemble of $r=10$ feedforward neural networks as a machine learning regressor of the fitness landscape. At each active learning loop, sequences are sampled and selected based on the model predictions and a reward function. First, the ensemble is queried with $10,000$ sampled sequences, and the Upper Confidence Bound reward function $J_i$ is calculated for each sequence $X_i$:

where $y_{ij}$ is the predicted value of the $i$th sequence with the $j$th model in the ensemble. The terms in Eq. (1) correspond to the mean and standard deviation of the predicted expression levels across an ensemble of $r$ neural networks for each sequence. The parameter $\alpha$ controls the balance between sequence exploration and exploitation. For sequence selection, to emulate scenarios with limited data acquisition capability, we fixed the batch size to the top $q=100$ sequences ranked by their reward function value. At each learning loop, the selected batch then goes into the evaluation step, and is added to the existing data for model retraining in the next loop.

To test the impact of the strategy employed for sequence sampling in active learning, we compared random sampling with a biologically inspired sampling based on directed evolution (DE) [44], whereby new sequences are generated by introducing mutations at $j$ positions starting from sequences in the previous active learning batch. The results (Fig. 2A) suggest that after four active learning loops, directed evolution sampling reaches better optima than random sampling across all levels of NK ruggedness. The final-batch genotype distribution on the NK0 landscape is shown in Fig. 2B, showing that active learning in tandem with directed evolution can effectively find optimal DNA sequences, in agreement with previous studies on protein fitness optimization [30].

We also tested the impact of exploration and exploitation in the reward function, as their balance is crucial to ensure a thorough landscape search while maintaining strong optimization performance [45]. We compared the phenotype distribution along the active learning loops in two exploration regimes using directed evolution for sequence sampling (Fig. 2C). A lower ratio resulted in a higher mean and a lower variance of phenotype, indicating a tendency to refine predictions around already identified high expression regions while limiting broader landscape exploration.

To compare active learning with traditional optimization approaches, we implemented three one-shot optimization strategies based on random screening (RS), strong-selection weak-mutation (SSWM) and gradient descent (GD). The results (Fig. 2D) across four rounds of active learning optimization suggest that active learning can outperform one-shot methods, particularly in complex landscapes with strong epistasis. We observed comparable optimization performance with SSWM in smooth landscapes ($K=0$), but under higher-order interactions, active learning can produce better optima than one-shot optimizers. Notably, sequences selected through active learning also exhibit a clear trend of performance improvement across training loops, likely due to refinements to the STE model through re-training at each loop. We also observed large batch-to-batch improvement of expression levels with active learning (Supplementary Figure S3). The success of active learning with directed evolution sampling can be attributed to its mutation from sequences that already exhibit satisfactory performance in each round of selection, which also produces STE models with improved accuracy (Supplementary Figures S4–S5). Since epistatic effects in regulatory sequences typically occur within short ranges and contiguous positions in a motif, we further explored active learning on a localized NK landscape, where interactions are restricted to neighboring positions only. The results (Supplementary Figure S6) suggest that the performance of both one-shot and active learning methods is similar to the unconstrained NK model (Fig. 2D).

Although the optimization results show that active learning can effectively ascend the expression landscape, we also observed a substantial gap between the achieved expression levels and the true global maximum; this is particularly evident in high ruggedness landscapes such as NK3, where the active learning optimum is 25 % lower than the global maximum (Fig. 2D). We explored whether an increased batch size could help improve the quality of the optimum. To this end, we ran active learning loops on the NK3 landscape with increasing batch sizes (Supplementary Figure S7). We observed a sizeable increase in the optimum when using batches of 1,000 sequences, but such gains were substantially less pronounced for larger batch sizes.

2.2. Optimization of promoter sequences in yeast

To test the utility of active learning in an experimentally characterized expression landscape, we explored the optimization of promoter sequences using a large STE dataset acquired in Saccharomyces cerevisiae [13]. These data include two sets of approximately $n=20$M and $n=30$M promoter sequences of fixed length ($L=80$nt), alongside expression readouts of a yfp fluorescent reporter in two different growth media. This dataset is ideal for our study because we can subsample the sequence variants to mimic different real world use cases, and test the ability of active learning to learn and transfer information across different experimental conditions (Fig. 3A). To query the landscape with sequences that were not screened in the original work, we employed a transformer-based STE regressor as a surrogate for the expression landscape; this regressor was developed in the original work [13] and achieved high accuracy on independent test sequences (Pearson $r>95$% in both growth media).

To mimic data scenarios encountered in applications, we considered a small initial batch of promoters in two relevant cases: one where initial promoters cover a broad range of protein expression levels, and another in which the initial promoters are enriched for low expression (Fig. 3B). This allowed testing the ability of active learning to traverse the landscape from substantially different initial starting points. We compared active learning to one-shot methods using the same total number of evaluated sequences. We employed an initial batch size $n=1000$ promoters, with $M=4$ active learning loops, and a batch size of $q=1000$ top performers according to the same reward function in Eq. (1) for an ensemble of 10 feedforward neural networks evaluated at 10,000 DE-sampled sequences; we adjusted the exploration/exploitation balance in the reward function to account for the longer sequence length compared to the previous NK landscape. In the case of uniform initial data, we observed (Fig. 3B) clear batch-to-batch improvement with active learning, which outperformed one-shot methods after two active learning loops. In the case of low expression initial data, one-shot methods struggled to escape local optima at intermediate expression levels, while active learning was able to traverse expression levels close to the measured maximum. The expression distributions in Supplementary Figure S8 also shows batch-to-batch improvements in expression across the active learning loops. Despite the comparable sequence diversity in both initial batches of sequences, we observed pronounced differences in the space of sequences explored by active learning (Supplementary Figure S9).

Examination of neural network predictions against ground truth along the optimization rounds suggests that in SSWM one-shot optimization, the fixed STE regressor tends to overshoot predictions as compared to the ground truth measurements computed with the surrogate transformer model (Supplementary Figure S10). In contrast, the STE regressor in the active learning loop progressively improves accuracy with each batch. The ability of active learning to iteratively improve model accuracy is particularly noticeable in the low expression scenario, whereby the narrow distribution of training labels leads to an initial regressor with poor generalization (Supplementary Figure S11). By iteratively supplementing the model with improved sequences, active learning can improve predictions and traverse the landscape toward better phenotypes.

Since the choice of regressor can affect optimization performance, we sought to employ convolutional neural networks as a STE model inside the active learning loop. This neural architecture has been shown to have improved inductive bias and to capture promoter motifs that correlate with expression [46]. The results in Supplementary Figure S12A suggest that convolutions can indeed improve the quality of the optimum, as compared to feedforward neural networks. Furthermore, we found that increasing the batch size can also improve the performance of the active learning loop (Supplementary Figure S12B), in agreement with our findings for the synthetic NK landscape (Supplementary Figure S7).

To test the robustness of the active learning approach, we introduced measurement noise to the initial batch of promoters, and thus forced the loop to start from poorer quality data. We gradually controlled the correlation between the initial batch of measurements and ground truth by perturbing expression levels with Gaussian noise with an increasing dispersion, as quantified by the coefficient of variation (Fig. 3C). The results show that, in line with expectation, the performance of active learning decreases for starting batches that are less correlated with the ground truth, which in turn increases the number of loops required to achieve a given expression level. However, after four loops the optimizer can effectively recover performance, achieving optima of comparable quality across all noise levels tested.

We sought to examine this robustness in more detail by simulating a data transfer scenario, whereby the loop is initialized with data acquired under different experimental conditions. This mimics use cases in which data from other laboratories is used to start an active learning pipeline. We initialized the sequence optimizer with promoters measured in growth medium B, and ran an active learning loop using the Transformer network pre-trained on medium A as a surrogate for new measurements. The results (Fig. 3D) show that the active learning loop can achieve similar expression levels for both initializations, likely because of the high correlation between initial batch sequences in the two growth conditions (Supplementary Figure S13). Furthermore, when initialized with sequences pre-optimized in medium B, the active learning loop was able to reach a further 7.4 % improvement in expression in the final batch in medium A.

2.3. Performance improvement with alternative sampling and selection

In this section, we explore further performance improvements via other strategies for sequence sampling and selection on the yeast promoter dataset. We first trialed alternatives to directed evolution sampling that aim to mimic other aspects of natural evolution. Specifically, we explored two alternative ways to introduce mutations in the sampled sequences (Fig. 4A, see Methods for details):

• •

  Genetic drift In this approach new query sequences are designed by assigning a probability of mutation to each site from a randomly picked sequence in the sequences from previous active learning loop. Each sequence was generated with 10 % probability of mutation from a random sequence of the batch in the previous loop.

• •

  Recombination We break down the sequences from the querying data of the latest learning loop, and recombine them randomly. In each loop, 10 breakpoints were randomly introduced, and the sequences from the previous loop were broken at these points and recombined.

To evaluate these sampling strategies, 10,000 sequences are sampled in each learning loop with the 3 biological sampling methods: DE, genetic drift and recombination, and the top 1000 are selected based on the reward function in Eq. (1). Fig. 4B shows the performance of sampling methods across four active learning loops. The results show that the sampling methods have a substantial impact on the active learning loop. While genetic drift performs slightly worse than DE, recombination outperforms DE, highlighting further potential of improving the active learning pipeline. In all three methods we observed a decreased but reasonable sequence diversity along the active learning loops (Supplementary Figure S9, Supplementary Figure S14).

We also observed motif enrichment in the batches of sequences sampled by directed evolution, genetic drift and recombination in the active learning pipeline. These motifs are known to bind to specific transcription factors [47]. For each AL loop, we calculated the average motif scores across the selected batch of 1000 sequences for the 244 relevant transcription factor motifs, and compared them with the average motif score in the initial batch, which has no bias towards any motif (Fig. 4C). Certain motifs were consistently enriched or depleted across all four batches and sampling methods, which suggests that active learning can capture biologically-relevant sequence features.

As another strategy to improve the active learning loop, we explored the use of introducing task-specific domain knowledge into the sequence selection step. To this end, we modified the reward function with a score that weights the presence of sequence motifs known to bind to specific transcription factors. Previous work has shown that such mechanistic information has strong correlations with expression levels [48] and helps with model generalization to new regions of the sequence space [49]. To this end, we modified the reward function by including a weight on a motif score (see Eq. (3) in Methods) mined from the YeTFaSCo database [47]. This strategy is motivated by previous work showing that most of the motifs are activators rather than repressors [48], so that maximization of the modified reward function can steer the search towards sequences enriched with relevant motifs and improved expression. The results (Fig. 4D) show that motif sampling improves active learning performance, indicating the effectiveness of including task-specific biological knowledge into the learning loop.

4.1. Protein expression datasets

4.2. Model evaluation

Model performance on the NK landscapes was scored with the coefficient of determination ($R^2$) between predicted expression and ground truth:

where $y_i$ and $f_i$ are the $i^{\text{th}}$ measurement and prediction on the test set, respectively, and $\overline{y}$ is the average expression measured. The $R^2$ score is exactly zero for the naive regressor, i.e. a model that predicts average expression in the test data for all variants in the test set, and negative $R^2$ values indicate a flawed model structure that is worse than predicting the average observation.

4.3. Active learning loop

4.4. One-shot optimization

For fair comparisons with active learning in Fig. 2, Fig. 3, one-shot optimizations were computed using MLP regressors trained on the same number of sequences as in each active learning loop. The MLP architecture was selected with 5-fold cross-validation across the same architectures in Table S1, trained with Adam optimizer and a constant learning rate of 0.001. Three methods were applied consistently across the NK landscapes and the promoter sequence landscapes, detailed next.

4.5. Motif enrichment and motif sampling

For the results in Fig. 4D, we included an additional motif score in the reward function:

where $r$ is the number of MLP predictors in the ensemble, the parameter $\alpha$ controls the balance between sequence exploration and exploitation, $P$ is the sum of the motif score, and $k$ is a penalty on the motif score; we employed parameters $\alpha =0.3$ and $k=0.66$. The motif score $P$ was computed from predicted binding probabilities between transcription factors (TF) and sequence motifs based on earlier work [49]. The 244 position frequency matrices (PFM) were obtained from the YeTFaSCo database [47], which provides the occurrences of each nucleotide at each position of 244 TF-specific motifs. Each PFM has dimensions of $4\times {\ell }_i$, where ${\ell }_i$ is the motif length. For each promoter sequence, we compute the binding probability for each motif $P_i$ by integrating over all sequence positions:

with $P_{i,j}$ being the probability of TF binding starting from the $j^{\text{th}}$ nucleotide in the 80-nucleotide sequence, derived from the corresponding PFM. The final value $P_i$ is the probability that the $i^{\text{th}}$ TF binds to at least one site within the entire sequence.

4.6. Data transfer across growth conditions

For our data transfer analysis (Fig. 3D), we employed the two pre-trained transformer models (medium A, medium B) from the original data source [13] as surrogates for ground truth data. For the data transfer across media, we first initialized the active learning loop with randomly sampled $n=$1,000 sequences measured in medium B. We then ran the active learning optimization on medium A using iterative collection of data from measurements in that medium. For data transfer with pre-optimization, we followed the same process but pre-optimized the $n=$1,000 initial sequences for medium B, and then ran the active learning loop for medium A. The pre-optimization process consists a separate round of four active learning loops on medium B, and the final batch of 1000 sequences with highest expression in medium B was employed to initialize the loop in medium A. The benchmark (Fig. 3D, blue) is the same as Fig. 3B (top).