Scientist Stories: Chenwei Zhang, Visualizing DNA hybridization trajectories with biophysics

Hello everyone. First, I would like to thank you all, MLCB organizers, to hold this amazing conference. My name is Chunwei Zhang. I'm a third-year PhD student at the University of British Columbia. Today I'm going to talk about DNA reaction visualization. DNA reactions are very important in applications of DNA computing, DNA storage, and other advanced DNA nanotechnologies. However, the interaction mechanisms are very elusive because they involve high-dimensional trajectories over secondary structure state space. So in today's presentation, I will introduce a novel method which integrates the deep graph embedding with biofreaks-informed features to map this high-dimensional trajectory data into low-dimensional space to visualize, therefore help us to understand the interaction mechanism. So first, what is a DNA reaction? A DNA reaction is the strand interaction where the strands folding or unfolding, binding or unbinding from the initial to the final states. And this process is stochastic. There are many types of DNA reactions. In our work, we are focused on DNA hybridization. This animation shows the gives an intuition about what the hybridization looks like. There are two unpaired complementary strands start binding by interstrand base pairs and eventually forming a fully paired double helix structure. Actually, this animation just shows a very simplified example, and it does not often happen in real. Instead, during the hybridization process, it's very possible to form hundreds of thousands intermediate states, which makes the visualization problem even more challenging. And this figure shows an intermediate secondary structure with a 3-stem hairpin structure, and the 3-stem hairpin have effect on the reaction rate, which makes the reaction mechanism even more complicated. So in today's talk, I will split the talk into 4 parts. First, I will give some brief background and then talk about the method. Next, I will show the results, and finally, I will point out our future direction. So before I'm talking about the reaction trajectory, let me first introduce a concept which is called the dot-princess notation. We are using this notation to represent a secondary structure. The notation is shown right here. So each dot refers to a base, and the open and closed parentheses refers to a base pair, and this plus sign separates two different strands. The reaction trajectory is a sequence of the secondary structures from the initial state to the final state. Along with the transition time from one state to the next. And we are interested in studying elementary steps, so the successive structure in the sequence only differs by one single base pair. In our project, we are using a kinetic simulator which is called MultiStrand to generate trajectories. Specifically, we studied an experimental hybridization reaction from Gao et al. And there are 25 bases in each strand. We found that the formation of 4-stem hairpins significantly slowed down the overall reaction. This big table shows the output from the multistrand for this reaction. The first column shows a reaction trajectory, and the second column shows the transition time, and the third column shows the free energy. The free energy is decreasing during the hybridization process. But this table only shows one single reaction trajectory, which may include hundreds of thousands of intermediate and elementary step secondary structures. So if we would like to visualize hundreds of thousands of trajectories, it will be more challenging. So one way to visualize these trajectories is using the coarse-grained method that groups the sector structures by a given policy and then laying out the trajectories. This plot shows a, a specified coarse-grained method. The x and y-axis refers to the number of base pairs in the final structures and not in the final structures, respectively. However, there are many limitations of this coarse-grained method. First, the single cell may contain states with different structures which makes the interpretation very difficult. Second, the layout trajectories are overlapped and is unable to distinguish, which is not good. Overall, the coarse-grained method seems not an ideal option for visualizing the interactions. So now we should question, what is a good visualization, or what is a good embedding for domain experts? We argue that a good embedding should have a good resolution for individual secondary structure. And also, we argue that a good embedding should preserve global structure, which means the state embedding should follow the trend of the course of a reaction. And a good embedding should also preserve local structure, which means that the neighboring states should have the similar structures which is only differ by a few base pairs, and the states that are far from each other should have very different structures. What's more, a good embedding should also lay out the trajectory smoothly, and there's no big jump between the elementary steps. Keeping these questions in our mind, we propose FATA, a new framework which integrates deep graph embeddings with bioFIX-informed features to visualize the gene interactions. And our method can create the fine-grained low-dimensional secondary structure state space with high resolution, and VEDA can also preserve good global and local structure, and the visualized trajectories are smooth and distinguishable, which could improve the interpretation of reaction mechanisms. And I will show the results later. So now let me switch over to the method. Our framework has 3 major parts. The first part is a scattering transform, which maps the graph to the real vector space. Each graph is a secondary structure where the nodes refer to the bases and the edges refer to the backbones and base pairs. The second part is a variational autoencoder, with the loss functions augmented by 3 biophysical features, the free energy, the minimal passage time, and the graph edit distance. And the third part is a fit model, which is a nonlinear dimensionality reduction method. So the input vectors are embedded by the autoencoder, and the embedding are then passed into the fit model to further reduce the dimensionality to 2D for the future plot. The top figure shows a detailed version of our model. In this model, the total loss is made up of 5 terms. In the first 2 terms are VAE losses, which are unsupervised, and the last 3 are supervised by FrIC's informed losses. The free energy loss is calculated by the mean squared error between the predicted energy term and the multistrand generated energy. The graph edit distance loss is calculated by taking the pairwise Euclidean distance between state i and j in the variational autoencoder latent space, subtract the graph edit distance between state i and j. Similarly, the minimal passage time loss is calculated by taking the pairwise Euclidean distance between state i and j in the VAE latent space, subtract the normalized minimal passage time from state i to j or from state j to i and multiply by an important weight vector, which is derived from the state probability. So now let's take a look at the results. This plot shows the embedded segmental structure state space. In this plot, each circle represents a segmental structure., and the color of the circle refers to the corresponding free energy. The initial and final states are highlighted by the green bulbs. From this plot, we could see that the state embedding follows a high-to-low energy trend from the initial state to the final state. This is consistent with the reaction process, indicating our embedding preserves global structure. We also investigated the neighboring states and found that they have very similar structures and energies, indicating that our embedding preserves local structure as well. We also noticed that there are two branches in the plot, and after investigation, we found that the two branches refer to two reaction pathways. The top branch refers to the strands binding from the outside to the inside, and the bottom branch refers to the strands binding from the inside to the outside. And laying out the reaction trajectories, which shows as colored curves, we find that the trajectories proceed nicely along the embeddings and there's no big jump occurred. And these two trajectories are separated from each other. From the trajectories, we also found a few dense regions. And we attributed these regions as energy basins, which include kinetic traps. So the reaction gets stuck in the basins and takes a huge amount of steps and time to get out of the basin towards the final structure. We are very interested in these energy basins. So in order to retrieve the kinetic trap information, we applied an unsupervised method called DBSCAN to to cluster the post-processed map, which filters out the low probability states. And there are 3 resulting clusters, and we identify the states with the minimal free energy as the kinetic trap in each cluster. The trap information is shown in each box, and for the trap in the blue cluster, we find that they have 4 stem hairpins on both strands. This is consistent with the previous computational analysis by Schwenk and all. In addition, we also noticed 2 more traps. The trap in the green cluster, which has only 1 3-stem hairpin in 1 strand, and we attribute this trap as a minor trap because it only has 1 3-stem hairpin. For the trap in the red cluster, we observe two three-stem hairpins in both strands, which is strong enough, so we attribute this trap as a major one. The three-stem hairpin is very different from the four-stem hairpin because they are in completely different position of the strand, so this three-stem hairpin is not an intermediate stage in the formation of the four-stem hairpin. Instead, it refers to a new reaction mechanism. And this is a new finding by only using our visualization method. And this finding can prompt domain experts to conduct further analysis of the specific traps which also slow down the overall reaction. Actually, we also did other case study for a completely different reaction. But for the sake of time, I won't talk about it here. So if you are interested, please feel free to check out our paper. We also compared VEDA with other common embedding methods where some methods are completely relied on the data itself, such as PCA and PHATE, and some methods are relied on part of biophysical features, such as GSAE. Because there's no ground truth for us to use to compare, so we designed a quantitative metric which called average distortion to measure the average trajectory step length in the normalized embedding map, and the lower value represents the better embedding results. So we find that our method achieved the lowest value, and what's more, we also Show the corresponding visualization plot by each method, and by comparison, we find that the VEDA has the most smooth trajectory layout. And this visualization result agrees with the distortion score, indicating our method has the best performance. So this result confirms that integrating bowfreaks-informed features with deep graph deep neural network improves the embeddings and reduces the trajectory distortion. By the way, we also tested T-SNE and UMAP, but they all give very bad results, so we don't show here. All right, so because of the promising results from VEDA, it guides us to further direction. So in future work, we plan to generalize VEDA For more complicated interactions, such as 3-way strand displacement, which plays a very important role in DNA computing, we also plan to explore VEDA with other types of datasets, such as RNA reaction and molecular dynamics trajectories, which is very important in drug discovery study. What's more, we also plan to further improve VEDA embedding so that it's able to separate the secondary structures in the strand level. So last, I would like to thank my supervisors, Ed and Kuan, and my coworkers, Jonah and Boye, and my collaborator, Eric, for their help with this project. All right, so that's all of my talk, and thank you for listening, and I'm open for any questions.