Scientist Stories: Jacob Schreiber, Dissecting cis and trans regulation of chromatin accessibility

Jacob Shriver. I'm a postdoc working with Anshul Kundaje at Stanford. I'm really excited to have this opportunity to show you some recent work we've been doing with a model that we're calling Dragonfruit. So most biological processes are driven by continuous sets of interactions between cis- and trans-regulatory factors. This is a dataset I'm about to show you that was originally collected by Mo and analyzed by Srivag from our lab. So I couldn't be here today because he's about to defend in a few hours. So on the left-hand side, we have a schematic showing the experimental design. In the middle, we're going to have a UMAP. This is going to be a reprogramming time course. The UMAP is going to be of the chromatin accessibility. And on the right-hand side, we're going to have gene expression values at landmark cell states. At day 0, we have just pure dermal fibroblasts. These are just pure sample gone from patients. Then they were infected by 3 Sendai virus vectors that contained the 4 Yamanaka factors. After day 2, it seems like a lot has happened, including the emergence of a population of cells that have very high expression of the 4 Yamanaka factors, indicating successful infection. After 2 more days, it looks like there's a transition into a population of cells that simply have high, around endogenous levels of expression of these 4 factors. Then we can see that more continues to happen as the days progress. Around day 8, it seems like things have stalled, and so we switch from a fibroblast media to an iPSC media. Things start to continue again, and around day 12, we get pre-iPSCs, which are cells which are similar to iPSCs but don't yet have the full differential characteristics of them. After letting that flesh out and waiting another 2 more weeks, we see the emergence of real pure iPSCs. So just by looking at the cells across the days that they've been collected, we can see the main trajectories that are going on. There's a no-factor uptake trajectory. These are fibroblast cells that just sat in culture waiting and nothing happened. The virus never got to them. Then we have the primary reprogramming trajectory. These are cells that began as fibroblasts, went through this continuous process, and ended up as iPSCs. But this is really interesting because these cells started with infection from the virus and the overexpression of the four Yamanaka factors, altering trans-regulatory state, protein concentrations. These factors then went and bound to their motifs according to the affinity of the motif, spacing with other motifs, and other characteristics, cis-regulation. By binding, these then alter endogenous gene expression, once again altering trans-regulatory state.. And this process just continues and continues until you reach equilibrium in the iPSCs. However, despite this process clearly being continuous, most modeling approaches don't actually treat it that way. When Saurabh first analyzed this data, he did so using the standard methods for the time, which basically involved clustering the cells and then looking at each cluster. So here he started off by clustering our data into 15 landmark states, not thinking that each one was a different cell type, just a landmark along the process. Then he aggregated cells according to their cluster. He ordered these clusters according to some notion of pseudotime, and then he inspected loci to see whether or not they exhibited dynamic activity with respect to that pseudotime. I'm going to make you all watch that again because it took me a long time to make this. But the locus that I was showing you is NANOG, which is one of the genes that's heavily involved in reprogramming. And we can see that during this trajectory, it starts to become more accessible. And so our question was, can we move away from this paradigm and actually model this continuous single-cell data at continuous or single-cell resolution? Basically, if our data looks like this, where each one of our cells is a different row and each one of our loci is a different column, can we predict and interpret each one of the pixels, each one of the loci at each one of the cells? And so we need a representation of our loci, which are just genomic sequence. We need a representation of our cells, which can be a very flexible, just vector representation that you think faithfully represents a cell. I'll talk more about that in a second. And we need a target to predict, and that's just the number of reads mapping at a locus. But if you've dealt with single-cell data before, you know it doesn't look like that. It looks incredibly sparse. And so what we actually do is we are trying to predict— we're trying to predict dynamic pseudobulbs, which is basically smoothing along a KNN graph. All right. So, what does the model look like? Well, for those who are familiar with it, we start off with a BPNet model that takes in nucleotide sequence, and it makes predictions for some sort of activity of interest. Because we are dealing with ATAC-seq data, which is generated by a TN5 enzyme, we have to account for the sequence specificity of the TN5 cutting enzyme itself. This is a really annoying characteristic of TN5, that it's supposed to be making cuts at things that are biologically driving accessibility, but it itself, you know, cares about the sequence. So this is something that was originally being proposed by Anna Sabina and Anusu Pamphuri in the ChromBPNet paper. And this is the ChromBPNet architecture, if you've heard of that before. So the final part is how do we incorporate cell representation?. And basically we just use multiple layers of neural network magic in order to pull out the parts of the features that matter. And here's the key of the model. What we are predicting are the parameters of the convolution layers of that BPNet model. This is key. We are not simply extracting a representation of the cells and concatenating it with a representation of the sequence. We are using the representation from the cells to predict the parameters of the convolution. We are processing sequence differently depending on which cell you are in. And then we add in the read depth in order to account for read depth differences. All right. So, the first key point about this is that quite literally, you end up with one predictive model for every single cell in your dataset. Everything you are doing at bulk resolution right now, calculating attributions, running ISM, doing in silico marginalizations, Any of that you can now seamlessly do at single-cell resolution. The second key point is that because this is just a vector representation of a cell, it can come from anything. I'm about to show you some results showing what happens when you use a representation that comes from the same ATAC-seq data you're about to predict. But if you have multimodal data, your representation can come from gene expression. And then you can do in situ knockouts where you say, if we knock out this gene, what happens to the chromatin accessibility landscape across an entire experiment? If you have spatial coordinates, you can simply concatenate those in. If you have your favorite batch correction method, you can run that and then feed those embeddings in to this model. So does this work? Quite simply, yes. That if we evaluate the single-cell resolution, what we're comparing against here is if we train one ChromoVP-NEF model for each one of the 15 clusters and make predictions for it, the same prediction for each one of the cells in that cluster. The more blue the cell is, the better Dragonfruit's predictions are over that baseline ChromBPNet model. And you can see that what's happening is basically that ChromBPNet does a pretty good job of modeling the center of clusters, which makes sense. If you're aggregating all the reads across a cluster, then you're going to do a pretty good job at getting the average. But Dragonfruit does better the further away from the centroid you are. All right, then we wanted to really zoom in and start to interrogate biology. So we looked at the promoter of the NANOG locus. NANOG again is the gene that— one of the genes that's heavily involved in reprogramming. And the key players here is that there's an OCTSOX heterodimeric site. There's an array of 3 weak KLF sites. Then there's a ZEB repressor. If we look at the attributions in fibroblasts, we see that nothing's happening because this gene isn't active in fibroblasts. If we look at the beginning of reprogramming, you can see that things start to pop out, which makes sense. Now, OXOX and KLF are expressed in cells, so they're free to bind to their sites. Then if we look at the end of reprogramming, you see that OXOX is still around, but that the activity around KLF has really diminished. With Dragonfruit, like I was mentioning, we're not restricted to simply looking at snapshots like this. If we want to see what's happening with this OXOX site, we can look at it continuously across the entire experiment. So we can just sum attribution values and color our UMAP according to that. To be clear, this isn't a global measure of aux aux. It's not like gene expression values or other methods that give you global measures of this motif. It's saying for this single locus, what is the utilization? What is the predicted utilization? If we look over at KLF, you can see there are different dynamics, that it dies off a lot faster, which makes sense because KLF is not endogenously expressed by any of these cells. It simply decays when the construct decays. Dragonfruit doesn't know what a repressor is, and yet it still learns that this ZEB site closes. It's a repressor of accessibility, and it turns off when the fibroblast program turns off. But we were interested by this. We noticed this kind of fuzz in the attribution space in both the Chrome BPNF models and the Dragonfruit models, and we wanted to investigate that further. And what we found is that because the array of KLF sites is weak, it is able to accommodate two weak ZIC sites. ZIC is a protein that only begins to be expressed in iPSC. So it turns on at the end of the reprogramming trajectory. And so we wondered, is it possible that what's happening is that these same nucleotides are transitioning from being utilized by KLF at the beginning of the trajectory to being utilized by ZIC at the end? We ran a bunch of computational experiments and I don't have a chance to really go into them, but we were able to confirm that. And so what we did is we basically created a local version of pseudotime based off this. We turned this into a clock. We said, can we order cells according to the progression of KLF utilization versus ZIK utilization? And so you can see this score with orange being more dominated by KLF and purple being more dominated by ZIK. When we compare this with a global measure of pseudotime, we see that it corrects a big mistake. That global pseudotime puts the iPSC cells in cyan before the pre-iPSC cells. This was a consistent error we were noticing, which simply isn't biologically possible. But when you use this biologically derived local clock, you see that it fixes the issue immediately. So we then ordered all of our cells by this version of pseudotime, and we looked at the attributions of a bunch of different of all the sites there. And you can see that what happens matches what you'd expect. Octin-Sox have higher utilization as you have both ectopic and endogenous expression of those proteins. And then it decays to basically endogenous levels of expression. KLF, because it's never endogenously expressed, simply decays as the construct does. But that this decay is mediated by the affinity of the KLF site. Stronger KLF sites have higher utilization at the beginning and diminish slower. And then there's that. All right. So this slide is really cool, and I really wanted to show you this. Basically, it— the nanogene is full of really cool stuff. But in the enhancer, there is this site which has a SOX binding site, which you might expect because SOX is one of the key players here. But then there's another SOX site nearby. And there's a third SOX site nearby. And then there's an aux site right next to it. And you're like, oh, that makes sense. It's an aux site. But then there's another oxo site right there. This is basically an oxox heterodimer that where the sox has a secondary ox and the ox has a secondary sox. And the spacing between them matters a lot because the sox in the main oxo site is actually quite weak compared to what you would expect, but it's stabilized by the fact that it's surrounded by these oxides and the spacing on the right-hand side between the oct and the socks is important because it, you know, it's positive, it's non-canonical. And then when you look on the other side, there's negative spacing between that oct and socks, which are kind of pushed together, making it a weak ox site. Not only is this interesting, it exhibits really interesting dynamic activity. If we look at dragonfruit attributions across this trajectory, you can see that that weak ox site diminishes in attributions quite rapidly, but that the weak sock site that's right next to it continues to have high attributions with what you'd expect a much higher affinity SOX site to, simply because it's being stabilized by the OX sites that are next to it. So I think that this is really a key point. We can't just look at these sites according to, is there a SOX binding site? We can't look at it binary. We can't just say like, what is the maximum affinity site? This recruits OctonSOX much higher than any individual OXOX binding site does. And this is crucial. To understanding how the genome is regulating genes and other biochemical phenomena. All right. So Dragonfruit explicitly models cis and trans-regulation. And by doing that, it's able to, you know, generalize well. That the accuracy is quite high compared to cluster-specific methods, suggesting it's really picking up on the continuous resolution of these experiments. These attributions reveal subtle and continuous changes in motif utilization, even in the presence of a much stronger motif that is overshadowing everything that's going on. That ox oxide in the promoter is really overshadowing everything in terms of accessibility. So, finally, just like a request, when dealing with genomics data, the more continuous you can do your modeling, the better. Remember, cell state isn't discrete. Motif affinity isn't discrete. Anything that you can model continuously is going to be closer to modeling the biological truth. So I want to thank my collaborators as well as the brains behind the operation.

Great.

So questions?

Really nice talk.

Thank you. I was wondering how you actually train the models. I mean, do you have any identifiabilities? Because the models is exponentially more complicated than the original.

Bulk-level models, right? I don't know if I'd say that it's exponentially more complicated. The time of a forward pass through this model is the same as the forward time through a Chrome BPNAP model because the vast majority of time gets spent running a convolution across a sequence. So dynamically generating that bottom part uses like less than 1% of the total time. It is true that it takes more iterations in order to train, but it takes fewer iterations to train one Dragonfruit model than to train 15 ChromoDP-Net models on the same dataset.

Okay, thank you. Nice talk. I was wondering if you can talk a little bit about what would be the difference in terms of what your model learns versus modeling the data at the single-cell resolution, something like SC-Basset, Or if you want to address sparsity, something like KNN, sorry, like metacell level of single-cell data so you don't have these like broader clusters.

Yeah, so there are kind of a few, there are a few questions in there. The first is, what is the effect of doing KNN smoothing? And the simple reality is that nobody does modeling at truly single-cell resolution. The data is simply too sparse in order for any of these methods to do it. If you use PCA, you're not doing it at single-cell resolution. If you're doing anything that involves neighbors, you're not doing things at single-cell resolution anymore. And almost all of our steps require that. And so it is true that the neighborhood size is going to be something that matters and it's going to be something that's dataset specific, but it's not unique to the Dragonfruit framework.

So my question was going to be around like how much did this work like out of kind of like out of the box? So you have the cell representation, dense layers, and that actually gives you the the filters themselves. Is that how this works?

Yeah, so there are a few iterations of this, and it depends on what types of things you think in there, that you can be predicting the full convolution, all the weights, and that's what we did initially. You can also just predict the bias, and that for a lot of cases, that's simply saying what proteins are active or not active in this cell state. And that we think that there are reasons why you'd want to predict the full convolution, because spacing can make a difference. The interactions can make a difference. But that for a lot of the earlier stuff, we just used the predicted bias term.

Okay. And in terms of like initializing weights, so I mean, we, all this stuff we think about, about like the scale of the variance for the, the noise that initializes weights and things that, that just kind of works with default settings.

I got a reasonably good model out of the box. I needed to do a lot of tuning in order to get a good model. But the main things, like aux aux being active across the main trajectory and that type of thing, popped out immediately.

Thank you for the nice talk. I think it's a great step towards using model attributions as a proper tool to featurize the signal. I was wondering if you reconstruct the single-sine matrix using the attributions as features and infer the pseudotime based on this new representation of the matrix capturing the dynamics through model attributions, would it be more accurate than the pseudotime that you are using, which is based on RNA-seq?

Oh, so the global pseudotime we were using was based off of accessibility. But I think that the idea of a global pseudotime is not— it's useful, but it's also an approximation. The genome, as Anshul always says, the genome is asynchronous. Different things are happening at different parts of the genome at the same time. You can't order all cells according to one number, but you can order local processes based on what's going on there. And so the idea of having a local attribution-based pseudotime, I think, is useful for modeling what's going on at a local site. But I think that's not going to be great at modeling something that's going on at a totally different chromosome, some other Okay.

Sorry, I think we better move on for the sake of time. I'm sure Jacob will be around to answer more questions. So let's thank him again. Next up, we have Karen Issiev from my lab. He'll be telling us about single-cell splicing analysis.

Yeah, I should be fine. The previous one was also fine. You see