Scientist Stories: Magnus Rattray, Transcription Dynamics in Whole Embryos and Single Cells

Welcome to a new joint seminar of LS Life and the NCT Data Science Community. We are back after a really long time now. It's my great pleasure to welcome our speaker, Magnus Redway. He's a professor of computational and systems biology at the University of Manchester. Also, a warm welcome to the host of today's event, Wolfgang Huber from EMBL. Wolfgang, the stage is yours.

Well, thank you very much everybody for organizing this, and of course to Magnus coming here. So Magnus did his bachelor in mathematics and physics and then a PhD in computer science, both at the University of Manchester, after which he moved to do a postdoc in neural computing at Aston University in Birmingham. In 1998, he became a lecturer and later senior lecturer at the Department of Computer Science in Manchester. From 2010 to 2012, he had a stint as a professor of statistical bioinformatics and machine learning in Sheffield. After which he moved back to be a full professor of computational and systems biology at the Faculty of Life Sciences in Manchester since 2012. He's also director of the Institute for Data Science and Artificial Intelligence in Manchester, a Turing Fellow at the Alan Turing Institute, and also the Turing University Lead for Manchester. He's also the director of the ELLIS Unit for Manchester. Magnus mentioned that they are organizing this unit, a summer school for machine learning in health and biology, probably next June. So that's something that you can keep your eyes open for. It's bound to be very good and interesting. As a Swede in the north of England, Magnus has an unshakable and sort of quite fine sense of humor, and he's also responsible for putting Gaussian processes into the world of biinformatics in Manchester into the world of computational biology. So the floor is yours, I'm really looking forward to your talk, Magnus.

Okay, thank you.

Let me just share my screen. Thanks for that very nice introduction, and I'm really sorry that I can't actually physically be in Heidelberg, but maybe at some point in the future.

So can everyone see the screen okay?

OK, perfect.

Yes, as Wolfgang mentioned, I'm responsible for bringing Gaussian processes into bioinformatics. And there'll be a little bit of that in the talk today. So I'm going to talk about—

let me just see if I can—

two sort of projects— well, two parts of the same project, really, today. And in the first part, I'm going to talk about how we can infer transcription and degradation rates from whole embryo total RNA-seq timecourse data using Gaussian process methods. And then in the second part, I'm going to talk a little bit more about single-cell transcription. And really, the talk is about using models to try and uncover some of the mechanisms regulating these processes. So in the first part, I'm going to try and disentangle the effect of degradation on gene expression. And in the second part, I'll talk about how models can be used to try and infer promoter dynamics.

So promoter bursting type dynamics.

So we're going to be looking at embryonic development in Drosophila. And we're going to look at the earliest stage where the maternal RNAs are being cleared out and being replaced by new zygotic gene activation. And we have a total RNA-seq time series, and in panel B there, you can see a PCA plot of the RNA-seq data. And the developmental times are shown in the left panel. And you can see a nice kind of clustering together of the replicates. So there's 3 replicates at every time point. And the PCA is really showing the time evolution through the dynamics quite nicely. And you can see actually that PC1 is really well correlated with time. And indeed, recently, people looking at embryonic development have sort of bypassed staging the embryos and just use pseudotime. And we can see that if we use pseudotime here, pseudotime would probably correspond to PC1, and we'd be able to kind of infer the time back from this high-dimensional data. So the data looks kind of nice. And then PC2, I guess, if you look at this sort of L-shaped PCA plot, I think the first arm of the L here is to do with the maternal RNAs clearing out, and then the second arm is the zygotic genes coming in.

So you can really see this kind of difference in the dynamics in this nice PCA here.

And then in panel D, I'm showing the gene expression for a selection of genes. And we can see examples where the expression is purely maternal, and the maternal RNA is just being degraded away. We see some examples where the maternal is going down, and then some new zygotic expression is coming in. We see examples that's purely zygotic new expression that peaks really, really early. And then we see some others where it kind of peaks later. So there's really nice different dynamics in this data. Now this is at the gene level. Later on, I'm going to show some of this at the transcript level. And it turns out some of these maternal zygotic up and down dynamics can really be separated into different transcripts where the maternal one is dying out and a new transcript is coming in. So So really, you need to look at this at the transcript level. And then because we have total RNA-seq, it means that we can actually quantify pre-mRNA levels using intronic reads. So in the bottom right panel here, I'm showing the proportion of intronic reads in the data.

And you can see that—

initially, we just have maternal mature RNAs, so there's really no intronic reads. And then as the zygotic gene activation starts, then there's really new RNAs being produced, new pre-mRNAs, so we're getting a nice kind of intronic signal that we're going to use in our modeling. So if we look at these intronic reads, we see that they really precede the expression of the exonic reads from the same genes, as expected.

So the pre-mRNA expression peak precedes the mature RNA expression peak, and we see that nice kind of general pattern.

These here are just the purely zygotic genes, so now we're focusing on genes which are switching on during development, and we're not looking at the maternal RNAs. And in our modelling, we're going to focus on these zygotic genes. So in order to model— so what we really want to do here is get a handle on how degradation is shaping the dynamics..

And so we can use a mathematical model to try and get a handle on that.

So both the pre-mRNAs and the mature mRNAs are produced through the Pol II elongating along the gene. And if you can think of this as sort of base pair per minute, units, and then the transcription rate is going to be the same for the introns and exons, assuming that the polymerase is moving along the gene at the same rate.

And then we're also assuming that splicing out of the intronic RNA is happening at the same rate across the gene, but maybe different rates across different genes.

And so the introns are being produced and then spliced out. So that's what the first equation is saying.

And then the second equation is saying that the mature, or the exons that contribute to mature transcripts are being produced, and then those are being degraded out through mRNA degradation. And those parameters might be different from gene to gene. And we're interested in inferring what those parameters are from data.

Now in Drosophila, there have been experiments where people put probes in introns and looked at the degradation of introns to model the splicing half-life.

And it's known that the half-lives are pretty short compared to the timescales we're looking at here. So the median is 2 minutes.

And so we think it's appropriate to just make a large splicing rate approximation. And if you do that, then you find that the level of the intronic signal just is equal to the transcriptional rate divided by this splicing rate. So then you can substitute that into the second equation, and you get a simple equation for the the production of mature mRNAs, which has a contribution from splicing. So it's a splicing rate times the level of pre-mRNAs. And then there's a degradation, which clears out the mRNAs. So the two terms are production term and the degradation term. Okay, so that's just showing that equation again. So we've got data about pre-mRNAs, and we've got data about mature transcripts. It's a time series, there's noise, there's replicates, and we want to do a nice statistical model that allows us to actually infer the parameters from data in a consistent way. So in order to do that, we need a way to model the pre-mRNA dynamics, P(t), and we need a way to infer the parameters, S and D, from data. And I'm going to argue that a good approach for doing that is using Gaussian processes, and it's an approach we've used several times in the past. And so I'll spend a little bit of time now introducing what is a Gaussian process and why is it useful for this type of inference task. Okay, so people may be familiar with a Gaussian distribution and in particular multivariate Gaussian distribution. A Gaussian distribution is parameterized by a mean function and by, sorry, a mean vector and a covariance matrix. Okay. And then if you have a mean and a covariance matrix, you can draw random samples from a Gaussian distribution. And on the left here, I'm showing one sample drawn from a 25-dimensional Gaussian distribution.

Sorry, just checking the chat there.

From a 25-dimensional vector. So, and it's drawn from a Gaussian distribution with a covariance matrix shown on the right-hand side. So white here means 1 and black means 0. And so this matrix has, high numbers along the diagonal and is 0 at the sort of extremes. And the index n is the index along this 25-dimensional vector. Okay, so n 1 up to n 25. So in the middle of the covariance matrix here, we see that indices that are close together are highly correlated And indices that are far apart are very weakly correlated or don't have any correlation. Now, if I draw one 25-dimensional vector from this random number generator, I'm going to get a sample that looks like the one on the left-hand side if I plot the indices along the x-axis. And you can see it looks a little bit like a function. But it is just a 25-dimensional vector. It's not a function at the moment. A Gaussian process is just the limiting case where you take that n and you make it a continuum. So you take n to infinity. The matrix on the right-hand side becomes a covariance function. So at the moment, it's covariance matrix. But if that n vector becomes infinite, it becomes smoothed out and it becomes what's called a covariance function or a kernel function, and the points then look like functions. And here's an example of that. Okay, so now it's the same, it's really— I've drawn these in the same way, but instead of drawing 25 points, I've drawn a few hundred points. And now that I draw a few hundred points, it really looks like a function. And the Gaussian process is formally where you go to infinity here, so you really have a distribution over functions. So I like to think of this as a random function generator, and it generates these wiggly curves, and the characteristics of these curves is determined by the covariance function. So in this particular case, the covariance function is shown at the top here. It's called the exponentiated quadratic or squared exponential or RBF.

It's kind of overloaded with all these terms because it's used so much.

And it has two parameters, an amplitude and a length scale. The amplitude is the height the typical height of the curves up and down in the y-axis, and the length scale is the frequency with which the curve is crossing the zero line in the x-axis here. So the length scale is how wiggly the curves are, and the amplitude is their kind of scale up and down. And this particular choice of covariance function gives you smooth, infinitely differentiable functions like these ones I'm showing here, which might be suitable if you're modeling, for instance, bulk RNA-seq data like the data we're looking at the first part of the talk. So these are random functions, but they're smooth and have these kind of typical characteristics. And we say that f here, this function of time is drawn from a GP, which is Gaussian process. And there's different covariance functions you can choose depending on how you want to model your data. OK, so the top left here is the one I just presented, the squared exponential. That gives smooth, irregular, but stationary functions. You could, if you liked, put periodicity in the model so you can draw smooth periodic functions at the top right. You could be modeling stochastic data. So we've used rough Gaussian processes to model microscopy data similar to the data in the second part of the talk, although I'm not using GPs there, but we have done in the past. And then you might have rough covariance like the one in the bottom left, or you might have something that's both a little bit periodic and a little bit rough. Like the one at the bottom right. So the covariance function really captures the modeling assumptions in Gaussian process inference. So you might choose them based on prior knowledge about your domain. And so they capture the smoothness assumptions, they might capture the regularity assumptions, the stationarity assumptions, all the kind of assumptions that you have or the modeling choices you have can be captured in these covariance functions. So how do we use Gaussian processes for inference? Well, Bayesian inference is very natural with Gaussian processes because all the Bayesian kind of framework can be carried out exactly because these are Gaussian distributions.

Because we can do integrals, over Gaussians, then that allows us to do inference in a very natural way with Gaussian processes.

So the way it works is before we see any data, we have a belief about the type of time series that we're expecting to see, and that forms our prior distribution over functions. So if I'm modeling some bulk RNA-seq data, it's averaged over many cells, And therefore I think it's going to be smooth. There might be some experimental noise, but the underlying thing that I'm measuring should be smooth because it's averaged over many, many cells. Therefore I make the assumption that this thing is going to be smooth in time, and that's my prior. Now I make some measurements.

Here I'm making 4 measurements with very little noise. So there's very little noise, it's very clean data.

So I've observed 4 things, but there are many functions consistent with those 4, even though there's very little noise in the data, there's still a lot of uncertainty because I don't know what's going on between these points. So given these 4 points, I have uncertainty that's captured by my posterior distribution over functions. And the posterior itself is also a Gaussian process. So there's a distribution of functions consistent with the data that is my posterior distribution over functions. And then if I see more data, I can see that this is a sine wave in this case. So the data really constrains the function.

So given the smoothness constraint, then there are very few functions actually consistent with these very noise-free observations..

And so I've become very certain within the data about what's going on outside the data.

Then I'm becoming more uncertain.

And I've not, although this is a periodic function, I've not used that in my prior.

If I had a periodic prior, I'd also have quite high certainty away from the data as well. I'd be able to extrapolate, but in this case, my uncertainty grows as I move further from the data.

So Gaussian processes provide a nice way to make inferences from time course data, also spatial data. They're very popular and they really allow us to capture the uncertainty in our inferences in a nice way.

So in this project, we've used GPs 3 times in the analysis of the data.

So the first way we used it is to filter the intronic time series and the transcript time series to find things that have evidence of temporal changes. So we did that using a simple likelihood ratio test where the two models are a Gaussian process on the left and just a constant model on the right. And so we can fit the two models and just use the likelihood from Gaussian process regression to score the left-hand one.

And if there's a high score, then we keep that in the set of genes that we are going to analyze.

And we do that both for the introns and the transcripts because we need a nice signal in both of those datasets in order to do further modeling. So this is quite an old approach from Freddy Kalitas and Neil Lawrence a long time ago, but actually we've updated it with a negative binomial likelihood for counts data, which we used in this paper recently.

So we use the GPCounts package to do this, which really takes into account the noise characteristics of counts data.

The second approach that we've used in this data analysis is to do some clustering. And the nice thing about GPs is they really allow you to combine time series replicates in a nice principled way. And so here's an example where we have 3 genes, which are the rows at the top here, and then we have 8 replicates, so that's the columns here.

So this isn't— this is a different data set.

It's also from Drosophila, but it's from a different study. And then orange here is the cluster profile. So the cluster profile has— captures variation because different genes in the cluster are slightly different versions of the cluster mean, and the columns capture the fact that replicates vary a little bit, and they can even vary in where the time points are. So if you have time series where you collected data at different times, which actually happened in this particular experiment, so you can see on the right-hand side there were only 4 measurements earlier, there were more. They did later experiments later on, presumably because they thought there was something interesting going on later on that they hadn't collected earlier. So, so kind of lots of batch effects here.

And the GP can model these kind of experimental batches, time series batches. And people have used this model quite a bit now actually in lots of different applications because it's a nice way to capture this kind of batch variation within clustering.

So we're using this clustering on our data.

And then the third thing, which is the main thing I want to talk about, is that we're using GPs to model production and degradation.

So, if f here is a draw from a Gaussian process, so it's one of these random functions shown in red, then if M is produced at a rate which is proportional to F and then degraded away, so F here would be our pre-mRNA levels and M would be our mature mRNA, then if D is low, um, the half-life is long, the function M kind of integrates the pre-mRNA, so, so it kind of grows, and then peaks later than the pre-mRNA and sort of remains high. If I have a higher degradation rate, then I have a shorter half-life. The whole M function comes down because of the higher degradation rate, so there's just less mRNA in that case. But also you can see the shape of the curve changes.

The peak is earlier, and the peak is closer to the pre-mRNA peak.

In this case. And we've used this previously, actually, in this study that I've got down at the bottom here. We used it to combine Pol II ChIP-seq time series data with RNA-seq data.

But in the current study, we're going to use it on intronic reads from the same RNA-seq experiment.

And f and m together form a joint two-dimensional Gaussian process.

In the literature, these things are called multi-output GPs, and it's a special kind called a latent force model.

So you can think of the f as a latent force driving m, and you can model m and f together in their correlated Gaussian process, and you can write down the covariance function and generate GP functions like the red and black one here exactly. So if we look at some synthetic data with the same time points and same kind of level of noise that we have in our data, we can see how good this GP model is for inferring the degradation rate. And in the top here, we simulate data with a low degradation rate, so quite a long half-life of 90 minutes.

And we can see, and then we do MCMC sampling to infer the parameters, the D and the S, the splicing rate and the degradation are inferred from the data.

And we get a posterior distribution that captures the uncertainty in our degradation rate. And we can see that it's quite nicely peaked across the true value in this synthetic data.

And, you know, there's quite a lot of uncertainty in the half-life, but the GP has the advantage that it allows us to actually quantify that uncertainty. So if you have uncertainty, it's good to be able to quantify it.

And GPs are quite good in this situation where the number of time points is not too high, so there's going to be uncertainty in this type of inference. If we have higher degradation rates, then you can see that the, again, the estimation is good. And you can see that the peak here, so the blue is the pre-mRNA and the red is the RNA. You can see that the peaks are closer together in the high degradation situation.

And notice that the scale on the mRNA is different.

Okay, so the, in the top one, the mRNA is going to a higher level because the degradation rate is lower.

So really it's the shape that's allowing the model to learn about the degradation rate here.

It's not the relative height of these curves, because the relative height can be affected also by the splicing rate S, which is inferred. So if we apply this to our data, so I'm showing an example of the mapped reads at a particular time point in the top left here.

The pre-mRNA, we're quantifying just with the purely intronic reads. So we're excluding any part of the pre-mRNA which overlaps an exon of any gene of any type, because it's total RNA-seq. So we've got to be careful.

And I'll show you why in a second.

And then the mRNAs, we're quantifying using a transcript-aware algorithm.

So I think we use Callisto here. To do that. And then we're getting time course.

And these are two example genes.

They're actually transcripts of these genes. And the first one is short half-life. And you can see that the peaks are closer together, and the inferred degradation rate is low—

sorry, high.

And in the bottom one, we have a longer half-life gene.. You can see that the pre-mRNA comes up, and then the mRNA size slowly integrates that signal, and it rises slowly, and it doesn't really peak in this time frame.

It peaks a bit later on, and you can see that the degradation is inferred to be lower. There's a bit of uncertainty.

We use MCMC, so we have the distribution over the degradation rates, but you can see that they're quite clearly separated into low and high in this case. And when we look at the genes that we can model that pass all our filters for modeling, I think we fit the model to about 260 genes.

And we find a great variety in half-lives in these genes. And if you look at the ontology terms, you can see that short half-life genes are often transcription factors, DNA binding, also cell-adhesion proteins, whereas the longer half-life ones, we find things like signalling receptor binding, which, which kind of makes sense, I think. So it's known that transcription factors often have quite a high turnover, because they're very dynamic components of gene regulatory networks. So they need to have that kind of dynamic range. There's a really nice figure in the supplementary of this paper about looking at different genes and the transcript level differences. And so here's an example on the left, where the gene has two transcripts. And one of them is both maternal and zygotic. And the other one is just purely zygotic. So it's only coming on And we only model the zygotic one in this particular case. And you can see where the reads are kind of appearing in the bottom left here for that gene.

And then in the right-hand case, there's 4 transcripts. 2 of them remain low throughout, but 1 is really switching down from maternal, and the other one is switching on as zygotic.

And you can see those examples there.

Also, you have to— you can see that you have to be careful here, right?

Because it's total RNA-seq, so you get signal for tRNAs, and you get signals for this pseudogene here, which perhaps is a real gene, I'm not sure. But I guess it's a gene where they don't know what protein is, and it's quite highly expressed. So it's quite important when you're doing your RNA-seq processing to really be mapping at the transcript level, and really keep a large set of genes in your transcript set that you're mapping to, because total RNA-seq is really capturing everything.

So if we look at the pre-mRNA clusters, then we see a variety of dynamics.

And then we can take one of those pre-mRNA clusters and look at the mature mRNA clusters associated with it..

And we can see that there's a range of dynamics. So the short half-life genes really are peaking quite sharply.

But then the longer half-life genes just continue to integrate the signal and continue to rise. So there's really a variety of dynamics based on half-lives. So these are the half-lives of the genes that were modeled within these subclusters. You kind of cluster pre-mRNA and then cluster mature mRNA associated with that pre-mRNA and see quite a bit of difference in the dynamics.

As I mentioned before, the degradation rates correlate with expression peak time differences.

And it's not a perfect correlation, but there's definitely a relationship where large delay genes, large delays between the peaks are associated with much longer half-lives. So short delays correspond to quick kind of turnover genes, as one would expect. And then we did some imaging to look at the characteristics of 11 genes which fall into different categories of half-life. So we did some SM-FISH to really look at these in detail.

And so one thing we looked at is the conformation. So how close together are the 3' and 5' ends of these mRNAs?

And it appears that the short half-life category are more compact in general.

So this kind of compaction can be associated with the degradation, degradation machinery. So sometimes the mRNA has to be compacted to be degraded efficiently.

Could also be associated with translation, but we didn't find a strong association between our degradation categories and any sort of translation signal in our data. So we think maybe it's more likely that these compaction results are to do with degradation. Now, P-bodies are these cytoplasmic granules that are thought to be associated with both translation and degradation in eukaryotes.

And they're marked by ME3A.

31B here. So we looked at whether there was an association between our RNAs and P-bodies. And we use something called a P-body co-localization index, which is basically if you look at the SM-FISH, and you have dyes for the mRNA and for ME31B, you look at whether they're lying on top of each other more often than you'd expect by chance. Given the total number of mRNAs and ME31B signals. And it looks like, yes, the short half-life category have a higher co-localization with P-bodies, which is suggestive of a role of degradation within these P-bodies. And then we also looked at some RIPC data, some publicly available data for ME31B binding. So this is an RNA binding factor. And we also found an association with short half-life and greater ME31B binding. Looks like there's some kind of role for that factor, and that factor is involved in these pBodies. And then we looked within pBodies to see whether there was a bias in the number of lone 5' end mRNAs and lone 3' end mRNAs, so where we only pick up one of those signals. And it looks like there's a bias towards having less 5' end signal than 3' ends. And this is as a proportion of lone ends in the total image. So it's kind of, you know, so we're controlling for the number of lone ends of these two types. So this is kind of indicative of the potential 5' to 3' degradation that's been shown to be a mechanism that can occur within these P-bodies, and is different from the sort of more standard 3' to 5' degradation that you'd expect elsewhere in the cytoplasm. Cytoplasm. So this suggests there might be a role for this XRN1 factor in degrading these mRNAs. And this was very striking in all 11 of the genes that we imaged, we get a strong signal for this. Okay, so that's the first part. Sorry, I've taken up most of the time with that. So I'm going to quickly go through the second part, which is about really getting down to the single-cell level and looking at promoter bursting with a different type of probabilistic model. So, the MS2 system allows live cell imaging.

So rather than SM-FISH, where you're looking at the fixed image, here we're really looking at dynamics. And the MS2 system has these fluorescent probes in the 5' UTR. And those 5' UTR probes are persistent throughout transcription.

So you get a really strong signal for transcription, which allows you to look at the dynamics of transcription.

But the challenge is that because elongation takes a while, as the Pol II goes down the gene, the signal is integrated across time.

And that sort of convolutes the bursting that we think is going on underneath the signal. So we have to use a model to deal with that convolved process. So hidden Markov models are a useful way to model something, a telegraph process, a process that's going on and off, which we use to model transcriptional bursting. And so here we have a latent state, which is the on or off of the promoter.

And then we have some data X.

And if our probes were measuring things instantaneously, we could just use a straightforward HMM and do inference. And when you do inference in an HMM, you have to use algorithms that scale with the size of the state space. So here's an example of the forward algorithm, which is a message passing algorithm that we use for inference in these type of HMM things. And the message here is that the computation time scales with the state space. Now, Nicholas Lammers and others introduced something called a compound state hidden Markov model to deal with the fact that the MS2 system is integrating the signal through elongation. And so when we're picking up the fluorescent signal, the fluorescent signal can come from polymerase at different points along elongation. And that means you actually have a window, a sort of convolution of the signal of the promoter through time. And so there's a window determined by the amount of time it takes for a gene to elongate. And the hidden marker state space then scales as 2 to the power of this window size. And for one of our genes of interest, USH1, the window size is 19 and the state space is half a million. So it's a really big state space for hidden marker model inference. So to overcome this, so we can really just use the code from that paper. So to overcome this, we had to come up with a more efficient way of doing inference in the forward algorithm.

And so we did that by a dynamic truncated state space approximation where during the forward algorithm, we used the fact that the algorithm is ranking profiles.

So here we have these binary profiles through time. Here the time window is 4. So the possible states are 000, 001, and so on.

Through all different combinations.

And through time, the potential number of binary strings is growing, but we're pruning them down with this truncated state space approximation by taking away the unlikely sequences dynamically through time. And that allows the complexity of the algorithm to be bounded. So using that allowed us to really apply these models to long genes. And then we can look at how good it is.

In synthetic data, we have an inferred promoter state.

The integrated signal is shown on the right, and you can see that it's really smoothed out by this convolution process. But the inference is pretty good, even though the convolution is over quite a long window in this case. We also looked at the effectiveness of parameter estimation as we allow this truncated state space to grow. So as we allow more memory in our approximation, we kind of approach the known parameters in synthetic data or the inferred parameters using the full model in smaller cases where we can actually fit the full model. And the difference really kicks in around window size 17. So if we get windows bigger than 17, the full model really kind of blows up and we can't use it above 20, whereas the truncated model just stays the same complexity depending on the size of the number of allowed states. And so these kind of models we use to get a handle on promoter bursting dynamics. So here's the same Drosophila system as previously, but now we're looking at the spatial distribution of individual nuclei distributed through the embryo.

So before we were averaging gene expression across the whole embryos, now we're really looking at single cells.

And we can actually fit this HMM model in individual nuclei in this embryo.. And we can infer these K-on, K-off transcription parameters from the data. And what we found in this system was that really it's K-on which is driving regulation in a couple of genes that we were interested in. So looking at these genes, we saw that most of the gene expression changes were affected by the K-on rate, so the rate at which the promoter is switching on. The K-off rate remains rather similar across different expression levels. So expression here is driven by a particular regulating protein that has a particular spatial distribution. That's why there are these different spatial distributions of expression. And furthermore, the experimentalists then did some promoter swap experiments on this gene. So they swapped in a promoter from a different gene, and they found that the dynamics was very similar. And therefore, the hypothesis is that this K-on parameter is really driven by the enhancers around this gene, which are regulating this promoter switch-on rate. In this particular system. So that's an example of how this model is really giving us insights into this single-cell dynamics. Okay, so that was an example of how we're using these probabilistic time series models, and we think it's a useful way often to try and uncover regulatory processes and In these systems, we can really look in vivo, right? So like in the first case, there are experimental techniques which attempt to measure degradation by switching off gene production, for instance, but we're really trying to do it in vivo. We're really interested in what's happening that's regulated by things in the cell that would be really affected by those type of perturbations. Gaussian processes are good when things are smooth and They're really flexible. They're really good for relatively small numbers of time points or small numbers of spatial positions, I think, because then you have high uncertainty. So it's useful to be able to quantify that uncertainty. And latent Markovian models— so I talked about HMMs. We also use SDEs, stochastic differential equations, are really useful for capturing single-cell dynamics. So in future, we want to model the spatial regulation of degradation. So some of these target genes maybe have quite sharp spatial response, which indicates some kind of role for spatial regulation of degradation that we can maybe do with spatial-temporal model. We're interested in going non-stationary, especially in transcriptional bursting case.. And of course, since there's now single nucleus RNA-ATAC-seq for this kind of really early embryonic system, it'd be interesting to take a look at that with these kind of models also. Just finally, just acknowledge my team.

So Hillary leads the experimental lab, and then I put the experimental people in orange and the computational people in blue.

And Jenny, who I think is in the audience, she's a hybrid, so she's in orange and blue, and she does it all.

So, and thanks to the funders, Wellcome Trust and UKRI. Thank you. Thank you very much, Marcus. That was a fantastic overview of a beautiful model class that is very powerful, not so very tractable, and sort of Elegant. We can take questions. I think people can either type them into the chat or just raise their hands. I can give it a go. So I was wondering about your clustering of time curves which you showed on the, uh, of the pre- and mRNA levels. Um, are there any wise— I mean Since you already have this sort of nice modeling approach elsewhere, are there any good ways of getting particularly informative or distinct meaningful cluster centroid curves? Or do you, I mean, do like which clustering approach do you use? And for instance, is there a way to get particularly marked curves that have good positions of extrema and so on?

So we were using the gpclust package for that clustering. So the gpmean function of the cluster gives you the shape of that cluster.

And it's using Dirichlet process prior to select the number of clusters. So it also selects the number of clusters. But the Dirichlet process prior does have parameter which kind of tunes the number of clusters. And often the distinctness of the clusters is kind of determined by that hyperparameter in a way. So it's a little bit difficult to know from application to application what's the best choice.

Should you have like 100 clusters? Which are really— I'm thinking of like when you have time courses, maybe you, you know, you probably want one which is just flat and one which is just going up and one which is just going down, one which is maybe going up and down.

So it's all these kind of basic intuitive shapes. We don't want that, okay? And that's what people using simple models want. We don't want this sort of simplistic thinking.

We want all the kind of curves that's driven by the data. And this, I think, well, there are clustering algorithms that do, I think, use prototypical shapes, right? Which can be useful if you have strong prior beliefs about what the shape should be. But I think in our case, we're really interested in quite subtle differences. So we're interested in these really subtle—

like I showed in the example, we have these different half-lives really sculpting the curves subtly. Okay. And really there aren't clusters, there's a continuum. So it's not even real clusters, right?

I don't really— clustering is a useful way of visualizing the data, but really it's a continuum.

So it's just a way of visualizing a continuum.

I don't really believe in clusters.

So the question from Vitaly. Yeah, I can also ask this. So my question is very much related to what you just talked about. If we want to understand what are the drivers of cell identity, an obvious driver is transcriptional regulation. So I'm wondering, based on your work looking at RNA degradation rates, how much of a factor in defining cell identity both across cell types, development time and space degradation rate plays compared to transcriptional regulation?

Yeah, that's a good question. I don't think we've quantified it. But it would be kind of possible to quantify now, at least for the ones where we've modelled it, we can look at, like, how much of the mature transcript dynamics is determined by degradation differences versus transcriptional differences. I think the transcriptional differences will still dominate because in that clustering you saw that there was a big variation in the pre-mRNA shapes. So degradation, I think, is kind of shaping things a little bit beyond that. But there's definitely genes where the degradation is leading to the signal integrating, so the gene just continues to accumulate, versus something that's very sharply changing.

So if you've got trans— so I suppose— I mean, that's very nice to see, but the question is whether it changes, whether degradation rates change, like rather than providing a way to integrate information versus actively changing?

Yeah, well, we—

there's not strong evidence in the whole embryo time series of that, but I think there's evidence in spatial data that it looks like that is the case spatially in this system. But we are working on that. So I think, yeah, that's a good question.

Is it regulated for one gene in the system in a useful way. I think the answer is going to be yes for some genes, but I haven't got that data yet.

Yeah, thanks for explaining.

Very interesting work.

Thanks.

Oh, there's a question. Yeah, thanks, Magnus. Great to see this update. I was wondering also about the clustering, maybe two related connected questions. The one is, What about the exceptions? You know, other interesting genes that just don't form to any cluster? You know, you want to look at the exceptions and whether that's useful. And the second, maybe related topic is, is there a way to annotate these clusters into sort of regulatory, you know, classes by whatever overlaying with the ATAC data or looking at sequence motifs? You know, that you can go one step further and leverage the statistical strength of having many, many genes that behave in the same way, that maybe share specific sequence patterns or epigenetic states.

Is that feasible? Yeah, I think absolutely. I think that the— one would imagine— so we haven't done that, but one would imagine that the pre-mRNA clustering should be more cleanly related to regulatory factors at the transcriptional level. So separating out the clusters for the pre-mRNA and the mature mRNA should make that modeling better, with the proviso that you get less genes with a strong intronic signal. So you're throwing away some genes, but the signal for those genes should be cleaner.

Like, um, and, um, yes, we should do that.

We haven't done that kind of looking at the genomic data yet. We are now looking at— because there's quite a recent dataset from Eileen Furlong from Heidelberg on single nucleus RNA-seq and ATAC-seq.

So also there's some ATAC-seq data available for this very early embryonic thing that we can look at as well. So we kind of integrated that single cell data to get that ATAC data, but we're just looking at that now. We haven't done it yet.

But yeah, there's more to do with these clusters, definitely. Thanks. Can I have one quick question? So, so I mean, you made a nice— you made a point that you can fit very nicely smooth curves with this. In particular, maybe in spatial applications, you will probably want something that is mostly smooth, or often something that's mostly smooth and sometimes has steps. Yeah.

How easy is it to extend to such models?

Yeah, so I think in spatial statistics, people look at these kind of boundary type models as well, right?

So you can kind of capture these things.

So I think that, and then really you want to also do non-stationary modeling in that case.

So really kind of learning the covariance function from data rather than having this strong prior.

That everything is stationary. So I think you can go a long way by doing things that are more non-stationary in that kind of spatial case. So yes, I think that's an area where people have done quite a lot of nice work.

Vitaly, you got another question? Yeah, yeah. So if there is still time, It's related to what you responded to Oli, that you see, you would expect that pre-mRNA is more informative for transcriptional regulation. I'm wondering if there is any evidence whether nuclear RNA, be it pre-mRNA or mature RNA, whether nuclear RNA is more informative than cytoplasm, like thinking about single-cell RNA-seq technologies.

Yeah. Well, in a way, though, the—

so in this stage of the embryo, we don't really have cells yet.

We have nuclei, and then we have a pool of mRNA.

So actually, the single-cell data is only single-nuclei data in this system. There is no single cell. So I guess that in order to capture proper transcription, you have to look at the nuclei in the system. It's nuclei in a soup, kind of, until the cells are kind of forming during this embryonic time course. So at the end of it, you've got kind of proper cells. But at the beginning, you've got this kind of soup of mRNAs. So yes, I think nuclear RNA-seq would be more informative, but if you're doing single-cell RNA-seq, that's all you can do in this system anyway. By the way, although we looked at that nuclear RNA-seq, there is a little bit of maternal in there as well. So some cytoplasmic RNAs are sneaking into the nuclei somehow.

I'm not sure how. Okay, I think we're just right for lunchtime. So thank you again, Magnus, very much for this beautiful, clear talk.

Thanks for inviting me. Hope to see you soon in person.

Okay, bye.

Thank you. Thanks, everybody.

All right, thank you very much.

Goodbye.