Systematic morphological profiling of human gene and allele function via Cell Painting

1The dramatic increase in human genome sequence data has created a significant bottleneck. The number of genes and variants known to be associated with most human diseases has increased dramatically. Unfortunately, the next step - understanding the function of each gene and the mechanism of each allele in the disease - typically remains non-systematic and labor-intensive. Most commonly, researchers painstakingly design, develop, and apply a disease-specific or biological process-specific assay.

2Over 30% of genes in the human genome are of unknown function and even annotated genes have additional functions yet to be uncovered. Furthermore, even when a gene’s normal functions are known, methods are lacking to predict the functional impact of the millions of genetic variants found in patients. These gaps must be filled in order to convert the promise of human genome sequence data into clinical treatments.

3Therefore, there is a widespread need for systematic approaches to functionally annotate genes and variants therein, regardless of the biological process or disease of interest. One general approach depends on guilt-by-association, linking unannotated genes to annotated ones based on properties such as protein-protein interaction data, sequence similarity, or, most convincingly, functional similarity. In the latter category are profiling techniques, where dozens to hundreds of measurements are made for each gene perturbation and the resulting profile is compared against profiles for annotated genes. Various data sources can be used for profiling; gene expression is one that can be performed in relatively high-throughput and it has been proven useful in predicting gene function. In fact, high-throughput mRNA profiles were recently used to cluster alleles found in lung adenocarcinoma based on their functional impact, a precursor to therapeutic strategy for variants of previously unknown significance.

4Images are a less mature data source for profiling but show tremendous promise. Morphological profiling data is complementary to transcriptional profiling data and is less expensive. Morphological profiling has succeeded across several applications, including grouping small-molecule perturbations based on their mechanism of action, and grouping genes based on morphological profiles derived from cells perturbed by RNA interference (RNAi). One limitation of RNAi for morphological profiling is that the number of measurements must be limited or else the resulting profiles are dominated by off-target effects, especially seed effects. Some computational solutions have shown some promise in overcoming this problem for gene expression profiling, but their utility is unproven for image-based profiling, and regardless RNAi does not permit analysis of gene variants, only knockdown. Modification of genes via CRISPR will require new libraries of reagents and is as yet untested in morphological profiling.

5In the proof-of-concept work presented here, we tested morphological profiling using overexpression in human cells as a general approach to annotate gene and allele function. We profiled a reference series of well-known genes, and a small number of variants thereof, by Cell Painting. In particular, we wondered whether the information content of this strategy would outweigh potential limitations (e.g., due to cellular context or expression level). We found that the approach successfully clustered genes and alleles based on functional similarity, revealed specific morphological changes even when present in only a subpopulation of heterogeneous cells, and uncovered novel functional connections between important biological pathways.

6To profile each exogenously expressed gene (or allele therein), we used our previously developed image-based profiling assay, called Cell Painting. This microscopy-based assay consists of six stains imaged in five channels and revealing eight cellular components: DNA, mitochondria, endoplasmic reticulum, Golgi, cytoplasmic RNA, nucleoli, actin, and plasma membrane (Figure 1A). In five replicates in 384-well plate format, we infected U-2 OS cells (human bone osteosarcoma cells), chosen for their flat morphology and previous validation in the assay, with an arrayed ‘reference’ expression library of 323 open reading frame (ORF) constructs of partially characterized functions (Supplementary file 1A), a subset of which have been previously described. Of these, we prioritized analysis of the 220 constructs that were most closely representative of the annotated full length transcripts (see Materials and methods). Morphological profiles were extracted using CellProfiler for image processing, yielding 1384 morphological features per cell, and Python/R scripts for data processing, including feature selection and dimensionality reduction (Figure 1B, and see Materials and methods). This computational pipeline yielded a 158-dimensional profile for each of 5 replicates for each gene or allele tested.

7Not all genes are likely to impact cellular morphology given the limitation of our experiment; using a single cell line at a single time point under a single set of conditions and stained with six fluorescent labels. We therefore first asked what fraction of these ORFs impacted morphology. Surprisingly, we found that 50% (110/220) of these ORF constructs induced reproducible morphological profiles distinct from negative control profiles (Figure 2A, and see Materials and methods). Next, we ruled out the possibility that position artifacts may have artificially inflated this result by taking an alternative pessimistic null distribution which takes well position into account (Figure 2—figure supplement 1). Therefore, we conclude that a single ‘generic’ morphological profiling assay can detect signal from a substantial proportion of genes in our reference set. We next turned to testing whether those signals are biologically meaningful and can lead to novel, unbiased discoveries about gene function.

8Given that technical replicates produce similar morphological profiles, we next evaluated whether similarities between profiles induced by different constructs are meaningful. We began with the simplest case: for a subset of genes in the experiment, a ‘wild-type’ sequence (see Materials and methods for important definitions) was captured in more than one ORF construct (23 pairs). These pairs either correspond to different physical cloning events and preparations but with highly similar full-length sequence (as defined in Methods; category a: nine pairs), or a substantive difference in their nucleotide sequence, for example, isoforms (category b: 14 pairs). We found that, as expected, the phenotypes of over-expressed wild-type ORFs of the same gene were more similar to each other, on average, than to randomly selected genes. Of the 23 pairs for which both wild-type ORFs yielded a phenotype distinguishable from negative controls, 22 (~96%) of the pairs’ profiles were correlated more than expected by chance (Figure 2B, the one pair not meeting that threshold was in category b), confirming that different constructs with biological similarity indeed produce similar morphological profiles.

9This result also confirms that the sequence differences seen in separately cloned wild-type constructs do not generally have a major functional impact, but we caution that any individual construct of interest may have an impactful mutation; thus the raw sequence data should be examined and testing alternate constructs for a gene may be recommended. Note that if, for example, only 50% of wild-type pairs showed high profile correlation, it would remain ambiguous whether it was caused by poor assay quality or by constructs’ sequence mismatches. But in this particular case the mentioned near perfect consistency rules out either of the two possibilities. We also note that the 23 pairs analyzed here are located in different well locations on each plate; this result therefore also rules out widespread artifacts, such as plate position effects or metadata errors.

10We suspected that the small number of engineered constitutively activating alleles for certain genes would, on average, yield a stronger phenotype than their wild-type counterparts. We indeed found that correlations between replicates of the constitutively activating allele were typically higher than correlations between replicates of the wild-type version of a gene (Supplementary file 1B; p-value=0.012, one-sided paired t-test).

11We hypothesized that genes in pathways known to affect cellular morphology (RAC1, KRAS, CDC42, RHOA, PAK1, and genes related to the Hippo pathway) would be more likely to yield a morphological phenotype distinguishable from negative controls than other genes in the analysis. Indeed, we found this to be true (Fisher’s test p-value=3.7 × 10−3) (Figure 2C). Reassured by this validation, we were curious which pathways would be most and least likely to yield detectable morphological phenotypes, recognizing that ‘pathways’ are neither separate nor well-defined entities. We found genes manually annotated as being in the Hippo, Hedgehog, cytoskeletal reorganization, and Mitogen-activated protein kinases (MAPK) pathways were more likely to result in a phenotype, whereas genes annotated as belonging to the JAK/STAT, hypoxia, and BMP pathways were among the least likely to yield a phenotype under the conditions tested (Figure 2—figure supplement 2 and Supplementary file 1C). Nevertheless, the majority of pathways could be interrogated by morphological profiling.

12Given the caveats and limitations of overexpressing genes (see Discussion), we next tested whether image-based profiling of expression constructs could capture relationships among genes known to be functionally related. Because a reliable and complete map of all gene-gene connections is not available, we evaluated the accuracy of our results via two approaches.

13First, we compared our data to protein-protein interaction data from BioGRID. This is imperfect ground truth for judging our predictions because two proteins might physically interact without producing the same morphological phenotype when overexpressed, and genes in the same pathway might regulate the same phenotype without any physical interaction. Nevertheless, we expect that the corresponding proteins of gene pairs with highest profile similarity are more likely than average to physically interact. Indeed, looking at wild-type versions of genes showing a detectable phenotype (the 73 genes represented in the 110 constructs), the ratio of verified gene connections among the top 5% correlated gene pairs (9%, 13 verified out of 143 possible combinations) is significantly higher than that of other gene pairs (5%, 128 verified out of 2485 possible; Fisher’s test p-value=0.04; Supplementary file 1D).

14Second, we manually annotated each gene for the pathway with which it is associated. This approach is based on expert opinion and thus imperfect knowledge of all genes’ function; furthermore many pathways interrelate, and genes in the same pathway are not expected to have identical phenotypes given that their functions are rarely identical (most notably, overexpression of some may activate while others suppress a biological pathway or process). Nonetheless, we expect pairs of genes whose morphological profiles correlate highly to be more likely than average to be annotated in the same pathway vs. different pathways. Using the same 73 genes as in the previous analysis, the ratio of gene connections with the same-pathway annotation in the top 5% most-correlated gene pairs was 20% (29 pairs out of 143), significantly higher than the ratio for the remaining pairs (6%, 139 pairs out of 2485; Fisher’s test p-value =; Supplementary file 1E).

15Having quantitatively established that morphological profiling is sensitive, robust, and captures known gene-gene relationships, we explored these relationships in a correlation matrix (Figure 3 bottom left and Figure 3—figure supplement 1). The overall structure, with multiple groupings along the diagonal, is consistent with the fact that the 110 constructs (73 unique genes) that showed a phenotype had been annotated as representing 19 different pathways. That is, we did not see large, homogeneous clusters, as would be expected if morphological profiling was sensitive to perturbation but not highly specific. This rules out uniform toxicity induced by a large number of genes, for example. Neither did we see only signal along the diagonal, which would have indicated no strong similarity between any gene pairs.

16We next created a dendrogram (Figure 3) and defined 25 clusters (see Materials and methods and Figure 3—figure supplement 2) to explore the similarities among genes. Pairs of wild-type ORFs almost always clustered adjacently, consistent with our quantitative analysis described above (Figure 2B). After retaining only one copy of replicate ORFs, we found that the majority of clusters (19 out of the 22 clusters containing more than one gene) were enriched for one or more Gene Ontology terms (Supplementary file 1F), indicating shared biological functions within each cluster.

17Using this dendrogram, we began by interrogating three clusters that conformed well to prior biological knowledge. First, we analyzed Cluster 20, containing the two canonical Hippo pathway members YAP1 and WWTR1 (more detail in Supplementary file 2 [PDFs A2–A20 and B2–B20 ], and in a later section of the text). Both are known to encode core transcriptional effectors of the Hippo pathway, and a negative regulator of these proteins, STK3 (also known as MST2), is the strongest anti-correlating gene for the cluster (Supplementary file 2 [PDF A20], panel c1).

18Second, we noted Cluster 21 is comprised of the two phosphatidylinositol 3-kinase signaling/Akt (PI3K) regulating genes, PIK3R1 and PTEN, both frequently mutated across 12 cancer types in The Cancer Genome Atlas (TCGA). These results are consistent with previous observations that certain isoforms of PIK3R1 reduce levels of activated Akt, a dominant negative effect. AKT3 is in a cluster anti-correlated to the Cluster 21 ((Supplementary file 2 [PDF A21, panel b1]).

19Third, we examined three clusters (19, 6 and 3) that included many MAPK-related genes. Cluster 19 is the largest example of a tight cluster of genes already known to be associated; it includes four activators in the RAS-RAF-MEK-ERK cascade: KRAS, RAF1 (CRAF), BRAF, and MOS. Notably, two constitutively active alleles of these genes, BRAFV600E and RAF1L613V, form a separate cluster (Cluster 6) adjacent to their wild-type counterparts. Furthermore, the constitutively active RAS alleles HRASG12V and KRASG12V are in the next-closest cluster (Cluster 3), which also contains MAP2K4 and MAP2K3 (known to be activated by Ras), as well as CDKN1A. By contrast, MAPKs that are known to be unrelated to the RAS-RAF-MEK-ERK cascade, such as MAPK14 in Cluster 5, are far away in the dendrogram.

20Overall, these results support the notion that connections between genes can be efficiently discovered using our approach.

21We hypothesized that the specific morphologic features that segregated each of the clusters would provide insight into gene function. Examining images (Supplementary file 2 [PDF -19A, panel 3]) or rank-ordered lists of features that distinguish individual profiles or clusters (Supplementary file 1G) is tedious and lacks sensitivity for all but the most obvious of phenotypes, confirming that quantitative morphological profiling is more sensitive than the human visual system.

22We therefore devised several strategies to enhance biological interpretability from these experiments and applied these in combination. First, we grouped features into meta-features based on their type of measurement, i.e., shape, texture, intensity, etc., and the cell constituents to which they are related, to create a Feature Grid (Figure 4A). Second, we performed unsupervised grouping of features by mapping the top 20 most-distinguishing features for each cluster onto a plane, creating a Feature Map (Figure 4B), in which highly correlated features are mapped nearby each other (see ‘Feature Interpretation’ in Methods for an explanation of individual feature names). In certain cases, these visualizations revealed the nature of the morphological phenotype (e.g., nuclear shape abnormalities Supplementary file 2 [PDF 7A]), but for others these approaches did not suffice to yield an obvious phenotypic conclusion (e.g., for Cluster 19, Figure 4A and B).

23Third, we hypothesized that leveraging the single-cell resolution of image-based profiling might be highly sensitive in enhancing interpretation, particularly for cases where only a subset of cells is distinctive from negative controls. To test this, for each given cluster of genes together with negative controls we identified 20 subpopulations using k-means clustering on single cell data. We calculated the abundance of cells in each of the 20 subpopulations to determine which are over/under-represented relative to controls for the given cluster (corresponding images are shown; Supplementary file 2 [PDFs 1B–25B]). For example, the MAPK pathway activators in Cluster 19 show increased prevalence of a subpopulation of cells with strongly asymmetric ER, mitochondria, and Golgi staining, indicating a cell polarization phenotype (Figure 4C, and Supplementary file 2 [PDF 19B], Categories one and two), for which there is evidence in the literature. This phenotype was not captured by manual inspection nor the first two approaches (e.g., Supplementary file 2 [PDF 19A, panels a2 and b2]).

24Encouraged by this, we supplemented the morphological map by compiling these and other visualizations into PDF files for each cluster, summarized in Figure 5 and provided in full as Supplementary file 2. We also noticed that certain subpopulations were similar across several clusters (Figure 3—figure supplement 3 shows sample cell images of each such subpopulation); we annotated their enrichment/de-enrichment on the dendrogram (Figure 3).

25Using these visualizations, we began by interrogating three adjacent and correlating clusters (Clusters 4, 7, and 11) contain wild-type and mutant alleles of CDC42, a gene encoding a Rho family GTPase with diverse roles in cell polarity, morphology, and migration. Cluster 4 contains the constitutively active mutant CDC42 Q61L as well as MAP3K2 and MAP3K9. The highly similar Cluster 7 contains the dominant negative alleles CDC42 T17N and RAC1 T17N, a related RAS superfamily member. That activating and inhibiting alleles would yield similar phenotypes when overexpressed is not surprising for CDC42. Cluster 7 also contains isoforms and alleles of AKT: specifically, AKT3 and the constitutively active E17K alleles of both AKT1 and AKT3. Akt is known to be essential for certain Cdc42-regulated functions and vice versa. Finally, the nearby Cluster 11 (which is discussed in more detail later) contains the wild-type form of CDC42 as well as TRAF2, a canonical NF-κB activator; these two are known to interact and share functions in actin remodeling. We also note that anti-correlating genes to these clusters (generally in Clusters 13 and 21) are consistent with existing knowledge, including (a) AKT family member AKT1S1 (a Proline rich AKT substrate, PRAS40, Supplementary file 2 [PDF 7A, panels b1 and c1]) (b) CDK2 (a known target of Akt), (c) PIK3R1 and PTEN in Cluster 21, described previously, which have known interactions with AKT. Thus, all of these connections have previously been identified.

26Subpopulation visualization revealed that Clusters 4, 7, and 11 are enriched in cells that are huge and binucleate (Figure 3, example images shown in Supplementary file 2 [PDF 4B]). Genes in all three clusters also show irregularities in DNA content, namely, an enrichment in cells with sub-2N DNA content, a decrease in cells with 2N DNA content, and, for most genes, a decrease in cells with S phase and 4N DNA content, indicating a significant amount of DNA fragmentation and thus apoptosis (DNA histograms in Supplementary file 2 [PDFs 4B, 7B, and 11B]). These phenotypes are consistent with these genes’ known role in the cell cycle and cell polarity.

27As a second test case, we examined Cluster 8, which contains PRKACA (the catalytic subunit α of protein kinase A, PKA) and two of its known substrates: GLI1 (a transcription factor mediating Hedgehog signaling), and RHOAQ63L (a Ras homolog gene family member). The highly similar Cluster 10 contains the wild-type RHOA, as well as ELK1 which is also linked to the Rho GTPase family and PKA.

28We investigated the morphological changes causing these genes to cluster. RhoA is a known regulator of cell morphology and cell rounding is a known related phenotype. We found that indeed all members of Clusters 8 and 10 significantly induce cell rounding (Supplementary file 1H). Although cell count is lower for genes Clusters 8 and 10, the degree varies greatly (from z-score −0.67 to −3.02, Supplementary file 2 [PDFs 8A and 10A, panel a1]), ruling out that simple sparseness of cells explains their high similarity in the assay. As well, the overall DNA content distribution of the cell populations appears relatively normal (Supplementary file 2 [PDFs 8B and 10B]). Subpopulation extraction provides a satisfying biological explanation for these clusters’ distinctive phenotype: the increased roundness and strong variation in intensity levels (per the Feature Grid) across the population stems from an increased proportion of telophase, anaphase, and apoptotic cells (Figure 3 and Supplementary file 2 [PDFs 8B and 10B]).

29We therefore conclude that the morphological map can link related genes to each other and that the morphological data can provide insight into their functions, particularly with the help of subpopulation visualization.

30We wondered whether novel relationships might emerge from our unbiased classification of gene and allele function based on morphologic profiling. We noticed that the known regulator of NF-κB signaling, TRAF2 (in Cluster 11, together with CDC42), yields a signature strongly anti-correlated to YAP1/WWTR1 (Cluster 20), which encode the transcriptional effectors of the Hippo pathway, YAP (Yes-associated protein) and TAZ (Transcriptional co-activator with a PDZ-domain). The Hippo pathway and NF-κB signaling are critical regulators of cell survival and differentiation, and dysregulation of these pathways is implicated in a number of cancers, but we found no evidence in the literature (in particular through BioGRID) of physical interaction between the proteins encoded by Cluster 11 genes and Cluster 20 genes. Confirming our approach, a functional connection between CDC42 (Cluster 11) and YAP1 (Cluster 20) has been identified: deletion of CDC42 phenocopies the loss of YAP1 in kidney-specific conditional knockouts in mice. Still, the NF-κB pathway (and in particular the Cluster 11 member TRAF2), has not been closely tied to YAP and TAZ in human cells (see Discussion).

31We first wanted to characterize Clusters 11 and 20 to confirm that relationships within each cluster are supported in the literature. Indeed we found evidence for most of the within-cluster connections. CDC42 and TRAF2 (Cluster 11) physically interact and share functions in actin remodeling. As described in a prior section YAP/TAZ (Cluster 20) are known to share functional similarities in the Hippo pathway, being regulated by, and also regulating, cytoskeletal dynamics. Consistent with these known functions, we found that a core effector of the Hippo pathway which functions to restrict YAP/TAZ nuclear activity, STK3 (which encodes the Mst2 kinase), has a morphological signature strongly anti-correlated to YAP1/WWTR1 (Supplementary file 2 [PDF 20A, panel c1]). We note that although STK3 and TRAF2 are both moderately anti-correlated with YAP/TAZ (Cluster 20), STK3 and TRAF2 are not themselves highly correlated, indicating each has a different subset of phenotypes that anti-correlate to YAP/TAZ. We also note that two clones that express another regulator of YAP activity, STK11, form Cluster 22 which falls nearby YAP1/WWTR1; a connection between STK11 and YAP has been identified (albeit with opposite directionality, identified via knockdown of STK11). Further, YAP1 is among the highest anti-correlating genes to REL (data not shown; REL is a singleton in the dendrogram and thus not in a cluster), whose protein product, c-Rel, has a known connection to TRAF2. These results reaffirm that the Cell Painting-based morphological signatures are a useful reporter of biologically meaningful connections among genes in these pathways.

32Given the striking inverse correlation between YAP1/WWTR1 and TRAF2, we sought to confirm a negative regulatory relationship between the Hippo and NF-κB pathways by multiple orthogonal methods.

33First, we explored the observed inverse morphological impact using the Cell Painting data. The morphological impact of genes in Cluster 11 and 20 is quite strong (median replicate correlation is at the 74th and 81st percentile, and average within-group correlations are 0.66 and 0.73). Subpopulation analysis showed that Cluster 20 (YAP1, WWTR1) is enriched for cells that are slightly large, slightly elongated, and have disjoint, bright mitochondria patterns, whereas Cluster 11 (TRAF2, CDC42) is de-enriched for those subpopulations and instead enriched for binucleate cells, very large cells, and small cells with asymmetric organelles (Figures 3, 6A and B).

34Second, given that YAP/TAZ are transcriptional regulators, we analyzed gene expression data. Using the same constructs as in our Cell Painting experiment, we found an anti-correlated relationship at the mRNA level, consistent with the anti-correlation we had seen in morphological space. To do this, we used Gene Set Enrichment Analysis and publicly available data, which includes data from four to nine different cell lines at one to four time points (https://clue.io). Time point refers to the duration of treating the cells with over-expression constructs until the time gene expression readouts are made. This analysis revealed that the NF-κB pathway is the pathway most enriched among genes whose overexpression results in down-regulation of known YAP1 targets, CTGF, CYR61, and BIRC5 (Benjamini and Hochberg (BH) adjusted p-value = in Supplementary file 1I, and Figure 6C), with TRAF2 being among the genes contributing to this enrichment (Supplementary file 1I). We also saw enrichment of NF-κB pathway members when testing a data-driven set of targets of YAP1/TAZ (Figure 6—figure supplement 1, see Materials and methods). In the inverse analysis, genes that alter the levels of TRAF2/REL common targets are weakly enriched in Hippo pathway members (Figure 6—figure supplement 2, see Materials and methods). This is consistent with the hypothesis that NF-κB members can downregulate YAP/TAZ targets but not strongly vice versa.

35Finally, we more directly confirmed negative crosstalk between NF-κB effectors and YAP/TAZ using a synthetic TEAD luciferase reporter that is YAP/TAZ responsive. Importantly, these confirmatory experiments used different cellular contexts and perturbation constructs versus the original Cell Painting data. Co-expression of the NF-κB pathway effectors TRAF2 or C-REL with YAP or TAZ led to significantly lower reporter activity than expression of YAP or TAZ alone (Figure 6D1 and D2). Intriguingly, mutants of YAP or TAZ that are insensitive to negative regulation by the Hippo pathway (YAP-5SA and TAZ-4SA;) remained sensitive to suppression of transcriptional activity by TRAF2 and C-REL, indicating that the negative relationship we identified may be independent of canonical upstream Hippo pathway signals (Figure 6D3 and D4).

36We conclude that connections among genes can be profitably analyzed using morphological profiling of overexpressed genes via the Cell Painting assay. In a single inexpensive experiment, we were able to rediscover a remarkable number of known biological connections among the genes tested. Further, we found that morphological data from the Cell Painting assay, together with novel subpopulation visualization methods, can be used to flesh out the functionality of particular genes and/or clusters of interest.

37By adopting a two-pronged approach, merging this Cell Painting morphological analysis with transcriptional data, we were able to identify an unexpected relationship in human cells between two major signaling pathways, Hippo and NF-κB, both under intense study recently for their involvement in cancer. Through validation of these clustered genes, we have identified that YAP/TAZ-directed transcription is negatively regulated by NF-κB pathway effectors and our data suggests a novel regulatory mechanism that is independent of upstream Hippo kinases.

38To date, there has been little evidence of the intersection between these important signaling pathways. Recent work examining osteoclast-osteoblast differentiation has suggested that Hippo pathway kinases, such as Mst2, may affect the NF-κB pathway through phosphorylation of IkB proteins, thereby promoting nuclear translocation of NF-κB transcription factors. TAZ was found to be a direct target of NF-κB transcription factors and its expression is regulated via NF-κB signaling. Our work, however, supports a possible additional mode of interaction, whereby regulators of NF-κB signaling directly regulate the function of Yap and Taz as transcriptional co-factors. Recent work has demonstrated, in Drosophila, that NF-κB activation via Toll receptor signaling negatively regulates the transcriptional activity of Yorkie, the homolog of YAP/TAZ, through activation of canonical hippo pathway kinases. The work described here identifies, for the first time in a mammalian system, that a negative regulatory relationship exists between NF-κB activation and YAP/TAZ transcriptional function. Furthermore, we have identified that this regulation of YAP/TAZ occurs in a manner that is independent of Hippo pathway-mediated phosphorylation events on YAP/TAZ, suggesting a more direct relationship between NF-κB and YAP/TAZ signaling.

39In this work, we tested quantitatively and explored qualitatively the connections among genes revealed by morphological profiling. Our underlying hypothesis was that functionally similar genes would generally yield morphologically similar cells when overexpressed, and indeed we found this to be the case. Still, some discussion of this point is warranted. Most commonly, gene overexpression will result in activation of the corresponding pathway via amplification of the endogenous gene’s function. However, it is important to note that the profiling strategy to discover functional relationships does not assume or require this. For example, overexpression could also disrupt a protein complex, producing a trans-dominant negative effect that results in precisely the opposite phenotypic effect. In still other cases, overexpression of a particular gene may not affect any of the normal functions of the gene (producing a false negative signal), or trigger a stress response (yielding a confounded profile), or produce a complicated response, due to feedback loops. Further, artifactual phenotypes could be seen, e.g., if overexpression yields a non-physiological interaction among proteins or toxic aggregates. Nevertheless, despite these caveats and complications, our results indicate that valuable information could be gleaned from the similarity and dissimilarity of the morphological perturbations induced by gene overexpression. Using overexpression avoids the complications of RNAi off-target effects (often due to seed effects), which were far more prevalent (impacting 90% of constructs in our recent study).

40In addition to functionally annotating genes, as demonstrated here, one particularly appealing application enables personalized medicine: it should be feasible to use morphological profiling to predict the functional impact of various disease alleles, particularly rare variants of unknown significance. This has recently been successful using mRNA profiles. Thus, an even more exciting prospect would be to combine mRNA profiles with morphological profiles to better predict groups of alleles of similar mechanism, and ultimately to predict effective therapeutics for each group of corresponding patients.

41We make all raw images, extracted cellular features, calculated profiles, and interpretive visualizations publicly available, providing an initial morphological map for several major signaling pathways, including several unexplored connections among genes for further study (see Supplementary file 2). Expanding this map to full genome scale could prove an enormously fruitful resource.

42The Reference Set of human cDNA clones utilized here has been previously described; ~90% of these constructs induce expression of the intended gene greater than two standard deviations above the control mean. Briefly, wild-type ORF constructs were obtained as Entry clones from the human ORFeome library version 8.1 (http://horfdb.dfci.harvard.edu) with additional templates generously provided by collaborating laboratories, and cloned into the pDONR223 Gateway Entry vector. In addition, here, to maximize coverage of cellular pathways, we included additional clones with minimal sequence deviations from the intended templates. Sanger sequencing of Entry clones verified the intended transcripts and, if applicable, the intended mutation. Entry constructs and associated sequencing data will be publicly available via www.addgene.org and may also be available via members of the ORFeome Collaboration (http://www.orfeomecollaboration.org/), including the Dana-Farber/Harvard Cancer Center (DF/HCC) DNA Resource Core DNA Repository (http://www.dfhcc.harvard.edu/core-facilities/dna- resource/) and the DNASU Plasmid Repository at ASU Biodesign Institute (http://dnasu.asu.edu/DNASU/Home.jsp). Clone requests must include the unique clone identifier numbers provided in the last column of Supplementary file 1A (e.g. ccsbBroadEn_12345 as an example for a specific entry clone and ccsbBroad304_12345 as an example for a specific expression clone). ORFs were transferred to the pLX304 lentiviral expression vector by LR (attL x attR) recombination.

43For simplicity, throughout this paper ‘wild-type’ refers to ORFs found in the original collection without a particular known mutation intentionally engineered. Due to natural human variation, and occasional cloning artifacts, there are often non-identical matches of such constructs to reference sequence; these differences are fully documented for each construct and sequence data will be publicly available through AddGene, in addition to the sequencing data for the original Entry clones for the genome-scale library.

44U-2 OS cells (human bone osteosarcoma cells), RRID:CVCL_0042, were obtained from ATCC and propagated in the William Hahn lab; they were not additionally authenticated prior to this experiment. The cell line tested negative for mycoplasma prior to this experiment. HEK293T cells, RRID:CVCL_0063, were obtained from ATCC. The cell line was validated by STR profiling (Genetica DNA Laboratories) and was negative for mycoplasma as measured by MycoAlert Mycoplasma Detection Kit (Lonza, Walkersville, MD).

45We followed our previously described protocol except for durations of some steps. Briefly, cells were plated in 384-well plates and transduced with lentiviral particles carrying ORF constructs the next day. Viral particles were removed 18–24 hr post-infection and cells cultured for 48 hr until staining and imaging (72 hr total post-transduction). The experiment was conducted in five replicates, each in a different plate. The number of replicates being five was decided based on prior experiments.

46The Cell Painting assay followed our previously published protocol. Briefly, eight different cell components and organelles were stained with fluorescent dyes: nucleus (Hoechst 33342), endoplasmic reticulum (concanavalin A/AlexaFluor488 conjugate), nucleoli and cytoplasmic RNA (SYTO14 green fluorescent nucleic acid stain), Golgi apparatus and plasma membrane (wheat germ agglutinin/AlexaFluor594 conjugate, WGA), F-actin (phalloidin/AlexaFluor594 conjugate) and mitochondria (MitoTracker Deep Red). WGA and MitoTracker were added to living cells, with the remaining stains carried out after cell fixation with 3.2% formaldehyde. Images from five fluorescent channels were captured at 20x magnification on an ImageXpress Micro epifluorescent microscope (Molecular Devices): DAPI (387/447 nm), GFP (472/520 nm), Cy3 (531/593 nm), Texas Red (562/624 nm), Cy5 (628/692 nm). Nine sites per well were acquired, with laser based autofocus using the DAPI channel at the first site of each well.

47The workflow for image processing and cellular feature extraction has been described elsewhere, but we describe it briefly here. CellProfiler software version 2.1.0 was used to correct the image channels for uneven illumination, and identify, segment, and measure the cells. An image quality workflow was applied to exclude saturated and/or out-of focus wells; six wells containing blurry images were excluded, retaining 1914 plate/well combinations in the experiment. Cellular morphological, intensity, textural and adjacency statistics were then measured for the cell, nuclei and cytoplasmic sub-compartments. The 1402 cellular features thus extracted were normalized as follows: For each feature, the median and median absolute deviation were calculated across all untreated cells within a plate; feature values for all the cells in the plate were then normalized by subtracting the median and dividing by the median absolute deviation (MAD) times 1.4826. Features having MAD = 0 in any plate were excluded, retaining 1384 features in all. The image data along with the extracted morphological features at the per-cell level were made publicly available in the Image Data Repository under DOI http://dx.doi.org/10.17867/10000105.

48The code repository for the profiling and all the subsequent analysis is publicly available at https://github.com/carpenterlab/2017_rohban_elife (with a copy archived at https://github.com/elifesciences-publications/2016_rohban_submitted). We will next explain details of each analysis step implemented in the code. Single cell measurements in each well and plate position are summarized into the profiles by taking their median and median absolute deviation (abbreviated as ‘MAD’ or ‘mad’ in some tables) over all the cells. Although this method does not explicitly capture population heterogeneity, no alternate method has yet been proven more effective. We also include the cell count in a sample as an additional feature. This results in a vector of 2769 elements describing the summarized morphology of cells in a sample. We then use the median polishing algorithm after obtaining the summarized profiles, to remove and correct for any plate position artifacts. For each feature, the algorithm de-trends the rows, i.e. by subtracting the row median from the corresponding feature of each profile in that particular row. Next, it de-trends the columns in a similar way using column medians. The row and column de-trending is repeated until convergence is reached in all the features. For the rest of the analysis we considered only the constructs which have more than 99% sequence identity to both the intended protein and gene transcript, to avoid testing uncharacterized mutations/truncations.

49Not all of the morphological features contain useful reproducible information. We first filter out features for which their replicate correlation across all samples (except the negative controls) is less than 0.30, retaining 2200 features. Subsequently, a feature selection method is used. Briefly, starting with features (measurements) that we identify as essential, a new feature that contributes the most information with respect to those that have been chosen, is added to the set. The contribution of each feature to the already-selected features is measured by the replicate correlation of the residue when the feature is regressed on the already selected features. This is repeated until the incremental information added drops below a threshold. The original method proposed in overfits in its regression step when the original data is very high dimensional. As a remedy, in the regression step we only use features that have a Pearson correlation of more than 0.50 with the selected features thus far. This prevents overfitting of regression when the dimensionality of selected features grows. We stop feature selection when the maximum replicate correlation of residue is less than 0.30.

50The feature selection method greatly removes redundancy, but because of the non-optimal ‘greedy’ strategy, some redundancy remains. Principal component analysis is then applied to keep 99% of variance in data, resulting in 158 principal components being selected.

51The features measured using CellProfiler follow a standard naming convention. Each feature name is made up of several tokens separated by underscores, in the following order:

52Prefix which could be either empty or ‘mad’. This means that the feature is calculated either by taking median (no prefix) or median absolute deviation (‘mad’ prefix) of the relevant measurement over all the cells in a sample.

53Cellular compartment in which the measurement related to the feature is made, i.e., ‘Cells’, ‘Cytoplasm’, or ‘Nuclei’. Note that features labeled ‘Nuclei’ are based on segmentation of nuclei using Hoechst staining, ‘Cells’ are based on segmentation of the cell edges using the RNA channel, and ‘Cytoplasm’ is the subtraction of the aforementioned compartments.

54Measurement type, which can be either ‘Intensity’, ‘Texture’, ‘RadialDistribution’, ‘AreaShape’, ‘Correlation_Correlation’, ‘Granularity’, and ‘Neighbors’. Note that ‘Correlation_Correlation’ measures, within a cellular compartment, the correlation between gray level intensities of corresponding pixel pairs across two channels (specified in the next tokens in the feature name). Note also that the relative positioning of a cell is measured in the ‘Neighbors’ category.

55Name(s) of channels in which the measurement is made, if appropriate (omitted for AreaShape and Neighbors).

56Feature name. The precise measurement name appears at the end. A description of each metric can be found in the CellProfiler manual (http://cellprofiler.org/manuals/current/)

57Our method to identify which genes produce a discernable profile involves first normalizing each profile to the negative controls, such that a treatment’s median replicate correlation becomes a surrogate for phenotype strength. In the case that a treatment does not show a phenotype different from the negative control, its replicates would center around the origin in the feature space. This would consequently decrease the median replicate correlation. On the other hand, a phenotype which is consistently observed in the replicates and is significantly different from the controls results in the replicates to concentrate in a region far from the origin in the feature space, and hence a high median replicate correlation value.

58The cutoff for ‘discernible’ is set based on the top fifth percentile of a null distribution. The null distribution is defined based on the correlations between non-replicates (that is, different constructs) in the experiment. Treatments whose replicate correlations are greater than the 95th percentile of the null distribution are considered as ‘hits’ that have a morphological phenotype that is highly reproducible (Figure 2A).

59At this point, for strong treatments, all profiles of the replicates are collapsed by taking the average of individual features. 110 out of the 112 selected ORFs were significantly different from the untreated profiles in the feature space. That is, their average Euclidean distances to the untreated profiles were higher than 95th percentile of untreated profile distances to themselves. This shows these two alternative notions of phenotype strength–replicate reproducibility and distance to negative control–are consistent. We restrict all the remaining analyses to the 110 ORFs.

60In this analysis, mutant alleles were removed and we considered only one wild-type allele for each gene with a detectable phenotype, retaining 73 genes. We calculated a threshold to identify significantly correlated gene pairs. We picked the threshold to minimize the probability of error in classifying wild-type clone pairs versus different-gene pairs. To do so, we found the value at which the probability density functions of the two groups intersect; this value (here, 0.43) can be proved to have the desired property. This approach results in about 5% of the gene pairs being categorized as highly correlated. We next formed a two by two contingency table, where the rows correspond to two groups of gene pairs, determined by whether they have high profile correlation or not. Similarly, the columns also correspond to two groups of gene pairs, determined by whether the corresponding proteins have been reported to interact in BioGRID (or alternatively have been annotated to be in the same pathway; Supplementary file 1C and 1D). This table was then used to perform a one-tailed Fisher’s exact test.

61A dendrogram was created based on the Pearson correlation distance and average linkage, using the hclust function in R (Figure 3).

62Gene clusters were formed by cutting the dendrogram at a fixed correlation level, 0.522, which was chosen using a stability-based measure. The measure is defined as follows: the local clustering stability is measured for a range of candidate cutoffs, from 0.43 (used earlier to test consistency to protein interaction data) to 0.70. The point with highest stability was chosen (Figure 3—figure supplement 2), and the stability measure was defined as the proportion of treatments whose clusters do not change if the cutoff is slightly changed by a small amount,.

63In order to extract cell categories (subpopulations) and subpopulation enrichment laid over the dendrogram in Figure 3, we applied k-means clustering on the normalized single cell data for each gene cluster and the control. Data normalization was carried out on a plate-wise basis by z-scoring each feature using the control samples as reference. In order to avoid curse of dimensionality, we restricted the dataset to the features obtained from the feature selection step mentioned earlier. We set k = 20 to be the number of subpopulations. The algorithm was run for at most 5000 iterations. Each cell was assigned to the subpopulation for which it has the shortest Euclidean distance to its center. Then, the number of cells belonging to each cell subpopulation was counted and the proportion in each subpopulation for genes in the cluster was compared against that of the control. If the change in proportion of a cell category was consistent across the genes in the cluster, the cell category is shown in the Supplementary file 2 (type B PDFs). To quantify this consistency, we used the inverse coefficient of variation of the change in a category proportion. If this quantity exceeded one, we called the change consistent and included the corresponding cell category in the PDFs. Images of cells which have highest similarity to the category center in the feature space are then used to interpret and give name to each cell category (Figure 3—figure supplement 3)

64For this purpose, we used a replicate of the original experiment but with L1000 gene-expression readouts, which is provided in the supplemental data; i.e. cell line, time point, and ORF constructs are the same. This data is different from the data used in creating GSEA plots, which entails multiple cell lines and time points. The mRNA levels are all normalized with respect to the negative control. For each replicate of the overexpression construct, we sort the expression levels of landmark genes and take the list of top and bottom 50 landmark genes. Then, to find targets of the gene related to the construct, we find the landmark genes among this list which has shown up at least in p% of replicates/clones of the gene. In particular, we set p to 33% for YAP1, 50% for WWTR1, TRAF2, and REL. Then, we simply take the intersection of predicted targets of YAP1 and WWTR1 (and similarly TRAF2 and REL, separately) to get their common targets. These targets are then used to produce Figure 6—figure supplements 1–2.

65In order to produce Figure 6C, we specified the three known targets of YAP/WWTR1 (CYR61, CTGF, and BIRC5) and queried for ORFs resulting in down-regulation of these genes. This scores each ORF (out of the 430 in the dataset) based on the observed change in mRNA level of the specified YAP/WWTR1 targets, across between four to nine different cell lines and between one to four time points. For each ORF, we then sought the summarized score which takes the mean of 4 largest scores across time point/cell line combinations. Finally, the ORFs were sorted based on the summarized score, and top 30 ORFs were tested for enrichment in different pathways (Supplementary file 1I). We used the ‘clusterProfiler’ package in R and the KEGG pathway enrichment analysis implemented in it for creating the GSEA plot.

66Wild-type and mutant sequences of WWTR1 (TAZ) (4SA: S66A, S89A, S117A, and S311A) and YAP1 (5SA: S61A, S109A, S127A, S164A, and S397A) were previously generated and cloned into the pCMV5 backbone; these constructs are distinct from those used in the original Cell Painting data set. TRAF2 and REL were cloned from the original constructs (using Broad ID# ccsbBroadEn_01710 and ID# ccsbBroadEn_11094, respectively), into pCMV5 expression vectors. These were sequenced and confirmed to BLAST against the appropriate Broad clone ID. The empty pCMV5 backbone was used as the control condition. The Tead luciferase reporter construct, 8xGTIIC-luciferase was a gift from Stefano Piccolo (Addgene plasmid # 34615).

67HEK293T cells, RRID:CVCL_0063, were transfected using Turbofect (ThermoFisher Scientific) according to manufacturer’s protocol. All cells were co-transfected with a β-galactosidase reporter plasmid (pCMV-LacZ from Clontech) as a transfection control. Cells were lysed 48 hr following transfection. Lysates were mixed with firefly luciferase (Promega) according to the manufacturer’s protocol and luminescence was measured using a luminometer (BioTek). Lysates were mixed with o-nitrophenyl-β-D-galactoside (ONPG) and β-galactosidase expression was determined spectrophotometrically by measurement of absorbance at 405 nm following ONPG cleavage. All luciferase readings were normalized to β-galactosidase expression for the sample. Statistical analysis was conducted using a two tailed unpaired Student’s t test. The data shown in Figure 6D are from triplicate samples within a single experiment and is representative of replicate experiments.