Scientist Stories: Alexander Sasse, Biophysical convolutional neural networks

Thanks a lot for giving me the opportunity to present this work in progress that I've been working on. I'll present to you a strategy, a new strategy, to model RNA-seq data from genomic sequence with convolutional neural networks. So I guess we've heard from Jacob before. Sequence-to-function models have been established as a very useful tool to model molecular phenotypes that we can measure, for example, with the ATAC-seq. Usually these models use a genomic sequence around the measured phenotype as input, send it through a convolutional neural network, and then make the prediction of the measured phenotypes. Then afterwards we can go actually back with these models, and that's nice because they're using convolutions. We can then see what kind of features they learn that are important to the models to make a prediction, and often these features actually resemble the sequences that are bound by the biological factors that cause this measured phenotype. Most of these models are trained in a multitask fashion where earlier layers learn sort of like a combination of features or like a useful representation that is then used by the last layer only to be separated into different, for example, cell types or different data modalities to make predictions for all of these from this representation. And so the goal of these models is really eventually that we can use them to perform experiments that we wouldn't be able to do at the bench. And so for example, with these models, we can then go in, go back, and then say, all right, we perform any type of variation to the sequence, send it through the model, and let it predict what would be the effect of this predictions. So these models are really meant to be variant effect predictors eventually. And so while they have performed pretty well on a whole bunch of variant effect prediction benchmarks, earlier this year we set out to basically test them and see how well are they performing on predicting differences between expression of different individuals. And so what we used was the Informa model. So at that point, like the the most elaborate, largest model that was out there. And we used whole genome sequencing data of 839 people that was paired with RNA-seq data. And so we sent these sequences that were like 200,000 base pairs long, which had about like 100 to 1,000, couple thousand variations to the reference through the model and then compared the predictions of the model with the observed data. And so here on the right, you see this for one gene, PDX11. You see the measured values on the x-axis and the predicted values on the y-axis, and you can see that the model is kind of able to identify clusters of people with high, medium, and low expression that results in like a pretty high correlation value. But when we did this for 6,000 expressed genes in frontal cortex, we observed that So there's like a 50-50 chance that the model either gets a positive or negative correlation. However, the interesting thing is that even though it's a negative correlation, a lot of these negative correlations were also significantly strong. And so when we investigated this further, we had this idea that this could be potentially that the model learns false regulatory grammar. That basically would give it the wrong direction of the effect. However, when we investigated this further, we saw that it's actually not like a misinterpretation of the regulatory grammar, but it's more like an uninformed bias towards variants that are close to the TSS. So basically, the model doesn't really know how to predict changes with differences between individuals. And so one of the reasons we think that the model isn't able to predict these differences is that it doesn't know anything about post-transcriptional regulation. So the former and also most of the other models are trained on large datasets that are mostly covering transcriptional processes. But as the study here on the left suggests, we can see enrichments of causal variants not only in regions that are associated with transcriptional processes, but also especially in regions that are associated with post-transcriptional processes, such as, for example, the 5' UTR, the 3' UTR, or even the coding region. And so while these datasets exist, for example, CLIP-seq for the location of RNA binding protein binding, or other sequencing techniques for measuring mRNA stability, there are very few of those compared to other types of datasets. And also when they exist, it's often just like in a single cell type. And so we're suggesting to use RNA-seq directly to train these models on. And historically, it has been a little bit challenging to model RNA-seq because of large-range interactions and also because of the different layers of regulation that the mRNA goes through that eventually leads to the abundance of mRNA. However, Just recently, there's been a bunch of papers that actually used a similar approach that I will present here to model RNA-seq. And so our idea is basically go back and reprocess RNA-seq first, similar to what the approaches that I just mentioned that just came out. And so how are we doing this? Basically from this very simple model here of the life cycle of an mRNA. We can determine that the transcription rate is sort of proportional to the amount of pre-mRNA, and we can determine that the degradation rate is proportional to the amount of pre-mRNA and mature mRNA. And so with this knowledge in hand, we can turn to RNA-seq, and we can, we can think about like how we can approximate these, these two values with basically the data that we have. And so We look at RNA-seq, intronic reads are likely to come from the premature mRNA, whereas exonic reads will be used as a proxy for mature mRNA. So once we use, once we start covering, mapping basically reads to intronic reads and exonic reads, we can use that information to then get 3 different data modalities from the same RNA-seq dataset. So that is transcription rate, expression and degradation rate. And most importantly here is that now that we can do this with every data set, we get this basically for every cell type and every condition that we have RNA-seq data for. And so with this new modeling strategy in hand, I want to turn now to a, to a data set that was recently published by the Immunological Genome Project where 14 immune cells were stimulated with 6 different interleukins. Interleukins are important signaling molecules that basically, after binding to common and specific receptors, induce a cascade of events that eventually lead to changes of gene expression, as you can see here for one T cell after IL-2 stimulation in this volcano pot. And so basically the data matrix that we want to model is this one here where you see the 14 different cell types up here and you see the stimulation of these cell types in colors underneath. And so we want to use a sequence-to-function CNN model to understand what are the cis-regulatory patterns that basically cause these changes between different stimulations. So yeah, again, we're not— while we're modeling the entire data matrix with our model, what we're actually looking at in terms of like making good predictions for is like how well are we doing in predicting differences in a cell type across different conditions. So this is shown here for like one gene as an example. So this is really what we're targeting at, is like this getting a good Pearson correlation between between the measured and the predicted counts for these genes in a single cell type across different conditions, just to keep that in mind. And so, yeah, let's get back to the modeling part. So when we started modeling this with a standard model, a standard multitask model, just modeling the expression values here shown in yellow for the different cell types under different conditions, and then compared it to, um, basically a multitask, standard multitask model that was just expanded by these new data modalities that I just presented to you. So basically, just like, it has basically 3 times as many outputs as the original, um, original model. However, when we then compared basically the performance of this, again, like the correlation across different interleukins So against this, the standard model that only had like the expression values, we saw that there was no difference at all. So like these new data modalities didn't help us at all in this setting with our predictions. So we thought about it a little bit more and went back to the biology, thinking about like how is— how this is the— what's the biology behind mRNA abundance.. And as a simple model, it's basically just the fraction between the transcription rate and the degradation rate. However, we know that the degradation rate is influenced by factors that bind to the mRNA, for example, RNA-binding proteins or microRNAs that bind these in a sequence-specific manner similar to transcription factors. And so it just felt natural to use mRNA as a second input to the model to model these effects. And so what we usually use is basically this yellow part here, the sequence around TSS of this gene here, GATA1, for example. However, we're proposing to also use this exonic sequence here that comes from the axons of the gene. You might think in this case, well, this is kind of redundant, sort of the same sequence is already covered. However, now the sequence kind of reordered and potential like context effects that are located on different exons could now be in close distance, and that could be consistent across different genes, and therefore important for the model to learn the sequence context. And so this is eventually the model that we're proposing, calling this like a biophysical model because it uses DNA and RNA as input. And importantly, the DNA input in this case 40 kilobases is only predicting the transcription rate, whereas the exonic, like the mRNA part, and the mRNA model is only predicting the degradation rates that we have from the RNA-seq data. And then the expression of the eventual expression that we measure, it's only predicted indirectly as the log ratio between the transcription rate and the degradation rate. And so now when we look at the performance, again, it's in purple is the standard model that I showed you before that only predicts RNA. Yeah, only predicts RNA abundance. And in green, you see the new biophysical model. And here I'm showing you the performance is again the correlation for every cell type across different interleukins. And you see that for every cell type, we see significant improvements, especially even for cell types where we already had pretty good performance, like the PDC cells here on the right. We not only can split this into like different cell types, but we can also, because we have all these, these 3 data modalities, we can also classify these genes a little bit differently. So we can look, not— so basically we can classify the genes that have significant changes in mRNA abundance, but we can also say if this gene, if this change in mRNA abundance is also accompanied by a change or significant change in degradation rate. So we basically get—

oops, sorry—

three different classes of genes. Genes that cause changes of expression through changing degradation rate, genes that cause changes of expression changing transcriptional rate, and genes that cause it by sort of a combination of both. And so if we look at the performance differences between these three classes, like this last class kind of sticks out, and that is the class where— which of genes that basically drive the expression change through changes of degradation rates.. And so this kind of suggests to us that the model really, like, this new modeling approach really picks up on, like, that these modeling choices that we made really help them help the model to understand postconceptional effects that then cause, like, this improvement in the expression, in the prediction, especially for those kind of genes. So yeah, in summary, I hope I've shown you that current state-of-the-art models cannot predict differences between individuals' gene expression levels, show you a way of how we can use RNA-seq data to get information about RNA processing rates, and then how we can use this altogether in like a combined biophysical model to make better predictions for gene expression. And in the future, it's very exciting because there's model has also the ability now we can split basically the, the, the model into transcriptional and post-transcriptional parts, and therefore also make better analyses of like what influences post-transcriptionally and transcriptionally expression. But not only that, but also the processing rates directly. And yeah, with that, I would like to thank my group, everybody, every member, Sarah especially, as a great mentor, um, the collaborators at Imgen, and the funding sources, and you for listening. Thank you very much.

Any questions?

Very interesting talk.

Going back to this question of whether models can predict the direction of change, I'm wondering if you've gone back to each of these 3 outputs and asked if your model is getting the direction right versus just the magnitude of effect rate.

Yeah, yeah, we could, we could try that. Yeah, it's a different dataset that I'm working with right now, but yeah, we could. Retrain it on— yeah, I think that's definitely possible. That's a nice idea.

Thank you. Hi, so I was looking at like the measurements that you were— the way that you're calculating these rates. So you assume that they're pretty much constant, is that correct? Because you're saying they're proportional to the amounts of mature and nascent RNA. Assumptions are made for steady state. Yeah, so like, do you think that steady state assumption is reasonable? Because I think there's I mean, especially nowadays, there's like a lot of recent methods that have come out that don't make that assumption because it's not really recapitulating like the real kinetics of RNA, especially when you're looking at, for example, like transcriptional boosting. So I was wondering if you've ever considered that.

Yeah, I mean, I'm not really sure. I don't really know how to answer that question because I'm not an expert in whether it's true or not, but yeah, so I don't really know if I can comment on that.

Sorry.

Thank you.

Can you go back to the slide about the two different models where you were showing the biophysical model? As I understood, the biophysics, the, yeah, this one is getting an additional input, right? So you're giving it also exonic mRNA that the other model is not getting. Did you try basically like ablation experiments where given exonic mRNA, you get rid of the DNA component, or you basically put both as inputs and you try to predict the output directly?

Here, I've done— I had done this before, and it didn't— looked like it was working. That was sort of like on the process towards what we have right now. I did this, but didn't go back and do it as an ablation experiment.

Most of the concepts also remind me of RNA velocity and basically all the RNA velocity-based literature in single-cell RNA-seq. So this seems more like an effort going from sequence modeling towards RNA velocity, and maybe one can do it the other way around using sequence features and adding them to existing RNA velocity models. How do you think about, like, using the existing velocity models and adding sequences and additional modality to these models?

Yeah, I'm not really sure if I I'm sorry. Yeah, I'm not really sure how I would do this. There's definitely the idea that you could model sort of like RNA velocity with sequence and see what kind of sequence patterns kind of influence RNA velocity, but I'm not sure how we'd kind of combine all of this together. Have you thought about thought incorporating alternative splicing into this model and, you know, the existence of different isoforms which interfere with each other and can have different rates of degradation? Yeah, we've thought about that, and I think we've— I don't know, haven't really looked into it yet, like how well it works. But yeah, the idea would be that you could actually account for this, right? You could give it basically different isoforms isoform sequences, and then based on, you know, you would have to have like a, there's some sort of processing that assigns counts to the different isoforms. But then yeah, with that in hand, you could then basically totally work with different isoforms in this model approach.

Great. Okay, let's thank Alexander again. And we have lunch from now until 1:30, so please make your way back here at that time.

Yeah, yeah, sure.

Yes. We can just use this. Oh, sorry about that. Oh, that totally messed up everything.

I'm sorry.

No, no, no, I think it's going to be fine.

Sorry, which one is your—

this one? Yep. You just want to try to pull through. Just wondering if the sequence— if they're visible.

Yeah, yeah.

When it goes to the—

uh, do they normally for Google Slides