Abstract
Article Figures and data Abstract Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Single-cell measurements of cellular characteristics have been instrumental in understanding the heterogeneous pathways that drive differentiation, cellular responses to signals, and human disease. Recent advances have allowed paired capture of protein abundance and transcriptomic state, but a lack of epigenetic information in these assays has left a missing link to gene regulation. Using the heterogeneous mixture of cells in human peripheral blood as a test case, we developed a novel scATAC-seq workflow that increases signal-to-noise and allows paired measurement of cell surface markers and chromatin accessibility: integrated cellular indexing of chromatin landscape and epitopes, called ICICLE-seq. We extended this approach using a droplet-based multiomics platform to develop a trimodal assay that simultaneously measures transcriptomics (scRNA-seq), epitopes, and chromatin accessibility (scATAC-seq) from thousands of single cells, which we term TEA-seq. Together, these multimodal single-cell assays provide a novel toolkit to identify type-specific gene regulation and expression grounded in phenotypically defined cell types. Introduction Peripheral blood mononuclear cells (PBMCs) purified using gradient centrifugation are a major source of clinically relevant cells for the study of human immune health and disease (Böyum, 1968). Like most other human tissues, PBMCs are a complex, heterogeneous mixture of cell types derived from common stem cell progenitors (Laurenti and Göttgens, 2018). Despite the genome being mostly invariant between different PBMC cell types, each immune cell type performs an important and distinct function. Understanding the genomic regulatory landscape that controls lineage specification, cellular maturation, activation state, and functional diversity in response to intra- and extracellular signals is key to understanding the immune system in both health and disease (Satpathy et al., 2019; Wang et al., 2020; Zheng et al., 2020). Recent improvements in single-cell genomic methods have enabled profiling of the regulatory chromatin landscape of complex cell-type mixtures. In particular, droplet-based single-nucleus or single-cell assays for transposase-accessible chromatin (snATAC-seq, scATAC-seq, dscATAC-seq, mtscATAC-seq) allow profiling of open chromatin at single-cell resolution (Buenrostro et al., 2015; Lareau et al., 2019). Promising new methods have combined scATAC-seq with simultaneous measurement of nuclear mRNAs (e.g., sci-CAR, Cao et al., 2018; SNARE-seq, Chen et al., 2019; SHARE-seq, Ma et al., 2020) or in combination with cell surface epitopes (ASAP-seq, Mimitou et al., 2020). However, a unifying approach for all three modalities that can be applied to highly specified functional immune cell types has yet to emerge in the landscape of single-cell methods. We systematically tested whole cell and nuclear purification and preparation methods for PBMCs to overcome limitations that restricted previous assays to measurement of only nuclear components (ATAC and nuclear RNAs) or proteins on the cell surface. We found that intact, permeabilized cells perform extremely well for scATAC-seq, exceeding conventional scATAC-seq on nuclei by some measures (Figure 1b). This insight enables a new protocol analogous to Cellular Indexing of Transcriptomes and Epitopes (CITE-seq; Stoeckius et al., 2017) to measure both surface protein abundance and chromatin accessibility: integrated cellular indexing of chromatin landscape and epitopes (ICICLE-seq, Figures 1a and 3). Finally, we demonstrate that our optimized permeable cell approach can be combined with a droplet-based multiomics platform to enable the simultaneous measurement of three different molecular compartments of the cell: mRNA (by scRNA-seq), protein (using oligo-tagged antibodies), and DNA (by scATAC-seq), which we term TEA-seq after Transcription, Epitopes, and Accessibility (Figure 4). These assays enable a new, more unified view into the molecular underpinnings of gene regulation and expression at the single-cell level. Figure 1 with 6 supplements see all Download asset Open asset Improvements to scATAC-seq methods to enable permeabilized cell profiling. (a) Schematic overview of major steps in snATAC, scATAC, and ICICLE-seq methods. (b) Comparison of quality control characteristics of ATAC-seq libraries generated from nuclei isolation and permeabilized cells, with and without fluorescence-activated cell sorting. Top panels show signal-to-noise as assessed by fraction of reads in peaks on the y-axis and quantity of unique fragments per cell barcode on the x-axis. Lower panels display fragment length distributions obtained from paired-end sequencing of ATAC libraries. Colored lines represent barcodes that pass QC filters; gray lines represent barcodes failing QC (non-cell barcodes). All libraries were equally downsampled to 200 million total sequenced reads for comparison. Colors in (b) are reused in remaining panels. (c) Total coverage of Tn5 footprints summed across all transcription start sites (TSS). Tn5 footprints are 29 bp regions comprising the 9 bp target-site duplication (TSD) and 10 bp on either side, which represent accessible chromatin for each transposition event. (d) Total coverage of TSD centers summed over a set of 100,000 genomic CTCF motifs found in previously published DNase-hypersensitive sites (Meuleman et al., 2020). TSD centers are obtained by shifting +4 and −5 bp from the 5’ and 3’ ends of uniquely aligned fragments, respectively. (e) Barplot representations of the fraction of total aligned reads in various QC categories. Fragments overlapping a previously published peak set for peripheral blood mononuclear cell dscATAC-seq (Lareau et al., 2019) are in the ‘Overlap Peaks’ category. Unique fragments are the remaining uniquely aligned fragments that do not overlap peak regions. ‘Waste’ reads were not aligned or were assigned to cell barcodes with fewer than 1000 total reads. (f) Violin plots showing distributions of QC metrics. Median (wide bar) and 25th and 75th quantiles (whiskers and narrow bars) are overlaid on violin plots. Median values are also in Table 1. Note that the y-axis of the first panel is on a logarithmic scale; remaining panels are linear. snATAC: single-nucleus assays for transposase-accessible chromatin; scATAC: single-cell assays for transposase-accessible chromatin; ICICLE-seq: integrated cellular indexing of chromatin landscape and epitopes. Figure 1—source data 1 Single cell metadata and QC metrics for scATAC-seq experiments . https://cdn.elifesciences.org/articles/63632/elife-63632-fig1-data1-v2.zip Download elife-63632-fig1-data1-v2.zip Figure 1—source data 2 Fragment size distribution data . https://cdn.elifesciences.org/articles/63632/elife-63632-fig1-data2-v2.zip Download elife-63632-fig1-data2-v2.zip Figure 1—source data 3 TSS Footprint pileups. https://cdn.elifesciences.org/articles/63632/elife-63632-fig1-data3-v2.zip Download elife-63632-fig1-data3-v2.zip Figure 1—source data 4 CTCF Tn5 target site duplication center pileups. https://cdn.elifesciences.org/articles/63632/elife-63632-fig1-data4-v2.zip Download elife-63632-fig1-data4-v2.zip Figure 1—source data 5 Fraction of reads per alignment category. https://cdn.elifesciences.org/articles/63632/elife-63632-fig1-data5-v2.zip Download elife-63632-fig1-data5-v2.zip Table 1 QC metrics summary for experiments displayed in Figure 1. Source typeFACS depletionN pass QCMedian N fragmentsMedian mitochondrial readsMedian uniqueMedian in TSSMedian in peaksNucleiUnsorted4719734443 | 0.6%5247 | 71.6%1332 | 24.9%2306 | 43.9%NucleiDead/debris552610,52659 | 0.6%7284 | 69.4%2186.5 | 30.1%3647 | 50.5%NucleiDead/debris/neutrophils776919,972136 | 0.9%11,528 | 59.1%4846 | 41.5%7503 | 64.8%Permeabilized cellsUnsorted63295541100 | 1.9%3308 | 59.8%871 | 26.4%1390 | 42.2%Dead/debris69566733.5120 | 2.0%3795.5 | 56.8%1219.5 | 32.6%1874 | 50.2%Dead/debris/neutrophils984914,069514 | 4.0%4756 | 34.3%2536 | 54.4%3650 | 76.8%For all metrics to the right of median N fragments, both the absolute number and a percentage are provided. Median % mitochondrial and median % unique were calculated as a fraction of total fragments; % in TSS and % in peaks were calculated as a fraction of unique fragments. TSS: transcription start sites; FACS: fluorescence-activated cell sorting. Results Optimization of single-nucleus and single-cell ATAC-seq of PBMCs As a baseline for optimization and cell surface retention, we performed snATAC-seq as recommended by 10x Genomics, with a protocol based on the Omni-ATAC workflow (Corces et al., 2017). This single-nucleus assay utilizes a combination of hypotonic lysis, detergents, and a saponin to isolate nuclei without retaining mitochondrial DNA. After performing snATAC-seq using this method, sequencing, and tabulating data quality metrics (Materials and methods), we identified two major populations of cell barcodes (Figure 1b, left panel): (1) a large number of barcodes, shown in gray, that have a low number of unique fragments and a low fraction of reads in peaks (FRIP). These barcodes contain little useful information but consume 80% of total sequenced reads (Figure 1e, non-cell barcodes) at a sequencing depth of 200 million reads per library (20,000 reads per expected barcode). (2) Barcodes with higher quality as measured by FRIP (red points) that contain enough information to attempt downstream analysis. The loss of 80% of sequenced reads to non-cell barcodes is costly. Previous studies of scRNA-seq data have shown that cellular lysis can release ambient RNA that increases the abundance of low-quality barcodes and contaminates droplets, yielding barcodes with both cellular and ambient RNAs that reduces the accuracy of the transcriptional readout (Marquina-Sanchez et al., 2020). We reasoned that nuclear isolation protocols may cause the release of ambient DNA, causing a similar effect in scATAC-seq datasets. Optimization of nuclear lysis protocols, especially changing to less stringent detergents, provided increased FRIP and decreased non-cell barcodes (Figure 1—figure supplement 1, Figure 1—figure supplement 2a). Hypotonic lysis conditions used in these protocols may also be a biophysical stressor to the native chromatin state, as previously observed (Lima et al., 2018). To reduce perturbation of chromatin and retain the cell surface for multimodal assays, we performed cell membrane permeabilization under isotonic conditions to allow access to the nuclear DNA without isolating nuclei through hypotonic lysis. The saponin digitonin was used to cause concentration-dependent selective permeabilization of cholesterol-containing membranes while leaving inner mitochondrial membranes intact, preventing high levels of Tn5 transposition in mitochondrial DNA (Adam et al., 1990; Colbeau et al., 1971). Digitonin has previously been used for ATAC-seq assays under hypotonic conditions in Fast-ATAC (Corces et al., 2016) and plate-based scATAC-seq (Chen et al., 2018b) protocols. Permeabilization of intact cells under isotonic conditions greatly reduced the number of non-cell barcodes and their contribution to sequencing libraries (Figure 1—figure supplement 2b). Removal of neutrophils greatly increases PBMC scATAC-seq quality We observed that PBMCs purified by leukapheresis rather than Ficoll gradient centrifugation had consistently higher FRIP scores and fewer non-cell barcodes (Figure 1—figure supplement 2c). A major difference between Ficoll-purified PBMCs and leukapheresis-purified PBMCs was the presence of residual neutrophils in Ficoll-purified samples. We tested removal of dead cells and debris with and without removal of neutrophils using fluorescence-activated cell sorting (FACS) from PBMC samples with high neutrophil content (Figure 1a, Figure 1—figure supplement 4). When applied to either nuclei (Figure 1b, left panels) or permeabilized cells (Figure 1b, right panels), there was a large increase in FRIP and reduction in non-cell barcodes in our scATAC-seq libraries (Figure 1e). Removal of neutrophils did not have an adverse effect on leukapheresis-purified PBMC data (Figure 1—figure supplement 2c, right panel). Depletion of neutrophils using anti-CD15 magnetic beads also improved data quality (Figure 1—figure supplement 2d, Figure 1—figure supplement 5), though not to the same extent as FACS-based depletion. Staining and flow cytometry using an eight-antibody panel (Supplementary file 1) on Ficoll- and leukapheresis-purified PBMCs showed minimal effect of the magnetic bead treatment on non-neutrophil cell-type abundance (Supplementary file 2). Comparing single-cell and single-nucleus ATAC characteristics We assessed the quantitative and qualitative differences between nuclei and permeabilized cell protocols with and without sorting by performing both protocols on a single set of input cells. To fairly compare data quality across methods, we equally downsampled raw data at the level of raw sequenced reads per well for each sample that we compare directly (see Materials and methods for details). We then utilized a uniform set of transcription start site (TSS) regions (TSS ± 2 kb) and a previously published PBMC peak set (Lareau et al., 2019, Materials and methods) as reference regions for computation of fraction of reads in TSS (FRITSS) and FRIP, respectively. Permeabilized cells yielded many more high-quality cell barcodes than nuclear preps using equal loading of cells or nuclei (15,000 loaded, expected 10,000 captured, Table 1). scATAC-seq libraries prepared from nuclei contained many more reads originating from nucleosomal DNA fragments (Figure 1b, lower panels), and non-cell barcodes from nuclei (gray lines) contained more of these fragments than cell barcodes. Thus, an overabundance of mononucleosomal fragments may indicate non-cell fragment contamination. Libraries from permeabilized cells consisted almost entirely of short fragments, suggesting that permeabilization under isotonic conditions did not loosen or release native chromatin structure at the time of tagmentation (Figure 1b, lower panels). Previous bulk ATAC-seq studies have shown that differing nuclear isolation protocols lead to varying amounts of mononucleosomal fragments (Li et al., 2019a). In agreement with in vitro experiments studying the effects of low salt on nucleosomal arrays (Allahverdi et al., 2015), this further suggests that hypotonic lysis leads to alteration of chromatin structure, raising the possibility of artifactual measurements of accessibility in nuclei-based ATAC-seq. To assess the effect that this difference has on the data obtained by each method, we overlaid Tn5 footprints near TSS ( Figure 1c) and CTCF transcription factor binding sites (TFBS, Figure 1d). The signal at TSS was retained in permeabilized cells, but positions flanking the TSS (occupied by neighboring nucleosomes) had reduced signal compared to isolated nuclei (examined in detail in Figure 1—figure supplement 6). At CTCF motifs, we observed nearly identical patterns of accessibility in both nuclei and permeabilized cells, suggesting that scATAC-seq signal at regulatory TFBS is retained in permeabilized cells. Neutrophil and dead cell removal improved the quality of nuclear scATAC-seq libraries, which yielded the highest number of unique fragments, reads in TSS, and reads in peaks, though at the cost of fewer high-quality barcodes captured when compared to permeabilized cells (Table 1). Overall, permeabilized intact cells obtained by FACS had the highest FRIP and FRITSS scores, fewest non-cell barcodes, and greatest cell capture efficiency with only a modest increase in mitochondrial reads (Table 1, Figure 1e, f). Improved label transfer and differential analysis We next examined the effect of methodological differences on downstream biological analyses (Figure 2). Removal of neutrophils greatly improved the ability to separate various cell types in uniform manifold approximation and projection (UMAP) projections of both nuclei and cells (Figure 2a, b). To provide ground truth for label transfer, we performed flow cytometry on an aliquot of the same PBMC sample used for scATAC-seq above. A panel of 25 antibodies (Supplementary file 3) was used to determine the proportion of each of the 12 cell types used to label the scATAC-seq cells in the PBMC sample (Figure 2—figure supplement 1 and Supplementary file 4). Label transfer was performed using the ArchR package (Granja et al., 2020) to generate gene scores that enabled label transfer from a reference scRNA-seq dataset using the method provided in the Seurat package (Stuart et al., 2019) (Materials and methods). Using these tools, removal of neutrophils improved label transfer scores, and permeabilized cells yielded more cells with high label transfer scores than nuclei-based approaches (Figure 2b, c). In addition, permeabilized cells provided labels most similar to the cell-type proportions identified by flow cytometry (Figure 2d), with identification of CD8 effector cells only observed in scATAC-seq with permeabilized cells. All methods yielded fewer CD16+ monocytes than observed by flow cytometry, suggesting that CD16+ monocytes may be lost during scATAC-seq using either nuclei or permeabilized cells, or that label transfer methods were not conducive to identifying this cell type (Figure 2d). After labeling cell types, we used ArchR to call peaks for each cell type and perform pairwise tests of differential accessibility between each pair of cell types (Figure 2—figure supplement 2a). We found many more differentially accessible sites in both cells and nuclei after removal of neutrophils. Differential accessibility was also used to identify differentially enriched TFBS motifs in each cell type (Figure 2—figure supplement 2b). Without neutrophil removal (nuclei unsorted, top panel), we were unable to identify significantly enriched motifs in B cells and NK cells that were readily apparent in data from clean nuclei or permeabilized cells (bottom two panels). Together, these results demonstrate that neutrophil removal and the use of permeabilized cells allow for identification of specific cell types and TFBS motifs that are involved in regulation of gene expression. Figure 2 with 2 supplements see all Download asset Open asset Improvements to 2D projection and label transfer for scATAC-seq data. (a) Uniform manifold approximation and projection (UMAP) projection plots for corresponding datasets in Figure 1. Points are colored based on a common scale of fraction of reads in peaks, bottom right. The number of cells in each panel is displayed in Figure 1b. (b) UMAP projection plots colored based on cell type obtained by label transfer from scRNA-seq (Materials and methods). Colors for cell types are below, to the right. (c) To visualize the number and quality of transferred labels, we ranked all cells based on the Seurat label transfer score obtained from label transfer results and plotted lines through the score (y-axis) vs. rank (x-axis) values. (d) Barplot showing the fraction of cells in each dataset that were assigned each cell-type label. The top row shows cell-type proportions for the same peripheral blood mononuclear cell sample obtained by 25-color immunotyping flow cytometry (Materials and methods, Supplementary file 3, and Figure 2—figure supplement 1). Root-mean-square deviation values were computed by comparison of labeled cell-type proportions to values derived from flow cytometry (Supplementary file 4). Colors for cell types are to the right of the barplot. Figure 2—source data 1 Single cell UMAP coordinates and labeling scores. https://cdn.elifesciences.org/articles/63632/elife-63632-fig2-data1-v2.zip Download elife-63632-fig2-data1-v2.zip Figure 2—source data 2 Fractions of cells assigned to each type by flow cytometry and scATAC-seq. https://cdn.elifesciences.org/articles/63632/elife-63632-fig2-data2-v2.zip Download elife-63632-fig2-data2-v2.zip Joint measurement of accessibility and epitopes with ICICLE-seq Under standard scATAC-seq protocols, removal of the cell membrane severs the connection between the cell surface and the chromatin state of cells. To test our ability to measure cell surface proteins and chromatin state simultaneously on permeabilized cells, we modified our optimized permeabilized cell scATAC-seq methodology to incorporate measurements using commercially available barcoded antibody reagents, and we term this new method ICICLE-seq (Figure 3 and Materials and methods). The ICICLE-seq protocol utilizes a custom Tn5 transposome complex with capture sequences compatible with the 10x Genomics 3’ scRNA-seq gel bead capture reaction for simultaneous capture of ATAC fragments and polyadenylated antibody barcode sequences (Figure 3—figure supplement 1 and Supplementary file 5). Antibody-derived tags (ADTs) from oligo-antibody conjugates and ATAC-seq fragments can then be selectively amplified by PCR to generate separate libraries for sequencing (Figure 3—figure supplement 1). Due to the nature of fragment capture in this system, we obtain both a cell barcode and a single-end scATAC-seq read from the two ends of the paired-end sequencing reaction. We performed ICICLE-seq on a leukapheresis-purified PBMC sample using a 46-antibody panel (Supplementary file 6) and were able to obtain 10,227 single cells with both scATAC-seq and ADT data from three capture wells that passed QC unique ATAC fragments FRIP ATAC QC had a median of ADT unique per cell (Figure 3—figure supplement 2b). UMAP projection and ATAC label transfer on ICICLE-seq data had resolution similar to scATAC-seq on intact permeabilized cells after dead cell and debris removal (Figure However, the data quality was by the capture method and the single-end readout of fragments. the ICICLE-seq results that permeabilized cells enabled simultaneous capture of chromatin accessibility in the and high-quality cell surface We were able to the ADT data to and identify cell types based on their cell surface (Figure UMAP based on ADT and allowed identification of (Figure based on of markers with (Figure Figure 3—figure supplement 2c). Thus, ICICLE-seq provided a key of for directly capture of functional cell types and chromatin Figure 3 with 2 supplements see all Download asset Open asset profiling of chromatin accessibility and cell surface epitopes. (a) Uniform manifold approximation and projection (UMAP) projection of integrated cellular indexing of chromatin landscape and epitopes cells based on single-cell assays for transposase-accessible chromatin (scATAC-seq) data. are colored based on fraction of reads in 10,227 cells QC are (b) UMAP projection of scATAC-seq as in are colored based on cell-type labels obtained by ArchR label transfer (Materials and methods). (c) UMAP projection of ICICLE-seq cells based on data. are colored to the total number of unique across all (d) UMAP projection based on ADT as in colored to cell-type labels derived from expression (Materials and methods). (e) of median ADT values for each in each cell type labeled in are in each row between and the for each Figure data 1 type labels and UMAP coordinates for ICICLE-seq cells. Download Figure data 2 ADT data for ICICLE-seq. Download measurement of epitopes, and accessibility with TEA-seq the release of a commercially available platform for simultaneous capture of and ATAC-seq from single we reasoned that permeabilized cells be used to perform simultaneous capture of three major molecular DNA can be captured using scATAC-seq, RNA be captured using and protein abundance can be captured using polyadenylated antibody barcodes, which we term TEA-seq after Epitopes, and After and optimization of key we were able to obtain libraries on the 10x Genomics ATAC platform that combined all three of these measurements for thousands of single cells (Figure using a panel of oligo-tagged antibodies Table 6). After data for cells into we identified cell barcodes that passed the QC for scATAC-seq and had unique ATAC-seq fragments unique ATAC fragments, Figure and had by scRNA-seq RNA median Figure and ADT ADT Figure QC metrics are provided in Figure supplement 1. Figure 4 with 4 supplements see all Download asset Open asset measurement of epitopes, and (a) for the major steps in epitopes, and accessibility (b) unique single-cell assays for transposase-accessible chromatin (scATAC-seq) fragments and scRNA-seq unique for each TEA-seq cell In and barcodes are displayed in QC are by (c) unique scATAC-seq fragments and tags for each cell (d) scATAC-seq QC unique scATAC-seq fragments and fraction of reads in peaks scores for each cell total cells are displayed with unique ATAC QC are by (e) Uniform manifold approximation and projection (UMAP) projections generated using each of the three modalities cells QC barcodes) are in f). (f) A UMAP projection generated using that all three of the measured showing the (x-axis) and (y-axis) values for each peak that was found to be with the gene or antibody in TEA-seq data. at the show the distribution of scores for RNA and protein shows not found to be in each method were assigned a score of showing between peaks (red and the gene based on protein expression and gene expression are by colored based on the score scale for both panels to the The bottom panel shows the gene All coordinates are from the genome as in for the gene and as in for the gene Figure data 1 Single cell quality metrics for TEA-seq samples. Download Figure data 2 type labels and UMAP coordinates for TEA-seq samples. Download Figure data 3 to gene and peak to protein link Download the of scRNA-seq we were able to perform cell-type label transfer from RNA to RNA using Seurat label transfer (Stuart et al., 2019) rather than transfer from RNA to ATAC (Materials and methods). each we were able to perform reduction and UMAP projection to separate of the between cells (Figure as for (Figure b). However, these do not the of multimodal To to of our we extended a method by and et for paired analysis et al., 2020) to allow an number of simultaneously measured modalities to to the (Materials and methods). After this we generate a UMAP with from all three simultaneously measured modalities (Figure which cell-type We found that expression is across all three modalities for some as though we found that across all modalities is not a of functional