Analysis of single-cell CRISPR perturbations indicates that enhancers predominantly act multiplicatively

Abstract

A single gene may have multiple enhancers, but how they work in concert to regulate transcription is poorly understood. To analyze enhancer interactions throughout the genome, we developed a generalized linear modeling framework, GLiMMIRS, for interrogating enhancer effects from single-cell CRISPR experiments. We applied GLiMMIRS to a published dataset and tested for interactions between 46,166 enhancer pairs and corresponding genes, including 264 "high-confidence" enhancer pairs. We found that enhancer effects combine multiplicatively but with limited evidence for further interactions. Only 31 enhancer pairs exhibited significant interactions (false discovery rate <0.1), none of which came from the high-confidence set, and 20 were driven by outlier expression values. Additional analyses of a second CRISPR dataset and in silico enhancer perturbations with Enformer both support a multiplicative model of enhancer effects without interactions. Altogether, our results indicate that enhancer interactions are uncommon or have small effects that are difficult to detect.

Keywords: CRISPR; data simulation; enhancers; gene expression; generalized linear models; genome-wide CRISPR screen; regulatory screen; single-cell sequencing; statistical modeling; transcriptional regulation.

Graphical abstract

Figure 1

Detecting enhancer effects with GLiMMIRS (A) Schematic of GLiMMIRS-base, a GLM to infer the effect of individual enhancers on the expression of target genes. (B) Schematic of GLiMMIRS-int, a GLM to infer joint perturbation effects of two enhancers on target genes and for detecting interaction effects between the enhancers. (C) Schematic of GLiMMIRS-sim, which simulates data from single-cell CRISPR perturbation experiments with RNA-seq readout. It can be used to generate reference values for evaluating the performance of GLiMMIRS-base and GLiMMIRS-int. (D) Scatterplots comparing reference coefficient values generated by GLiMMIRS-sim to coefficient estimates from applying GLiMMIRS-base to the simulated data for $n=1000$ target sites (y axis) using two different values of $X_{perturb}$: (1) a perturbation probability (magenta), calculated as a function of gRNA efficiencies, and (2) an indicator variable based on the presence/absence of a targeting gRNA for the putative enhancer of interest (lavender). Pearson’s correlations (r) are denoted on plots (see also Table S2). (E) Quantile-quantile (QQ) plot of observed versus expected −log₁₀p values. The baseline values (red) are the results from GLiMMIRS-base ($n=609$). The Gasperini et al. values (blue) are previously published results ($n=664$). Mismatch gene (green) is a negative control in which GLiMMIRS-base was applied to randomly selected genes, rather than genes near the enhancer ($n=609$). Shuffled guides (purple) is a negative control in which the vector containing guide perturbation probabilities for each cell was shuffled ($n=609$).

Figure 2

GLiMMIRS-int detects interactions between pairs of enhancers in simulated data (A) Boxplots showing the number of cells containing gRNAs targeting both enhancers belonging to a testable pair in the Gasperini et al. dataset (gray; $n=46,166$ pairs) or in the simulated data with positive enhancer effects (colored; $n=500$ pairs per value of $\lambda$). To consider an enhancer pair to be testable, we required both enhancers to be located within 1 Mb of a common target gene and for the enhancers to be jointly perturbed in at least 10 cells. (B) Histogram of the number of unique gRNAs per cell for data simulated with different values of $\lambda$ and positive interaction effects (colored; $n=50,000$ cells per value of $\lambda$) and the Gasperini et al. dataset (gray; $n=207,324$ cells). (C) Power to detect interaction effects (y axis) in simulated datasets with varying MOIs ($\lambda$) and effect sizes (x axis) using GLiMMIRS-int. See also Table S5. (D) Histogram of significant (FDR corrected p value <0.1) effect sizes estimated by GLiMMIRS-base for $n=560$ individual enhancers in the Gasperini et al. dataset.

Figure 3

Enhancers act multiplicatively to control gene expression, but analysis of CRISPR perturbations provide limited evidence for interactions (A) Distribution of $\Delta AIC$, the difference in AIC between the best-fitting model and the lesser model for $n=264$ high-confidence enhancer pairs and corresponding target genes from Gasperini et al. In every case evaluated, the multiplicative model fit better than the additive model. (B) QQ plot of interaction coefficient p values for $n=264$ high-confidence enhancer pairs, where each individual enhancer had significant effects on the target gene expression. None of the interaction coefficients were significant (Benjamini-Hochberg FDR <0.1). (C) QQ plot of $n=46,166$ enhancer pairs in the entire testable set, where each constituent enhancer did not necessarily have a significant effect on gene expression. Significant interaction coefficients (FDR <0.1) are blue if one of the jointly perturbed cells was an outlier by Cook’s distance and red otherwise. (D) Volcano plot of interaction coefficients for the 46,166 enhancer pairs in the entire testable set. (E) An enhancer pair with a significant negative interaction on the expression of TIMM13. The y axis is TIMM13 expression estimated in cells lacking perturbations to either enhancer (None), cells with perturbations of one enhancer (E1, E2), and cells with joint perturbations of both enhancers (E1 + E2). The blue rectangle is the expected expression in joint perturbation condition under the null model of multiplicative enhancer effects (90% confidence interval [CI] estimated from 100 bootstraps). Whiskers are 90% CIs of expression estimates (from 100 bootstraps). (F) A gene (SNX27) and enhancer pair with a significant positive interaction. (G) Posterior estimate of frequency of enhancer pairs, with interactions estimated from the dataset of 264 high-confidence enhancer pairs. Each line corresponds to a different power to detect interactions. The x axis is the prior belief in enhancer interaction frequency. (H) Posterior estimate of frequency of enhancer pairs with interactions estimated from the 46,166 enhancer pairs in the entire testable set.

Figure 4

No evidence for enhancer interactions from an independent CRISPRi dataset and in silico perturbations (A) Volcano plot of interaction coefficients and −log10 p values estimated by GLiMMIRS-int from the Morris et al. CRISPRi perturbation dataset. All $n=36$ possible pairs of enhancers for the nine candidate enhancers for PTPRC were tested. Both significant interactions were driven by jointly perturbed cells with outlier expression levels. (B) Schematic of in silico perturbation strategy. We input 196-kb sequences into Enformer, where each input sequence contained both candidate enhancers and the target gene. We generated predictions from unperturbed WT sequences, sequences with the first enhancer (EnhA) scrambled, sequences with the second enhancer (EnhB) scrambled, and sequences with both enhancers scrambled (EnhA&EnhB). (C) We compared the Enformer-predicted change in expression of the double mutant (EnhA&EnhB) to the expected expression under a multiplicative model with no enhancer interactions. We analyzed 1,372 enhancer pairs and genes where Enformer-predicted WT expression was >10 and quantified the change in expression as the log ratio of mutant expression to WT expression. The expected expression for the double mutant was computed from the Enformer-predicted expression of the single (EnhA and EnhB) mutants under a multiplicative model of activity. The shading of points corresponds to the Enformer-predicted WT expression level. Pearson's r is indicated.