A Pooled Cell Painting CRISPR Screening Platform Enables de novo Inference of Gene Function by Self-supervised Deep Learning

Development and optimization of CellPaint-POSH

We made several modifications to optimize the Cell Painting assay ¹⁶ for pooled optical screening (workflow described in Fig. 1a, detailed plasmid maps for lentiviral delivered Cas9 and single-guide RNA (sgRNA) shown in Fig. S1a-b). First, MitoTracker used in Cell Painting is a stress-inducing live cell stain, and the sustained fluorescence of MitoTracker into the sequencing stages prohibited reliable in situ base calling (data not shown). We thus developed an RNA-based label for mitochondria that we term “Mitoprobe.” Briefly, using the structure of the human mitochondrial ribosome from the Protein Data Bank (3J9M) ³⁰, we identified the most likely RNA sequences to bind to the human mitochondrial ribosome’s 12s and 16s ribosomal RNA (rRNA) (Fig. 1b). We then optimized on solvent-accessible surface area (SASA) ³¹, probe length, GC content, absence of repetitive bases, and intra-sequence distance (Fig. S1c), and conjugated the 5’ end with Cy5. Simultaneous co-staining of fixed adenocarcinomic human alveolar basal epithelial (A549) cells with Mitoprobe and Tom20 antibody gives concordant staining patterns (Fig. 1c, Fig. S1d), demonstrating in a POSH-compatible mitochondrial RNA probe that can be applied to fixed cells (Table 1), streamlining all cellular stains into a single step in the workflow. The RNA-FISH based mitochondrial label is washed out during ISS chemistry and heating steps, avoiding optical interference with 4-color sequencing. Second, RNAse inhibitor Ribolock was implemented during all staining phases to prevent degradation of mitochondrial rRNA, sgRNA spacer containing transcripts and other cellular RNAs. Third, we observed that subjecting cells to Cell Painting prior to reverse transcription (RT) led to degradation of RNAs and ISS signals (data not shown). Moving RT prior to Cell Painting preserved RNAs and led to successful ISS (Fig. 1d).

• Download figure
• Open in new tab

Figure 1. Overview and development of POSH platform.

a. a general overview of the Pooled Optical Screening in Human cells platform (POSH). A pool of cells is generated that contain both constitutively-active Cas9, as well as lentivirally-delivered sgRNAs targeting a select number of genes. Cells are fixed and the sgRNA sequences are amplified via rolling circle amplification (RCA). The CellPaint imaging assay is then conducted on the cells followed by in situ sequencing by synthesis to match image data with transfected sgRNA and gene knockout. b. the protein structure of the human mitochondrial ribosome (transparent), 16s, 12s, and Val highlighted in teal, green, and orange, respectively. The exposed regions of the 12s and 16s rRNAs are highlighted in red. c. 60X water-immersion confocal fluorescence images of A549s treated with Hoechst, mitoprobe, and anti-Tom20 antibody, reveals close overlap between mitoprobe and Tom20 without nuclear staining. d. example raw images showing the multiple CellPaint tiles, as well as the several sequencing-by-synthesis (SBS) tiles paired to the phenotypic data.

The finalized 5-stain morphological assay panel includes Hoechst, concanavalin A (ConA), wheat germ agglutinin (WGA), phalloidin, and mitoprobe, which stain the nucleus, endoplasmic reticulum (ER), membranes, actin, and mitochondria respectively. This leaves one vacant channel for flexible use, such as a hypothesis-specific biomarker. After imaging cells, ISS was conducted to identify the sgRNA within each cell (Fig. 1d, Methods).

We created a ML enabled data pipeline for high-throughput processing of phenotypic imaging and sequencing by synthesis (SBS) acquisitions in order to create a [morphological feature] * [genetic perturbation] matrix for downstream analysis (Fig. 2a). The process utilized Hoechst staining from all acquisitions to create registration tables, ensuring all morphological phenotyping and SBS acquisitions were spatially aligned and removing the need for manual coarse alignment of field of views during image acquisition, simplifying the lab process. Our approach also includes an improved methodology for base-calling by training a 3-layer fully convolutional neural network (FCN) (Fig. 2a), which increases the deconvoluted cell recovery rate from 66.6% to 78.8% (Fig. 2b), and the determined sgRNA sequences were then matched to the original library design, with ∼2.0, ∼4.0, and ∼4.0 median amplicons detected per amplicon-presenting cell across three separate screens, respectively (Fig. 2c). Samples were sequenced beyond the minimum number of nucleotides needed to call a specific sgRNA within the library. We targeted a minimum Hamming distance of 3 between any two sgRNAs, enabling sgRNA correction in the event of single-nucleotide mismatch. Across the 3 screens, this recovered an additional 5.9%, 5.3%, and 34.8% of cells respectively (Fig. 2d), making the ISS workflow more robust against occasional low-quality cycles or mistakes in base calling. We developed a simple metric for measuring the fidelity of POSH dots using the signal-to-total-ratio (STR, Methods), and found the quality of ISS to be consistently high across cycles of sequencing and across experiments in general (Fig. S2a-d). A single-cell dataset is generated by cropping cell image tiles centered on each nucleus, masked by its corresponding cell mask and associated with a sgRNA identity based on the mapped barcode locations (Fig. 2a). Morphological features are determined using both classical feature extraction and self-supervised learning described below.

• Download figure
• Open in new tab

Figure 2. Computational pipeline and QC of POSH platform.

a. a general overview of the computational pipeline for converting raw phenotyping and SBS images into usable tiles and feature matrices. Sequencing by synthesis (SBS) and readout images are registered using Hoechst staining, followed by amplicon base calling and alignment, in parallel with illumination processing, cell segmentation, and tiling. Multiple analysis methods can then be implemented such as Deep Learning-based methods or direct featurization, followed by embeddings and gene correlation calling. b. % cells recovered with a valid sgRNA barcode from Feldman et al¹² (grey), our study using classical blob detector (orange), and our data using ML (green). c. number of valid amplicons per cell across the three screens. Most cells contain at least one valid amplicon. d. improvement of cell count based on Hamming correction of miscalled sgRNAs.

Proof-of-Concept (PoC) with Morphological Regulators

To demonstrate that high-content imaging using Cell Painting captures rich morphological information and can be used to classify gene functions, we set out to perform a small scale CellPaint-POSH experiment against 124 genes with known morphological impact ³². A sgRNA library was synthesized to target six key pathways: mitochondrial translation (25 genes), proteasome (27 genes), actin/kinesin (23 genes), unfolded protein response (10 genes), golgi-ER retrograde processes (13 genes), and microtubule/dynein (19 genes), along with 7 miscellaneous genes and 124 nontargeting/intergenic control sgRNAs (Fig. 3a, Supp Data 1), at a representation of 10 sgRNAs per gene for robust analysis. A549 cells were processed with CellPaint-POSH protocol with the aforementioned staining panel. We first assess the fitness effect to validate gene knockouts (KOs). As expected, DepMap-determined common essential genes ³³ like COPA, KIF11, and HSPA5, as well as several of the proteasome genes, had lower representation than other sgRNAs (Fig. S3a); this provides validation that the KOs are effective and that sgRNAs are correctly sequenced and decoded.

• Download figure
• Open in new tab

Figure 3. Design and classical image feature based analysis of 124-gene PoC screen.

a. summary of genes targeted in the screen in order to create distinct morphological effects. b. comparison of morphology network gene-gene edges to those from StringDB. Edges defined as strong by StringDB show similarly high edge strengths in our morphological network. c. AUC-ROC of CellPaint POSH screen using StringDB network as ground truth, showing strong overlap (AUC=0.77). d. UMAP projection of CellStats features showing genes (filtered by AUC > 0.55) clustered by their biological functions. Colors represent Leiden communities, node size represents the median similarity of aggregate sgRNA embeddings targeting the same gene. e. volcano plot to identify genes that show significant variation in nucleus size compared to intergenic control sgRNAs. Insets show example tile images of the nuc_size_hi and nuc_size_lo cells that led to the overall disruption in the distribution of nuclear sizes for their respective knockout populations. Intergenic control in orange. f. volcano plot of features comparing MRPL57 KOs to intergenic controls. g. histogram of the cytoplasmic 95th-percentile intensity feature in intergenic-sgRNA-presenting and MRPL57-sgRNA-presenting cells. h. example images of CellPaint on intergenic and MRPL57 KO cells. i. volcano plot of features comparing ACTB knockouts to intergenic controls. j. distribution of mean cytoplasmic phalloidin intensity of control and ACTB KO cells. k. example images of the cells with high cytoplasmic phalloidin staining leading to the skewed distribution shown in (j). l. volcano plot of features comparing GOSR2 knockouts to intergenic controls. m. distribution of perinuclear conA intensity of intergenic and GOSR2 KO cells. n. example images of the cells with high perinuclear conA intensity leading to the skewed distribution shown in (m).

In total, 163,090 cells were analyzed using a classical image featurization engine similar to CellProfiler ³⁴ that extracted 1,301 morphological features in four broad categories: (i) localized pixel intensity statistics, (ii) geometric features of cell segmentation masks, (iii) features characterizing textures that emerge from the different types of staining, and (iv) correlations of intensities across multiple channels (Supp Data 2, Methods). We term this classical featurization engine “CellStats.” Each unique sgRNA was assigned a ’morphology signature’ based on the mean features of all cells containing that sgRNA (Methods). Cosine similarity was determined between these averaged feature vectors to assign similarity scores among all sgRNAs. Similarities between sgRNAs targeting the same gene were specifically compared (Fig. S3b), with a high similarity indicating effective sgRNA cutting and strength of a KO’s morphological phenotype. sgRNA signatures were then further aggregated to the gene level. To systematically assess the accuracy of the biological network emerging from our morphological analysis against prior literature, we compared it to StringDB ³⁵, which constructs a gene-gene interaction network in which interactions are weighted by a combination of literature scrubbing, gene-gene interaction databases, co-expression analysis, and organism transfer, and used this as a ground truth comparator. We first compared low, medium, and strong edges in our network to those formed by StringDB and found a very significant increase in correlation between networks at higher StringDB cutoff score (Fig. 3b). We additionally determined the area under the curve of the receiver operating characteristic (AUC ROC) using StringDB edge > 0.95 to define true positives, resulting in a 0.77 AUC, indicating strong overlap with the established network for the same genes (Fig. 3c). We visualize the gene embeddings using Uniform Manifold Approximation and Projection (UMAP) ³⁶ grouped by Leiden community detection³⁷ to further illustrate the core clusters, notably recreating functional clusters corresponding to the golgi / ER, proteasome, mitochondrial translation (2 distinct clusters), and microtubule / dynein classes (Fig. 3d). Comparison with StringDB and visual inspection of gene clustering both demonstrate that CellPaint-POSH can reconstruct known gene-gene interaction networks.

We took two approaches to further investigate the morphological features driving the gene clustering. First, we identify genes affecting a given morphological feature, such as nuclei size, by computing each KO’s Z-score for this feature (Fig. 3e). Genetic drivers of enlarged nuclei size included gene KOs that yielded strikingly large nuclei (e.g. KIF11), or multiple nuclei within the same cell (e.g. KIF23), which are consistent with nuclear size, mitotic arrest and spindle polarity phenotypes reported on the same genes in Funk et al¹³. We also observed a significantly smaller nuclear size in TLN1 KO cells, likely linked to the inability of the cells to form adhesion with their substrate or neighboring cells due to lack of talin-1 protein. Other strong modifiers of nucleus size tended to fall within the Microtubule / Dynein subgroup or encode for members of the actin-related protein (Arp) family, which form key chromatin-remodeling complexes ³⁸.

In the second analysis approach, which we term “differential morphology analysis,” we examine all CellStats features modulated by a particular genetic perturbation, relative to the control cells. Differential morphological features are similar in concept to differentially expressed gene analysis commonly used in comparing two sets of transcriptomics data. As an example, we selected one gene that represented a single point in one of the tightly-connected clusters within the network (MRPL57, Fig. 3d). We conducted Kolmogorov-Smirnov (K-S) tests between MRPL57 cell features and intergenic-targeting sgRNA controls, and found that perinuclear mitochondrial intensity, variance, and texture were the most significantly different features of MRPL57 KO (Fig. 3f-g), which matches the mitochondrial dysfunction anticipated with the KO. Representative plots of MRPL57 KO cells driving these mitochondrial scores depict the low-mitochondria phenotype (Fig. 3h). We next examined cytoskeletal protein ACTB. Interestingly, phalloidin scoring indicated a significant increase in actin presence within the cells; this provides evidence for compensatory up-regulation of other actin isoforms upon loss of β-actin ³⁹ (Fig. 3i-k). Lastly we examine Golgi SNAP Receptor Complex Member 2 (GOSR2), whose loss of function likely leads to progressive myoclonus epilepsy in patients ⁴⁰. Examining the differential morphological features upon GOSR2 KO (Fig. 3l), its unique morphological phenotype was driven largely by increased intensity, variance, and texture of Golgi/ER specific ConA staining (Fig. 3m-n), matching its likely biological mechanism ⁴¹, and providing a disease phenotype that could be adapted into a disease-modifier screen.

Mapping MoA with Self-Supervised Vision Transformer Learning

Encouraged by our morphology-driven PoC, and inspired by previous efforts that use Cell Painting information to cluster compounds by their MoA ^24,27, we next attempted to assess the performance of CellPaint-POSH at phenotypic clustering of genes by their annotated mechanism of action, irrespective of previous reported morphological phenotypes. We curated a list of 300 genes whose gene products are targeted by tool compounds with well annotated MoAs ^42–44 (Supp Data 3). We built self-supervised vision transformer (DINO-ViT) models ⁴⁵ that can extract meaningful image representation with no labels, and compared between the following imaging representation techniques: 1) classical morphological featurization (“CellStats”); 2) self-supervised DINO-ViT model, trained on ImageNet ⁴⁶ containing ∼1.2 million natural image data: e.g. animals, vehicles, and tools (“ImageNet-dino”); 3) DINO-ViT embeddings trained on ∼1.5M single cell Cell Painting images from the 300-gene MoA experiment (“CP-DINO 300”, Methods, full cell tile image dataset provided in Supp Data 4). All three featurization methods are capable of morphologically classifying many genetic perturbations from non-targeting control cells (Fig. 4a-b), including ImageNet-dino which is not trained on cellular morphology data. This is consistent with previous reports of transfer learning ⁴⁷. Nevertheless, CP-DINO 300 trained on bioimaging data yielded a more informative embedding that has higher median prediction accuracy than the other two models (Fig. S4a-b), and correctly classified more perturbations with better accuracy (Fig. 4c). CP-DINO 300 also recovered more known biological relationships from StringDB as measured by cosine similarity of the aggregate gene KO embeddings (Methods) than the other two models (Fig. 4d).

• Download figure
• Open in new tab

Figure 4. DINO trained embedding outperforms Cellstats and ImageNet-trained embedding.

a-b. Comparison of feature embedding methodologies based on median AUC of binary classification of KO from WT for each genetic perturbation. a. CP-DINO 300 vs. CellStats b. CP-DINO 300 vs. ImageNet-dino. c. Number of genetic perturbations classified at AUC > 0.55 by different representation models. d. True positive StringDB edges predicted (at 5% false positive rate) performance using gene-gene similarity calculated from different representation models. e. UMAP projection of DINO-ViT features showing genes (filtered by AUC > 0.55) clustered by their cellular component localization. Colors represent Leiden communities, node size represents the median similarity of aggregate sgRNA embeddings targeting the same gene. Cell tile images and top neighbors (cosine similarity in parenthesis) of select KOs are shown.

Similar to CellStats, CP-DINO 300 representation of the Cell Painting assay allows reconstruction of gene networks, and the identification of new pathway components in a hypothesis-free way using Leiden community detection (Fig. 4e) ³⁷. Specifically, we were able to reconstruct the genetic modifiers of glycoprotein biosynthesis, mitochondrial translation, actin/cytoskeletal organization, lipid metabolism, PI3K/Akt activation, mTORC1 signaling, PRC2 complex, where morphological similarity associates genes from the same pathways / protein complexes. For example, genes with the highest similarities to SUZ12 included the other two PRC2 complex members EZH2 and EED ⁴⁸ (Fig. 4e); genes similar to mTORC1 activator RHEB included other pathways activators mTOR and PDPK1 ⁴⁹. Interestingly, the network was capable of clustering key components of the lipogenesis pathways without the use of an explicit lipid stain, including the core fatty acid synthesis enzymes (ACLY: ATP Citrate Lyase, ACACA: Acetyl-CoA Carboxylase Alpha, SCD: Stearoyl-CoA Desaturase, and FASN: Fatty Acid Synthase), and upstream AKT signaling regulators (PIK3C3: PI3K Regulatory Subunit 3, PIK3R4: PI3K Regulatory Subunit 4) that all contribute to lipogenesis ⁵⁰ (Fig. 4e). Representative images of gene KOs with the most distinguished phenotypes (highest binary predictive scores) from each cluster are shown (Fig. 4e).

Deep Learning Enables Novel Discovery via Druggable Genome Pooled Optical Screen

We next sought to scale CellPaint POSH to the druggable genome to demonstrate the throughput needed for discovery screens, and to further assess the generalizability of Cell Painting across more diverse sets of genes and pathways. To this end, a library of 1640 genes was assembled based on the Tier 1 list of the “Druggable Genome” library ¹⁷, spiked in with a subset of morphological genes and known mTORC1 pathway genes for controls (Supp Data 5). In addition to Cell Painting, anti-phospho-S6 (pS6) antibody with AlexaFluor 750-conjugated secondary antibody was used in the 6th channel as an established biomarker ⁵¹ used extensively in mTORC1 studies including genome-wide screens ^22,52,53. In order to improve throughput and reduce labor-intensive manual work, we built and deployed fully automated liquid handling and microscopy for the ISS workflow (Methods, video demonstration in Supp Data 6). We formally assessed sgRNA library abundances as determined by ISS with that measured by genomic DNA next-gen sequencing (NGS), by comparing the enrichment and depletion of sgRNAs targeting fitness genes. We again re-discovered known essential genes ³³, and the measured fitness effects were comparable between POSH and NGS (Fig. S5a-c), indicating that there is little detection bias in POSH ISS.

We explored if the added diversity in the 1640-Gene druggable genome POSH dataset would lead to better self-supervised image representation. To that end we trained CP-DINO 1640, and found that while it performs similarly to CP-DINO 300 at binary classification of genetic KOs vs. negative controls in the 1640 dataset (Fig. 5a-b), CP-DINO 1640 captures more semantically meaningful structure in the data as demonstrated by its more accurate predictions of StringDB gene-gene interactions than CP-DINO 300, ImageNet-dino, or CellStats (Fig. 5c).

• Download figure
• Open in new tab

Figure 5. Morphological analysis of druggable genome discovery screen.

a. Comparison of feature embedding methodologies (CP-DINO 1640 vs. CP-DINO 300) based on median AUC of binary classification of KO from WT for each genetic perturbation. b. Number of genetic perturbations classified at AUC > 0.55 by different representation models. c. True positive StringDB edges predicted (at 5% false positive rate) from 1640 druggable genome dataset using gene-gene similarity calculated from different representation models. d. True positive StringDB edges predicted (at 5% false positive rate) from a held-out 124 gene PoC dataset using gene-gene similarity calculated from different representation models. e. True positive StringDB edges predicted (at 5% false positive rate) from 300 gene MoA dataset using gene-gene similarity calculated from different representation models. f. UMAP projection of DINO-ViT features showing genes (filtered by AUC > 0.54) clustered by their biological functions. Colors represent Leiden communities, node size represents the median similarity of aggregate sgRNA embeddings targeting the same gene. Cell tile images and top neighbors (cosine similarity in parenthesis) of select KOs are shown.

We next evaluated all image representation models on held out data: the 124-gene morphology PoC perturbation dataset. While all 3 self-supervised DINO-ViT models generalized out of experiment and outperform CellStats at predicting StringDB networks (Fig. 5d), CP-DINO 1640 in particular surpassed prediction accuracy from CP-DINO 300 and ImageNet-dino (Fig. 5d), and faithfully reconstructs gene networks (Fig. S5) by Leiden clustering. The superior performance of CP-DINO 1640 is unlikely a result of trivial memorization, as the 1640-genes druggable genome library and 300-genes MoA library share similar numbers of overlapping genes with the 124 PoC library (30 and 26 genes respectively). Lastly, when evaluated on the 300-MoA dataset, both CP-1640 and CP-300 performed similarly in StringDB prediction (Fig. 5e), likely because the latter already learned the full variance in the corresponding dataset.

We constructed gene networks from the druggable genome dataset using CP-DINO 1640 model and Leiden community detection (Fig. 5f). Overall, many genetic pathways and known protein complexes emerged from network analysis of imaging features from the unbiased morphological profiling. As an example, a cluster of morphologically-related genes formed by MAP4K4, MAP3K2, PIK3CA, EGFR and CREBBP / EP300, all known regulators of EGF signaling signaling ^54–57. As another example, top related genes to TGFBR1 included other known TGFb genes such as TGFBR2, TGFBR3 and SMAD3 (Fig. 6d)⁵⁸. Among the top genes similar to NCSTN is PSENEN, both of which encode for the gamma-secretase complex components nicastrin and presenilin enhancer 2, respectively ⁵⁹. A cluster with strongest gene-gene similarities include ACTR10, POLD1, TUT1, CDK7, TUBB, and TOP2A. Interestingly, while their KO phenotypes are similar, these genes perform different biological functions (e.g. TUT1 encodes a nucleotidyl transferase that functions as both a terminal uridylyltransferase and a nuclear poly(A) polymerase, while PLK4 is involved in centriole duplication). We note that most of these are reported essential genes ³³ and confirmed by sgRNA abundances in our experiment (Fig. S5d). Indeed, PLK4 is a known essential regulator of cell cycle ⁶⁰, and it is possible that nucleotidyl transferase TUT1 may also play an important role in cell division. Consistent with this hypothesis, CellStats analysis suggested that TUT1 KO induces large nuclei size (Fig. 3e). In addition, Luke et al. recently reported both TUT1 and PLK4 KOs produce larger cell area, nucleus area and nuclear DNA integrated intensity, consistent with the phenotype of cell cycle regulator genes such as CDC proteins ¹³. Both Deep Learning (DL) and CellStats analysis nominate the potential role of TUT1 in cell cycle regulation.

• Download figure
• Open in new tab

Figure 6. Combining morphology and pS6 biomarker analysis for mTORC1 pathway discovery and validation.

a. Comparison of feature embedding methodologies (CP-DINO 1640 with and without pS6 information) based on median AUC of binary classification of KO from WT for each genetic perturbation. b. Top neighbor genes to mTOR and TSC2 (mTORC1 activator and inhibitor respectively) predicted with and without pS6 are shown. Numbers indicate cosine similarity. c. the pipeline for conducting single-cell based analysis via Kolmogorov-Smirnov 2-sided, 2-tailed statistical test with benjamini-hochberg correction, in contrast to methods used for flow-based screens such as MAUDE (Fig. S7a). d. results of single-cell based, one-feature analysis of cytoplasmic phospho-S6. Teal genes are well established mTOR upregulators, orange genes are well established mTOR downregulators, black genes are potentially novel pS6 modifiers. e. representative pS6 field with segmentation. f. example pS6 and mitoprobe staining images from control, known, and potentially novel regulators of mTORC1. g-h. differential feature analysis on AURKAIP1 and HSD17B10.

Another major pathway that emerged from gene network analysis is mTORC1, where 3 known mTORC1 complex proteins form a tight cluster: RPTOR, mLST8, mTOR (Fig. 5f) ⁶¹. Other members of the core mTORC1 cluster include well established activators of the PI3K/AKT/mTOR pathway such as CSNK2A1 and PDPK1 ^62,63. In contrast, mTORC2 core complex component, RICTOR, clusters away from this group, likely because of its primary role in a distinct mTOR complex with different functions and morphological consequences on the cell. In addition, mTORC1 inhibitors TSC1 and TSC2, whose loss of function (LoF) mutations lead to the eponymous genetic disorder ⁶⁴, cluster away from mTORC1 activators as expected (Fig. 5e).

Given the clustering of mTORC1 regulators, we asked if the inclusion of the pathway biomarker pS6 may improve the discovery of relevant genetic perturbations. We thus trained CP-DINO models with Cell Painting + pS6 channel information, and found that while both CP-DINO 1640 and the new CP-pS6-dino 1640 have similar sensitivity in predicting genetic perturbations, pS6 notably improves the detection of many mTORC1 inhibitors (TSC1, TSC2, DDIT4, TBC1D7) (Fig. 6a). Indeed, while mTORC1 activators (defined as top genes morphologically similar to mTOR) are similarly predicted with or without pS6 (Fig. 6b), interestingly mTORC1 downregulators (defined as genes morphologically similar to TSC2) include relevant hits such as TBC1D7 and DDIT4 only when pS6 information is taken into model building ^65,66. In order to further validate morphologically discovered mTORC1 regulators against the ground truth biomarker, we systematically analyzed genes that influence pS6 intensity with two methods: 1) a binning approach, grouping cells with similar levels of pS6 intensity changes, simulating enriched bins of cells produced by FACS-based CRISPR screens (Fig. S7a); 2) a distribution based approach, taking full advantage of the single cell resolution of our data (Fig. 6c, Methods). Both analyses resulted in comparable hits, effectively identifying many of the established regulators of the mTORC1 pathway ^22,52 (Fig. 6d, Fig. S7b). Most top mTORC1 activators as defined by morphological similarity to mTOR, including mTOR, PDPK1, RPTOR, mLST8, CSKN2A1 mentioned above (Fig. 6b), are also top regulators of the pS6 biomarker (Fig. 6d). Interestingly, mTORC1 hits discovered by the pS6 assay can be further stratified into different subclusters by generic morphological analysis (Fig. 5e). For example, certain pS6 hits GFM1, HSD17B10, MRPS5, MRPL21, and AURKAIP1 form a morphological cluster that is separate from other pS6 hits and core mTORC1 regulators (e.g. mTOR, PDPK1, and RPTOR) (Fig. 5e), suggesting that the former group of genes may regulate mTOR biology through a distinct mechanism. Based on their biological functions ^67–69, we hypothesize that mitochondrial translation and metabolism may be involved in mTORC1 regulation. While the role of mTORC1 in regulating mitochondrial function has been reported⁷⁰, the reverse causation is not fully understood. We performed differential feature analysis on these genes, and showed that genes such as AURKAIP1 and HSD17B10 not only lower pS6 levels in the cell but also have a significant impact on mitochondrial staining (Fig. 6e-f). Furthermore, CellStats based differential feature analysis shows that pS6 and mitochondria-related morphological features are among the top changes upon AURKAIP1, HSD17B10, and MRPS5 KOs (Fig. 6g-h, Fig. S7c). Consistent with our findings, HSD17B10 was recently discovered as a novel mTORC1 regulator in a genome-wide screen against pS6 ²². In addition, mitochondrial translation defects caused by the deletion of another mitochondrial elongation factor mtEF4 in C. elegans and mouse have been shown to regulate mTOR and cytoplasmic translation ^71,72 as a compensatory mechanism to cope with mitochondrial stress. We hypothesize that other mitochondrial genes like AURKAIP1 and HSD17B10 may have an effect on mTORC1 in our experiment. Interestingly, mutations in HSD17B10 have been linked to recurrent seizures and intellectual disability ^73,74, two symptoms also commonly observed in mTORopathies ⁷⁵.

Pooled morphology CRISPR screening for hypothesis-free exploration

Pooled CRISPR screens have emerged as a powerful technique to interrogate multiple genetic interventions in parallel; they are considerably more cost effective and scalable than arrayed screens, and are able to greatly ameliorate batch effects. CellPaint-POSH combines the benefits of pooled CRISPR screening with the rich cell state information provided by Cell Painting. While biomarker based screens can enumerate pathway regulators along a unidimensional axis (e.g. pS6 intensity), morphology analysis can identify rich structures spanning multiple biologies. Indeed, across 3 experiments, morphological data clustered genes that affect the same cellular components (e.g. proteasome, golgi/ER, mitochondria inner membrane), molecular functions (e.g. protein glycosylation, cytoskeletal organization, fatty acid synthetic process, chromatin modification) or biological pathways (e.g. EGF, mTOR, TGFb), and was able to do so without prior hypothesis or explicit biomarkers/assays. Morphological guilt-by-association analysis also helps generate novel hypotheses, for example linking uridylyltransferase to cell cycle regulation. In the same example, we note that functionally distinct genes (e.g. PLK4 and TUT1) may have similar morphological phenotypes when deleted, suggesting shared functions for genes that are not known to be related.

We also expand the CellPaint panel with an additional channel, and show that the resulting “CellPaint + 1” design enables the use of a biomarker (e.g. pS6) that provides built-in validation of morphological discoveries from the same experiment. We demonstrate that the rich morphology readout can stratify genes that manifest similarly in biomarker studies (e.g. pS6) into sub-clusters representing molecular mechanisms (e.g. core mTORC1 v.s. mitochondria translation), demonstrating the benefit of coupling high dimensional data with ML analysis. We envision that our hypothesis-free approach can be a useful tool for disease modeling and drug target discovery, particularly when the disease biological pathway or biomarker is not completely understood, similar to recent applications of perturb-seq in disease research ⁷⁶. Techniques that further increase the dimensionality of cellular imaging, such as multiplexed protein / RNA detection ^77,78, may provide finer grained and more diverse phenotypic readouts and reveal additional biological insights.

Pooled optical screening does have some important technical limitations. First, pooled optical screening, like any other high-content image-based assay, can be sensitive to imaging artifacts such as poor focus and imaging channel bleedthrough, requiring careful image QC and staining panel design. Second, ISS signals can be low in certain cell-types with lower expression levels of sgRNA spacer containing transcripts, requiring higher cell input to compensate for lower data recovery (data not shown). Third, pooled optical screening requires instance (per cell) segmentation, which is difficult in cell-types that have overlapping structures (e.g. neurons with extensive and long dendrites, hepatocytes that grow in 3D, data not shown). We have shown that our use of deep-learning assisted cell segmentation helps ameliorate some of these challenges, but further development along those lines will improve the analysis for these more challenging cell types. Fourth, while we demonstrate that high-dimensional screening provides strong phenotypes that cluster effectively, the interpretation of these phenotypes can be more challenging than in transcriptomics. Leveraging both DL analysis for better image representation and classical image analysis for interpretability may be a good combination strategy, as a way to explain why TUT1 may be related to cell cycle regulators for example.

During the final preparation of this manuscript, a similar method (PERISCOPE) emerged that combined Cell Painting and pooled optical screening to build an unbiased morphological atlas ⁷⁹. Our work is differentiated in a number of ways that enable the broad applicability and scaling of our approach. On the experimental side, we have addressed the spectral overlap challenges between the CellPaint panel and the ISS barcodes by using a novel, structural biology inspired RNA-FISH-probe for mitochondrial staining, avoiding the need for custom-conjugation probes. We also developed a closed-loop automation system for in situ sequencing, which was successfully utilized in our druggable-genome discovery screen (Supp Data 6); this capability is critical for scaling throughput and ensuring process and data robustness. Computationally, our use of ML enabled significant improvements to many critical steps during POSH analysis, including base calling, cell segmentation and image representation. For base calling, our approach using a 3-layer convolutional network followed by global registration and transformation of coordinates can be executed without manual alignment of field-of-view images at acquisition time and without manual parameter tuning, yielding improved sensitivity and a simplified lab workflow. Moreover, the use of ML for cell segmentation and guide calling considerably enhanced our ability to utilize more of the cells and hence greatly increased the throughput of our platform. Finally, we note that while we have focused our presentation on our work in A549 cells, we have also successfully applied CellPaint-POSH in HepG2s, as well as in iPSC-derived cortical neurons, motor neurons and stellate cells (data not shown); we believe it can be adapted to most cell types that could be grown and transduced in vitro.

ML-based analytical methods

Morphological analysis is heavily influenced by the features extracted, and we rigorously compared classical image feature extraction with self-supervised DL methods. Classical image analysis uses expert engineered features such as intensity, texture, shape and size of cells and their internal compartments. While such an approach is simple and interpretable, it incorporates human bias and does not capture the full extent of biological variability in the dataset. Deep learning is an alternative approach that is less susceptible to human bias, but existing models are trained on natural image datasets and may not align well with the statistical patterns in bioimaging data. Recently, weakly supervised DL and self-supervised auto-encoder methods have been used to learn embedding directly from cellular images ^80,81. Our work leverages DINO and Vision Transformer (ViT) ⁴⁵, which has been shown to achieve state-of-the-art performance in learning representations from natural images. We demonstrated that self-supervised models with ViT architecture trained on single cell feature embeddings provide better predictive accuracy and improved gene interaction prediction when compared with classical single cell featurization methods and ImageNet pretrained self-supervised models.

The information richness of a high throughput and high dimensionality CellPaint-POSH dataset, combined with the power of deep learning featurization, allows a diverse range of analyses to be conducted based on a single experiment, greatly accelerating biological discovery. For global gene network analysis, we borrowed heavily from standard single-cell RNAseq methods and found that “guilt by association” was one of the most effective ways of identifying known and novel components of pathways. Our work demonstrates that this approach helps identify both known and novel biology via a simple and uniform analysis. As a future direction, other dimensionality reduction techniques that define and explore multi-dimensional latent spaces, similar to singular value decomposition or non-negative matrix factorization used in single-cell RNAseq analysis ^82,83, may further improve the computational resolution of mapping gene-function relationships and help reveal pleiotropic gene function (such as in the case of TUT1).

Importantly, we also showed that the single cell feature representation models trained from one dataset generalized well to an unseen dataset, allowing a pre-trained model to be easily deployable in interpreting future experiments without expensive and time-consuming de-vno training every time. We note that when comparing CP-DINO 1640 with CP-DINO 300, the former captures more semantic structure in the data, even though it is trained with a similarly sized dataset as the latter (both ∼1-1.5million cell tile images). We reason that this is due to the greater degree of phenotypic diversity present in the 1640 dataset thanks to the five-fold more genetic perturbations. We predict that increasing the in-house training data size and (even more importantly) phenotypic variance -e.g. scaling to genome-wide POSH experiments or even experiments done in different cell types - may further improve the DL model and make it more generalizable. In addition, publicly available non-POSH datasets such as the JUMP Cell Painting collection ⁸⁴ can be incorporated into the training data. Notably, the self-supervised ViT architecture requires no labels, and thus can be trained on bioimaging datasets without upfront annotations, circumventing a significant bottleneck and enabling model training on a broad range of internal and external data sets. Thus, over time, we can build richer foundation models that provide an unbiased distillation of the phenotypic landscape of cell morphology. These models can be used to interrogate new gene-phenotype relationships and deliver new insights regarding diverse pathways and biologies, paving the way towards novel biological discoveries.

sgRNA Design, Synthesis and Ordering

Oligo libraries were designed that targeted 1) morphological phenotypes (124 genes, 10 sgRNAs per gene), 2) the druggable genome (1640 genes, 4 genes per sgRNAs) and 3) mechanism of action study (300 genes, 4 sgRNAs per gene). The gene list was then put through a sgRNA selection pipeline that pulls sgRNAs (∼20bps long) from Brunello ⁸⁵ and Toronto KnockOut version 3.0 ⁸⁶ sgRNA repositories. sgRNAs targeting the same region were removed by filtering on a Levenshtein distance of less than 3 and a max subsequence match length of greater than 15. sgRNAs without a BsmBI restriction enzyme cutting site used for golden gate cloning were discarded. The sgRNAs were then sorted by their on-target score per gene as reported if available ⁸⁵. A greedy approach was used to iteratively pick N sgRNAs per gene in the ascending order of the number of sgRNA designs available per genes. In order to allow for correction of at least 1 base-pair error throughout the library, we also set a hamming distance cutoff requirement of 3 for selecting any subsequent sgRNA and select the sgRNAs iteratively until we have sgRNAs for all genes in our input design list. Non-targeting and intergenic sgRNAs were then placed in each library as controls. The last step of the pipeline appends the oligo sequence flanking the sgRNA that allows cloning into our backbone, while including the dialout primers needed for library specific amplification. Libraries were duplicated or triplicated before submitting to a vendor for synthesis, based on the size of the initial sgRNA list. This helped to prevent jack potting and propagation of errors while also keeping the relative ratio of controls to targeting sgRNAs the same. We then ordered the library from Twist Bioscience.

Library Cloning

Library pools received from Twist Bioscience were spun down at 11,000 rcf for 2 mins then resuspended to 5-10ng/mL in IDTE pH8.0 1xTE solution (11-05-01-13, IDT) . We incubated the freshly resuspended library pool for 5 minutes at room temperature then vortexed for 20-30 seconds, then spun down for 2 minutes at 11,000 rcf.

All of our oligo libraries were cloned into a customized lentiviral CROPseq-like vector (plasmid backbone) of which some features include puromycin resistance, TagBFP, padlock sites and Esp3I/BsmBI sites for cloning (Supp Fig 1b). Per library to be cloned, the backbone was digested using 2 ug of plasmid backbone, 2 uL of FastDigest Esp3I (BsmBI) enzyme (ThermoFisher, FD0454), 0.5 uL 100 mMM DTT (Bioworld, 3483-12-3), 2 uL of 10X FastDigest buffer, and molecular biology (molbio) grade water (46-000-CI, Corning) added to a final volume of 50 uL. Solution was then incubated in a thermocycler at 37C for 2 hrs. Then 1 uL of fast alkaline phosphatase (ThermoFisher, EF0652) was added and the sample was incubated at 37C for 45 mins then heat inactivated at 75C for 5 mins. We then followed the E-Gel CloneWell II Agarose Gel protocol to gel-purify the digested backbone. Sample was eluted and DNA was quantified using a Nanodrop. Typical recovery of DNA was 50-70% of the starting solution. Digested vector was then stored at -20C.

A 200uL dialout PCR solution was made in molbio grade water (46-000-CI, Corning) to a final concentration of 1x Q5 PCR (M0492S, NEB), oligo pool (2.5ng final), 1x evagreen (Biotium, 31000T), 0.25 uM dialout primers (see Table 1). qPCR was run on the sample according to the protocol with 1x 98C for 10 seconds > 16x 56C for 30 seconds, 72C for 30 seconds, and 72C for 30 seconds > 1x 72C for 120 seconds. qPCR was run until the last sample entered the exponential phase, up to a maximum of 16 cycles. PCR product was cleaned up using SPRIbeads (Beckman Coulter, B23318) at 0.65x ratio (ex: 32.5uL SPRI beads and 50uL of sample). DNA was eluted typically in 50 uL elution buffer after two 80% ethanol washes. We then quantified the DNA sample using a Qubit 4 Fluorometer (invitrogen, Q33238) and the high sensitivity 1x dsDNA kit (invitrogen, Q33231).

Golden gate assembly was done using the NEB Golden Gate Assembly Kit (BsmBI-v2) E1602 protocol. In short, the NEB reaction was prepared in a PCR strip tube to a final amount of 200 ng of pre-digested CROP-seq backbone, 80 ng of the insert (sgRNA library), 2.5 uL T4 DNA Ligase Buffer, 2.5 uL of NEB Golden Gate Enzyme Mix (BsmBI-v2) Mix (NEB, E1602L), and molbio grade water added to a total of 50uL reaction solution. PCR strip tube with solution was then placed in a pre-warmed thermocycler with a temperature profile of 42C for 60min > 60C for 5min > hold at 4C.

The libraries were then cleaned and concentrated using SPRI using a 1X template:bead ratio (see SPRI cleaning above). Samples were then eluted by adding 15 uL of molbio grade water to the dried beads and keeping ∼12 uL, being sure not to carry through any beads. Next the samples were transformed into electrocompetent cells (Biosearch Technologies, 60242-2). One 0.1cm cuvette (Bio-Rad, 1652089) for each transformation and clean DNA in water were prepared on ice. 1 mL of recovery media (Biosearch Technologies, 80026-1) was transferred to each culture tube in advance. The electrocompetent cells were then thawed on ice. The entire volume of the clean plasmid library (∼10-12 uL) was added to the entire volume of cells (∼50 uL), then mixed gently with a pipette tip 10 times while avoiding bubbles. Entire contents of the DNA+Cell solution was then transferred to the cuvette. The cuvette was then gently tapped to drop cell solution between the metal plates while continuing to avoid bubbles. The solution was then electroporated with the following settings on a Bio-Rad GenePulser Xcell machine (bio-rad, 1652660): 10 uF, 600 Ohms, 1800 Volts. Immediately after the electrical pulse, 1000 uL of Recovery Media was added to the cuvette, then cells were resuspended by pipetting up and down, and the solution was transferred to the 15 mL culture tube (Falcon, 352057). They were recovered at 37 C for 1 hour in an incubator shaking at 300 rpm. Samples were then transferred to 100 mL of LB broth (Corning, 46-050-CM) + Carbenicillin (C2135, Teknova) in 500 mL culture flask (corning, 431401) and incubated at 32.5 C at 160 rpm overnight (16 hrs). Next we plated dilutions to estimate total transformation efficiency as a quick validation. A 0.5 mL aliquot from the culture flask was pipetted into a 1.5 mL eppendorf tube and then diluted twice in recovery media by 1:100 so that the final dilution from the 2 mL recovery solution (before it went into the large 500 mL flask) was 1:10,000 and 1:100,000 respectively. Each dilution was then seeded into carbenicillin agar plates (Teknova, L1011) by adding the full volume of the tube to the plate and then distributing the solution by shaking with seeding beads (Cole-Parmer, UX-01850-33). Plates were then incubated overnight (∼16 hrs) at 37 C. The colonies were counted the next day and coverage was calculated by: Calculated coverage was compared to expected coverage for any major deviations. The remaining was grown in 100 mL cultures in a 500 mL flask at 30 C overnight at 300 rpm and harvested 16 hours post-transformation. Midiprep was then performed on the samples using NucleoBond Xtra EF plasmid purification, with a final elution volume of 1500 uL. Samples were stored at -20C for long term storage. Dialout and plasmid library distributions were verified by deep sequencing.

Lentivirus Production

Lenti-X HEK293T cells (632180, Takara) were thawed and cultured in 10% FBS (Gemini, 100-800) in high glucose DMEM (gibco, 10569-010) for at least 1 week before transfection. Cells were passaged every 3 to 4 days when they reached 70-80% confluence. On day 0 of transfection, cells were passaged into 6 well plates (Falcon, 353046) at 1 million cells per well and 2mLs per well of media. On Day 1, cells were transfected using an edited version of the lipofectamine 3000 reagents and protocol (ThermoFisher, L3000075). In short, two solutions were first made. The first with 167 uL per well of OptiMEM (gibco, 51985034) and 14.4 uL of L-3000 reagent. The second solution was made with the following values per well: 167 uL of OptiMEM, 1056 ng of psPAX2 (addgene, plasmid #12260), 704 ng of pMD2.G (addgene, plasmid #12259), 1408 ng of plasmid library, and 12 uL of P3000 reagent. For the second mix it was important to add all of the plasmids first before adding P3000, to prevent the plasmids from crashing out of solution. After mixing the two solutions independently by quickly vortexing, we then added the first solution to the second solution, mixed well and then let incubate for 10 minutes. We then added 333 uL of solution drop wise to each well. On day 2 cells were treated with a final concentration of 2mM viralboost (Alstem, VB100) solution in media by adding dropwise to the wells. On day 3, the virus was harvested by collecting the supernatant and straining it through a 0.45um syringe filter (millipore sigma, SLHV004SL). The virus solution was aliquoted into cryovials and immediately frozen at -80C for later use. Virus was titered on A549s at a wide range of concentrations to determine appropriate volumes to achieve an MOI of approximately 0.25 for screens. Titer was measured via flow by percent BFP positive cells and then the following equation was used to scale the transduction for the screen:

A549-Cas9 line generation

CAG-Blast-Cas9 plasmid was purchased from Horizon Discovery (Horizon Discovery, CAS10141). A cas9 lentivirus was created as described above. A549s in culture were then transduced with the cas9 lentivirus. Using a selection titer strategy with blasticidin, we selected cells at an MOI of 1. The A549-Cas9 line was then expanded and banked in 10% DMSO (Sigma, D2650-100ML), 10% FBS (Gemini, 100-800), in low glucose DMEM (gibco, 11885-084).

Cell Culture and Lentivirus Transduction

A549-Cas9 line was thawed into T225 flasks (Thermo Scientific, 159934) and then expanded and selected for Cas9 positive cells with 10 ug/mL blasticidin (gibco, A11139-03) in 10% FBS (Gemini, 100-800) in low glucose DMEM (gibco, 11885-084) for 6 days. They were passaged at 70-80% confluence every 3 or 4 days using TrypLE (gibco, 12604039). On day 6 post thaw, cells were infected with lentivirus plasmid pool at an MOI of 0.3 for 24 hours. Cells were then passaged into multiple T225s for expansion and selected with 1ug/mL puromycin (gibco, A11138-03). Blasticidin selection was discontinued at this point. Transductions and selections were designed to ensure that cell coverage never fell below 300x cells per sgRNA. Cell pellets were collected for next generation sequencing at day 1, 4, 7, and 10 post transduction by spinning down passaged cells at 300g for 5 minutes, aspirating supernatant, and immediately storing at -80C. On day 7 post transduction, cells were seeded into 6-well glass bottom, black plates (Cellvis, P06-1.5H-N) and taken off of puromycin selection. Cells were cultured in 10% FBS in DMEM medium until day 10 post transduction.

Next Generation Sequencing Sample Processing

NucleoSpin Blood kit (Macherey-Nagel, 740951.50) was used to isolate the genomic DNA. We ran indexing PCR on our NGS samples by creating the following mix per sample: 1x Kapa Hifi readymix HS (Roche, 7958935001), 0.5uM common primer (Table 1), 0.5uM indexing primer (Table 1), 5% DMSO (Sigma, D2650-100ML), 1-10ug of DNA sample, filled up to 100uL with molbio grade water. The thermal profile for the PCR is 95C for 5 min > 95C for 30sec, 56C for 30sec, 72C for 20 sec, 24 cycles > 72C C for 10min > hold at 12C. Samples were split into at least 4 wells of a 96wp PCR plate to minimize jackpotting. DNA was cleaned and quantified using the protocol described in the PCR cleanup section above. Samples were sequenced with the MiSeq V2 with PE50 V2 kit with custom primer according to manufacturer’s recommendations.

Demultiplexing was done by running the BCL files through a BCL-to-fastq pipeline. Fastqs were then used to align sgRNAs to our sgRNA library using an exact match approach of up to 18 base pairs of the 20 total base pairs of the sgRNA. Log2 enrichment was calculated by first getting the cell fraction as: Fold Change was then calculated as: Log2 of Fold Change (FC) was then calculated. Lastly, the FoldChange was then normalized to the intergenic by doing the following: Plots were made using normalized log2 FC.

MitoProbe

Coordinates for the human mitochondrial ribosome from the Protein DataBase (PDB) (3J9M) were used. We used Visual Molecular Dynamics (VMD) for initial inspection of the PDB. For all other calculations, we used standard python libraries and MDTraj.

Within 3J9M, for all RNA sequences between the lengths of 15 to 25 (inclusive), we computed the SASA for all atoms in the sequence using MDTraj. We then computed the mean, median, and standard deviation of the SASA for all RNA sequences. We then filtered the sequences to only keep sequences that had a fully continuously observed sequence in the PDB. This yielded roughly 4379 sequences though this also contains overlapping sequences (i.e length 15 sequences are in length 16 to 25 length sequences). We sorted these sequences by their mean and median atomic SASA values. We sorted the list by its SASA values and then filtered the sequences by their GC content, absence of repetitive bases, and intra-sequence distances. Out of the final 8 sequences, we picked the first 2 sequences (lengths 20 and 23 base pairs) with the highest SASA (Figure 1) and a large physical separation within the PDB. These selected sequences occupy distinct positions on the mt-ribosome. Oligonucleotides of these sequences were ordered with all uracil (U) bases in the sequence replaced with thymine (T), coupled to 5’ dyes and purified by HPLC (Table 1).

Cell-Paint POSH protocol

Druggable Genome Screen Improvements

Computational Methods

Data Capture and Storage

The images acquired on the microscope are exported in raw tiff format and the corresponding metadata for each image as specified in Table 5 are recorded.

CellPaint Image Processing

The raw fluorescence field of view images acquired from the microscope have uneven illumination due to camera and vignetting artifacts. In order to create a flat-field illumination image, the raw images are scaled by an illumination correction function as estimated from the mean over field of view images followed by gaussian smoothing as described ⁸⁷. The illumination corrected images are then min-max percentile (min percentile = 0.1, max percentile = 99.9) scaled based on per-channel intensity percentiles computed from plate-level intensity histograms and converted to 8-bit unsigned-integer type (min value = 0, max value = 255) for subsequent processing.The field of view images are segmented using CellPose ⁸⁸ to obtain single cell and nuclei boundaries. The field of view images are then cropped around the centroids of each of the segmented nuclei and masked by the corresponding cell segmentation mask to create tiles with a single cell in context.

in-situ sequencing by synthesis Image Processing

The sequencing by synthesis procedure described generates a dataset in which the full plate is imaged several times, with 4-color stationary dots showing variable fluorescent signatures corresponding to A/C/T/G that need to be converted to sequencing base calls.

Single Cell sgRNA identity assignment

To assign sgRNA barcode identities to cropped single cell images obtained from CellPaint, we transform the barcode locations, sgRNA sequence and their corresponding gene KO identities from the first cycle of SBS to the CellPaint acquisition space using the same odometry transformation methodology described in the registration section above. The barcodes in CellPaint acquisition space are then assigned to single cell contexts using the cell boundaries and then mapped to corresponding the cropped tiles to obtain a labeled perturbation dataset of single cell images.

Explicitly Engineered Image Featurization (CellStats)

To localize the signal, we capitalize on the two segmented areas and derive features that are computed in the nucleus, the cytoplasm, and also the perinuclear regions by iteratively dilating the segmentation mask corresponding to the nucleus. First, we start by capturing characteristics of the pixel intensity empirical distributions corresponding to each localization mask. In particular, we extract percentiles from the empirical distributions at discrete steps, along with their means, variances, and standard deviations. Additionally, we compute pairwise Pearson and Spearman correlations among the different channels in our assays. Second, focusing on the localization masks of the nucleus and the cytoplasm, we capture characteristics of each cell’s geometry. Key attributes of cell geometry include the area of the nucleus and the cytoplasm and its convexity ⁹⁰. Lastly, we extract features that have been heavily utilized by the computer vision community to characterize the textures emerging in the nucleus and the cytoplasm from the selected types of staining. In particular, statistics derived from suitable wavelet transforms and region covariance descriptors are employed towards capturing such information. Detailed lists of the features in each of the aforementioned categories are provided in Supp Data 2.

Self-Supervised Learning of Single-Cell Image Representations

Self-supervised learning pertains to methodologies that are capable of learning latent space representations over input distributions without having an explicit label assigned to them. Recently, self-supervised learning methods [DINO, SimCLR, CLIP] have been shown to improve generalization and quality of learned representations as compared to supervised learning methods ^45,91,92. In this work, we use DINO ⁴⁵ a self-supervised learning technique that uses knowledge distillation between student and teacher networks parameterized by θs and θ_t respectively to learn rich representations from images. Both networks are trained to output a probability vector over K dimensions computed as: where τ_s is the temperature parameter that controls the sharpness of the output distribution. The student network is trained by minimizing the cross-entropy loss over outputs (as given by Eq below) with the teacher network using stochastic gradient descent. where H(a, b) = − a log b, x_g1 and x_g2 are augmented global views and V is a set of augmented local crop views of the input images. The teacher network is updated using an exponential moving average of past weights of the student network as given by: Our motivation to choose DINO as our self-supervised framework stems from their superior performance in k-Nearest Neighbor benchmarks on ImageNet dataset as reported ⁴⁵. We also use a Vision Transformer architecture ⁹³ to parameterize both the student and the teacher network parameters.

Feature Normalization

In order to adjust for intensity variations and batch effects introduced by assay and imaging systems, we normalize each feature embedding by the feature embeddings of non-targeting controls in each well. We employ Robust Center-Scale Normalization where for each well, we subtract all the embeddings by the median of non-targeting controls in the well and scale by the median absolute deviation of the non-targeting controls.

Feature Aggregation

For extrinsic evaluation of the usefulness of embeddings and understanding biological relationships from embeddings, we first filter the embeddings based on their phenotype score computed by a linear binary classification model between nontargeting and every targeting gene KO. We then compute a single sgRNA KO representation by taking the median over all the embeddings corresponding to that sgRNA KO. The gene-KO level aggregate representations are obtained by taking a mean over all the representations of sgRNAs targeting the gene.

Differential Morphology Analysis and Gene Enrichment Analysis

In several cases, more focused analysis of a single knockout or single feature was necessary. In the case of gene-wise differential morphology analysis, the Z-scored cell populations containing the given knockout, as well as the non-targeting/intergenic control cells within the same pool, were subset from the rest of the pool. Between these two samples, a two-sample Komogorov-Smirnov (K-S) test was applied for each feature with the Z score of the sub-population relative to controls. Due to the large number of features, as well as their redundancy, clear trends in significantly-changed features were annotated; more specific feature results were determined for all genes within the Druggable Genome Screen and created as a report file. A similar approach was applied in the case of feature-wise gene enrichment analysis, in which, for the given feature, each subset of cells receiving the same genetic knockout was compared to non-targeting/intergenic via a two-sample two-tailed K-S test.

Comparison to FACS Binned 1-Dimensional Analysis

We sought to compare how POSH-based single-cell data on a single feature would compare to a simulated binned FACS screen using the same data. For the single cell resolution based analysis, K-S testing and mean Z were determined as described above for mean cytoplasmic pS6 intensities. To simulate a FACS screen using a pseudo-binning approach, we used the well-normalized Z score matrix and determined the 0-15%, 15-30%, 70-85%, and 85-100% percentile bins, and assigned all cells within the screen to these respective bins. We then assumed perfect sequencing and assignment of the sgRNAs within cells that had been assigned to the bins. MAUDE analysis ⁹⁴ was then conducted to predict the likely distribution shift of cytoplasmic phospho-S6 within each knockout subset, leading to predicted mean Z and FDR for each gene.

Evaluation of Single-Cell Embeddings by Linear Classifical Model

We evaluate the intrinsic representational power of the different embedding methodologies by measuring their performance on downstream predictive tasks. We evaluate the embeddings by their performance on multiple binary classification tasks for discriminating between non targeting controls and samples corresponding to a gene KO. For each gene KO we train a Elastic Net Logistic Regression model with the task of binary classification between nontargeting and gene KO samples with a leave-one-well-out cross validation schema for all the wells in the experiment. The binary prediction area under the receiver operating characteristic curve (AU-ROC) is computed on the validation samples for each gene in every validation split and the median over AU-ROC is obtained as the scoring function for each embedding. We report the number of gene KOs with an AU-ROC > 0.55 and the mean AU-ROC for the subset of genes on prediction tasks for each embedding methodology as an objective metric for evaluation and comparison.

Evaluation of Biological Meaningfulness of Aggregate Embeddings

UMAP Visualization of Embeddings

We visualize the aggregate gene-KO embeddings using a UMAP projection with cosine similarity as the metric. The size of each gene-KO node in the UMAP is adjusted based on the similarity of sgRNAs targeting the gene. The gene-KO nodes were colored using Leiden community detection algorithm ³⁷ on the UMAP projection.