Systematic morphological profiling of human gene and allele function reveals Hippo-NF-κB pathway connectivity

We hypothesized that human genes and disease-associated alleles might be systematically functionally annotated using morphological profiling of cDNA constructs, via a microscopy-based Cell Painting assay. Indeed, 50% of the 220 tested genes yielded detectable morphological profiles, which grouped into biologically meaningful gene clusters consistent with known functional annotation (e.g., the RAS-RAF-MEK-ERK cascade). We used novel subpopulation-based visualization methods to interpret the morphological changes for specific clusters. This unbiased morphologic map of gene function revealed TRAF2/C-REL negative regulation of YAP 1/WWTR1-responsive pathways. We confirmed this discovery of functional connectivity between the NF-κB pathway and Hippo pathway effectors at the transcriptional level, thereby expanding knowledge of these two signaling pathways that critically regulate tumor initiation and progression. We make the images and raw data publicly available, providing an initial morphological map of major biological pathways for future study.


Introduction
The 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 (Amberger et al. 2015) . Unfortunately, the next step understanding the function of each gene and the mechanism of each allele in the disease typically remains nonsystematic and laborintensive. Most commonly, researchers painstakingly design, develop, and apply a diseasespecific or biological processspecific assay.
Over 30% of genes in the human genome are of unknown function (Leonetti et al. 2016) 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.
Therefore, 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 guiltbyassociation, linking unannotated genes to annotated ones based on properties such as proteinprotein interaction data, sequence similarity, or, most convincingly, functional similarity (Shehu, Barbará, and Molloy 2016) . 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 highthroughput and it has been proven useful in predicting gene function (Lamb et al. 2006) . In fact, highthroughput 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 (Berger et al. 2016) .
Images are a less mature data source for profiling but show tremendous promise. Morphological profiling data is complementary to transcriptional profiling data (Wawer et al. 2014) and is less expensive. Morphological profiling has succeeded across several applications, including grouping smallmolecule perturbations based on their mechanism of action (Caicedo, Carpenter 2016; BougenZhukov et al. 2016) , and grouping genes based on morphological profiles derived from cells perturbed by RNA interference (RNAi) (Mukherji et al. 2006; Boutros and Ahringer 2008; Fuchs et al. 2010; Pau et al. 2013 . One limitation of RNAi for morphological profiling is that the number of measurements must be limited or else the resulting profiles are dominated by offtarget effects, especially seed effects (Singh et al. 2015) . Some computational solutions have shown some promise in overcoming this problem for gene expression profiling (Schmich et al. 2015) , but their utility is unproven for imagebased 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.
In the proofofconcept 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 wellknown 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.

Morphological profiles from Cell Painting of expression constructs are sensitive and reproducible
To profile each exogenously expressed gene (or allele therein), we used our previously developed imagebased profiling assay, called Cell Painting (Gustafsdottir et al. 2013; Bray et al. 2016) . This microscopybased 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 (Fig. 1A). In five replicates in 384well plate format, we infected U2 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 (Supp . Table 1), a subset of which have been previously described (E. Kim et al. 2016) . Of these, we prioritized analysis of the 220 constructs that were most closely representative of the annotated full length transcripts (see Methods). Morphological profiles were extracted using CellProfiler for image processing, yielding 1,384 morphological features per cell, and Python/R scripts for data processing, including feature selection and dimensionality reduction (Fig. 1B, and see Methods). This computational pipeline yielded a 158dimensional profile for each of 5 replicates for each gene or allele tested.  Fig. 3). Replicates are merged to produce a profile for each gene which is then compared against others in the experiment to look for similarities and differences.
Not all genes are likely to impact cellular morphology given the limitation of our experimentusing 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 ( Fig. 2A, and see 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 (Supp. Fig. 1). Therefore, we conclude that a single "generic" morphological profiling assay can detect signal from a substantial proportion of human genes. We next turned to testing whether those signals are biologically meaningful and can lead to novel, unbiased discoveries about gene function. Constructs yielding a reproducible phenotype ought to have a median correlation among replicates that is higher than the 95th percentile of correlations seen for pairs of different constructs; this is true for 51% (112 out of 220) of the constructs (as shown). Additionally, we removed two constructs that passed that filter but whose profiles were highly similar to negative control profiles (not shown), leaving 110 constructs (50%) for further analysis. (B) Of wildtype ORF pairs that both yielded a distinguishable phenotype, 96% showed significant correlation to each other. Correlations between the 23 pairs of constructs that are clones of the same gene (although with potential sequence variation or possibly different isoforms) were almost always much higher than correlations between pairs of constructs related to different genes. (C) Genes in pathways thought to regulate morphology were more likely to yield detectable phenotypes vs. the remainder of genes in the experiment (87% vs. 48%).

Morphological profiling is robust, showing expected relationships
Given 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 "wildtype" sequence (see 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 fulllength sequence (as defined in Methods; category a: 9 pairs), or a substantive difference in their nucleotide sequence, for example, isoforms (category b: 14 pairs). We found that, as expected, the phenotypes of overexpressed wildtype 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 wildtype ORFs yielded a phenotype distinguishable from negative controls, 22 (~96%) of the pairs' profiles were correlated more than expected by chance (Fig. 2B, the one pair not meeting that threshold was in category b), confirming that different constructs with biological similarity indeed produce similar morphological profiles.
This result also confirms that the sequence differences seen in separately cloned wildtype 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. 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.
We suspected that the small number of engineered constitutively activating alleles for certain genes would, on average, yield a stronger phenotype than their wildtype counterparts. We indeed found that correlations between replicates of the constitutively activating allele were typically higher than correlations between replicates of the wildtype version of a gene (Supp. We hypothesized that genes in pathways known to affect cellular morphology (RAC1, KRAS, CDC42, RHOA, PAK1, and Hippo) 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 pvalue = 3.7 × 10 3 ) (Fig. 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 welldefined entities. We found genes manually annotated as being in the Hippo, Hedgehog, cytoskeletal reorganization, and Mitogenactivated 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 (Supp. Fig. 2). Nevertheless, the majority of pathways could be interrogated by morphological profiling.

Morphological signature similarity captures known genegene relationships
Given the caveats and limitations of overexpressing genes (see Discussion), we next tested whether imagebased profiling of expression constructs could capture relationships among genes known to be functionally related. Because a reliable and complete map of all genegene connections is not available, we evaluated the accuracy of our results via two approaches.
First, we compared our data to proteinprotein interaction data from BioGRID (Stark et al. 2006) . 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 wildtype 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 2,485 possible; Fisher's test pvalue = 0.04; Supp. Table 3).
Second, 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 samepathway annotation in the top 5% mostcorrelated gene pairs was 20% (29 pairs out of 143), significantly higher than the ratio for the remaining pairs (6%, 139 pairs out of 2,485; Fisher's test pvalue = ; Supp. Table 4). .53 0 7 1 9

An initial morphological map of gene function
Having quantitatively established that morphological profiling is sensitive, robust, and captures known genegene relationships, we explored these relationships in a correlation matrix ( Fig. 3 bottom left and Supp. Fig. 3). 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. We next created a dendrogram (Fig. 3) and defined 25 clusters (see Methods and Supp. Fig. 4) to explore the similarities among genes. Pairs of wildtype ORFs almost always clustered adjacently, consistent with our quantitative analysis described above (Fig. 2B). 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 (Supp. Table 5), indicating shared biological functions within each cluster.
Using 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 Cluster 20A and 20B PDFs, and in a later section of the text). Both are known to encode core transcriptional effectors of the Hippo pathway (Johnson and Halder 2014) , and a negative regulator of these proteins, STK3 (also known as MST2), is the strongest anticorrelating gene for the cluster (Cluster 20A PDF, panel c1).
Second, we noted Cluster 21 is comprised of the two phosphatidylinositol 3kinase signaling/Akt (PI3K) regulating genes, PIK3R1 and PTEN, both frequently mutated across 12 cancer types in The Cancer Genome Atlas (TCGA) (Kandoth et al. 2013) . These results are consistent with previous observations that certain isoforms of PIK3R1 reduce levels of activated Akt, a dominant negative effect (Abell et al. 2005) AKT3 itself is in the top five strongest anticorrelating genes for Cluster 21 (Cluster 21A PDF, panel c1). (bottom left inset, full resolution is available as Supp. Fig. 3). Only the 110 genes/alleles with a detectable morphological phenotype were included. The rows and columns are ordered based on a hierarchical clustering algorithm such that each blue submatrix on the diagonal shows a cluster of genes resulting in similar phenotypes. The correlations were then used to create a dendrogram (main panel) where the radius of the subtree containing a cluster shows the strength of correlation. The 25 clusters containing at least two constructs are printed on the dendrogram in arbitrary colored fonts, while gene names colored gray and marked by asterisks are those that do not correlate as strongly with their nearest neighbors (i.e., they are singletons or fall below the threshold used to cut the dendrogram for clustering). Each colored arc corresponds to a cell subpopulation as noted in the legend. Line thickness indicates the strength of enrichment of the subpopulation in the cluster samples compared to the negative control. Solid vs. dashed lines indicate the over vs. underrepresentation of the corresponding subpopulation in a cluster, respectively. Note that the number next to each cluster in the dendrogram is referenced in the main text and corresponds to the numbered supplemental data file for each cluster.
Third, we examined three Clusters (19, 6 and 3) that included many MAPKrelated genes. Cluster 19 is the largest example of a tight cluster of genes already known to be associated; it includes four activators in the RASRAFMEKERK cascade: KRAS, RAF1 (CRAF), BRAF, and MOS. Notably, two constitutively active alleles of these genes, BRAF V600E (H. Davies et al. 2002) and RAF1 L613V (Wu et al. 2011) , form a separate cluster (Cluster 6) adjacent to their wildtype counterparts. Furthermore, the constitutively active RAS alleles HRAS G12V and KRAS G12V (McCOY, Bargmann, and Weinberg 1984) are in the nextclosest cluster (Cluster 3), which also contains MAP2K4 and MAP2K3 (known to be activated by Ras (Shin et al. 2005) ), as well as CDKN1A (Jalili et al. 2012) . By contrast, MAPKs that are known to be unrelated to the RASRAFMEKERK cascade, such as MAPK14 in Cluster 5, are far away in the dendrogram.
Overall, these results support the notion that connections between genes can be efficiently discovered using our approach.

Visualization approaches to assist interpretation of morphological signatures
We hypothesized that the specific morphologic features that segregated each of the clusters would provide insight into gene function. Examining images (Cluster 19A PDF, panel 3) or rankordered lists of features that distinguish individual profiles or clusters (Supp. Table 6) 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.
We therefore devised several strategies to enhance biological interpretability from these experiments and applied these in combination. First, we grouped features into metafeatures 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 (Fig. 4A). Second, we performed unsupervised grouping of features by mapping the top 20 mostdistinguishing features for each cluster onto a plane, creating a Feature Map (Fig. 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 distinguishing Cluster 7A PDF), but for others these approaches did not suffice to yield an obvious phenotypic conclusion (e.g., for Cluster 19, Fig. 4A and 4B).
Third, we hypothesized that leveraging the singlecell resolution of imagebased 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 kmeans clustering on single cell data. We calculated the abundance of cells in each of the 20 subpopulations to determine which are over/underrepresented relative to controls for the given cluster (corresponding images are shown; Cluster 1B, 2B, … 25B PDFs. 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 (Fig. 4C, and Cluster 19B PDF, Categories 1 and 2), for which there is evidence in the literature (Šamaj, Baluška, and Hirt 2004; Elsum, Martin, and Humbert 2013; Godde et al. 2014) . This phenotype was not captured by manual inspection nor the first two approaches (e.g., Cluster 19A PDF, panels a2 and b2).
Encouraged by this, we supplemented the morphological map by compiling these and other visualizations into PDF files for each cluster, summarized in Fig. 5 and provided in full as Supplemental Data files for each cluster. We also noticed that certain subpopulations were similar across several clusters (Supp. Fig. 5 shows sample cell images of each such subpopulation); we annotated their enrichment/deenrichment on the dendrogram (Fig. 3). showing the greatest difference between the cluster and negative controls are shown, based on largest absolute value of zscores (full resolution version is available in Cluster 19A PDF). They are mapped in 2D space such that features that are highly correlated with each other across all genes' profiles are placed close together and thus can be interpreted together. Blue/red colored names indicate positive/negative sign of the zscore (i.e., blue indicates that the cluster shows higher values than controls). According to this map, the average intensity of AGP, RNA and Mito shows high variation for cells within samples in Cluster 19 (e.g., large mad_Cytoplasm_Intensity_MeanIntensity_AGP, where the prefix "mad" refers to median absolute deviation, a robust form of standard deviation). (C) Sample images of a subpopulation of cells enriched and deenriched for all genes in Cluster 19. Cells with asymmetric organelle distribution are highly overrepresented for genes in the cluster, and cells with more even distribution of organelles are less abundant. Scale bars are 39.36 long. Pixel intensities are multiplied by 5 for display. m μ clusters, there are two corresponding Supplemental PDF files. Left: Cluster A PDF files (e.g., "Cluster_1A.pdf") provide an overview of data about the cluster. Panel a1 lists the genes/alleles in the cluster as well as expert annotations regarding related pathways and the cell count (as a zscore) for each gene/allele. Panel b1 contains the average correlation of the cluster to other clusters, indicating uniqueness of the cluster's morphological phenotype. Panel c1 lists the top five negatively correlated gene/alleles to the cluster. Panel a2 shows the Feature Grid summarizing categories of morphological features distinguishing the cluster from the negative control. Panel b2 shows the Feature Map displaying the names of the top 20 morphological features distinguishing the cluster from the negative control, positioned based on similarity. Panel c2 shows a correlation matrix for just those genes/alleles in the cluster. Panel 3 contains sample images of fields of view of cells expressing each gene/allele in the cluster, along with images of the control for comparison. Right: Cluster B PDF files contain multiple plots aiming to illustrate the phenotype based on singlecell data, including cell subpopulation enrichment/suppression in the cluster. First, a histogram of singlecell DNA content is shown for all cells from all genes/allele treatments in the cluster, indicating the overall cell cycle distribution. Next, bar plots show (for the cluster overall and for each gene in the cluster) which of 20 subpopulations of cells are enriched and suppressed relative to negative controls. Finally, each subsequent page of the PDF is devoted to the subpopulations whose representation differs from negative controls in a statistically significant way, whether enriched or suppressed (subpopulations which are very small in both the cluster and negative control samples are omitted). For each subpopulation, a bar plot shows the top 10 mostdistinguishing feature names (versus negative control cells). Then, sample images are shown of individual representative cells from each subpopulation.
Using these visualizations, we began by interrogating three adjacent and correlating clusters (Clusters 4, 7, and 11) contain wildtype and mutant alleles of CDC42, a gene encoding a Rho family GTPase with diverse roles in cell polarity, morphology, and migration (Melendez, Grogg, andZheng 2011; Martin 2015) . Cluster 4 contains the constitutively active mutant CDC42 Q61L (Nobes and Hall 1999) as well as MAP3K2 and MAP3K9. The highly similar Cluster 7 contains the dominant negative alleles CDC42 T17N (Nobes and Hall 1999) and RAC1 T17N (S. Zhang et al. 1995) , a related RAS superfamily member. That activating and inhibiting alleles would yield similar phenotypes when overexpressed is not surprising for CDC42 (Melendez, Grogg, and Zheng 2011) . Cluster 7 also contains isoforms and alleles of AKT: specifically, AKT3 and the constitutively active E17K alleles of both AKT1 and AKT3 (M. S. Kim et al. 2008; M. A. Davies et al. 2008 . Akt is known to be essential for certain Cdc42regulated functions (Higuchi et al. 2001) and vice versa (Stengel and Zheng 2012) . Finally, the nearby Cluster 11 (which is discussed in more detail later) contains the wildtype form of CDC42 as well as TRAF2, a canonical NF B activator; these two are known to interact and share functions in actin remodeling (Marivin et al. 2014) . We also note that anticorrelating 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 (Kovacina et al. 2003; Wiza et al. 2014 , Cluster 7A PDF panels b1 and c1) (b) CDK2 (a known target of Akt (Maddika et al. 2008) ), (c) PIK3R1 and PTEN in Cluster 21, described previously, which have known interactions with AKT (Cheung and Mills 2016; Hemmings andRestuccia 2015) . Thus, all of these connections have previously been identified.
Subpopulation visualization revealed that Clusters 4, 7, and 11 are enriched in cells that are huge and binucleate (Fig. 3, example images shown in Cluster 4B PDF). Genes in all three clusters also show irregularities in DNA content, namely, an enrichment in cells with sub2N 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 Cluster 4B, 7B, and 11B PDFs). These phenotypes are consistent with these genes' known role in the cell cycle and cell polarity (Chircop 2014) .
As 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) (Asaoka 2012) , and RHOA Q63L (a Ras homolog gene family member) (Lang et al. 1996; RolliDerkinderen et al. 2005 . The highly similar Cluster 10 contains the wildtype RHOA, as well as ELK1 which is also linked to the Rho GTPase family and PKA (Bachmann et al. 2013; Murai andTreisman 2002) .
We 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 (Oishi et al. 2012) . We found that indeed all members of Clusters 8 and 10 significantly induce cell rounding (Supp. Table  7). Although cell count is lower for genes Clusters 8 and 10, the degree varies greatly (from zscore 0.67 to 3.02, Cluster 8A and 10A PDFs, 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 (Cluster 8B and 10B PDFs). 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 ( Fig. 3 and Cluster 8B and 10B PDFs).
We 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.

An unexpected relationship between the Hippo pathway and regulators of NF B signaling (Clusters 11, 20, and 22)
We 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) (Grech et al. 2004; Tada et al. 2001 , yields a signature strongly anticorrelated to YAP1/WWTR1 (Cluster 20), which encode the transcriptional effectors of the Hippo pathway, YAP (Yesassociated protein) and TAZ (Transcriptional coactivator with a PDZdomain). 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 (Varelas 2014; Hoesel and Schmid 2013; Tornatore et al. 2012) , 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 kidneyspecific conditional knockouts in mice (Reginensi et al. 2013) . 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).
We 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 withincluster connections. CDC42 and TRAF2 (Cluster 11) physically interact and share functions in actin remodeling (Marivin et al. 2014). 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 core effector of the Hippo pathway which functions to restrict YAP/TAZ nuclear activity, STK3 (which encodes the Mst2 kinase) (Meng, Moroishi, and Guan 2016) , has a morphological signature strongly anticorrelated to YAP1/WWTR1 (Cluster 20A PDF, panel c1). We also note that two clones of another protein known to influence 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 (Mohseni et al. 2014) ). Further, YAP1 is among the highest anticorrelating genes to REL (data not shown; REL is a singleton in the dendrogram and thus not in a cluster), whose protein product, cRel, has a known connection to TRAF2 (Jin et al. 2015) . These results reaffirm that the Cell Paintingbased morphological signatures are a useful reporter of biologically meaningful connections among genes in these pathways.
Given 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.
First, 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 withingroup 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 deenriched for those subpopulations and instead enriched for binucleate cells, very large cells, and small cells with asymmetric organelles (Fig. 3 and 6A and 6B). Second, 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 anticorrelated relationship at the mRNA level, consistent with the anticorrelation we had seen in morphological space. To do this, we used Gene Set Enrichment Analysis (Subramanian et al. 2005) and publicly available data, which includes data from four to nine different cell lines at one to four time points (http://lincscloud.org). This analysis revealed that the NF B pathway is the pathway most enriched among genes whose overexpression results in downregulation of known YAP1 targets, CTGF, CYR61, and BIRC5 (Zhao et al. 2008) (adjusted pvalue = , Fig. 6C), with TRAF2 being 0 2 1 8 among the genes contributing to this enrichment (Supp. Table 8). We also saw enrichment of NF B pathway members when testing a datadriven set of targets of YAP1/TAZ (Supp. Fig. 7, see Methods). In the inverse analysis, genes that alter the levels of TRAF2/REL common targets are weakly enriched in Hippo pathway members (Supp. Fig. 8, see Methods). This is consistent with the hypothesis that NF B members can downregulate YAP/TAZ targets but not strongly vice versa.
Finally, we more directly confirmed negative crosstalk between NF B effectors and YAP/TAZ using a synthetic TEAD luciferase reporter that is YAP/TAZ responsive (Dupont et al. 2011) . Importantly, these confirmatory experiments used different cellular contexts and perturbation constructs versus the original Cell Painting data. Coexpression of the NF B pathway effectors TRAF2 or CREL with YAP or TAZ led to significantly lower reporter activity than expression of YAP or TAZ alone ( Fig. 6D1 and 6D2). Intriguingly, mutants of YAP or TAZ that are insensitive to negative regulation by the Hippo pathway (YAP5SA and TAZ4SA; (Zhao et al. 2008) ) remained sensitive to suppression of transcriptional activity by TRAF2 and CREL, indicating that the negative relationship we identified may be independent of canonical upstream Hippo pathway signals ( Fig. 6D3 and 6D4).

Discussion
We 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.
By adopting a twopronged 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/TAZdirected transcription is negatively regulated by NF B pathway effectors and our data suggests a novel regulatory mechanism that is independent of upstream Hippo kinases.
To date, there has been little evidence of the intersection between these important signaling pathways. Recent work examining osteoclastosteoblast 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 (Lee et al. 2015) . TAZ was found to be a direct target of NF B transcription factors and its expression is regulated via NF B signaling (Cho et al. 2010) . 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 cofactors. 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 (B. Liu et al. 2016) . 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 pathwaymediated phosphorylation events on YAP/TAZ, suggesting a more direct relationship between NF B and YAP/TAZ signaling.
In 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 transdominant negative effect that results in precisely the opposite phenotypic effect (Veitia 2007) . 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 nonphysiological 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 offtarget effects (often due to seed effects), which were far more prevalent (impacting 90% of constructs in our recent study (Singh et al. 2015) ).
In 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 (Berger et al. 2016) . 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.
We 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 Supp. Data files). Expanding this map to full genome scale could prove an enormously fruitful resource.

cDNA constructs used for expression
The Reference Set of human cDNA clones utilized here has been previously described (E. Kim et al. 2016) ;~90% of these constructs induce expression of the intended gene greater than 2 standard deviations above the control mean. Briefly, wildtype 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 For simplicity, throughout this paper "wildtype" 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 nonidentical 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 genomescale library (Yang et al. 2011) .

Lentiviral transduction for morphological profiling
We followed our previously described protocol (E. Kim et al. 2016; Berger et al. 2016) except for durations of some steps. Briefly, cells were plated in 384well plates and transduced with lentiviral particles carrying ORF constructs the next day. Viral particles were removed 18-24 hr postinfection and cells cultured for 48 hr until staining and imaging (72 hr total posttransduction).

Cell staining and imaging
The Cell Painting assay followed our previously published protocol (Bray et al. 2016

Image processing and feature extraction
The workflow for image processing and cellular feature extraction has been described elsewhere (Bray et al. 2016

Profiling and data preprocessing
The code repository for the profiling and all the subsequent analysis will be made publicly availabe at https://github.com/carpenterlab/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 (Ljosa et al. 2013) . We also include the cell count in a sample as an additional feature. This results in a vector of 2,769 elements describing the summarized morphology of cells in a sample. We then use the median polishing algorithm immediately after obtaining the summarized profiles, to remove and correct for any plate position artifacts. For each feature, the algorithm detrends the rows, i.e. by subtracting the row median from the corresponding feature of each profile in that particular row. Next, it detrends the columns in a similar way using column medians. The row and column detrending 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.
Not 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 2,200 features. Subsequently, a feature selection method is used (Fischer et al. 2015) . 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 alreadyselected 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 (Fischer et al. 2015) 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.
The feature selection method greatly removes redundancy, but because of the nonoptimal "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.

Feature interpretation
The 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: • Prefix 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. • Cellular 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. • Measurement 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. • Name(s) of channels in which the measurement is made, if appropriate (omitted for AreaShape and Neighbors). • Feature name. The precise measurement name appears at the end. A description of each metric can be found in the CellProfiler manual (http://d1zymp9ayga15t.cloudfront.net/CPmanual/index.html)

Identifying ORF constructs that are distinguishable from negative controls
Our 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. The cutoff for "discernible" is set based on the top 5th percentile of a null distribution. The null distribution is defined based on the correlations between nonreplicates (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 ( Fig. 2A).
At 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.

Comparison of morphological connections between genes to proteinprotein interaction data and pathway annotations
In this analysis, mutant alleles were removed and we considered only one wildtype 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 wildtype clone pairs versus differentgene pairs. To do so, we found the value at which the probability density functions of the two groups intersect (Supp. Fig. 6); this value (here, 0.43) can be proved to have the desired property (Duda, Hart, and Stork 2012) . This approach results in about 5% of the gene pairs being categorized as highly correlated. We next formed a 2 by 2 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; Supp. Table 3 and 4). This table was then used to perform a onetailed Fisher's exact test.

Creation of a dendrogram relating genes to each other, and agglomerative clustering by cutting the dendrogram
A dendrogram was created based on the Pearson correlation distance and average linkage, using the hclust function in R (Fig. 3).
Gene clusters were formed by cutting the dendrogram at a fixed correlation level, 0.522, which was chosen using a stabilitybased 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 (Supp. Fig. 4), 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, . 002 ε = .

Subpopulation extraction
In order to extract cell categories (subpopulations) and subpopulation enrichment laid over the dendrogram in Fig. 3, we applied kmeans clustering on the normalized single cell data for each gene cluster and the control. Data normalization was carried out on a platewise basis by zscoring 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 supplemental PDFs. To quantify this consistency, we used the inverse coefficient of variation of the change in a category proportion. If this quantity exceeded 1, 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 (Supp. Fig. 5)

Identifying targets of a gene using a datadriven approach
For this purpose, we used a replicate of the original experiment but with L1000 geneexpression 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 Supp. Fig. 7 and 8.

Gene Set Enrichment Analysis
In order to produce Fig. 6C, we specified the three known targets of YAP/WWTR1 (CYR61, CTGF, and BIRC5) and queried for ORFs resulting in downregulation 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 (Supp .  Table 8). We used the "clusterProfiler" package in R and the KEGG pathway enrichment analysis implemented in it for creating the GSEA plot (Yu et al. 2012) .

Luciferase Reporter assay
Wildtype 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, 8xGTIICluciferase was a gift from Stefano Piccolo (Addgene plasmid # 34615).
HEK293 cells were transfected using Turbofect (ThermoFisher Scientific) according to manufacturer's protocol. All cells were cotransfected with a βgalactosidase reporter plasmid (pCMVLacZ from Clontech) as a transfection control. Cells were lysed 48 hours 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 onitrophenylβDgalactoside (ONPG) and βgalactosidase expression was determined spectrophotometrically by measurement of absorbance at 405nm 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 Fig. 6D are from triplicate samples within a single experiment and is representative of replicate experiments.

Supplementary Data/Figures
Supplemental data: We will upload the L1000 gene expression data to GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70138) and the image data along with the extracted morphological features at the percell level to the Image Data Repository ( http://idrdemo.openmicroscopy.org/about ).
As supplemental data files, we provide: • Two PDFformatted files (A and B) containing details about each of the 25 clusters (50 total). • Files containing the percell and perwell features extracted by CellProfiler.
• The CellProfiler pipeline used to quantify the images (.cppipe format).
Supp. Figure 1: Position artifacts do not contribute to the hit rate seen in the experiment. We were concerned that position artifacts may result in overestimating the replicate correlations because replicates of the same treatment are assigned to the same well location across different plates in the experiment (that is, it was infeasible to scramble well locations). We ruled out this possibility by taking an alternative pessimistic null distribution which takes well position into account. In contrast to Figure 2A, which shows a 51% hit rate, a more pessimistic alternative null distribution is shown here (left), calculated based on the replicate correlation of pairs of negative controls in the same position only. We consider this less reliable because the number of such pairs is small (26) and we excluded edge wells; nevertheless the hit rate increases slightly, to 60%.

Supp. Figure 2: Strength of morphological phenotypes, according to annotated pathway.
Morphological phenotype strength is calculated as the average replicate correlation for genes that experts manually annotated genes as belonging to each pathway. The number of genes tested in each category is shown in parentheses after the pathway name (two wildtype clones of the same gene are only counted once, but a mutant allele is counted separately). The red line shows the threshold beyond which an individual gene's profile would be considered to yield a distinguishable phenotype. Error bars indicate the deviation in the replicate correlation among the genes associated with the pathway.

Supp. Figure 3: Correlation among the 110 genes/alleles with a detectable morphological phenotype.
The rows and columns are ordered based on a hierarchical clustering algorithm such that each blue submatrix on the diagonal shows a cluster of genes resulting in similar phenotypes. The scale bar depicts Pearson correlation.
Supp. Figure 4: Smoothed stability score across different cutoffs, in order to choose a threshold for cutting the dendrogram to form clusters. The maximum occurs at threshold = 0.522. Smoothing is done by taking the moving average of order 0.02. The stability score is defined as the proportion of treatments whose clusters are not affected if the cutoff is increased or decreased by a small amount ( ). 002 ε = .
Supp. Figure 5: Common cell subpopulations seen across more than one cluster. These names are used to annotate clusters of genes in Fig. 3. Example images shown are taken from individual clusters. Scale bar is 63 and image intensities are log normalized. References to size and shape in the subpopulation m μ legends refer to both the nucleus and cell borders, unless otherwise noted.
Supp. Figure 6: Selection of the threshold to determine highly correlated gene pairs. Correlation probability densities of pairs of different genes (labeled 'nonclone') and pairs of wildtype clones cross at 0.43, making 0.43 the optimal threshold to determine highly correlated gene pairs. Figure 7: Gene Set Enrichment Analysis (GSEA) reveals that overexpression constructs sorted based on their similarity to YAP1/WWTR1 overexpression (in terms of impact on particular mRNA targets), are enriched for regulators of the NF B pathway (Enrichment Score pvalue = 0.0019). mRNA targets common to both YAP1 and WWTR overexpression include INPP4B, MAP7, LAMA3, STMN1, and TRAM2, which are positively regulated, and SPP1, IER3, RAB31, and GPR56, which are negatively regulated.

Supplementary Tables
Supp.