Titrating gene expression using libraries of systematically attenuated CRISPR guide RNAs

(a) Representative plots illustrating gating strategy to select cells for analysis. (b) Comparison of relative activities obtained from two replicate transductions. Relative activity was defined as the fold-knockdown of each mismatched variant (GFP_{sgRNA[non-targeting]}/GFP_{sgRNA[variant]}) divided by the fold-knockdown of the perfectly-matched sgRNA. The background fluorescence of a GFP^– strain was subtracted from all GFP values prior to other calculations. n = 57 sgRNAs; r² = squared Pearson correlation coefficient. (c) KDE plots of GFP distributions 10 days after transducing K562 GFP⁺ cells with the perfectly-matched sgRNA, a non-targeting sgRNA, and each of the 57 singly-mismatched variants. Fluorescence of GFP^– K562 cells is shown in gray. Although most GFP distributions are unimodal, some are broadened compared to those with the perfectly matched sgRNA or the negative control sgRNA. This heterogeneity could be a consequence of the random integration of the GFP locus, cell-to-cell differences in expression of the dCas9-KRAB effector in our polyclonal cell line, the amplification of gene expression bursts by long GFP half-lives, or a combination of these factors. Two replicate transductions were evaluated for each sgRNA (see panel b); data from one replicate are shown here.

(a, b) Comparison of growth phenotypes (γ) of all sgRNAs derived from replicates of the (a) K562 (n = 119,201 sgRNAs) and (b) Jurkat screens (n = 119,229 sgRNAs). Marginal distributions are normalized to element count in each category. r² = squared Pearson correlation coefficient for targeting sgRNAs (mismatched and original). (c) Comparison of γ of perfectly matched sgRNAs from the K562 screen in this work and a previously published K562 screen¹⁴ (average of two replicate screens). n = 4,830 sgRNAs; r² = squared Pearson correlation coefficient. (d) Comparison of γ of perfectly matched sgRNAs in K562 and Jurkat cells reveals substantial differences, likely reflecting cell-type specific gene essentiality (average of two replicate screens). n = 4,892 sgRNAs; r² = squared Pearson correlation coefficient. (e) Comparison of mismatched sgRNA relative activities in K562 and Jurkat cells, classified by the difference in γ of the corresponding original guide. n = 15,103 (left) and 26,409 (right) sgRNAs; r² = squared Pearson correlation coefficient. (f) Distribution of mismatched sgRNA relative activities for sgRNAs with 1 mismatch (left) or 2 mismatches (right). (g) Distribution of mismatched sgRNA relative activities stratified by sgRNA GC content, grouped by mismatches located in positions –19 to –13 (PAM-distal region), positions –12 to –9 (intermediate region), and positions –8 to –1 (PAM-proximal/seed region). n = 282-7,592 sgRNAs. (h) Distribution of mismatched sgRNA relative activities stratified by the identity of the 2 bases flanking the mismatch, grouped by mismatches located in the three regions as in g. n = 155-2,031 sgRNAs. (i) Distribution of mismatched sgRNA relative activities stratified based on whether or not the invariant first G of the sgRNA (position –20) matches the genome, grouped by mismatches located in the three regions as in g. n = 4,267-11,524 sgRNAs. (j) Comparison of mean CRISPRi relative activities from large-scale screen and cutting frequency determination (CFD) scores²⁷. Values are compared for identical combinations of mismatch type and mismatch position; mean relative activities were calculated by averaging relative activities for all mismatched sgRNAs with a given combination. n = 228 mismatch type/position combinations; r² = squared Pearson correlation coefficient. (k) Distribution of sgRNA series by number of sgRNAs with intermediate activity (0.1 < relative activity < 0.9), using only sgRNAs with a single mismatch (top) or all mismatched sgRNAs (bottom). Lines in violin plots in panels g, h, i denote distribution quartiles.

(a) Comparison of growth phenotypes measured in replicate screens after 4, 6, or 8 days of growth from t₀. Data from Day 4 were used for all subsequent analyses. n = 35,830 sgRNAs; r² = squared Pearson correlation coefficient. (b) Comparison of relative % knockdown (quantified via RT-qPCR) and mean relative growth phenotype for 10 intermediate-activity constant region variants paired with two targeting sequences against DPH2. Data represent the mean of technical triplicates. (c) Relative activities of constant regions paired with all 30 targeting sequences, ranked by the average strength of each constant region and displayed as rolling means with a window size of 50. (d) Distribution of all pairwise correlations of constant region relative activities within and between gene targets. n = 30 and 1,350 for intra-gene and inter-gene comparisons, respectively; indicated p-values are derived from a two-tailed Student’s t-test; dashed lines in violin plots indicate the distribution quartiles. (e) Relative activity of each indicated target sequence:constant region pair vs. the mean relative activity of the respective constant region for all targets. Growth phenotypes (γ) with the unmodified constant region are indicated in the figure legends. Lines represent rolling means of individual data points.

(a) Graph of the CNN model architecture. (b) Example of 5-fold cross-validation using only the training dataset, further analyzed in the subsequent two panels of this figure. A similar scheme was used to optimize hyperparameters for the CNN model, albeit with 3-fold cross-validation to allow for larger training sets in each split. (c) Model loss, measured as root mean squared error, for training and test data over 30 training epochs. Each line represents one of 5 splits diagrammed in panel b. The final models used for our predictions were trained for 8 epochs, as additional cycles only reduced training loss without significant improvement in validation loss (i.e., the model becomes overfit). (d) Stability of the model with different input data. For each split in panel b, 20 independent CNN models were trained for 8 epochs on the same data. The root mean squared error on the test set for each model is plotted as a blue dot. Box plots indicate the interquartile range of each distribution. (e) Model loss for the final CNN ensemble. Each line represents one of 20 models trained for 8 epochs on the entire training set. (f) Explained variance of validation sgRNA relative activities for each individual model (black), and for the mean prediction of all 20 models (red). n = 5,241 sgRNAs evaluated for each model; r² = squared Pearson correlation coefficient. (g) Validation error stratified by mismatch position. (h) Validation error stratified by mismatch type. (i) Comparison of CNN prediction error (difference between measured and predicted activity) and off-target specificity score for all sgRNAs in the validation set. Off-target specificity scores were calculated using CRISPRi relative activities as described in the Methods. n = 5,241 sgRNAs; r = Pearson correlation coefficient. (j) Partitioning of sgRNAs into bins based on relative activity in the large-scale K562 screen. (k) Confusion matrix showing the fraction of sgRNAs in each actual (measured) activity bin that were assigned to each predicted bin by the CNN model. Each row sums to 1. (l) Statistics indicating the requisite number of randomly sampled sgRNAs from each activity bin to have a given probability of selecting at least one sgRNA with true activity in that bin. Simulations are based on the probabilities outlined in the confusion matrix (panel e). (m) Similar to panel l, with random sampling from bin 2 (relative activity 0.37-0.63) to yield at least one sgRNA with intermediate activity (0.1-0.9). We tested several sampling schemes (e.g. drawing from bin 1, 2, 3, or combinations of these), and found this method to empirically give the highest success rate for selecting sgRNAs with intermediate activities.

(a) Comparison of measured relative growth phenotypes from the large-scale screen and predicted activities assigned by the elastic net linear model. Marginal histograms show distributions of relative activities along the corresponding axes. n = 5,241 sgRNAs; r² = squared Pearson correlation coefficient. (b) Comparison of measured relative activity (relative knockdown) in the GFP experiment and predicted relative sgRNA activity. n = 57 sgRNAs; r² = squared Pearson correlation coefficient. (c) Comparison of predicted relative activities from the linear model and the neural network, based on the validation set of singly-mismatched sgRNAs. n = 5,241 sgRNAs; r² = squared Pearson correlation coefficient. (d) Regression coefficients assigned to each feature in the linear model. 228 features (gray, blue) describe the position and type of mismatch; 42 features (gold) carry other information about the sgRNA and genomic context surrounding the protospacer. These features are detailed in subsequent panels. (e) Linear coefficients for features of the sgRNA and targeted locus. TSS; transcription start site. (f) Linear coefficients for features covering positions in the distal, intermediate, and seed regions of the targeting sequence (highlighted blue in panel d).

(a) Composition of the compact library, in terms of previously measured relative activities in the large-scale screen (dark purple), or predicted relative activities assigned by the CNN model ensemble (light purple). Perfectly matched sgRNAs, which by definition have relative activities of 1.0, comprise 20% of the library but were not included in the histogram. (b) Distribution of mismatch positions and types for singly-mismatched sgRNAs in the compact library, for previously measured (dark purple) and CNN-imputed (light purple) sgRNAs. (c) Heatmap showing the distribution of mutated positions for doubly-mismatched sgRNAs in the compact library. (d) Comparison of growth phenotypes measured in each K562 replicate screen 4- and 7-days post-transduction. Data from Day 7 were used for all subsequent analyses. n = 25,518 sgRNAs; r² = squared Pearson correlation coefficient. (e) Comparison of growth phenotypes measured in each HeLa replicate screen 6- and 8-days post-transduction. Data from Day 8 was used for all subsequent analyses. n = 25,518 sgRNAs; r² = squared Pearson correlation coefficient. (f) Comparison of growth phenotypes of original (perfectly matched) sgRNAs in HeLa and K562 cells (γ, expressed as the average of two replicate screens). n = 4,810 sgRNAs; r² = squared Pearson correlation coefficient. (g) Measured vs. predicted relative activities of CNN-imputed sgRNAs in K562 cells (left) and HeLa cells (right). A small number of points beyond the y-axis limits were excluded to more clearly display the bulk of the distribution. n = 6,147 sgRNAs; r² = squared Pearson correlation coefficient. (h) Comparison of sgRNA composition and model error for the large-scale and compact libraries. The CNN-imputed guides had substantially higher predicted activities than those for the large-scale validation set; higher predicted activity was generally associated with higher model error for the validation (red) and imputed (blue) sgRNA sets, consistent with the discrepancy in model performance on each set. (i) Distribution of the number of intermediate-activity mismatched sgRNAs targeting each gene in the compact library. The number of genes with at least 2 intermediate activity sgRNAs is indicated above each histogram; sgRNA activities were quantified for 1907 and 1442 genes in K562 and HeLa cells, respectively. Note that here activities are aggregated by gene as opposed to by series, as was done in Supplementary Fig. 2i. (j) Comparison of phenotypes measured in replicate screens after 12 days of growth in the drug screen. n = 25,518 sgRNAs; r² = squared Pearson correlation coefficient. (k) Comparison of vehicle- (γ) and lovastatin-treatment (τ) growth phenotypes for all sgRNAs in the compact library. Knockdown of HMG-CoA reductase (HMGCR) greatly sensitizes cells to lovastatin, compared to knockdown of other genes such as tubulin (TUBB). n = 25,518 sgRNAs.

(a) Schematic of Perturb-seq strategy to capture single-cell transcriptomes with matched sgRNA identities. (b) Summary of sequencing and perturbation assignment statistics. (c) Distribution of number of cells captured per perturbation. Median: 122 cells per perturbation; 5^th to 95^th percentile: 66 – 277 cells per perturbation. n = 19,587 cells. (d, e) Comparison of (d) growth phenotypes (γ) and (e) relative activities measured in the large-scale mismatched sgRNA screen and in the Perturb-seq experiment. Differences are likely due to the different timescales and the different vectors used. n = 128 sgRNAs; r² = squared Pearson correlation coefficient.

(a) Distribution of target gene expression levels, quantified as target gene UMI count normalized to total UMI count per cell. Cell numbers for each perturbation are listed in Supplementary Table 14. Box plots inside violin plots denote quartile ranges (box), median (center mark), and 1.5 × interquartile range (whiskers). (b) Mean target gene expression levels for target genes with low basal expression levels.

Expression is quantified as raw target gene UMI count. Cell numbers for each perturbation are listed in Supplementary Table 14. Box plots inside violin plots denote quartile ranges (box), median (center mark), and 1.5 × interquartile range (whiskers).

(a) Distributions of total UMI counts in cells with the perfectly matched sgRNA against the indicated genes. Cell numbers for each perturbation are listed in Supplementary Table 14. Box plots inside violin plots denote quartile ranges (box), median (center mark), and 1.5 × interquartile range (whiskers). (b) Left: Comparison of median UMI count per cell and relative growth phenotype in cells with sgRNAs targeting BCR, GATA1, or POLR2H or control cells. Right: Comparison of median UMI count per cell and target gene expression. (c) Cell cycle scores (Methods) for populations of cells with individual sgRNAs. (d) Fraction of cells in indicated cell cycle phase for populations with a negative control sgRNA or sgRNAs targeting CAD. (e) Magnitudes of gene expression change of populations with perfectly matched sgRNAs targeting indicated genes. Magnitude of gene expression change is calculated as sum of z-scores of genes differentially expressed in the series (FDR-corrected p < 0.05 with any sgRNA in the series, two-sided Kolmogorov-Smirnov test, Methods), with z-scores of each gene in individual cells signed by the average direction of change in the population. Cell numbers and violin plots are as in a. (f) Comparison of magnitude of gene expression change to growth phenotype (γ) for all perfectly matched sgRNAs in the experiment. (g) Comparison of relative growth phenotype and magnitude of gene expression change for all individual sgRNAs, as in Fig. 6f but without increased transparency for individual series. (h) Comparison of magnitude of gene expression and target gene knockdown, as in Fig. 6g but without increased transparency for individual series. (i) Comparison of relative growth phenotype and target gene expression, as in Fig. 6f. (j) Comparison of measured growth phenotype (γ, not normalized to strongest sgRNA) and target gene expression, as in Fig. 6f.

(a) Clustered correlation heatmap of perturbations. Gene expression profiles for genes with mean UMI count > 0.25 in the entire population were z-normalized to expression values in cells with negative control sgRNAs and then averaged for populations with the same sgRNA. Crosswise Pearson correlations of all averaged transcriptomes were clustered by the Ward variance minimization algorithm implemented in scipy. Cell numbers for each perturbation are listed in Supplementary Table 14. (b) UMAP projection, distribution of cells with indicated sgRNAs, target gene expression (rolling mean over 50 cells), and magnitudes of transcriptional changes for all differentially expressed genes and selected ISR regulon genes (rolling mean over 50 cells) for cells with knockdown of ATP5E or control cells. n = 2,781 cells total (negative control: 2,084 cells; ATP5E (0.070): 101 cells; ATP5E (0.554) 136 cells; ATP5E (0.914): 137 cells; ATP5E (1.000) 175 cells; ATP5E (1.185) 148 cells). See Methods for details.