Software

Motivation

Design of complementary DNA probes is not an easy task, and presents a huge problem space. For a probe designed to hybridize to a sequence of 40 nucleotides, there are 4⁴⁰ unique combinations. That's 1.2×10²⁴ combinations, which is more than the number of seconds the universe is old! This is even before considering that we need to design two DNA sequences, the probe and sink, with contrasting hybridization abilities. Brute force laboratory testing to achieve optimal probe and sink combinations is simply not possible. Existing computational tools for probe design have several limitations. We aimed to use software/modeling techniques to design tools to increase probe design efficacy and verify the probe sequences we had previously designed.

Previous Methods

Previous work has developed a computational approach to choosing probe and sink combinations¹. The authors utilize thermodynamic models to estimate the probe's fluorescence when exposed to a template strand of DNA. Using these parameters, they design a simulation which starts with a genetic algorithm to explore a broad solution space. Once a range of ideal local coordinates are chosen, a hill climbing algorithm performs a final optimization on the probe and sink combination. While this approach is helpful in exploring solutions, there are a few weaknesses we wanted to address:

1. The genetic algorithm prohibits the hill climbing algorithm from optimizing probes at a range of genomic coordinates.
2. The algorithm doesn't allow for mismatches in the probe-template interaction which drastically reduces the solution space.
3. The simulation doesn't include fail-safe statistics to inform users of potentially false positive probe designs.

Gaps in Knowledge

It has been shown that the location of mismatches in DNA sequence is extremely influential on its binding energy. Recent research suggests that DNA is least affected by mismatches at the ends of DNA strands and at the center of DNA³. However, there are few DNA modeling techniques that account for this information. We developed our modeling strategy with awareness of the current modeling shortcomings. We wanted our model to be scalable, accurate, and reflective of physical trends.

Machine Learning

Machine learning is a modeling technique that uses algorithms to learn trends not explicitly told to the program. We feel that these techniques are well suited for predicting probe binding energy since the problem represents a complex physical process that may benefit from a learning algorithm.

Machine learning is predicated on data. Without good data, quality machine learning cannot happen. In order to train a machine learning algorithm to predict free binding energies of DNA, we needed to have a dataset of DNA binding energies and corresponding information about the nature of that interaction such as nucleotide matching and sequence similarity. In the following paragraphs, we present data exploration on the learning dataset in order to demonstrate where the model may perform well and where the training data is lacking.

For our training data, we used a dataset with 695 samples². 340 of those samples are complementary, while 355 samples are non-complementary. The variance in sequence complementation allows us to train a machine learning model capable of predicting imperfect complements. In this dataset, the environmental sample temperatures where hybridization occurred assume a bimodal distribution with numerous experiments conducted at 25 or 37 degrees (Figure 1). Sequence lengths largely exist between 5-20 base pairs (Figure 2). Predictions within these parameter spaces will be more reliable due to higher sampling and more extensive training.

We hypothesized that capturing sequence parameters about both probe and template DNA could be insightful to machine learning models. Thermodynamic information could be encoded into a lower-dimensional space through sequence information, and machine learning algorithms could thereby learn to be very accurate through a relatively modest dataset. We wanted to include information relating to the content of purines and pyrimidines, as well as the continuity of the sequence.

To this end, our parameters included relative GC content, AT content, GA mismatches, GT mismatches, all other mismatch combinations, relative Hamming distance, and Ratcliff-Obershelp distance. All parameters relating to relative content aimed to encode hydrogen bonding information. Hamming distance is an absolute measurement of how similar two strings are⁴. We used this to compare the reverse complement to the template to understand how different the two sequences were. We also wanted spatially encoded information to understand where mismatches occurred. Ratcliff-Obershelp distance investigates the relative length of the longest unbroken substring between two sequences⁵. We believed Ratcliff-Obershelp distance could help the model understand where breaks in sequence matching likely occurred.

Hamming Distance

Ratcliff-Obershelp Distance

Data is often messy. Considering the wide range of possible errors that exist in the training data we acquired, we needed to find ways to create a robust processing scheme. Quantile transformation is a technique which maps potentially skewed distributions onto a normal distribution. This is a robust preprocessing technique that helps models to learn with outlier data¹⁴. However, it is also subject to frequent misuse, specifically when applied to the entire data set before making the train-test split. We made sure to meticulously exclude any data leakage from our training pipelines to ensure uncontaminated data. This ensures our model is independently validated and that we can be more confident of the results we achieved.

In order to learn the processed data, we utilized the random forest machine learning model. A random forest is a robust, ensemble machine learning model which aggregates and averages the predictions of many decision trees⁸. Decision trees are a simple learning model which make predictions based on random splits in the dataset. Since our training parameters were all scaled, we felt the random forest was a good architecture to account for class imbalance and the variance in outcome target values.

In order to create a generalized model able to predict oligo sequences longer than ~20 base pairs, we needed to scale our target data. Instead of directly predicting the free energy, we predict a free energy value proportional to the sequence length. We call this target value the Nucleotide Efficiency Coefficient (NEC). It is worth noting that we trained a separate machine learning model using scaled learning parameters, but we also included the length of the sequence. SHAP analysis, a game theoretic approach of evaluating parameter importance^{6, 7}, of this model revealed that sequence length was the most important parameter (Figure 3). Therefore we believe that NEC introduces some error, but is still the most fair way of evaluating probe sequences.

There are other programs capable of predicting the free binding energy of DNA duplexes, such as MutliRNAFold, UNAFold, and others. Those programs predict free binding energy with a Pearson correlation coefficient (r) ranging from 0.67 to 0.69 among the exact complement samples of our dataset². Tulpan et. al. was able to achieve an r value of 0.92². With our generalized scaled model, using mismatches and no sequence length information, we achieved an r value of .84. While this r value is not as high as Tulpan et. al. we proceeded with our model because it is more generalizable to sequence and is less stringent in DNA complement requirements.

One of the reasons we believe our model performs well is because we include information about the difference in substrings. Many previous models assume that where mismatches exist in the sequence is independent of free energy. Our model captures the SNP location dependence described in source³. This is demonstrated in Figure 4.

Simulation

We aimed to design a simulation around our machine learning program that was able to take advantage of the NEC. Our algorithm is a hill climbing algorithm which simultaneously optimizes a probe and sink combination to a given SNP.

The probe and sink are fundamentally different probe sequences with different goals. However, their combined goal can be identified as producing the greatest difference in fluorescence when being exposed to the mutant vs. the wild type template. To that end, an optimal solution should account for both the probe and sink sequences when determining an optimal combination. We designed a loss function to capture this goal.

Loss Function

We aim to minimize the loss function. However, this presents a challenge: since the NEC is a negative number, why doesn't the loss function blow up (binding energy for probe or sink driven to zero) where three of the parameters are zero and only one sequence binds to the reference frame?

The loss function exploits the codependency of P_w to P_m and S_w to S_m. Since both parameters refer to the same sequence and refer their NEC to very similar templates, the values are similar. Therefore the contrasting nature of the loss function prohibits any single probe from dominating.

Now, we can start to implement a proper simulation. We start by taking in two reference frames (wild and mutant) and two probes (probe and sink). We mutate the probe and sink randomly a number of times specified by the user. Then, by iteratively mutating and observing the effect of the NEC, we attempt to converge the solutions to the lowest minima. Many samples in the simulation may occupy the same solution. We therefore have a natural false positive safety mechanism built into the simulation. Predicted solutions with a higher degree of convergence are more likely to be valid solutions.

In order to end the solution, convergence must be defined. Simply, we say convergence is achieved if all samples, when randomly mutated, do not yield a more positive result compared to their prior sequence. In other words, each of the samples has gotten stuck at the bottom of a mutation free energy well.

Of course, achieving convergence for large simulations is extremely time consuming. Therefore, we integrate a more lenient and practical implementation into the simulation. Fundamentally, we implement the threshold parameter. The threshold parameter is a number from 0-1, which indicates that once that proportion of samples are not positively mutated in a single step, the simulation can stop. We also offer an alternative parameter called patience to avoid false positive convergence results. Patience indicates the number of consecutive steps the threshold must be reached in order to end the simulation.

In order to encourage more convergence and to increase the solution space explored, we also took inspiration from genetic algorithms. In genetic algorithms, there is a small chance unfavorable solutions are chosen for the next generation⁹. Likewise, in our simulation, there is a small chance an unfavorable mutation is chosen. We believe that this may help our simulation explore a discontinuous solution space. This could help our simulation reach non-obvious, but nonetheless optimal solutions.

Results

We have open sourced our simulation as Command Line Interface (CLI). We have provided documentation and a GitHub so that the community may take advantage of our work. Our CLI is capable of taking user probes and genes, running the simulation, ordering results, visualizing convergence, and doing reverse machine learning analysis to verify results. Preliminary results for F2RL3 can be found in our F2RL3 Probe Design Results table. The link to our CLI can be found at the Michigan iGEM 2023 GitHub repository.

Demo

Finally, we have also produced a demonstration video to guide users through the CLI tool. This video guide is available on our 'Michigan iGEM' Youtube page. The guide assumes minimal knowledge of command line interfaces and walks users through every step of the simulation. Users will be able to download the CLI, run simulations, and perform proper analysis of the generated solutions through our tutorial.