Systematic comparison of high-throughput single-cell RNA-seq methods for immune cell profiling

Abstract

Background: Elucidation of immune populations with single-cell RNA-seq has greatly benefited the field of immunology by deepening the characterization of immune heterogeneity and leading to the discovery of new subtypes. However, single-cell methods inherently suffer from limitations in the recovery of complete transcriptomes due to the prevalence of cellular and transcriptional dropout events. This issue is often compounded by limited sample availability and limited prior knowledge of heterogeneity, which can confound data interpretation.

Results: Here, we systematically benchmarked seven high-throughput single-cell RNA-seq methods. We prepared 21 libraries under identical conditions of a defined mixture of two human and two murine lymphocyte cell lines, simulating heterogeneity across immune-cell types and cell sizes. We evaluated methods by their cell recovery rate, library efficiency, sensitivity, and ability to recover expression signatures for each cell type. We observed higher mRNA detection sensitivity with the 10x Genomics 5' v1 and 3' v3 methods. We demonstrate that these methods have fewer dropout events, which facilitates the identification of differentially-expressed genes and improves the concordance of single-cell profiles to immune bulk RNA-seq signatures.

Conclusion: Overall, our characterization of immune cell mixtures provides useful metrics, which can guide selection of a high-throughput single-cell RNA-seq method for profiling more complex immune-cell heterogeneity usually found in vivo.

Keywords: High throughput sequencing; Immune-cell profiling; Single cell; Single-cell RNA-seq; Transcriptomics.

Fig. 1

Overview of high-throughput single-cell benchmarking experiments. Experiments were performed using four immune cell lines to benchmark cell recovery, transcript detection sensitivity, concordance to bulk RNA-seq and differentially-expressed gene identification

Fig. 2

Library-pool and cell-capture efficiencies: a Cell capture efficiency was measured by the number of cell identifiers (CIDs) above the inflection point of the rank ordered reads/CID plot (knee plot) relative to the number of cells loaded on the instrument. Horizontal lines indicate theoretical capture efficiency based on bead/cell loading concentrations or manufacturer’s guidelines. b Library pool efficiency was measured by the number of reads in CIDs above the inflection point

Fig. 3

Transcript detection sensitivity: a Distributions of unique molecular identifiers (UMIs) and b genes detected in EL4 cells by sample are plotted. c Numbers of UMIs or d genes detected versus numbers of reads per cell for each cell type are plotted. e Accumulated average numbers of genes detected from aggregated data of subsamples up to 50 cells are plotted. f Dropout modeling (dropout rate versus FPKM of bulk sequencing) for EL4 cells by method are shown. A left-shifted curve indicates higher sensitivity, that is, fewer dropouts at lower expression levels. Sensitivity of methods for EL4 cells ranked in the following order: 10x 3′ v3 > 10x 5′ v1 > 10x 3′ v2 > ddSEQ > Drop-seq > ICELL8 3′ DE > ICELL8 3′ DE-UMI. Cells with high mitochondrial expression rates were excluded from this calculation

Fig. 4

Correlation to bulk RNA-seq: a Pearson correlation (r) of cell identifiers (CIDs) to bulk RNA-seq data using highly-expressed variable genes. Only r values above 0.2 were included in plot. b Average Pearson correlation using all genes for aggregated data of 50 subsamples of up to 50 cells are plotted

Fig. 5

Differentially-expressed (DE) gene detection: a Fold change (FC) versus false discovery rate (FDR) calculated using a hurdle model (MAST) for mouse genes in EL4 vs IVA12 cells. Shown is a representative subsample of mouse cells (n=199) using the 10x 3′ v2 method demonstrating the criteria for declaring DE genes (FDR < 10^− 4); DE genes are highlighted in red. b Number of significant DE genes calculated using MAST between EL4 and IVA12 cells by method. Error bars represent the 95% confidence interval. The total number of significant DE genes are plotted in red, the number of DE genes with > 1.5-fold difference in expression in bulk RNA-seq (5868 genes) are plotted in cyan. c Median bulk RNA-seq expression (FPKM) of all significant DE genes (red) or DE genes with > 1.5-fold difference (cyan). Error bars represent 95% confidence interval

Fig. 6

Cell recovery by cell identifier (CID) thresholding: a Example of using the transposed log-log empirical cumulative density plot of the total counts of each CID to identify cell-containing droplets . Common thresholding points, the ‘knee’ and the ‘inflection’ are indicated with arrows. The knee is the point at which the signed curvature is minimized, the inflection is the point at which the first derivative is minimized. b The fraction of cells above the knee or inflection are plotted. c Fraction of cells below mitochondrial rate threshold (listed in Supplement Table 4) relative to knee point. Samples are colored by cell sample mixture listed in Supplement Table 2

Fig. 7

Dropout rates by cell type: a Distribution of reads across cell types is plotted by method. b Dropout rate models for cell types are shown. c 10x 3′ v3 cells were binned by number of unique molecular identifiers (UMIs) and distributions of nUMIs for each cell type in each bin are plotted. d Gene Detection 50 (GD₅₀) rates, expression level at 0.5 probability of the dropout model, are plotted for each cell type in 10x 3′ v3 experiments by bin. e Dropout models in each bin for EL4, IVA12 and Jurkat cells are plotted along with the model for TALL-104 cells in bin 1

Fig. 8

mRNA capture variation across peripheral blood mononuclear cells (PBMCs). a Single-cell data generated with the 10x 3′ v3 and 10x 5′ v1 chemistries were projected onto an annotated PBMC CITE-Seq reference dataset. b Violin plots of log normalized expression of common immune-cell markers. c nGenes and d nUMIs detected for each cell class from each method. CTL-cytotoxic, TCM – T central memory, TEM – T effector memory, Treg – T regulatory cells, dnT – double negative T, gdT – gamma delta T, MAIT - Mucosal associated invariant T, pDC - Plasmacytoid dendritic cell, ASDC - AXL+ dendritic cell, cDC – classical dendritic cell, MC – monocyte, Eryth-Erythrocyte, ILC - Innate lymphoid cell, HSPC - Hematopoietic stem and progenitor cell