Scientist Stories: Shushan Toneyan, Cis Regulatory Interactions from Large Scale Neural Networks

Hi everybody, I'm Shushant Anayam, a PhD student at the CULab at Cold Spring Harbor Laboratory, and today I will talk about interpretability of large-scale deep neural networks. Large-scale cis-regulatory elements, or CREAs, are DNA sequence elements that regulate gene expression.

Transcription factors and cofactors bind to these elements and lead to up- or downregulation of transcriptional activity. Although progress has been made in cataloging these CREAs, little is known about the rules by which enhancing and silencing silencing CRISPs interact with each other in order to shape gene expression levels. Experimentally, enhancer CRISPs can be mapped to their promoters using CRISPR interference or CRISPRi, where putative regulatory elements or pairs of elements are targeted for inactivation and gene expression changes are monitored.

However, these experiments are incredibly labor-intensive, and also higher-order interactions and redundancies within these CRISPs can complicate perturbation studies. Recently, the trend in deep neural networks predicting gene expression levels from the input DNA sequence is to expand the input sequence size from tens of kilobases to hundreds of kilobases. And the motivation for this is to predominantly incorporate more distal regulatory elements. For instance, Informer is trained to predict a large set of ENCODE epigenetic and CAGE coverage tracks from an input 200-kilobase sequence. Very briefly, CAGE quantifies transcription start site activity, and can be taken as a proxy for gene expression levels.

At the scale of 200 kilobase input, Informer includes the majority of mapped enhancer-gene pairs.

But because of the black box nature of this model, the authors of the paper had to validate that the model is making predictions for biologically meaningful reasons using interpretability methods. For example, they used in silico mutagenesis of a 200 base pair region around the disease-associated variant and revealed an affected transcription start site. However, it is difficult to understand the importance or significance of such motif-level, short-level explanations at 200-kilobase scale.

Another common approach for validating models is to replicate experiments in silico and compare to experimental results.

This has been extended to CRISPRi experiments as well. For example, here the authors have used attribution methods to assign importance to parts of the input and then compared the salient regions of the input to the experimentally validated enhancer-gene pairs. However, this method is of course limited by the availability of experimental data and therefore cannot tell us what the model has learned beyond known biology. Because of these two limitations, and to try to overcome this, we propose CREAM, a new interpretability framework for large-scale sequence-to-function deep neural networks. CREAM is an in silico sequence perturbation toolkit at the query scale And unlike the wet lab, CREAM performs experiments in silico through the lens of the deep neural network, allowing us to push the boundaries of the number and scale of different perturbation experiments we can perform.

But since this combinatorial space of the potential experiments is incredibly large, we have curated a set of different experiments that each addresses a specific biological question about the Cree grammar.

And we use the model as an in-silico laboratory, keeping in mind that the model shortcomings and biases will affect our conclusions.

Markers.

And then we apply the same experiment to a large set of sequences to find global patterns of gene regulation. To demonstrate the utility of Cram today, I will focus on K562 cell line and cage predictions as an output, and specifically targeting TSS activity as an outcome. Using Cram in our manuscript, we focus on a large set of questions, but today in the interest of time, I will only cover these four. And to start, we're going to address the question of how these regulators affect TSS activity in the model predictions. So to address this, we selected around 4,000 sequences centered at active genes. For each sequence, we dinucleotide shuffled the nucleotides everywhere in the sequence except for a 5-kilobase tile in the center of the sequence containing the TSS.

Then we compared the wild-type predictions versus the perturbed predictions.

We observed 3 scenarios depicted here. On the right, on top, the signal decreases after shuffling the sequence context, and that's shown in green. Therefore, the shuffling destroyed some enhancing sequences in the context, and that led to the decrease in the activity.

So refer to these sequences as enhancing contexts.

In the middle panel, the TSS activity virtually does not change after the shuffling, and therefore these sequences are robust to their sequence perturbations, and so we call this neutral context. And at the bottom, the predicted TSS signal increases after the context shuffling, which is indicated in the red color, suggesting that there were silencers in the context that got destroyed after the shuffling.

So refer to these sequences as silencing context.

And we have used these three sequence categories in subsequent experiments to understand the differences between these. So next, we asked how compatible different TSSs are with different sequence contexts. So to address this question, we tested the compatibility between the TSS and different contexts by pasting the central 5-kilobase TSS bin in various genomic sequences. For example, we took the CSS tile from enhanced, or tiles that we call enhanced TSSs, or TSSs from enhancing contexts, and pasted this into other enhancing contexts. So the results for this experiment are shown in green, and the normalized measure that we take is the mutant over wild type activity. On average, we observe lower than 0.5 activity fold change, meaning that the TSSs in this new enhancing context have more than twice lower activity compared to their original contexts. And this indicates that the original enhancers were kind of tuned to activate their specific TSSs. Potentially suggesting some sort of compatibility rules between enhancers and TSSs. And as a control, we also inserted the same TSSs into neutral and silencing contexts and observed, as expected, observed a drop in the TSS activity. We also repeated the same experiments with TSSs from neutral contexts, and again, as expected, observed virtually no difference in the activity. And just to note here, the majority of these TSSs are, or 70% of these are from housekeeping genes, and this is kind of a neat result to relate to biology as well. And then finally, we tested the TSS from silencing contexts. And in this case, when we look at the red box, we can see that the majority, or on average, the activity of these TSSs remains the same around 1, meaning that the silencers do not display the same kind of compatibility rules as enhancers and have more general silencing effect on the TSSs. And again, we as a control put these TSSs in enhancing and neutral contexts, and again, as expected, we saw higher activity in this because they kind of escape the silencing effect of their contexts. Okay, so moving on from this coarse-grained whole-context perturbations to higher-resolution, 5-kilobase-style perturbations, we can, you know, ask the question whether individual CREEs, which we define here as 5-kilobase sequence elements, are sufficient to explain the previously observed context effects.

So whether one tile of 5 kilobases can explain the effect of the whole context.

And to test this within a sequence of interest, we shuffled everything except for the central 5-kilobase tile, and then we systematically reinserted the original 5-kilobase tiles at their original positions and measured the difference between the wild-type sequence predictions and the perturbed sequence predictions. And on the right, the distribution of these effects is plotted grouped by sequence context. In enhancing context sequences, we mostly see positive effect CREEs, while in neutral and silencing contexts, there is a range of positive and negative effect CREEs. And surprisingly, we don't really see any individual CREEs that are reaching a normalized CRE effect of 1, meaning that none of them can individually on their own explain the whole effect of the context. And therefore, essentially, TSSs rely on multiple CREs to get regulated, I guess. And then, so in this experiment, we showed that there are multiple CREs acting on the same TSS. And so the next final question that we ask is, what set of CREs maximally enhance or silence TSS activity? And since it is unclear if the perturbation of one CRE affects the behavior of the other CREs in the sequence, we devised a greedy search strategy to find the optimal set of enhancing trees. In the greedy search for enhancers, we iteratively remove or shuffle the most enhancing tiles that drive TSS activity, and hence the removal of these tiles should lead to the largest drop in TSS activity because we're searching for the enhancing tiles. At the first iteration, we search for the tile perturbation that maximally drops the TSS activity, and then we keep this tile shuffled and record the normalized prediction of the TSS.

At the next iteration, we keep the CRE shuffled and run another round of search for the most enhancing CRE within the remaining tiles.

We continue the search, and for each iteration plot the normalized TSS signal. And here are the average results of the greedy search for enhancers, again for those 3 sequence categories that I defined earlier. In enhancing context sequences in green, on average perturbing 4 to 5 CREs is necessary to reduce the TSS signal by more than 80%. Interestingly, even the silencing and neutral contexts contain enhanced increase, but they don't have the same kind of effect on their respective TSSs, so the activity never actually drops to zero for those sequences. We repeated the same greedy search, but for silence increase in this case, using the same exact set of sequences, but now optimizing for—

oh, that didn't render right—

but optimizing for increasing the TSS activity, gradually removing the most silencing tiles. And interestingly, we can see that silencers have smaller effect sizes on average and work in larger groups. And again, in enhancing a neutral context, we see different silencing tiles as well. So together, it seems that Informer is integrating the net effect of numerous silencers and enhancers that exist in the same context. And so in the last set of experiments, we performed the greedy search before because we suspected that perturbing one CRE could affect the behavior of the other CRIs. So in this experiment, we're going to ask whether this is indeed the case or CRIs are acting in an independent additive manner. So in this comparison, we started with a hypothetical scenario where each child does act independently, and to calculate the model predictions in this hypothetical scenario, we added the independent effect sizes of CRIs in the same order as in greedy search. And to calculate the independent effect size of CRIs, we calculated the first iteration's effect sizes as the wild type minus mutant, and then subtracted these from the wild type in each step of the hypothetical additive search. In this example, toy diagram, the hypothetical additive prediction is the same as the greedy search results. And indeed, when we look at the average of the experimental results, we can see that Informer on average is integrating things in an additive manner, but when we stratify the results, we can see interesting cases of redundancy and cooperativity. So very briefly, in case of cooperative interactions, the simultaneous presence of both CRISPs is necessary to ensure the TSS activity. In this toy diagram, when both are present, the prediction is 1. When either one is removed individually, we observe a dramatic reduction in TSS activity, and this corresponds to a large individual tile effect size. For instance, the prediction of 0.2 means that the shuffle effect of that tile is 1 0.2, which equals to 0.8.

Using these values, we can compute the hypothetical additive effect size by summing individual effect sizes and subtracting from the wild type prediction, which is in this toy diagram just 1.

While actual predictions don't ever go below 0, the hypothetical additive model under cooperative interactions can lead to negative values, as you can observe in the experiments with Informer. And on the other hand, redundant enhancers by definition can compensate for each other. So in this case, shuffling either enhancer can lead to a small drop in the TSS activity, which translates to small effect sizes. And during the grid search, as the redundant tiles are shuffled, the predictions drop to zero, while the additive model never actually reaches zero.

And this is why the predictions are different essentially from each other.

So this scenario can happen in individual or even pairwise Creep perturbations in CRISPRi experiments, complicating the results of the characterization of CREEs. Okay, and with that, I wanted to essentially conclude saying that multiplex-in-silico perturbations are powerful tools to interrogate these large-scale deep neural networks, and CREAM is such a tool that provides a suite of different perturbation experiments to try to address different query grammar questions. Today I've talked about how we can interpret informers specifically to ask questions about how cooperative, redundant, or additive interactions can lead to a certain prediction. And then further observations that we have made in the paper include TSS, quidistance dependence effects, as well as the sigmoidal nature of gene expression.

And finally, we hope that similar methods can be used to guide wet lab experiments to essentially improve the efficiency of these.

And some of the future steps in this ongoing project are that we have already kind of tried to start implementing. In addition to the additive model, we also hope to incorporate multiplicative model in this final comparison. We have to characterize biochemical features and motif enrichment in the identified CREAs, which will hopefully shed some light on the compatibility rules between TSSs and CREAs. We also hope to perform a similar analysis at scales smaller than 5 kilobases, even going down to motif level. The complication here is that the search space for this expands dramatically, obviously, but we hope to use a nested experiment set in order to narrow down the motifs of interest within a tree. And finally, we hope to extend this work to further models such as Borzoi and PuffinD and other cell lines. And with that, I would like to thank the CoLab members and encourage you to check two posters from our lab, number 42 and 45, which are on rethinking initializations for neural nets and the new attribution method that uses interpretable surrogate models.

And thank you all for your attention.

Thanks for the great talk. I was wondering if you broke down any of your results by the class of genes that you're looking at. So, for example, like housekeeping TSSs might behave differently than developmental genes or very tissue-specific genes.

Great question. So, we have looked at the— if you remember the 3 categories of the context that I defined in the beginning of the talk, within those we have seen that the neutral context sequences are enriched for housekeeping genes.

So, like 70% of those, around 70%, are in that category. Which is kind of, I guess, what you would expect because those are supposed to be expressed all the time independent of the regulation.

Yeah. Yeah, that's pretty much the extent to which we have done that so far.

I also had a question, and then David. So related to that question, I guess in FORMER or any sequence-based model, the prediction ability on different sequences is very different. Some genes are very well predicted and others aren't. And I guess, that kind of uncertainty, can you incorporate that into your analysis, or do you think it affects your analysis?

Yeah, I guess, could I rephrase the question to, if the sequence is predicted with high accuracy, is that going to affect the—

yeah, I think it would definitely affect it. We haven't really looked specifically at, for example, cell lines that are predicted very poorly.

But this particular one on average is predicted very, you know, like the correlation is pretty nice with the experimental data. What we have done though is we've looked at the training test validation split and ensured that there's no differences in the conclusions essentially if we only look at one of the folds.

Yeah.

So the compatibility finding is really interesting to me. And I'm curious if you, between TSS and sequence context, I'm curious if you have any understanding of the mechanistic basis of that. Is it binding of like the same transcription factor about the TSS and in the context of CREES around it or something else?

Yeah, that's definitely an area of active research in our lab currently. Yeah.

We're hoping that motif-level analysis tools will make us, you know, like, understand that better.

Hi.

So I'm curious about, you know, you did this in silico mutagenesis experiments to kind of try and understand how they affect the outputs. Have you tried looking at some of the intermediate activations to gain, you know, a mechanistic in the computational sense understanding of how your model makes these predictions?

Yeah, we haven't really so far compared it with other methods, but we do hope to, especially at the motif level.

We are planning to compare our results to results from attribution methods, for example.

And we have done a very, very brief analysis of inserting, like, parts of the sequence that are identified by motif attribution methods back into the 5-kilobase files and seeing how those importance scores translate to changes in the model prediction. And it doesn't seem to be a very direct translation, essentially.

But yeah, I will not tell too much more because that is something that we're going to expand on. Thank you.

There is one last question, and maybe we can switch the presentations already while you Thanks for great presentation.

I have one question. So Informer has a lot of cell types at the output. So I wonder, did you explore if an enhancer is enhancing effect in some of the cell types and silencing effect in the other cell types?

We haven't really done that in this current version of version of—

going to extend the— and very quickly, it doesn't like choke some of the files. I haven't really seen any cases of like complete— but that's definitely something that we're still analyzing. Thank you for the question.