Scientist Stories: Reuben Saunders; Genome-scale Perturb-seq

So I'm Ruben, I'm a grad student in Jonathan's lab. And today I'm gonna be telling you about some work that we've been doing recently, most of which is now published, but then also specifically highlighting some sort of open analytical questions that I'm hoping will be interesting to the people in the MIA community. All right, not working.

Um, it's fine. Uh, at the risk of sounding like a, a broken record, the development of CRISPR-Cas9 has really been transformative for mammalian experimental genetics. Uh, you know, it's obviously a very powerful tool potentially for medical or therapeutic applications, but for the context of this talk, We really love CRISPR because it for the first time gives us easy access to the cell in terms of, you know, potent manipulations of gene function. And CRISPR tools have really been applied in two ways in experimental genetics. One is in the form of pooled CRISPR screens, which were developed in part at the Broad Institute by Eric Lander's lab. And these are experiments where you take a a pool of cells, and you usually infect them with a complex library of lentiviruses that express sgRNAs. And this usually leads to a situation where every cell gets one sgRNA, and then that sgRNA directs Cas9 to conduct some effector activity, which could be, for example, cutting a specific spot of the genome or recruiting transcriptional repressors or activators to a specific spot in the genome. And then this pool of cells is subjected to some kind of selection process, for example, growth or some sort of cell sorting process. And then finally, after the selection, the cells are— the sgRNAs are extracted from the cells and deep sequenced. And then you can calculate sort of an enrichment score for every sgRNA. And then the other big utility of CRISPR is the ability to finally easily, relatively easily, create isogenic cell lines. So this could be a context where you can take a cell type and you can introduce a specific mutation, validate that you made that specific mutation, and then you can analyze those cells by any of the sort of modern powerful methods in molecular biology or in cell biology, like RNA-seq or microscopy or proteomics. That said, there has been a technological gap which has been exposed by the really rapid advances both in these sort of breadth-first genetic screens and depth-first genetic screens. And that trade-off is that these breadth-first screens, they can assay many, many genes, but they tend to get quite simple phenotypes. And importantly, you have to be very explicit about what phenotypes you're going to screen on before you do your screen, because you just get like sort of a one-dimensional scalar score for every gene, and you really, you get what you screen for. You won't, it's your, some of the opportunities for serendipity are a little bit reduced.

And then the other—

on the other side, in these sort of depth-first methods, even now it's still pretty labor-intensive to create isogenic matched cell lines. It can take a couple months of time from a researcher. And so that really limits the scale of these experiments. So pretty often you can get these really complex high-dimensional phenotypes, but you tend to get them for a relatively small number of genes. And you really have to sort of prestate what genotypes you're interested in before you begin to step down your endeavor, also potentially diminishing the opportunity for serendipitous discovery. So the sort of the conceptual framework that I'm going to be talking about today really relies on a new category of methods that wrap—

that map many genotypes to rich high-dimensional phenotypes. And there are really some conceptually very different kinds of analyses that you can do when you can map many genotypes to many phenotypes— to many high-dimensional phenotypes or phenotypic elements.

And these could, for example, be like, you know, the RNAs in the cell or proteins or even like morphological features in the context of microscopy. And so, for example, one kind of analysis that you can do when you are measuring the impact of many genotypes on many phenotypes is you can understand how similar the phenotypes are of different genotypes. So for example, you can, you know, use some similarity metric across these phenotypes, and that can often give some insight into how functionally related the two genes whose phenotypes you're sort of messing with might be. The other, and one that I'm really actually quite excited about, is that you can directly probe how similar the phenotype elements are. So for example, if you have two RNAs and these two RNAs tend to covary across many perturbations, that might indicate something about their functional relationship. And this could also be true for two proteins or even two morphological features, the cell size and the nuclear size that's scaling together across many perturbations. Finally, when you have many genotypes and many phenotypes, you can begin to take steps toward thinking about the system of a cell as a sort of an availability of states in some high-dimensional space where perturbations, you know, can sort of drive the cells in some different position in this high-dimensional space.

And without getting sort of too get too technical about manifold claims or something like that.

I think we probably all kind of think that if you can understand all of the different states that can exist across a large variety of perturbations, you might be able to sort of learn the underlying structure that makes some states possible and other states impossible, which might give some fundamental insight into cellular function. I'm just trying to represent this concept as like an available genotype-phenotype surface of some sort. So over the past few years, a number of groups have developed methods that can potentially overcome this trade-off and record high-dimensional phenotypes for a lot of genotypes. And the method that I'm gonna talk about today was a method developed in part here at the Broad in Aviv's lab, in part in our lab and by a number of other groups. They were all kind of pushing toward the same idea, which is that you can take a pool of cells, infect them with a pool of sgRNAs, each giving the cell one genetic perturbation. And then you can do single-cell sequencing on this pool of cells. And then for every cell, you can recover both the identity of the guide RNA that was present in that cell, which is sort of analogous to the genotype of that cell. And then a single-cell transcriptome, which is analogous to a phenotype.

Although this method, you know, the power of this method was immediate— was immediately apparent when we started this current project, due to technical limitations, the scale of these experiments had been limited to only a couple hundred perturbations, only a couple hundred genes.

And as we kept doing relatively small-scale perturb-seq experiments targeting a couple hundred genes that we preselected from some previous screen or just from some literature review or something like that, we began to realize that there might be some emergent value to doing very large perturb-seq experiments to potentially create all genotypes in a single cell. And before I get started on the talk, I want to mention that this has been a really fantastic collaborative project in our lab with Joseph, an amazing MD-PhD student who's now back at UCSF completing his MD, Angie, a really awesome technician, and then Tom, who used to be a postdoc in the lab and now runs an amazing group at MSK and also Daphne. So just to get a bit into the experimental overview, What we did is we took a pool of cells that expressed a CRISPRi effector, dCas9-KRAB, which represses transcription. And we transduced them with a library of guide RNAs targeting essentially all expressed genes in a single cell type, in these K562 CML cells. And then also smaller sublibraries, one that was replicated in K562s and one in a very different cell type. RP1s. And then we did droplet-based single-cell sequencing by 10x, and what we ended up with was a gigantic matrix of millions of cells with known perturbation identities and then 20,000 genes whose expression we're measuring.

And really what I'm going to tell you about today is our efforts to sort of extract biological meaning from this gigantic matrix.

So there are a couple analyses that we can do. One is we can do some clustering analysis to try to summarize the expression programs and perturbation phenotypes. One is we can cluster perturbations by the similarity of the transcriptional response they induce to try to define gene function in an unbiased way. And then I'm going to talk a little bit less about this today, but one thing that we almost didn't expect is that there are sort of some other analyses that we could do beyond just this matrix per se, such as measuring splicing or total RNA content of the cells or cell cycle position, which we call composite phenotypes that we could sort of unexpectedly go back to our data and do computational screens on. So one thing that I said, which you hopefully flagged and wanted to call me on is that I said that we're using the sgRNAs as a reflection of genotype. And that's not strictly true. Really, all we're measuring when we measure the sgRNAs in a proteomic experiment is the presence of that sgRNA in a cell. And so one of the advantages of using the CRISPRi technology is that we are directly knocking down genes at the RNA level. And then we're using an RNA-seq readout. So we can actually confirm the linkage between our sort of predicted phenotype as measured by the fact that we detect that sgRNA inside the cell and the real— sorry, predicted genotype and the real genotype by actually directly measuring the on-target knockdown activity of each sgRNA. And so here I'm showing this in the two cell types. In K562s, we achieved median 85% knockdown across the about 10,000 genes that we targeted. The RPE1s, we used a slightly better dCas9 construct that absolutely we wish we had used for the KFASX2 experiments as well, but it wasn't available at the time.

And for the RPE1s, we achieved a median 92% knockdown of the genes that we were targeting.

This sort of let us establish that up to the limit of like 80%, 85% knockdown, measuring the presence of an sgRNA in a cell really does indicate that we made the sort of genotypic change that we think we made.

Good question.

How do you define knockdown?

Knockdown? Oh, yeah. So the number of counts of the target gene in cells that are targeted by that target gene versus the number of counts of that target gene in cells with guide RNAs that don't target any genes.

So any drop is counted as a knockdown, or you need to dip it all the way to zero?

Sorry. So one key thing to mention is that this is an average knockdown. So we're essentially averaging the expression of all droplets that contain a given sgRNA for that sgRNA targeting the gene, and then also for control sgRNAs. So this is a pseudo-bulk knockdown measurement. For genes that are pretty high expression, we can also look at like single-cell knockdown. And, you know, with the CRISPRi method, single-cell knockdown tends to be really quite homogeneous.

But obviously, you know, there are many genes that we're measuring here where the mean expression is less than 1 count per droplet.

And so for those, this is definitely a—

a pseudo-pulse knockdown.

Thank you.

Yeah, great question.

So, you know, there are many analyses that we could do with this data. One of the first that we wanted to do is to identify how many of our perturbations are actually causing any evidence of statistical differential expression in the single-cell RNA-seq. And so We tried to be extremely conservative with the statistics that we used throughout our analysis and to pretty much only use nonparametric physical tests, even at the cost of a little bit of physical power. And so with a very simple Mann-Whitney test for differential expression, we found that about 30% of the perturbations that we did targeting these 10,000 expressed genes caused a detectable transcriptional phenotype of more than 10 differentially expressed genes. And just as a bit of a bootstrapping benchmark, 2% of the non-targeting control guides that we used caused the phenotype at this 10-gene cutoff. One thing that we noted immediately is that perturbations targeting essential genes, so genes whose expression is required for the proliferation of K562s in this in this culture system were much more likely to cause a transcriptional response, especially a large transcriptional response, than expressed genes whose knockdown does not impact the growth of the cells.

You know, not totally surprising in some ways, kind of a nice benchmarking of the results of the experiment.

And so you can see, you know, we were able to robustly measure perturbations that induced, you know, 50 differentially expressed genes and also perturbations that induce 10,000 differentially expressed genes.

Although what exactly that means is something that we can talk about potentially during the discussion.

That's quite surprising, I guess, because on the previous slide you showed that many of the perturbations induced successfully knocked down the gene that they were targeting.

But I guess not that many relative—

maybe like half of those didn't have any wider—

Yeah, I mean, it kind of makes sense. You know, not every gene expressed in a cell type is required for the function of that cell in the specific culture conditions that we're expressing it in.

You know, there could be potentially for a lot of genes, you know, not that much sort of purifying selection to maintain their expression at a specific on or off level in a given cellular state. They're just kind of there and not not really doing anything. Or maybe, you know, if we were to expose the cells to some other perturbation and then do this experiment, we would find an effect there. But then in the sort of, you know, rich media with FBS condition that we were doing the experiment in, it wasn't detectable.

But yes, that's absolutely a very important takeaway of this experiment. You know, about 20% of the genes that are expressed in a cell type tend to be required for growth of that cell type in rich media conditions. And here, you know, 30% of the genes that are expressed cause any detectable transcriptional response when you perturb them.

So there are lots of genes that don't impact growth but do cause a detectable transcriptional response, but most don't impact either growth or a transcriptional response to the level that we can measure here in this one cell type.

So there— yeah, so once we had that, there were a number of sort of more matrix bioanalyses that we tried to do to sort of make sense of the experiment, make sense of these data. So one of the simplest was to cluster perturbation to define gene function. And so the way that we did this is basically to calculate the mean transcriptional profile for each perturbation versus internal non-targeted control cells. So we basically do this by averaging the technical z-normalized transcriptional changes relative to the distribution of gene expression in non-targeted control cells.

And then we compared the transcriptional profile for each perturbation versus the transcriptional profile for all other perturbations.

And one very quick benchmark of this was to take a look at genes or at perturbations, the gene that they target are both subunits in a sort of gold standard, you know, quaternary structure protein complex. Like in this case, we use the quorum dataset of protein complexes. And the idea here is that if you target two genes that are both sort of essential subunits of a protein complex, then targeting either of those two genes should basically cause the same transcriptional response. And we can, you know, assess that by similarity metrics.

So in this case, we're using Pearson correlation.

And so you can see that in fact, you know, gene perturbations that target two members of the same quorum complex tend to give much more highly correlated transcriptional responses than sort of random pairs of perturbations. So when you actually do this, and this in this case I'm showing this analysis limited to genes that cause more than 50 differentially expressed genes because those tend to sort of cluster better, you end up with this highly structured similarity matrix. Here I'm kind of showing it as a cluster map and What's really fun to do is you can look along the diagonal of this cluster map and you can kind of keep zooming in and really see a tremendous amount of structure in these data. And then when you finally get to the level where you can actually see individual gene names, if you're a researcher scrolling through this on your laptop, it's really gratifying to see that a tremendous amount of known biology, of known cell biology, is apparent when you zoom in far enough. So along the diagonal, you basically end up with dozens to hundreds of sort of core cellular processes, genes that are known to be involved in core cellular processes, all appearing clustered together, even if there's no obvious reason to think that targeting this particular process would cause a specific transcriptional response. Clearly the transcriptional responses that we're measuring are accurate enough fingerprints.

They can put together the components of the COP9 signalosome or the STAGA complex.

Question.

So how do you understand that in the sense that— OK, so let's say that I take— if I take any of the member of the genes in the integrator complex, and if I CRISPRiA. Why do you see the same transcriptional phenotype? Is it like, how do you understand it? Is it because that elicits the same sort of feedback mechanism in the cell?

Yeah, so I actually want to move away from the integrator complex because I'm going to get back to that in a second, and I think you'll understand why I moved away from it.

But let's take the stargardt complex or something like that.

So yeah, I think in a lot of these cases, You know, you literally have this protein structure that's only stable if you have all of the pieces present.

You just— it may be competitively degraded by the proteasome because it has a lot of exposed hydrophobic surfaces if you don't have all the pieces present. And so, you know, targeting any one piece basically leads to the whole function of the complex being removed from the cell. And then what we're measuring here is the cellular response at this point I don't want to make an adaptive claim about that response, but just the fingerprint response to the depletion of that complex.

Yeah, so I think in a lot of cases it's literally causing the same effect, which is like the degradation of that complex, or at least the absence of function of that complex from the cell. And just because this is the Broad Institute, I thought I should show some sort of dimensionality reduction embedding. So this is a minimal distortion embedding of all of the different perturbations that have more than 50 differentially expressed genes, where the embedding makes an effort to put perturbations that induce similar phenotypes nearby to each other. It's kind of conceptually analogous to UMAP or t-SNE.

I like the representation because it gives a global sense of the complexity of the data. A bit better than heat map.

And I should say we have a website which has an interactive version of this that you can zoom in on or highlight, and it can be a fun way to explore the relationships between genes with all, you know, abstractions and issues with dimensionality embedding as an obvious caveat. So the reason why I was cautioning you about the integrator complex is that the integrator complex had a really fascinating and unexpected behavior where that the subunits of the integrator complex, which have basically been defined by proteomics, they clustered as 3 basically completely distinct modules. And this was really quite a unique behavior. You know, in almost every situation, all of the subunits of a complex kind of clustered together, but that was not the case for the integrator. For the integrator, we had one complex that contained the subunits 10, 13, 14, and also this uncharacterized gene C7R26. There was this other cluster that had 2, 5, 7, and 8. And then this other cluster that contained 3, 4, 9, 11—

I guess not 11 in this version of the clustering—

plus some other genes that were known to be involved in snRNA biogenesis. And there was really remarkably low intercluster correlation. So one thing that we checked is that we just wanted to confirm that this uncharacterized gene, C7RP26, we wanted to know whether it really was a core component of the integrator complex. And so we were able to co-IP it with the integrator complex, at least with this specific 10, 13, 14, 15 component. And also we actually purified integrator 10, 13, 14, and also CUS26 over size exclusion column. We actually, we've renamed CUS26 as INTEF15, and you know, one kind of cool thing about this is that people have had difficulty, you know, purifying this part of the integrator complex before, and that's probably because they were trying to purify it without this component, INTEF15/CUS26, And so from this genetic prediction based on RNA-seq phenotype, we were able to specifically identify the protein structure that this uncharacterized gene was in.

Now to give it off to the structural biologists and the biochemists and have them run with it, develop some sort of mechanism. And so here I'm showing again this pretty low inter-cluster correlation relative to the intra-cluster correlation.

And one thing that was pretty interesting to do was to compare these modules that we define by RNA-seq phenotype with the structure of the integrator complex, parts of which have recently been identified with cryo-EM. And so we could see that these components here in this tight cluster were in what the structural biologists had called the backbone builder module of Integrator, whereas the 43911 MEDIC cluster was in this endonuclease module of the Integrator, which actually contained the known active site of the Integrator complex. And then there was no structure of the 10-13-14 C7O26 module. So we actually used AlphaFold to make a prediction of these 4 complexes together. Take that as a caveat.

AlphaFold prediction looks pretty darn good, but it is absolutely a model.

And AlphaFold predicts that this structure forms quite nicely, looking a little bit like Barney. Dinosaur, and maybe some biochemical evidence that it suggests that it may sit kind of here in the complex. But, you know, the reason we find this exciting is that, you know, here we have this biochemical complex which, you know, structural biologists love, like looking at a protein and saying that it has specific domains. But here we can really say that in a totally unbiased way, totally unbiased sense, these different domains of integrator have different functions. In the sense that disrupting these different domains causes like very different transcriptional responses. So I think it really motivates sort of future work to understand why these different genetic modules of the integrator induce such different responses.

They must really have different functions. And I actually think, you know, the function of this 10, 13, 14, 15 module really is quite distinct from, from the other proposed functions of integrator that have been present in the literature.

So we can certainly cluster perturbations to define gene function. We also tried to summarize expression programs. So we use these CRISPRi perturbations to drive cells into hundreds or thousands of extremely diverse transcriptional states. And kind of our hope is that we can analyze these states together to really learn the core transcriptional modules that are available to a cell, at least with like a specific cell type like these K562 cells. And our dream is really to sort of decompose each transcriptional response into these core modules. And so here I'm showing an extremely naive example, like the response to perturbation P1, we can say it's like alpha times this generic core Myc response. Plus beta times this generic core HSF1 response.

These are transcription factors, plus gamma times this generic GATA1 response, which is a cell identity transcription factor in hematologic cells. And there are many, many mathematical methods that we could potentially use for this decomposition.

There are also many different forms of the data that we could try to apply this decomposition to. We could try to apply it to the single-cell transcriptomes, but those have the issue that they're quite noisy. We don't sequence these data that deeply just based on cost. We could apply it to the pseudo-bulk transcriptional changes, so the average transcriptional changes in response to each perturbation.

And we could do things like PCA or ICA or various nonlinear methods that we could use to try to do this decomposition.

And so the reason that I'm including this section of this talk at this seminar is that we have not found a good way to do this decomposition. Way that simultaneously decomposes the responses in a mathematically elegant way, but also yields human-interpretable components of these core modules. You mentioned ICA or this type of forecasting.

So what is the result with the line ICL?

In every case, you can explain most of the variants, but you can't interpret— at least I have not been able to find some unifying feature of the components that helps me present them as like the MYC module or the HSF1 module or something like that. But you can certainly explain the variant with all of these methods. And so, you know, just as an extremely simple example, we could take the suitable transcriptional responses and we can do PCA. And we can explain, you know, 80% of the variance with 20 principal components, which I do think is like, is fundamentally interesting from a question of like how many important dimensions in gene expression space are there in like a single cell type. But then when we actually began to do like manual inspections of the PCs or sort of more sophisticated inspections of the PCs by like gene set analysis or something like that, we really find them hard to interpret. At least in this case, all of the PCs, the top loading in positive and negative directions is basically like the same genes that are just like remixed, the same 30 genes that are just remixed across all of these different PCs. And so here I'm just showing that for this gene CDKN1A or CXCLH, although it's a little bit cut off for people who are here in the audience, that this is not the MYC component and the ATF1 component and the GATA1 component. —which again, I'm not saying that we'd necessarily expect it to be because obviously there are a lot of pretty strong limitations to a method like PCA, and we don't expect this to be a linear combination of modules—but nonetheless, we have not come up with a really elegant way to do this decomposition. And so, because decomposition has been challenging, we've pursued approaches that use clustering instead. Which is a little bit less limiting in a sense than a decomposition method. So the core idea here, you know, was really first formulated by Mike Eisen and David Botstein in the '90s, is that genes that are co-regulated across many different cellular states may form a functional module. So, you know, Mike did this on yeast expression data from microarray back in the back in the '90s. And our hope is that because we have, you know, this RNA-seq data across these, you know, 10,000 diverse perturbations, and crucially because these are sort of internally controlled measurements that actually have very low noise, maybe we can really identify the core modules based on this analysis. That said, one interesting feature of these data is that the transcriptional responses to perturbations are not that clumpy, and so clustering is not necessarily so easy to do. So I'm going to compare that to perturbation with the kind of like, yeah, you know, the same analysis on basically the transpose of matrix. So if you take the perturbation-perturbation correlation and you make a cluster map, you can see that really this incredible amount of clustering along the diagonal. And then if you take this and embed it to the neighbor-adjacent, you can see very distinct clumps of perturbations that behave— that cause very, very similar transcriptional responses, and so therefore have very high levels of inter-cluster correlation relative to the intra-cluster correlation. And when we did the same analysis on the transpose, we find that there is less clumpiness of the correlations. So there is absolutely some structure along the diagonal, but there is much less, except for this one large cluster, which represents the ribosomal proteins that basically covary.

And then this other cluster, which represents the mitochondrial genes that basically covary.—

if we then do this embedding, there's certainly structure, but there are only a couple of clusters that really clearly sort of stick out to this very obvious analysis. And although these are all real and meaningful, they're not sort of representative of like most genes. Instead, there are specific subsets of genes that have specific features that make them so clumpy. And so I just want to give one example of these clusters. Culture of lincRNAs. So when we actually manually inspected this, people who have spent time looking at RNA-seq datasets, this is like the weirdest set of gene lists. Like whenever you see a gene list like this that have all these antisense transcripts or these ACE transcripts or something like that, you know something odd is taking place biologically. And what's going on here is that these are basically like weird RNAs that are transcribed and then efficiently degraded by the nuclear exosomes. And so when you deplete subunits of the exosome, all of these weird RNAs accumulate. And then this is what they'll read map to when you actually do your alignments. And so the reason that these genes are so highly correlated with each other and therefore form this clump is that knockdown of any of the subunits of the exosome causes this super specific transcriptional response where these genes all get activated by like 30 V4s or something like that relative to non-targeting control cells. So We can do clustering and it absolutely is, it's a very meaningful analysis, but there are limitations of clustering. For example, clustering doesn't allow genes to be in multiple clusters most of the time. So if we're doing clustering, we find like a cluster that represents the integrated stress response, which is one of the core cellular responses to stress and a cluster that we call as the unfolded protein response, there are absolutely genes that seem to really— and here I'm just showing correlation of each transcript in the experiment with the centroid of the UPR module and the centroid of the ISR module. And you can see that there are some genes that— and then they're colored by whether the clustering algorithm which we use in this case, which is like HTBScan, defines these genes as really being in these modules. And you can see that there are absolutely genes that are very much UPR genes. And there are genes that are very much ISR genes, but there are all these genes in between that seem to be activated both under conditions that activate the UPR and under conditions that activate the ISR. And, you know, one limitation of the clustering approach is that it just ignores these genes because they don't sort of like behave in this very discrete and unique way. So, oh yeah.

So if I were to speculate, and I would like to see what's your response to this, could Could be because CRISPR perturbations are very different from the type of innate stimuli that cells receive in the environment that they're evolved to actually show a response to. Absolutely.

Yeah, I think that's extremely important. So I'm going to try to get to that. It's a little bit of a difficult claim to make, except to say that clustering doesn't work as well here as it does in like, you know, single-cell organism RNA-seq data or something like that. And that it seems like most of the states that we're achieving through our CRISPRi perturbations are super specific to the perturbation that we're causing, but they may not necessarily reflect like a specific negative feedback loop or a specific adaptive module.

And so that can make them very difficult to sort of classify or decompose into sensible quantities. So I believe that. But I'm not totally sure how to move forward with that. If you do this PCA analysis and you compare two cell types, it's excellent at identifying the genes that distinguish those two cell types. And those are often core transcription factors or something like that. But within the cell type with this set of CRISPR perturbations, which are essentially all the different cellular states that you can achieve with single-cell perturbation, with single gene perturbations, Yeah, it's a little bit more difficult to get at.

And yeah, so the approach that we took, which we're not in any way arguing is optimal, but did nonetheless yield significant insight, was basically to try to find a space where the neighbors are a little bit closer together. And this has some similarities and methods like weighted gene network correlation analysis or something like that. So in this case, basically we use the graph Laplacian, and here I am showing Euclidean distance in this space constructed from the graph Laplacian—

sorry, the eigenvectors of the graph Laplacian— and then taking that and embedding the neighbors who are adjacent.

And so this is moving cells closer, moving transcripts closer to their neighbors, and then doing clustering.

—which obviously will improve the quality of the clustering, but may also increase the false positive rate of clustering. And so when we do this, we do end up with nicer clusters, which we can collect with some method like identity-based clustering.

And we ended up with 37 different transcriptional response clusters. About half of which we could look at them and identify some clear unifying biological feature of that clustering, whether it's like the, you know, this cluster that is really like the activation of proteasome subunits or myeloid differentiation or erythroid differentiation or ribosome biogenesis. But, and these are actually all the different sort of perturbation clusters. And so we can see, which perturbation clusters activate which or repress which sort of transcriptional response clusters.

But more than half of the transcriptional response clusters, we couldn't find any unifying biological feature. It doesn't mean they're not real.

It just means that they escape obvious human or gene set analysis.

Yeah, so that's why we're sort of like still searching for the best way to understand these highly specific transcriptional responses.

That said, we have found that once we've defined these modules in this sort of completely unbiased clustering-based manner, we can then go and look back at the perturbations and see which perturbations activate each of these transcriptional modules. And so one cool feature of this analysis for us, since our —our lab often works in cellular stress responses—is we could sort of in a completely data-driven way re-identify the distinct unfolded protein responses and the integrative stress responses and find the drivers of the UPR and the ISR. So for example, we found that this gene HspA5, which is really essential for protein folding in the ER, strongly activates the UPR.

Whereas this gene HspE1, which is a mitochondrial chaperone, when knocked down strongly activates the integrated stress response.

And K562s, they're kind of part of the way down a hematologic differentiation hierarchy where they're kind of stuck because they're leukemia cells. And we found perturbations that drive them toward an erythroid phase and also perturbations that drive them to a more myeloid phase. And one thing that we're quite excited to do in the future is do these experiments in other cell types like pluripotent stem cells or HSCs where that occupies sort of highly biologically meaningful states and then find the perturbation, find all those different directions which they go. Since this kind of differentiation phenotype often includes a large set of genes that all change together, and so this module is very easy to detect with this kind of analysis. So I want to get back a little bit to what you were asking. So the genes that best distinguish perturbations are not necessarily the genes that are in these modules. So something that, you know, we did pretty often as we were analyzing these data was to sort of train classifiers to learn gene sets that best distinguish perturbations, since we have all these sort of labeled cells and we can try to find which genes when we look at them, like, can best classify cells as, like, labeled with this perturbation versus labeled with other perturbation, and then check that against the true list of the guide RNAs. And something surprising is that there's the sort of usual suspect list of so-called housekeeping transcripts that almost always do a great job of distinguishing cells based on which perturbation they have. And this is really surprising because these housekeeping genes are like the controls in like GPCR experiments.

These are the things that are supposed to not change across perturbations, yet they're doing an excellent job of classifying cells between perturbations.

And so here I'm just showing like the genes that best distinguish the subunits of the integrator complex. And you can see that— sorry, the label may have disappeared a little bit— but you can see that a lot of these are like, you know, ribosomal proteins or vimentin or something like that, which are not necessarily—

like, vimentin is like a component of the cytoskeleton. These are not the pieces that you would expect would change in response to perturbations.

You know, for example, the ribosomal proteins, they're proteins that form a stoichiometric complex, like the large and small subunits of the ribosome.

And so it would be very surprising if changes in the levels of these RP transcripts were doing a good job of like distinguishing between different perturbations.

That said, for reasons that I would say we don't necessarily understand, these housekeeping genes do an excellent job of distinguishing between most perturbations. So here I'm showing, you know, an embedding similar to one that I was showing earlier that's done on 2,300 very—

like, the most variable genes we're measuring in our data. I'm sorry, it's cut off here, but for people in the live audience, you can see like the RPU one showing 1,100 perturbations that cause more than 50 differentially expressed genes. And here I'm looking at just the 89 radical protein genes.

So these are genes that if you're analyzing your data, it's like RPL4 or RPS19 or something like that. And it's worth noting that most of these genes are in a global sense very highly correlated with each other. They're much more highly correlated with each other than they are with other random genes. In the genome. And yet there's still enough dimensionality in the expression of these 89 RP genes that we can reconstruct like, you know, 50 different clusters of perturbations, 35 of which are essentially identical to the clusters that we can reconstruct when we use all genes in the transcriptome.

And we do not understand this phenomenon.

That said, For people who are doing like tissue atlas experiments or something like that, something that we've seen is that these same housekeeping genes also do an excellent job of distinguishing between cell types and tissue atlas experiments. I don't understand why.

And people often ignore these kinds of transcripts when they're doing like single-cell RNA-seq atlases of tissues.

Yet nonetheless, there's a tremendous amount of signal.

That's encoded in these expressions.

So the random forest classifier, does it show— do those housekeeping genes show the same predictive value for all perturbations? Have you looked at the accuracy of predicting individual perturbations?

So I mean, the precise gene set is a little bit different between different perturbations. But as I'm kind of showing here, you know, across all of these different extremely diverse categories of perturbations, these RPs, which are just a small subset of the housekeeping genes, do an excellent job of prediction.

As I said, there are absolutely some cases where they don't do as well as sort of more specific gene sets.

I'm just wondering whether it could be— I mean, there's a lot of details here, like For example, if you're training this on only the perturbations that show the very strong phenotypes, that's correlated with probably core functional genes.

And many of them are probably part of key complexes like the ribosomal complex. Yeah, yeah, yeah, that's absolutely true. So they do well both for genes that cause very strong responses and genes that cause less strong responses.

But I'm absolutely not making the claim that they do better than more specific subsets of genes in all cases, just that based on the fact that these are like qPCR controls, I would like to say that they don't have any predictive value.

Right. Yeah.

So to kind of sum this up, again, in a way that's hopefully more thought-provoking than satisfying, some perturbations cause transcriptional responses that reflect really obvious transcriptional programs.

Or feedback loops, you know, things that may be adaptive transcriptional responses.

But at least in this cell type, in this culture context, most strong perturbations cause modest but, you know, specific fingerprint changes in the expressions of many genes, including these housekeeping genes. And these changes really are sort of the most distinctive changes that are downstream of each of these perturbations.

And so this is why if you read our paper carefully, we basically say that what we're doing is we're clustering to summarize expression programs. We don't want to make a stronger claim about what we're actually doing in this analysis than that at this point, even though our long-term goal would really be to link this back to transcription factors. So one other analysis that we can do— the analysis that I've been showing up to this point are mostly these pseudo-bulk analyses.

But this is a single-cell RNA-seq dataset. And so one fun form of the analysis is we can look at how homogeneous the response, the single-cell responses to perturbations are. And so the specific way that we did this is with these actually PCA-based leverage scores. So it's, you know, you take the first 20 components of a PCA of this 2 million by 20,000 dimensional matrix. And then you compare the loading, you compare like the values of these 20 components for cells with a given perturbation. It could be like 200 cells that are targeting WDR3 or something like that versus non-targeting control cells.

And you actually, we tend to normalize these for the distribution of leverage scores and non-targeting control cells, make it a little bit easier to think about and to have entered at 0.

So we find that most perturbations tend to cause pretty homogeneous transcriptional responses, but some perturbations tend to cause much more variable transcriptional responses. And we actually just classify this with like the standard deviation of leverage score. And so something that stood out to us is that when we look at the number of differentially expressed genes, which correlates almost perfectly with the sort of average leverage score,— and then we compare that to the standard deviation of the leverage score, there are some genes that really stand out as having a very high number, a very large amount of sort of single-cell variability relative to the mean transcriptional change. And this I'm just kind of including because it's fun, and also it was really a brilliant analysis from my colleague, Joshua Elkington, that I really did not believe was going to work despite some significant sort of somewhat similar analysis from the Broad. So some of these genes that cause these highly heterogeneous responses, what Joseph realized is that they weren't actually causing like a mean transcriptional change. Instead, they were increasing the odds that cells acquire karyotypic abnormalities during replication. And the way that he realized this is that, you know, normally you look at like gene expression changes sorted by like, you know, ensemble gene ID or something like that. But then when you resort them by like position in the genome, you find that if you look at like a single cell, so each of these is like the xenormalized gene expression of a single cell, you would find like one chromosome where like all the genes go down by 50% or something like that, or all the genes go up by a factor of 2. The interesting thing about this is that if you look at two individual cells, they have different karyotypic abnormalities. And so they have, you know, there's no mean transcriptional change in response to that perturbation, but each of them has a very distinct sort of individual single-cell transcriptional response. And so here I'm just showing —this is like gene expression in this set of non-targeted control cells versus this gene PPK, which is like a checkpoint kinase, which is known that when you inhibit it or when it goes away, you accumulate a lot of karyotypic abnormalities. And so you can really see that every one of these cells has a very different transcriptional profile when you look at it across the genome. But every one of these cells has also gained or lost basically at least one chromosome and sometimes more.

And so one really cool feature of this that we would have never imagined when we were actually doing the experiment is that we could go through the data, you know, months after the experiment was done and do basically a retrospective genome-wide screen for chromosome instability.

And, you know, we actually did this both in RPE1s and K562s, which are quite different from a cell biological standpoint, but we found genes that significantly increased the odds of developing karyotypic abnormality in both RPE1s and K562s. And these genes tend to be genes involved in the cell cycle or in DNA repair.

And so we're still working on how to do broad analysis of the single-cell aspect of this data versus the pseudo-bulk analysis. In part, the single-cell data are just—

they're a lot slower to work with in most cases than the pseudo-bulk analysis. But this is, we thought, was an example of like an exceptionally unexpected single-cell analysis that, you know, led to some surprising biological insight.

And also within, in this case, like the first screen for the genes that when you knock them down cause karyotypic abnormalities, which I'm not really sure how you would screen for that in a different way. So in the last few minutes, I just want to talk quickly about the transcriptional response to mitochondrial stress. I think this is particularly relevant for people who are doing single-cell sequencing analysis themselves.

So the mitochondria, they have really diverse functions in like respiration, metabolism, iron-sulfur biogenesis, and also they're actually quite important in apoptosis. And, you know, the endosymbiosis took place at the origins of eukaryotes where an alpha proteobacterium was taken up by the sort of progenitor eukaryotic cell. And over time, most of the genes that were encoded on the alphaproteobacteria genome were transferred to the nuclear genome.

But 13 of these protein-coding genes remain in the mitochondrial genome.

So the nucleus in like a human cell encodes more than 1,000 mitochondrial proteins. And these proteins are synthesized from the cytoplasm by the cytosolic ribosomes.

And then they're actually imported into the mitochondria across either one or two mitochondrial membranes.

And, you know, for sort of reasons related to our particular interests, we were very interested in understanding the coordination of the nuclear genome with the state of the mitochondria. So how does the nuclear genome sense the state of the mitochondria and then respond to specific mitochondrial stressors? And And kind of surprisingly to us, we found that there was basically a generic nuclear transcriptional response to all mitochondrial perturbations.

So basically whether you were targeting like protein import into the mitochondria or complex I respiration or complex II respiration or ATP synthase, which is complex V, or the mitochondrial ribosome or something like that, there was basically a generic nuclear response to that perturbation. Which I'm sort of illustrating here as like just a high level of Pearson correlation between the nuclear transcriptional responses for these mitochondrial perturbations. And really, this response was almost entirely explained by the integrated stress response, which is caused by EIF2 alpha phosphorylation.

And that's actually sent by kinases that directly monitor the status of the mitochondria. And this was surprising because there are non-mitochondrial perturbations that also cause strong ISR activation.

And on the nuclear genome, we couldn't even identify which perturbation targeted mitochondrial genes versus which perturbations target other processes that when disrupted also activate the ISR.

And something that was very surprising to us is that we began to actually inspect the mitochondrial genes themselves. And we found that although there was this sort of very generic response in the nuclear genome to mitochondrial perturbations, these 13 mitochondrial genes had in many cases exquisitely specific responses to the disruption of different processes in the mitochondria. So for example, we could easily see that the disruption of like mitochondrial DNA polymerase has a different response from the disruption of the large subunit mitochondrial ribosome versus complex I or II or III or IV, just kind of a set of representative mitochondrial perturbations that are then clustered through, you know, through a set of data in a data-driven way.

And there's essentially no literature precedent for these kinds of mitochondrial gene changes in response to mitochondrial perturbations.

For mitochondrial genes, the polyadenylation is not exactly the same as the polyadenylation of nuclear genes.

And so we wanted to make sure that this was a real transcriptional change that we were measuring and not some weird feature of the fact that we're using poly-A selected mRNAs and that we're actually measuring like changes in total gene level and not changes in like the specific levels of polyadenylated transcripts.

And so we actually did a bulk RNA-seq experiment that targeted a subset of perturbations and found—

and there was no poly-A selection in this bulk seq experiment.

And we found it's a correlation of 0.8 between the single-cell changes and the bulk changes. So these changes to mitochondrial gene expression are absolutely real.

And, you know, we really were struck at how much better mitochondrial genes were doing at classifying mitochondrial perturbation by the function of the gene they targeted than the nuclear genes were.

So here I'm just showing, you know, essentially all mitochondrial targeted perturbations that induce a significant differential expression response. And you can see that there's very, very little structure when you look at the correlation, pairwise correlation on all nuclear genes. But when you look specifically at mitochondrial genes, there is this really complex pattern of clustering with really a tremendous amount of structure. And when you analyze it, it really nicely lines up with known mitochondrial gene functions.

And in fact, we were able to make some predictions for uncharacterized genes, which specific mitochondrial process they were involved in, which we then identified. And we then proved that those were correct with other methods.

Is clustering transcriptional phenotypes, but only in the space of like 13 mitochondrial genes?

That's right. Yeah. So these 13 Pearson correlation of changes in these 13 genes versus Pearson correlation of changes of 22,000 nuclear genes.

Thank you. Yeah.

So there's much more signal with just these 13. Than there is with 2,000.

And I'm running a little bit out of time, but I'm essentially almost done.

What I do want to mention is that this is really surprising from a cell biological standpoint because unlike nuclear genes, the mitochondrial genome is polycistronic. There is one transcript that encodes 12 of the 13 mitochondria encoded genes. And they're transcribed essentially on a single transcript going all the way around. And then there is this post-transcriptional cleavage and polyventilation that leads to all of the individual genes. And then for God knows what reason, one of the mitochondrial genes is on the antisense version of the transcript and is transcribed separately.

And so what this means is that so far as we know, there is no opportunity for transcriptional control of mitochondrial genes, only for post-transcriptional control of mitochondrial genes through selective stabilization and degradation. And the complexity of that sort of post-transcriptional response, there's just no literature present for it. And we're very curious to understand what precise mitochondrial mechanisms could actually be contributing to that post-transcriptional complexity.

Just to sum up this part, we sort of define a new aspect of the architecture of mitochondrial stress, which is that, you know, essentially all disruptions to mitochondrial activity, at least in the cell types we looked at, cause ISR activation, but that there is specific mitochondrial transcriptional signatures to each of the different specific kinds of perturbation, whether that's like disrupting complex V or the large subunit or the integrated subcomponents. And something that is fun to think about, but we don't at this point have any evidence for or against, is whether there's autonomous regulation of these mitochondrial genes.

So there are many copies of the mitochondrial genome in each of the different mitochondria. And if you have independent local conditions in each mitochondria—

this one has an issue with complex V, or this one has an issue with protein folding—

Maybe you can more specifically tailor your response to local conditions if you can respond locally. Although again, that's really speculation at this point.

So just to sum up, this kind of genotype-phenotype mapping is really powerful at scale because there are certain kinds of analyses that really only work when you have this sort of large and comprehensive dataset.

More perturbation leads to quadratic increase in the number of pairwise comparisons and therefore, in general, much more power for the kind of clustering that we do. Rich phenotypes can drive serendipitous discovery, such as the ability to go back and look at genes that drive chromosomal instability.

I didn't get into too much here, but because this is an internally controlled experiment and we can essentially treat every cell as an independent experiment, in our analysis.

We have very high-end statistics, and therefore actually we have quite a lot of like capital P power in our analysis, but that there's still really tremendous room for analytical development to help make sense of these, these really complex datasets. Yeah, and just to sum it up, you know, perturb-seq is really great at identifying genes with similar functions based on transcriptional programs, and also at identifying new unexpected processes such as the independent regulation of mitochondrial genes. But there's a lot of need for new analytical methods to actually understand the causal mechanisms that link perturbations to the transcriptional changes. And so with that, I want to really thank Jonathan, who was an excellent mentor for this project, my fantastic colleague Joseph Elkington, Tom and Angie, and funding. And yeah, thank you so much for your attention.