Profiling the genetic determinants of chromatin accessibility with scalable single-cell CRISPR screens

Abstract

CRISPR screens have been used to connect genetic perturbations with changes in gene expression and phenotypes. Here we describe a CRISPR-based, single-cell combinatorial indexing assay for transposase-accessible chromatin (CRISPR-sciATAC) to link genetic perturbations to genome-wide chromatin accessibility in a large number of cells. In human myelogenous leukemia cells, we apply CRISPR-sciATAC to target 105 chromatin-related genes, generating chromatin accessibility data for ~30,000 single cells. We correlate the loss of specific chromatin remodelers with changes in accessibility globally and at the binding sites of individual transcription factors (TFs). For example, we show that loss of the H3K27 methyltransferase EZH2 increases accessibility at heterochromatic regions involved in embryonic development and triggers expression of genes in the HOXA and HOXD clusters. At a subset of regulatory sites, we also analyze changes in nucleosome spacing following the loss of chromatin remodelers. CRISPR-sciATAC is a high-throughput, single-cell method for studying the effect of genetic perturbations on chromatin in normal and disease states.

Fig. 1 |. CRISPR screens with single-cell combinatorial indexing assay of transposable and accessible chromatin sequencing (CRISPR-sciATAC) enables the joint capture of chromatin accessibility profiles and CRISPR perturbations.

(a) CRISPR-sciATAC workflow with initial barcoding, nuclei pooling and re-splitting, and then second round barcoding. (b) Comparison of bulk ATAC-seq chromatin accessibility profiles from K562 cells using Tn5 and TnY transposases and aggregated CRISPR-sciATAC single cell profiles from 11,104 cells. (c) Guide RNA (gRNA) reads mapping to human or mouse CRISPR libraries (n = 1986 cells). (d) ATAC reads mapping to human or mouse genomes (n = 721 cells). For display purposes, we removed one cell that had >10-fold the average number of ATAC reads. (e) Concordance between the percent of ATAC and gRNA reads mapping to the human and mouse genomes and human and mouse gRNA libraries, respectively, for each cell (n = 496 cells). (f) ATAC-seq fragment size distribution from K562 cells of bulk ATAC-seq data, aggregated CRISPR-sciATAC single cell profiles from 11,104 cells and one representative single cell from CRISPR-sciATAC. (g) Number of CRISPR gRNAs detected per cell. (h) Proportion of cells with 1, 2, or more than 2 gRNAs.

Fig. 2 |. CRISPR-sciATAC reveals changes in accessibility at HOX genes following loss of EZH2.

(a) Heatmap of chromatin accessibility Z-scores at histone and DNA modifications for different CRISPR perturbations (n = 3 gRNA per gene). We converted the fraction of accessible regions for each modification into Z-scores (using all cells in the screen). For visualization, we show the average Z-score for all cells receiving a particular gRNA. (b) Distances in the histone and DNA modifications accessibility profiles shown in panel a between gRNAs targeting different genes and gRNAs targeting the same gene. The distance metric is 1-(Pearson correlation of the Z-scores). (c) Pearson correlation between averaged accessibility Z-scores at histone and DNA modifications of the indicated number of single cells and the average profile of 400 single cells, for cells with either EZH2-targeting or non-targeting (NT) gRNAs. (d) UMAP representation of chromatin accessibility Z-scores at histone and DNA modifications from single cells receiving either EZH2 or NT gRNAs. Also shown is the same UMAP representation with single cells colored by TFBS accessibility enrichment scores for CBX2, CBX8, EZH2, POLR2B, and SIRT6. (e) H3K27me3 ChIP-seq coverage at the HOXA-D loci (top). Changes in accessibility at the HOXA-D loci in cells transduced with EZH2-targeting or NT gRNAs (bottom). *** denotes p < 0.001. (f) CRISPR-sciATAC fragments mapping to the HOXA locus in cells transduced with EZH2-targeting or NT gRNAs (n = 510 cells per condition). The sum of all ATAC fragments over the entire HOXA locus in cells transduced with EZH2-targeting and NT gRNAs is shown on the right. K562 H3K27me3 ChIP-seq coverage is shown at the bottom. (g) Gene expression (qPCR) of EZH2, HOXA3, HOXA5, HOXA11A, HOXA13 and HOXD9 for cells transduced with either EZH2-targeting or NT gRNAs. HOX gene expression in cells targeted by the two more effective EZH2-targeting gRNAs (g1 and g3, as defined by decrease in EZH2 expression compared to non-targeting gRNAs) is greater than in cells targeted by the less effective gRNA (g2) (Student’s t-test, p < 0.05 for HOXA3, HOXA5, HOXA11 and HOXA13, p = 0.09 for HOXD9).

Fig. 3 |. A CRISPR-sciATAC screen targeting 17 chromatin remodeling complexes uncovers widespread disruptions in accessibility upon SWI-SNF disruption.

(a) Chromatin remodeling complex subunits and cofactors targeted in the CRISPR library. (b) Heatmap of chromatin accessibility Z-scores at transcription factor binding sites (TFBSs) for the different chromatin remodeling complexes targeted in the screen. We converted the fraction of accessible regions for each TFBS into Z-scores (using all cells in the screen). For visualization, we first average over all cells for a particular target gene and then average over all genes in the complex. The histograms (left) show the distribution of Z-scores for each complex. The FACT complex is not shown due to a low number of single cells (n = 75 cells). (c) UMAP representation of the genes perturbed in the screen based on the TFBS differential accessibility Z-score profiles. Subunits of the SWI-SNF pBAF complex are labeled with filled circles and gene names. (d) The number of transcription factor binding sites with significant differential accessibility for cells that receive a specific gene-targeting CRISPR perturbation, as compared to cells that receive a non-targeting (NT) control gRNA (FDR q ≤ 0.1). SWI/SNF components and co-factors are highlighted in red. (e) The percent of ATAC fragments in enhancers and promoters in cells transduced with ARID1A-targeting and NT gRNAs. Each point is a single cell. K562 enhancer and promoter genome segmentation is from ENCODE (see Methods). (f) CRISPR-targeted chromatin complex genes with significant differential accessibility at enhancers and/or promoters. (g) Volcano plots showing significant changes in accessibility at TFBSs in cells transduced with ARID1A (left), SMARCA5 (middle) and RCOR1 (right) -targeting gRNAs. Standardized Z-scores are averaged over single cells. Points in red represent TFBSs with a significant change in accessibility (FDR q ≤ 0.1 and |Z-score| > 0.25).

Fig. 4 |. Nucleosome dynamics around transcription factor binding sites (TFBSs) following CRISPR targeting of chromatin remodelers.

(a) Schematic depicting the computational approach to identify changes in nucleosome positions around TFBSs. (b) The difference in nucleosomal distances in gene-targeted cells and nucleosomal distances in non-targeting cells (“peak shift”) across 7 TFBS following CRISPR targeting of chromatin remodelers (top). Bubble-plot of the peak shifts for individual TFBS (bottom). The color of the bubble corresponds to the peak shift (nt) and the size of the bubble represents the empirical p-value calculated by a label permutation test. (c) The number of nucleosome expansion and compaction events around TFBSs following CRISPR targeting of chromatin remodelers. (d) Coverage profiles of mononucleosomal fragments around AP1 binding sites in cells transduced with ARID1A-targeting (blue) and non-targeting (NT) (grey) gRNAs (top) and in cells transduced with EP400-targeting (blue) and NT (grey) gRNAs (bottom). Dashed lines represent the most highly covered base in each peak. (e) Peak shifts in TFBSs located in enhancers and promoters. Each point is a CRISPR-perturbed gene (average of all gRNAs for that gene). (f) Peak shifts in TFBSs located in enhancers and promoters in SFMBT1-targeted cells (left). Coverage profiles of mononucleosomal fragments in cells transduced with SFMBT1-targeting (blue) and NT (grey) gRNAs around AP1 binding sites in promoters (top) and in enhancers (bottom). (g) Peak shifts in TFBSs located in enhancers and promoter in SMARCB1 targeted cells (left). Coverage profiles of mononucleosomal fragments in cells transduced with SMARCB1-targeting (blue) and NT (grey) gRNAs around RAD21 binding sites in promoters (top) and in enhancers (bottom). For panels d, f, and g, the shaded regions represent s.e.m. (n = 3 gRNAs).