Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Jonathan Weissman (weissman@wi.mit.edu).

Materials availability

Plasmids and CRISPRi sgRNA libraries generated in this study have been deposited to Addgene.

Data and code availability

• Raw sequencing data are deposited on SRA under BioProject PRJNA831566. key resources table An interactive data browser including processed, downloadable single-cell and pseudobulk populations is available at http://gwps.wi.mit.edu.

• Our codebase for Perturb-seq analysis is available at https://github.com/thomasmaxwellnorman/Perturbseq_GI and https://github.com/josephreplogle/guide_calling.

• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Cell culture and lentiviral production

K562 cells were grown in RPMI-1640 with 25 mM HEPES, 2.0 g/l NaHCO3, and 0.3 g/l L-glutamine supplemented with 10% FBS, 2 mM glutamine, 100 units/ml penicillin, and 100 μg/ml streptomycin. hTERT-immortalized RPE1 cells (ATCC, CRL-4000) were grown in DMEM:F12 supplemented with 10% FBS, 0.01 mg/ml hygromycin B, 100 units/ml penicillin, and 100 μg/ml streptomycin. HEK293T cells were used for generation of lentivirus, and grown in DMEM supplemented with 10% FBS, 100 units/ml penicillin and 100 μg/ml streptomycin. Lentivirus was produced by co-transfecting HEK293T cells with transfer plasmids and standard packaging vectors using TransIT-LTI Transfection Reagent (Mirus, MIR 2306).

Cell line generation

CRISPRi K562 cells expressing dCas9-BFP-KRAB (KOX1-derived) were obtained from Gilbert et al. (2014). CRISPRi RPE1 cells expressing dCas9-BFP-KRAB (KOX1-derived) were obtained from Jost et al. (2017) and only used for growth screens. CRISPRi RPE1 cells were generated by stably transducing RPE1 cells (ATCC, CRL-4000) with lentivirus expressing ZIM3 KRAB-dCas9-P2A-BFP from a UCOE-SFFV promoter (pJB108) and sorting for BFP+ cells stably expressing the construct using fluorescence activated cell sorting. Cell lines were verified by monitoring BFP fluorescence over several generations to confirm stable integration and confirming knockdown of select surface markers by flow cytometry.

Library design and cloning

A distinct set of genes was targeted for each of the three large-scale Perturb-seq experiments. For the K562 day 8 genome-scale experiment, we targeted (i) genes expressed in K562 cells (ii) transcription factors as detailed in Lambert et al. (2018) (iii) Cancer Dependency Map common essential genes as defined in 20Q1 (iv) non-targeting control sgRNAs accounting for 5% of the total library. To define expressed genes in K562 cells, we used a combination of bulk RNA-seq data from ENCODE (https://www.encodeproject.org/files/ENCFF717EVE/) and 10x Genomics 3’ single-cell RNA-seq data (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146194), selecting a set of genes accounting for ~99% of aligned reads in both datasets. For the K562 day 6 essential-scale experiment, we targeted (i) Cancer Dependency Map common essential genes as defined in 20Q1 (ii) non-targeting control sgRNAs accounting for 5% of the total library. For the RPE1 day 7 essential-scale experiment, we targeted (i) 20Q1 Cancer Dependency Map common essential genes (https://depmap.org/portal/download/) (ii) a number of hand-selected genes with interesting phenotypes in the K562 genome-wide Perturb-seq dataset (iii) non-targeting control sgRNAs accounting for 5% of the total library. To define control perturbations, we randomly sampled non-targeting control perturbations from Horlbeck et al. (2016). A small number of genes were lost in this pipeline due to changes in gene annotation between datasets.

To minimize library size while maximizing knockdown, multiplexed CRISPRi libraries were constructed which targeted each gene with two unique sgRNAs expressed from tandem U6 expression cassettes in a single lentiviral vector, as previously described in Replogle et al. (2020). The Horlbeck et al. (2016) CRISPRi sgRNA libraries were used as a source of sgRNAs targeting each gene, with the optimal sgRNA pair targeting each gene selected based on a balance of empirical data with computational predictions. For strong essential genes (defined by a p-value<0.001 and gamma<-0.2 in the Horlbeck et al. [2016] CRISPRi growth screen), sgRNAs were ranked by growth. Then, for genes that produced a significant phenotype in previous CRISPRi screens, sgRNAs were ranked by a discriminant score multiplying the negative log₁₀ p-value by the effect size. Finally, for genes without any empirical evidence, sgRNAs were ranked according to the Horlbeck et al. (2016) hCRISPRi v2.1 algorithm. The full sgRNA content of the K562 day 8 genome-scale library, K562 day 6 essential-wide library, and RPE1 day 7 essential-wide library can be found in Table S1.

We adapted the protocol previously described in Replogle et al. (2020) to clone libraries with capture sequences for 3’ direct capture Perturb-seq. Briefly, an sgRNA lentiviral expression vector (pRS275/pJR101) was derived from the parental pJR85 (Addgene #140095), modified to incorporate a GFP fluorescent marker to avoid spectral overlap with BFP+ CRISPRi constructs and a UCOE element upstream of the EF1alpha promoter to prevent silencing (Table S4). A two-step restriction enzyme digestion and ligation cloning of oligos into pRS275/pJR101 was performed to maintain coupling of sgRNAs targeting the same gene. Oligos encoding the targeting regions of dual-sgRNA pairs were synthesized as an oligonucleotide pool (Twist Bioscences) with the structure: 5’- PCR adapter - CCACCTTGTTG – targeting region A - gtttcagagcgagacgtgcctgcaggatacgtctcagaaacatg – targeting region B - GTTTAAGAGCTAAGCTG - PCR adapter-3’. When ordering oligos, the representation of essential genes was increased to compensate for growth phenotypes (see below). Oligo pools were amplified, digested with BstXI/BlpI, and ligated into pRS275/pJR101. To add an sgRNA constant region and U6 promoter to the vector, pJR89 (Addgene #140096) was BsmBI-digested and ligated into the intermediate library.

K562 and RPE1 growth screens

Pooled sgRNA growth screens in K562 cells were used to quantify growth phenotypes of sgRNA pairs targeting expressed genes. CRISPRi K562 cells expressing dCas9-BFP-KRAB were transduced with lentiviral particles encoding the dual-sgRNA library by spinfection (1000g) with polybrene (8 ug/ml) to obtain an infection rate of~25%-35%. Screens were performed in biological replicate with the aim of maintaining 1000 cells per library element for the duration of the screen. Between day 2 and day 6 post-transduction, cells were selected for lentiviral infection using 1 ug/mL puromycin, replenished every 24 hours. On day 7 post-transduction, an aliquot of cells was harvested as an initial time point) The rest of the cell population was passaged for 10 more days and collected at final time point.

Pooled sgRNA growth screens in RPE1 cells were used to quantify growth phenotypes of sgRNA pairs targeting common essential genes. The CRISPRi RPE1 cell line expressing dCas9-BFP-KRAB was used for growth screens which took place before the publication of the next-generation ZIM3 KRAB domain. Cells were transduced in biological replicate with lentiviral particles encoding the dual-sgRNA library by replating cells into virus-laden media with polybrene (8 ug/ml) to obtain an infection rate of ~45%. Because RPE1 cells are puromycin resistant, we performed the screen without selection for sgRNA-infected cells, nonetheless maintaining an infection rate-corrected 1000 cells per library element for the duration of the screen. On day 6 post-transduction, an aliquot of cells was harvested as an final time point for direct comparison to the abundances in the plasmid library.

For both K562 and RPE1 growth screens, DNA libraries of the initial and final samples were prepared for deep sequencing by genomic DNA isolation and PCR amplification of dual-sgRNA amplicons as described previously (Nuñez et al., 2021; Replogle et al., 2020). First, a NucleoSpin Blood XL kit (Macherey–Nagel) was used to extract genomic DNA (gDNA) from cells. Then, isolated gDNA or plasmid DNA was amplified by 22 cycles (gDNA) or 13 cycles (plasmid DNA) of PCR using NEBNext Ultra II Q5 PCR MasterMix (NEB), appending Illumina adaptors and sample indices (oJR234 forward primer: 5’-AATGATACGGCGACCACCGAGATCTACACCGCGGTCTGTA TCCCTTGGAGAACCACCT-3’; index primers 5’-CAAGCAGAAGACGGCATACGAGATnnnnnGCGGCCGGCTGTTTCCA GCTTAGCTCTTAAA-3’). Amplicons were isolated by a 0.5–0.65X SPRI bead selection (SPRIselect Beckman Coulter #B23318). Sequencing was performed on a NovaSeq 6000 (Illumina) using a 19 bp read 1, 19 bp read 2, and 5 bp index read 1 with custom sequencing primers oJR326 (custom read 1, 5’- CGCGGTCTGTATCCCTTGGAGAACCACCTTGTTGG-3’), oJR328 (custom read 2, 5’- GCGGCCGGC TGTTTCCAGCTTAGCTCTTAAAC-3’), and oJR327 (custom index read 1, 5’- GTTTAAGAGCTAAGCTGGAAACAGCCGGCCGC-3’).

Perturb-seq experiments

The selection of time points for our experiments is based on a combination of previously published CRISPR screens, Perturb-seq experiments, and the goals of our experiment. The constraints on the design of CRISPR growth screens differ significantly from Perturb-seq experiments. In growth screens, an amplification of signal occurs over time as cells drop out of the population, so experiments often compare representation between an early time point (~3–5 days post-transduction) and a much later final timepoint (~14–28 days post-transduction). In contrast, in Perturb-seq screens, the phenotype is measured directly from the perturbed cells that are sampled on the day of scRNA-seq. Relatively earlier timepoints may then be advantageous because:

1. libraries remain more balanced, especially when studying perturbations targeting essential genes that are prone to dropping out over time at the representation levels typically used in scRNA-seq experiments; and

2. more ‘‘direct’’ phenotypic consequences of the genetic perturbation are observed (i.e., the transcriptomes reflect the cellular response to perturbation rather than later indirect consequences like cell death).

In contrast, the possible advantages of sampling cells at later timepoints are:

1. time is required to allow for CRISPR machinery to be expressed, the genetic perturbation to occur (in this case CRISPRi), and finally protein depletion to occur; and

2. for some perturbations, longer time points might be required in order to observe a phenotype (e.g., perturbations that result in buildup of cellular metabolites).

In designing our experiment, we wanted to ensure that we would sample the phenotypic consequences of perturbing essential genes which would quickly deplete from our library. We thus chose to sample at two time points in K562 cells as the phenotypic effects of different genetic perturbations can manifest at variable time points based on technical and biological factors. While each gene depletes and causes a cellular phenotype based on unique characteristics, the majority of genes are widely accepted to have phenotypes in between approximately day 6 to day 8 of the screen.

To perform our K562 day 8 genome-scale Perturb-seq experiment, library lentivirus was packaged into lentivirus in 293T cells and empirically measured in K562 cells to obtain viral titers. CRISPRi K562 cells were transduced via spinfection (1000g) with polybrene (8 ug/ml) with the target of obtaining an infection rate of ~10%. Cells were maintained at a viability of >90%, a coverage of 1000 cells per library element, and a density of 250,000 to 1,000,000 cells/ml for the course of the experiment. Three days post-transduction, an infection rate of 14% was measured, and cells were sorted to near purity by FACS (FACSAria2, BD Biosciences), using GFP as a marker for sgRNA vector transduction. Eight days post infection, the cells were measured to be 97% GFP+ (LSR2, BD Biosciences), >90% viable, and at a concentration of ~800,000 cells/ml (Countess II, ThermoFisher). Cells were prepared for single-cell RNA-sequencing by resuspension in 1X PBS with 0.04% BSA as detailed in the 10x Genomics Single Cell Protocols Cell Preparation Guide (10x Genomics, CG00053 Rev C). Cells were then separated into droplet emulsions using the Chromium Controller (10x Genomics) with Chromium Single-Cell 3° Gel Beads v3 (10x Genomics, PN-1000075 and PN-1000153) across 273 ‘‘lanes’’/’’GEM groups’’ following the 10x Genomics Chromium Single Cell 3ʹ Reagent Kits v3 User Guide with Feature Barcode technology for CRISPR Screening (CG000184 Rev C) with the goal of recovering ~15,000 cells per GEM group before filtering. Because the formation of droplet emulsions occurred in batches of 8 GEM groups over several hours, fresh populations of cells were obtained every hour to prevent alterations in single-cell transcriptomes.

To perform our K562 day 6 essential-scale Perturb-seq experiment, library lentivirus was packaged into lentivirus in 293T cells and empirically measured in K562 cells to obtain viral titers. CRISPRi K562 cells were transduced via spinfection (1000g) with polybrene (8 ug/ml) with the target of obtaining an infection rate of ~10% with maintenance of cells as described above. Three days post-transduction, an infection rate of 15% was measured, and cells were sorted to near purity by FACS (FACSAria2, BD Biosciences), using GFP as a marker for sgRNA vector transduction. Six days post infection, the cells were measured to be 93% GFP+ (LSR2, BD Biosciences), >90% viable, and at a concentration of ~600,000 cells/ml (Countess II, ThermoFisher). Cells were prepared for single-cell RNA-sequencing by resuspension in 1X PBS with 0.04% BSA as detailed in the 10x Genomics Single Cell Protocols Cell Preparation Guide (10x Genomics, CG00053 Rev C). Cells were then separated into droplet emulsions using the Chromium Controller (10x Genomics) with Chromium Single-Cell 3′ Gel Beads v3 (10x Genomics, PN-1000075 and PN-1000153) across 48 ‘‘lanes’’/’’GEM groups’’ following the 10x Genomics Chromium Single Cell 3ʹ Reagent Kits v3 User Guide with Feature Barcode technology for CRISPR Screening (CG000184 Rev C) with the goal of recovering ~15,000 cells per GEM group before filtering.

To perform our RPE1 day 7 essential-scale Perturb-seq experiment, library lentivirus was packaged into lentivirus in 293T cells and empirically measured in RPE1 cells to obtain viral titers. CRISPRi RPE1 cells expressing ZIM3 KRAB-dCas9-P2A-BFP were transduced via replating into virus-laden media with polybrene (8 ug/ml) with the target of obtaining an infection rate of ~10%. Three days post-transduction, an infection rate of 7% was measured, and cells were sorted to near purity by FACS (FACSAria2, BD Biosciences), using GFP as a marker for sgRNA vector transduction. Seven days post infection, the cells were measured to be 86% GFP+ (LSR2, BD Biosciences) and >95% viable (Countess II, ThermoFisher). After trypsin dissociation, cells were prepared for single-cell RNA-sequencing by resuspension in 1X PBS with 0.04% BSA as detailed in the 10x Genomics Single Cell Protocols Cell Preparation Guide (10x Genomics, CG00053 Rev C). Cells were then separated into droplet emulsions using the Chromium Controller (10x Genomics) with Chromium Single-Cell 3′ Gel Beads v3 (10x Genomics, PN-1000075 and PN-1000153) across 56 ‘‘lanes’’/’’GEM groups’’ following the 10x Genomics Chromium Single Cell 3ʹ Reagent Kits v3 User Guide with Feature Barcode technology for CRISPR Screening (CG000184 Rev C) with the goal of recovering ~15,000 cells per GEM group before filtering.

Perturb-seq library preparation and sequencing

For preparation of gene expression and sgRNA libraries, samples were processed according to 10x Genomics Chromium Single Cell 3ʹ Reagent Kits v3 User Guide with Feature Barcode technology for CRISPR Screening (CG000184 Rev C). To allow for parallel library preparation, samples were arranged in 96-well plates with magnetic selections conducted on an Alpaqua Catalyst 96 plate (#A000550). For sequencing, mRNA and sgRNA libraries were pooled to avoid index collisions at a 10:1 ratio. Libraries were sequenced on both (i) a NovaSeq 6000 (Illumina) according to the 10x Genomics User Guide and (ii) the Ultima Genomics ultra-high throughput sequencing platform.

For sequencing on the Ultima Genomics (UG) platform, final 10x libraries were converted using conversion primers that anneal to the R1 and R2 regions of the 10X library and contain a UG-specific adapter sequence overhang and sample index. 8 PCR cycles were used for conversion. Converted libraries were bead purified and quantified. After pooling libraries, pools were seeded and clonally amplified on UG sequencing beads and sequenced on a UG prototype Sequencer. Single reads were generated from the 10x 3’ libraries, reading the 10X cell barcode, unique molecular identifier (UMI), and 3’ end of the cDNA transcript. A specific sequencing protocol including a high volume of dT nucleotides was used to accommodate the high nucleotide consumption in the poly (dT) stretch of the cDNA. Following sequencing, the single reads were quality-trimmed and split into two sub-sequences corresponding to Read1 (10X cell barcode and UMI) and Read2 (cDNA), which were used as input to Cell Ranger for alignment.

rRNA analyses

K562s expressing Zim3-dCas9-2A-BFP were spinfected in biological duplicate (targeting sgRNAs) or quadruplicate (non-targeting sgRNAs) with lentivirus expressing GFP and an sgRNA. Two days after spinfection, the cells were sorted for GFP+ on a BD ARIA II. Sort purity was generally >95%. After the sort, cells were maintained in media supplemented with 4 ug/ml puromycin for four days and then recovered for two days. Cells were counted, collected by centrifugation, and harvested by vigorous vortexing in Tri Reagent (ThermoFisher AM9738).

RNA was extracted with chloroform according to the manufacturer’s instructions, quantified by nanodrop, and snap frozen. Small samples were diluted to 200 ng/ul and run on Bioanalyzer RNA nano chips (Agilent 5067-1511) according to the manufacturer’s instructions. Runs were aligned to the 18s peak and signal intensity was normalized to total RNA area.

Integrator co-depletion

K562s expressing Zim3-dCas9-2A-BFP were spinfected with lentivirus expressing GFP and an sgRNA. Two days after spinfection, the cells were sorted for GFP+ on a BD ARIA II. Sort purity was generally >95%. After the sort, cells were maintained in media supplemented with 4 ug/ml puromycin for four days and then recovered for two days. Cells were counted, washed twice with DPBS, and collected as pellets. The pellets were resuspended in SDS lysis buffer (100 mM Tris pH 8.0, 1% SDS), thermomixed at 95°/1500 RPM for thirty minutes, aliquoted, and snap-frozen.

Integrator co-immunoprecipitation

Human expression plasmids encoding codon-optimized INTS10 or His8-INTS10 were synthesized (Twist Bioscience) and transfected into HEK 293T/17 cells (ATCC CRL-11268) with FuGene HD (Promega E2311) according to the manufacturer’s protocol. Two days later, the cells were washed twice with DPBS and harvested with IP lysis buffer (25 mM Tris-HCl pH 7.4, 150 mM NaCl, 1 mM EDTA, 1% NP-40, 5% glycerol; ThermoFisher 87787) supplemented with protease inhibitors (ThermoFisher A32965). Lysates were nutated at 4° for 30 mins, clarified by centrifugation at 12,000xg for 10 minutes, and snap-frozen. Concentrations were measured with the BCA assay (ThermoFisher 23225).

Lysates were thawed on ice, supplemented with imidazole to 10 mM, and nutated at 4° for 30 minutes with cobalt magnetic beads (ThermoFisher 10103D) pre-equilibrated in IP lysis buffer + 10 mM imidazole. The beads were separated on a magnet, washed twice with lysis buffer + 10 mM imidazole, and eluted with lysis buffer + 300 mM imidazole.

Integrator purification

Human expression plasmids encoding codon-optimized HIS-INTS10, INTS13, INTS14, and C7orf26 were synthesized (Twist Bioscience) and co-transfected with ExpiFectamine 293 (ThermoFisher A14524) into Expi293 cells (ThermoFisher A14527) maintained in Expi293 medium (ThermoFisher A1435101) according to the manufacturer’s instructions. The cells were harvested after four days and snap frozen.

The pellets were resuspended in CHAPS Lysis Buffer (50 mM HEPES pH 8.0, 300 mM NaCl, 0.2% CHAPS, 10% glycerol, 1 mM TCEP, 1 mM EDTA, 0.5 mM PMSF, 1x protease inhibitors, 0.002% benzonase) and stirred at 4° for 30 minutes. The lysates were clarified at 120,000xg for 30 minutes, supplemented with 15 mM imidazole, and nutated for an hour with Ni-NTA agarose beads (ThermoFisher 25215) pre-equilibrated in CHAPS lysis buffer + 15 mM imidazole. The beads were loaded into a gravity column, washed with >10 volumes of wash buffer (50 mM HEPES pH 8.0, 300 mM NaCl, 10% glycerol, 1 mM TCEP, 1 mM EDTA, 15 mM imidazole), and eluted with wash buffer supplemented with 250 mM imidazole. The eluate was concentrated and buffer exchanged into SEC buffer (50 mM HEPES pH 8.0, 150 mM KCl, 10% glycerol, 1 mM EDTA) by ultrafiltration, and snap frozen.

The eluate was thawed on ice, passed through a 0.2 μM PES filter, and loaded onto an Superdex 200 Increase 10/300 GL column pre-equilibrated with SEC buffer. Fractions were collected and flash frozen.

Drosophila Integrator biochemistry

Integrator PRO-seq

Pro-seq was conducted largely according to published protocols with slight modifications (Reimer et al., 2021). K562s expressing dCas9-BFP-KRAB were spinfected with lentivirus expressing GFP and an sgRNA. Two days after spinfection, the cells were sorted for GFP+ on a BD ARIA II. Sort purity was generally >95%. After the sort, cells were maintained in media supplemented with 4 ug/ml puromycin for three days and then recovered for two days.

Cells were counted, harvested by centrifugation, and washed with cold DPBS. All subsequent steps took place at 4° or on ice. All solutions were made with RNase-free reagents and were 0.2 μm filtered and chilled before use. 12 million cells were pelleted by centrifugation, resuspended in 250ul of buffer W (10 mM Tris pH 8.0, 10 mM KCl, 250 mM sucrose, 5 mM MgCl2, 1 mM EGTA, 0.5 mM DTT, 10% glycerol, 1x protease inhibitor [ThermoFisher A32965], and 0.02% v/v SUPERase-In RNase inhibitor [AM2694], strained, and transferred to conical tubes that had been coated with 1% BSA in PBS overnight. The cells were permeabilized by dilution in 10 ml of buffer P (buffer W + 0.1% v/v Igepal CA-630 + 0.05% v/v Tween-20) and incubated for 5 minutes. The permeabilized cells were harvested by centrifugation at 400xg for 5 minutes, resuspended in 10 mL of buffer W, harvested by centrifugation at 400xg for 5 minutes, and resuspended in 250 ul buffer F (50 mM Tris pH 8.0, 40% v/v glycerol, 5 mM MgCl2, 1.1 mM EDTA, 0.5 mM DTT, 0.02% v/v SUPERase-In RNase inhibitor). ≥97% permeabilization efficiency was confirmed on a NucleoCounter NC-202 and permeabilized cells were snap frozen.

PRO-seq libraries were generated and sequenced by the Nascent Transcriptomics Core at Harvard Medical School according to their standard protocol. PRO-seq data were aligned and quantified using STAR (version 2.7.9a) with parameters alignEndsType=Local, outFilterMultimapNmax=20, outFilterScoreMinOverLread=0.3, and outFilterMatchNminOverLread=0.3. For comparison of gene-level expression profiles, gene counts were normalized for sequencing depth (reads per million), log-transformed, and subset to well-expressed genes (n=758 genes with >3000 rpm). Then, Spearman’s correlation was used to compare the similarity of expression profiles.

SDS-PAGE and western blotting

Samples were mixed with sample loading buffer (Licor 928-40004) supplemented with DTT and incubated at 95° for 5 minutes. SDS-PAGE was performed with pre-cast 4–12% gradient gels (ThermoFisher NW04127BOX) in MOPS (ThermoFisher B000102) according to the manufacturer’s instructions.

For Coomassie staining, gels were washed thoroughly in water, incubated with ReadyBlue Protein Gel Stain (Sigma RSB-1L) overnight, and destained in water. For western blots, proteins were transferred to nitrocellulose membranes by semi-dry transfer (Biorad 1704158) according to the manufacturer’s instructions. Membranes were rinsed in water and stained with Revert 700 Total Protein Stain according to the manufacturer’s instructions. The membranes were then rinsed in TBS, rocked with Everyblot Blocking Buffer (Biorad 12010020) at room temperature for > 30 minutes, and rocked with primary antibody overnight at 4°. The membranes were washed with TBST, and rocked with IR800CW-labeled secondary antibodies for 30–60 minutes, washed with TBST, and imaged on a Licor Odyssey CLx.

CD11b cell surface staining

K562s expressing dCas9-BFP-KRAB were co-spinfected with lentiviruses expressing GFP-sgKDM1A and mCherry-sgPTPN1. Eight days after spinfection, the cells were counted and harvested by centrifugation. Cells were washed with PBE buffer (DPBS + 0.5% BSA + 2 mM EDTA) and resuspended with α-CD11b-AF647 antibody diluted 1:50 in PBE. Cells were incubated at 4° in the dark for 30 minutes, washed twice with PBE, and analyzed on a BD LSRFortessa. The populations were gated from a single sample as sgKDM1A (GFP+,mCherry-), sgPTPN1 (GFP-, mCherry+), and sgKDM1A/sgPTPN1 (GFP+, mCherry+). Unstained K562s expressing either GFP or mCherry were used as single color compensation controls. AF647 was compensated with UltraComp eBeads Plus (Thermo 01-3333-42) labeled with α-CD11b-AF647.

Internally controlled growth assays

K562s expressing Zim3-dCas9-2A-BFP were co-spinfected in triplicate with lentiviruses expressing GFP-sgKDM1A and mCherry-sgPTPN1, or with lentivirus expressing GFP and a non-targeting sgRNA. Every two days, cells were analyzed for BFP, GFP, and mCherry on an Attune flow cytometer. Enrichment was calculated as sgKDM1A (BFP+, GFP+, mCherry-), sgPTPN1 (BFP+, GFP-, mCherry+), and sgKDM1A/sgPTPN1 (BFP+, GFP+, mCherry+) vs uninfected (BFP+, GFP-, mCherry-).

Bulk RNA-seq

K562s expressing dCas9-BFP-KRAB were spinfected in biological duplicate with lentivirus expressing GFP and an sgRNA. Two days after spinfection, the cells were sorted for GFP+ on a BD ARIA II. Sort purity was generally >95%. After the sort, cells were maintained in media supplemented with 4 ug/ml puromycin for four days and then recovered for two days. Cells were counted, collected by centrifugation, and harvested by vigorous vortexing in Qiazol (Qiagen 79306).

Total RNA was extracted with miRNeasy Mini columns (Qiagen 217004) according to the manufacturer’s instructions and sequencing libraries were prepared with TruSeq Stranded Total RNA Library Prep Human/Mouse/Rat kits (Illumina 20020596) according to the manufacturer’s instructions. Libraries were sequenced 2x150 on a NovaSeq (Illumina).

Bulk RNA-seq data were aligned and quantified using STAR (version 2.7.9a) with parameters alignEndsType=Local and outFilter-MultimapNmax=20. For comparison of gene-level expression profiles, gene counts were corrected for sequencing depth (reads per million), and the log₂ fold-change for each gene was calculated relative to within-replicate non-targeting control expression. The two replicates for each genetic perturbation were averaged in order to produce the final data.

Seahorse experiment

K562s expressing dCas9-BFP-KRAB were spinfected with lentivirus expressing GFP and an sgRNA. Two days after spinfection, the cells were sorted for GFP+ on a BD ARIA II. Sort purity was generally >95%. After the sort, cells were maintained in media supplemented with 4 ug/ml puromycin for four days and then recovered for three days. On the 9th day post spinfection, seahorse assay were plates were treated with Cell-Tak (Corning 354240) according to the manufacturer’s instructions. Cells were counted, collected by centrifugation, and resuspended in supplemented Seahorse XF RPMI (Agilent 103576-100). 150,000 cells were added to the Seahorse assay plate and attached via centrifugation at 200xg for 1 minute with no brake. After 30 minutes of recovery at 37°, the cells were subjected to a Mito Stress Test on a Seahorse XFe96 analyzer according to the manufacturer’s instructions.

Alignment, cell calling, and guide assignment

Cell Ranger 4.0.0 software (10x Genomics) was used for alignment of scRNA-seq reads to the transcriptome, alignment of sgRNA reads to the library, collapsing reads to UMI counts, and cell calling. The 10x Genomics GRCh38 version 2020-A genome build was used as a reference transcriptome. For specific applications discussed below, STARsolo (STAR version 2.7.9a) was used extract transcript features, including intronic and exonic alignments and alignment of reads to transposable elements.

Reads from the sgRNA libraries were mapped with Cell Ranger. To account for differences in sequencing depths across GEM groups from the same experiment, reads were downsampled to produce a more even distribution of the number of reads per cell across gemgroups, with a threshold of 1000 reads per cell in the K562 day 8 experiment, 800 reads per cell in the K562 day 6 experiment, and 3000 reads per cell in the RPE1 experiment. Guide calling was performed with a Poisson-Gaussian mixture model as previously described. For each guide, the mixture model was fit 100 times, selecting the maximum likelihood model from among the fits. After guide calling, each cell was categorized according to its guide identities as representing a single genetic perturbation or a multiplet (which may arise from lentiviral recombination or multiple cell encapsulation during droplet generation). Only cells bearing two guides targeting the same gene or a single guide were used for downstream analysis.

Downstream analyses were performed in Python, using a combination of numpy, scipy, Pandas, scikit-learn, pomegranate, infercnvpy, pygenometracks, scanpy and seaborn libraries.

Filtering and internal normalization of gene expression measurements

Our internal normalization approach is similar to the one described in Adamson et al. (2016). First, we identified ‘‘core’’ control sgRNAs. That is, within each experiment there are tens to hundreds of possible negative control sgRNAs that were synthesized to have similar base compositions to targeting sgRNAs (Horlbeck et al., 2016). Some of these by chance induce detectable phenotypes. We constructed a minimal set of control sgRNAs that are largely indistinguishable from each other using the following procedure: (i) We take all cells bearing all possible non-targeting sgRNAs and represent them by the vector of genes with mean >1 UMI count per cell. (ii) We z-normalize the expression of these genes: i.e. we subtract the mean and divide by the standard deviation. (iii) For each gene, we test for equality of distribution using the Anderson-Darling test (scipy.stats.anderson_ksamp) between all possible pairs of non-targeting control sgRNAs. (E.g. In the genome-scale dataset, there are 585 possible control sgRNAs and therefore $\left(\begin{array}{l}585\hfill \\ 2\hfill \end{array}\right)=170280$ pairwise comparisons.) (iv) We adjust the resulting p-values for multiple hypothesis testing using the Benjamini-Hochberg procedure. (v) For each potential control sgRNA, we compute the average number of differentially expressed genes relative to all other potential control sgRNAs. (vi) We set a dataset-dependent threshold on the number of differentially expressed genes (8 in the genome-scale dataset and 30 in the ‘‘K562 essentials’’ and ‘‘RPE1 essentials’’ datasets, which were more deeply sequenced and so had more genes passing the expression threshold) and kept all potential control sgRNAs that fell below the threshold. For example, in the genome-scale dataset this resulted in 514 control sgRNAs.

Next, we filtered cells based on quality metrics. We first computed scale factors to adjust for variable sequencing depths across gemgroups: we examined all core control cells (which make up ~4% of all cells), computed factors that equalized the mean UMI counts within these cells across gemgroups, and then applied these factors to all cells in the gemgroup to produce adjusted UMI counts. We then applied two quality filters, ensuring that cells passed a minimum adjusted UMI content filter (genome-scale dataset: 2000 UMIs, K562 essentials/RPE1 essentials: 3000 UMIs) and a maximum mitochondrial RNA filter (genome-scale dataset: <25%, K562 essentials: <20%, RPE1 essentials: <11%). (Mitochondrial RNA content is the fraction of total UMIs derived from mitochondrially-encoded genes.) These filter parameters were chosen by plotting adjusted UMI content vs. mitochondrial RNA content and manually setting thresholds that removed the low-quality cells.

Finally, we computed a normalized gene expression matrix for cells passing the quality filters via two normalization steps: (i) UMI count normalization: We scale expression within all cells so that their total UMI counts equal the median UMI count of core control cells within the experiment. (ii) Relative z-normalization: Within each gemgroup, for each gene, we compute the mean and standard deviation of expression within control cells and use these to z-normalize expression. In other words, if x is the expression of a given gene, it is represented by the score $z=\left(x-{\mu }_{control}\right)/{\sigma }_{control}$, where the mean and standard deviation are separately computed within each gemgroup.) The resulting scores should therefore be interpreted as ‘‘fraction of transcriptional effort’’ due to the UMI count normalization, with a scale set relative to control cells. Put simply, an expression score of +2 thus represents a gene expressed at a level 2 standard deviations above the mean in control cells.

Examining effects of normalization on batch effects

As described in the main text, we observed batch effects in the data (Data S2 [Figure i]). This variation appeared to track mostly with sets of 8 samples that went through the 10x Chromium instrument and library prep together, though the precise origin is unclear. To construct this figure, we normalized the data in two ways. Raw data normalization: (i) To adjust for variable sequencing depth, scale cellular UMI counts by factors chosen so that so that core control cells have the same total UMI counts across all gemgroups. (ii) Construct a gemgroup mean expression profile of all genes with mean >2 UMI counts per cell by averaging counts following normalization in previous step across all cells in the gemgroup. (iii) Scale the gemgroup mean expression profiles by dividing by their mean across all gemgroups. An expression value of 1 is then the mean across all cells across all gemgroups. Internal z-normalization. Normalize expression as described in previous section.

Data S2 (Figure i) compares the two normalization schemes. In both cases the ranges of the plot are chosen using seaborn’s robust option (which sets the min and max to the 2^nd and 98^th percentile of the data). Genes are clustered based on the raw data normalization and are in the same order in both panels. The gemgroups are presented in order based on how samples were multiplexed while performing the experiment as indicated by the color groupings at the top.

Energy distance test for identifying perturbations that induce altered transcriptional states

To compare distributions of expression states, we used tests derived from energy statistics, which allow for testing of equality of distributions when data are high-dimensional. In short, each cell is represented by a vector composed of its top 20 principal component scores, and we compare whether the distribution of these 20-dimensional vectors is equal or not between unperturbed control cells and cells bearing each perturbation. When these distributions differ, we can infer that the perturbation is causing some change either in the structure or distribution of transcriptional states within the perturbed cells.

To construct the distributions to compare, we first applied a series of filtering steps: (i) we removed cells that did not pass the UMI or mitochondrial RNA filters described in Internal normalization of gene expression measurements; (ii) as features, we took the z-normalized expression of all genes with mean expression >0.5 UMIs per cell; (iii) to dampen the effects of a handful of strongly induced outlier genes, we clipped any measure with a z-score greater than 10 to 10 (this only affects a handful of genes); (iv) finally, we applied principal components analysis (using sklearn’s PCA implementation, which will use randomized algorithms for datasets of this scale) and kept only the top 20 principal components. The test should therefore be interpreted as assessing gross changes in cellular transcriptional state.

To construct a null distribution, we randomly subsampled 5,000 control cells bearing non-targeting sgRNAs. (Subsampling was necessary for performance reasons.) For each perturbation, we then compute an estimator of the energy distance:

where each x_i is one of the control cells and each y_j is one of the perturbed cells.

In the limit of infinite data, the energy distance will be 0 between identical distributions and positive between non-identical distributions. We assess statistical significance in practice using a permutation test by permuting the labels of control and perturbed cells 10,000 times and estimating how frequently a larger energy distance would be observed by chance. The specific implementation is based on the python package torch-two-sample, modified to use numba for improved performance.

Gene-level differential expression testing using the Anderson-Darling and Mann-Whitney tests

Because of (i) biological differences in expression characteristics across different genes, (ii) the batch effects described above, (iii) incomplete penetrance of some perturbations, and (iv) heterogeneity of some gene expression programs, we opted to use non-parametric statistical tests rather than tests based on specific distributional assumptions about gene expression. Specifically, we z-normalize gene expression relative to control cells as described (Internal normalization of gene expression measurements) and for each gene test whether the distribution of normalized expression is identical between control cells bearing non-targeting sgRNAs and cells bearing each perturbation. We used two tests implemented in scipy: the Anderson-Darling test (scipy.stats.anderson_ ksamp), which is broadly sensitive to changes in distribution, and the Mann-Whitney U test (scipy.stats.mannwhitneyu), which tests whether one distribution is stochastically greater than another. For the Anderson-Darling test we extended the range of p values beyond those available in scipy’s implementation by computing the p-value for many values of the test statistic using R’s kSamples package and interpolating any intermediate values using scipy.interpolate.interp1d. For the Mann-Whitney test we used the asymptotic p values and excluded any perturbation with fewer than 10 cells. p-values in both cases were adjusted for multiple hypothesis testing using the Benjamini-Hochberg procedure to produce the final results.

Functional analyses of strong perturbations

We conducted functional analyses to gain insight into the types of gene perturbations that induced strong transcriptional phenotypes. ‘‘Strong’’ perturbations were defined by three criteria: (i) at least 50 differentially expressed genes at a significance of p < 0.05 by Anderson-Darling test following Benjamini-Hochberg correction; (ii) at least 25 cells that passed our quality filters; and (3) an on-target knockdown, if measured, of at least 30% (i.e. the target of perturbation was either knocked down by at least 30% or was not detected, a broad attempt to remove non-functional perturbations). ‘‘Weak’’ perturbations met the same criteria but had fewer than 5 differentially expressed genes. Strong and weak perturbations were largely similar in terms of representation and knockdown efficacy, as described in Figure S2. To look for classes of functional behaviors among strong and weak perturbations, we used the gseapy implementation of the Enrichr algorithm to compute gene set enrichment p-values within the KEGG2021 pathway gene set (with the set of all targeted genes in the experiment as the background list). To determine whether strong and weak perturbations fell in different subcellular locations, we used location annotations from Itzhak et al. (2016).

Global analysis and clustering of strong perturbations

The analysis presented covers 1973 perturbations that met three criteria: (i) at least 50 differentially expressed genes at a significance of p < 0.05 by Anderson-Darling test following Benjamini-Hochberg correction; (ii) at least 25 cells that passed our quality filters; and (3) an on-target knockdown, if measured, of at least 30% (i.e. the target of perturbation was either knocked down by at least 30% or was not detected, a broad attempt to remove non-functional perturbations). As features, we used a union of two sets of genes: (i) the top 10 differentially expressed genes for all perturbations (ordered by the value of the Anderson-Darling test statistic) and (ii) all genes of mean >0.25 UMIs per cell with variance in the top 30% of the dataset. We represented perturbations by their mean normalized expression profile across these 2319 highly variable genes. To prevent the direct targets of knockdown influencing results, the target gene value was replaced by 0 for the corresponding perturbation. For example, RPS5 gene normalized expression was set to 0 in the expression profile of the RPS5 perturbation, which is equal to the mean in control cells by construction.

Because clearly related perturbations sometimes showed variable absolute phenotypic strengths, we used correlation as a metric to compare profiles, since it is scale-invariant. We conducted two global assessments of the ability of these expression profile correlations to recall known biological relationships. First, curated complexes were obtained from the 03.09.2018 CORUM3.0 database (Giurgiu et al., 2019). We identified all complexes that had at least 66% of genes represented within the 1973 perturbations (based on matching gene symbols between the datasets), leading to 327 complexes. Each represented complex was then split into a series of links (e.g. if a complex contained genes A, B, and C, then it would be split into links A-B, B-C, and A-C). The figure plots the distribution of expression profile correlations of these annotated links versus the distribution of all possible links among the 1973 targeted genes. A similar analysis was then conducted using predicted protein links from the v11.5 of STRING (Szklarczyk et al., 2019) (9606.protein.links.v11.5.txt.gz) after mapping STRING protein IDs to gene names (using the ‘‘preferred_name’’ field in 9606.protein.info.v11.5.txt.gz). Among the 1973 genes in the figure there are 1,945,378 possible pairwise links between genes, 243,558 of which have scores within STRING. We binned these represented links into 6 equally spaced bins based on observed expression profile correlation. The figure shows kernel density estimates of the STRING score distribution within each bin made using seaborn’s violinplot (with cut set to 0 so that density estimates do not extend past observed data).

To identify clusters of related perturbations, we manually computed correlation distances between all pairs of expression profiles, and used HDBSCAN (metric=’precomputed’, min_cluster_size=4, min_samples=1, cluster_selection_method=’eom’) to identify 63 clusters. This procedure is intrinsically conservative due to the choice of metric and clustering algorithm, so many perturbations are not assigned to any cluster—our emphasis was on identifying the strongest signals rather than the most comprehensive. We then annotated the possible function of these clusters using a combination of manual lookup of related genes and automated annotation using CORUM complexes and STRING clusters. (STRING clusters are derived from 9606.clusters.info.v11.5.txt.gz and are labeled according to the ‘‘best_described_by’’ field.) We only assigned automated annotations when a cluster contained 75% or more of the members of a CORUM complex or STRING cluster. Aggregated information about clusters is provided in Table S3 which has the following fields:

Minimum distortion embedding of strong perturbations

The visualization in Figure 2D is a minimum distortion embedding (MDE) of the 1973 strong perturbations created using pymde v0.1.13. pyMDE solves MDE problems based on minimizing Euclidean distances. To adapt it to correlation distances, we first z-normalized each expression profile. (I.e. If a perturbation is represented by the vector x of 2319 highly variable genes, we computed $\widehat{\mathit{x}}=\left(\mathit{x}-\mathit{x}\right)/{\sigma }_x$.) Because of the polarization identity, minimizing the Euclidean distance between these normalized profiles is equivalent to minimizing the (square root of) the correlation distance between the unnormalized profiles:

where the second equality follows from the z-normalization and the scale-invariance of correlation.

We created two embeddings. First, we used pymde to embed the dataset into 20 dimensions. This ‘‘high-dimensional embedding’’ serves as an imputation step, as it distorts the geometry of the dataset so that clusters of related genes that may be driven by weaker overall correlations are allowed to form. Proximity within this embedding was used to identify genes, CORUM complexes, and STRING clusters that were near to the HDBSCAN clusters called on the raw data (which were well-preserved in the embedding) as described in Table S3. To construct the embedding, we initialized pymde using the spectral embedding of the dataset (using sklearn’s SpectralEmbedding with n_components=20, affinity=’nearest_neighbors’, n_neighbors=7, eigen_solver=’arpack’)), and then ran pymde’s ‘‘preserve_neighbors’’ function with embedding_dim=20, n_neighbors=7, and repulsive_fraction=5. pymde was run until convergence with a final average distortion of 0.0979 and final residual norm of 9.4e-06.

To produce the embedding in Figure 2D, we ran pymde with the same parameters but with the embedding dimension set to 2 (final average distortion 0.105, final residual norm 3.2e-06). The bold cluster labels in the figure correspond to the manual annotations mentioned in the previous section. A handful of changes were manually incorporated: (1) the cytochrome c-ubiquinol cluster was not detected by HDBSCAN, and was manually annotated (2) 4 clusters involving protein post-translational modifications (ubiquitination, sumoylation, acetylation, neddylation) were annotated with a single label of ‘‘post-translational modifications’’ (3) components of eIF3 split across two clusters that are next to each other in the embedding and are labeled as a single cluster (4) all clusters of unknown function were not labeled but are included in Table S3. The complex labels come from CORUM or STRING. A complex/cluster label was placed if and only if 75% of the members of a complex or cluster were close to each other in the 20-dimensional pyMDE embedding (‘‘close’’ meaning at or below the 5^th percentile of all pairwise distances). Redundant/duplicated clusters were manually deduplicated. The locations of the labels on the figure were then adjusted for readability. Label location is a decent proxy for, but not an entirely accurate representation of, cluster and complex locations.

Clustering of gene expression programs

We next turned to identifying conserved gene expression programs using a similar pipeline applied to the transpose of the expression matrix from the previous sections. Initial HDBSCAN clustering on the raw data did not yield very many clusters, which we attributed to the broad range in gene expression program size and dynamic range. To attempt to equalize for these factors, we performed the clustering on a high-dimensional embedding of the data. Each gene was represented by its expression across the 1973 perturbations in Figure 2 and we masked the targets of knockdown as described there to avoid target gene knockdown influencing clustering. We used pyMDE with the same normalization as above to encourage genes with correlated expression to be placed nearby to each other (20 dimensions, n_neighbors=7 and repulsive_fraction=5, final average distortion 0.145, final residual norm 5.1e-06). We then identified clusters using HDBSCAN applied to the embedding (metric=’euclidean’, min_cluster_size=10, min_samples=10, cluster_selection_method=’leaf’), producing 38 clusters. We performed similar analyses as in Figure 2 to annotate known CORUM complexes and STRING clusters. Cluster identities were then manually annotated using a combination of these automated annotations, manual gene searches, and gene set enrichment analyses conducted using Enrichr (Xie et al., 2021). These clusters are summarized in Table S3, which includes the fields:

To produce Figure 3B, we averaged expression within the 64 perturbation clusters from Figure 2 across the genes within the 38 gene expression clusters. Each element in the heat map therefore represents an average over both multiple related perturbations and multiple related genes. The labels were manually selected to highlight interesting features.

Screens of gene expression programs

To demonstrate the ability of Perturb-seq to conduct screens on aggregate phenotypes, we conducted two analyses to identify perturbations that strongly induced interesting expression programs. In the first comparison, we compared expression of genes associated with erythroid differentiation (gene expression cluster 15) to those associated with myeloid differentiation (gene expression cluster 21) in Figure 4D. We scored expression programs by taking the mean normalized gene expression of all genes in the associated clusters. We computed scores for all perturbations in the genome-scale dataset that were detected in at least 25 cells, and then z-normalized these scores to make scales comparable. The figure has labels on the 15 most outlying genes across the two programs. We then conducted an identical analysis comparing expression of an unfolded protein response cluster (cluster 2) to an integrated stress response cluster (cluster 12) in Figure 4C.

Analysis of composite phenotypes: total RNA, fraction mtRNA, fraction TE, RNA splicing

Composite phenotypes integrate data from across the transcriptome to describe global cellular features. While some derive from simple metrics, others rely on extracting information from the transcriptome beyond gene expression levels. As an example of a simple metric, the number of UMIs aligned to the transcriptome (GRCh38 version 2020-A) of single-cells was used to represent the total cellular RNA content. To produce Figure 4H, we averaged the total RNA content of cells for each perturbation, and z-scored the RNA content with respect to non-targeting controls. Similarly, to calculate the fraction of mitochondrial RNA per cell (fraction mtRNA), the sum of the expression levels of the 13 mitochondrial genome protein-coding genes (MT-ND6, MT-ND1, MT-ND2, MT-ATP8, MT-ND4L, MT-ND5, MT-ND3, MT-CO1, MT-CO2, MT-ND4, MT-ATP6, MT-CO3, MT-CYB) was divided by the total cellular RNA content for each individual cell.

The scTE (He et al., 2021) processing pipeline was used to quantify the expression of transposable elements in single cells. As transposable elements tend to be present in many degenerate copies throughout the genome, scTE allocates TE reads to TE metagenes rather than specific genomic positions. Reads were aligned to the genome using STARsolo (STAR version 2.7.9a) with the flags ‘ –outSAMattributes NH HI AS nM CR CY UR UY –soloFeatures Gene GeneFull SJ Velocyto –readFiles Command zcat –outFilterMultimapNmax 100 –winAnchorMultimapNmax 100 –outMultimapperOrder Random –runRNGseed 777 –outSAMmultNmax 1’ to allow multimapping. To avoid incompatibilities in cell calling between STARsolo and Cell Ranger, the output of Cell Ranger cell calling was used to define the STARsolo cell barcode whitelist using ‘–soloCBwhitelist’. Next, aligned reads were allocated to genes and TEs using scTE. The flag ‘-o nointron’ was used to prevent the quantification of TEs in gene introns. From single-cell transcriptomes, quantification of all TEs in the classes LINE, SINE, LTR, DNA, and Retroposon based on RepeatMasker were extracted. To calculate the fraction of repetitive and transposable element RNA per cell (fraction TE), the sum of the expression level of these TEs was divided by the total cellular RNA content for each individual cell. To produce Figure 4G, we averaged the fraction TE of cells for each perturbation.

The alignments from STARsolo described above were also be used to quantify RNA splicing. Due to (i) the sparsity of single-cell data (ii) the relationship between the fraction of spliced reads and gene expression levels, gene-wise RNA splicing was quantified at the pseudobulk level. From the STARsolo Velocyto output, the levels of spliced and unspliced reads for each gene were averaged across all cells bearing each perturbation, ignoring ambiguous reads. Then, for each gene, the fraction of unspliced reads was divided by the mean fraction of unspliced reads for that gene across all non-targeting control perturbations. The results in Figure 3F display a common set of genes across all perturbations.

Leverage scores for quantifying perturbation penetrance and variability

Our use of non-parametric tests in differential expression testing is in part to accommodate perturbations that may be incompletely penetrant or heterogeneous in effect. To attempt to quantify these features, we developed a scalar single-cell score to summarize how outlying each cell’s transcriptional state was relative to control cells. We used an approach based on leverage scores, which measure how outlying the rows or columns of a matrix are and which form the basis for many randomized algorithms (Ma et al., 2013). Specifically, we: (i) Construct an expression matrix consisting of all cells that pass the quality filters described in Internal normalization of gene expression measurements, and all genes with mean expression >0.25 UMI counts per cell. (ii) To dampen the effects of a handful of strongly induced outlier genes, we clipped any measure with a z-score greater than 10 to 10 (this only affects a handful of genes) (iii) To avoid the influence of gemgroup-level batch effects and variable sequencing depth, we then compute leverage scores separately within each gemgroup. Row leverage scores, corresponding to the cell axis of the expression matrix, are calculated as the squared norm of the top 20 left singular vectors within each gemgroup (computed via the truncated SVD routine scipy.sparse.svds with the arpack solver with k=20). We then normalize these scores so that the sum over all cells in the gemgroup is 1 (i.e. compute the leverage sampling probability distribution). (iv) Finally, to integrate leverage scores across gemgroups, we then take logs, and rescale by z-normalizing relative to the scores of control cells (subtracting their mean and dividing by their standard deviation). All the leverage scores presented in the figures are the leverage scores after this normalization procedure.

In Figure S6 we conducted various analyses to validate leverage scores as measures of phenotype and to use them to study penetrance of perturbations. We considered all perturbations that passed the following criteria: (i) >5 differentially expressed genes by Anderson-Darling test; (ii) detected in at least 25 cells that passed our quality filters; and (iii) the gene targeted by the perturbation was either undetectable in the expression data or knocked down by at least 30% if detected (i.e., we removed perturbations that appear non-functional). For these perturbations we compared the mean leverage scores to the number of differentially expressed genes found using the Anderson-Darling test. To assess reproducibility we then subset to perturbations that were present in both the genome-scale dataset and the K562 essentials dataset, including non-targeting controls.

In the analyses of perturbations targeting Mediator and the small subunit of the ribosome we only included perturbations that were (i) present in both the genome-scale dataset and the ‘‘K562 essentials’’ dataset and (ii) targeted the principal ‘‘P1’’ transcript identified by the FANTOM consortium. (A handful of genes also had perturbations targeting the P2 transcript that did not generally have effects.) Knockdown was computed as the ratio of mean (unnormalized) expression of the target gene within perturbed cells vs. that in cells with non-targeting sgRNAs. The plots show kernel density estimates of the distributions of the leverage scores of all cells with these perturbations constructed using seaborn’s violinplot (with cut set to 0 so that estimated distributions do not extend beyond the range of the data). The gray bars represent the 10%-90% quantiles of cells with non-targeting control sgRNAs for comparison.

Finally, in Figure 5B we used leverage scores to search for perturbations that had highly variable phenotypes. We considered all perturbations using the same criteria as in fig. S6. We used the standard deviation of the leverage scores as a metric for variability, as diagrammed in Figure 5A. The two examples in this figure are derived from actual data. The 20 labeled genes are the most outlying from the lowess local regression between the standard deviation of leverage scores and the log of the number of differentially expressed genes detected by the Anderson-Darling test (computed using statsmodels.nonparametric.smoothers_lowess.lowess).

Analysis of chromosomal instability and cell cycle

The framework described in inferCNV (Patel et al., 2014) and implemented as infercnvpy (https://github.com/icbi-lab/infercnvpy) was used to detect evidence of chromosomal copy number changes. The raw single-cell gene expression matrix was first filtered to remove lowly expressed genes (<0.05 UMIs per cell across the population), normalized for total UMI content (using scanpy.pp.normalize_total with target_sum=1e6, exclude_highly_expressed=False, and max_fraction=0.05), and log-scaled (using scanpy.pp.log1p). Then we used infercnvpy to compute rolling average gene expression changes for windows of 100 genes with a dynamic threshold of 1.5 standard deviations for noise filtering (using infercnvpy.tl.infercnv). For each single-cell, this generated a vector of CIN values across the genome. To label cells with likely karyotypic abnormalities, unstable karyotypic cells were heuristically defined as having ≥1 chromosome with evidence of changes in chromosomal copy number (nonzero CIN values) for >80% of the chromosomal length. For Figure 5G, the CIN score of genetic perturbations was calculated as the mean single-cell sum of squared CIN values, z-normalized relative to non-targeting control perturbations.

To show cell cycle effects in Figures 5D and 5E, we chose to use a dimension reduction approach. Previous approaches to cell cycle analyses in Perturb-seq have largely focused on supervised classification of cells into canonical cell cycle states. However, we found that these approaches did not allow for aberrant cell cycle states sometimes generated by genetic perturbations. As a summary of single-cell cell cycle states, we performed a Uniform Manifold Approximation and Projection (UMAP) dimension reduction based on the expression n=199 known cell cycle genes [obtained from Seurat (Satija et al., 2015) and (Adamson et al., 2016)]. From total UMI content normalized, log-scaled expression data, a neighborhood graph was computed (using scanpy.pp.neighbors with n_neighbors=30, method=’umap’, metric=’correlation’, and n_pcs=20) followed by UMAP embedding (using scanpy.tl.umap with default parameters). This UMAP revealed cells in canonical cell cycle stages when compared with other methods, but also naturally separated dying cells and putatively quiescent cells. Gates were drawn by manual inspection to approximately separate cells likely to be in S, G2/M, and G1/G0 cell cycle phases.

Analysis of transcriptional responses to mitochondrial stress

To compare the functional specificity of nuclear and mitochondrial genome responses to mitochondrial stress, we clustered mitochondrial perturbations based on their transcriptional phenotypes. The analysis presented in Figures 6A and 6C (K562 day 8 data) covers 268 mitochondrial perturbations that met three criteria: (i) at least 50 differentially expressed genes at a significance of p < 0.05 by Anderson-Darling test following Benjamini-Hochberg correction; (ii) at least 30 cells; and (iii) an on-target knockdown of at least 60%. As features, we used 1715 genes encoded in the nuclear genome expressed at >1 UMI per cell or the 13 protein-coding mitochondrial-encoded genes. Analogously, the analysis presented in Figures S12A and S12E (RPE1 day 7 data) covers 140 mitochondrial perturbations that met three criteria: (i) at least 20 differentially expressed genes at a significance of p < 0.05 by Anderson-Darling test following Benjamini-Hochberg correction; (ii) at least 30; and (iii) an on-target knockdown of at least 60%. As features, we used 2017 genes encoded in the nuclear genome expressed at >1 UMI per cell or the 13 protein-coding mitochondrial-encoded genes.

To cluster these data, we used correlation as a metric to compare z-normalized expression profiles, since it is scale-invariant. We clustered using HDBSCAN (metric=’correlation’, min_cluster_size=3, min_samples=1, cluster_selection_method=’eom’, alpha=1.0). As discussed above, this procedure is intrinsically conservative due to the choice of metric and clustering algorithm, so many perturbations are not assigned to any cluster. We used the hierarchical clustering output of HDBSCAN to manually identify and label groups of perturbations in the figures.

To compare the heterogeneity of mitochondrial genome expression, we used two approaches. First, we used a data-driven approach comparing the 38 gene expression programs described in Figure 4B across perturbations. For each program, we scored its expression within each genetic perturbation, and calculated a standard deviation of these scores across the different perturbations in the K562 day 8 experiment. The mitochondrial genome program discovered in this way consists of the 13 protein-coding mitochondrial genes, plus two MT-RNR2-like pseudogenes encoded in the nuclear genome which are likely to multimap with mitochondrial-encoded transcripts. In Figure S12C, we show a histogram of the variability of expression programs across perturbations. Second, we used a hypothesis-driven approach that compared the variation of mitochondrial genome responses by the protein localization of perturbations. We used data from the Human Protein Atlas to assign locations to different perturbations (https://www.proteinatlas.org/about/download table subcellular_location.tsv.zip; some locations were collapsed into supersets), excluding any dual-localized proteins. For all perturbations with the same localization and at least 50 differentially expressed genes at a significance of p < 0.05 by Anderson-Darling test following Benjamini-Hochberg correction, we calculated the variance of z-normalized gene expression profiles for each mitochondrially encoded gene. The data in Figure 7B (K562 day 8) and S12B (RPE1) represents the average variance across mitochondrially-encoded genes for each localization.

To quantitatively compare the specificity of the mitochondrial and nuclear transcriptional responses, we employed random forest classifiers. The analysis presented in Figure S12F (K562 day 8 data) covers perturbations that met three criteria: (i) cause at least 50 differentially expressed genes at a significance of p < 0.05 by Anderson-Darling test following Benjamini-Hochberg correction; (ii) at least 30 cells; and (iii) belong to either the small or large mitochondrial ribosomal subunit (identified by gene names beginning with MRP), ATP synthase (identified by gene names beginning with ATP5), proteostatic factors (including PAM16, DNAJC19, TFAM, HSPA9, TOMM40, TOMM20, PMPCB, LONP1, and HSPE1). As features, we used 2039 genes encoded in the nuclear genome expressed at >0.5 UMI per cell with variance in the upper 25% of the dataset, or the 13 protein-coding mitochondrial encoded genes. To train the random forest classifier to separate cells by perturbed complex (as described above), we used the scikit-learn implementation with extremely randomized trees with 1,000 trees in the forest and 100 features. The reported accuracy is the balanced accuracy defined as the defined as the average of recall obtained on each class. To visualize the data, we used a UMAP (metric=’correlation’, n_neighbors=10) generated on the nuclear or mitochondrial genome.