Multiplexed droplet single-cell RNA-sequencing using natural genetic variation.

Droplet single-cell RNA-sequencing (dscRNA-seq) has enabled rapid, massively parallel profiling of transcriptomes. However, assessing differential expression across multiple individuals has been hampered by inefficient sample processing and technical batch effects. Here we describe a computational tool, demuxlet, that harnesses natural genetic variation to determine the sample identity of each droplet containing a single cell (singlet) and detect droplets containing two cells (doublets). These capabilities enable multiplexed dscRNA-seq experiments in which cells from unrelated individuals are pooled and captured at higher throughput than in standard workflows. Using simulated data, we show that 50 single-nucleotide polymorphisms (SNPs) per cell are sufficient to assign 97% of singlets and identify 92% of doublets in pools of up to 64 individuals. Given genotyping data for each of eight pooled samples, demuxlet correctly recovers the sample identity of >99% of singlets and identifies doublets at rates consistent with previous estimates. We apply demuxlet to assess cell-type-specific changes in gene expression in 8 pooled lupus patient samples treated with interferon (IFN)-β and perform eQTL analysis on 23 pooled samples.

Droplet single cell RNA-sequencing (dscRNA-seq) has increased substantially the throughput of single cell capture and library preparation[1, 10], enabling the simultaneous profiling of thousands of cells. Improvements in biochemistry[11, 12] and microfluidics[13, 14] continue to increase the number of cells and transcripts profiled per experiment. But for differential expression and population genetics studies, sequencing thousands of cells each from many individuals would better capture inter-individual variability than sequencing more cells from a few individuals. However, in standard workflows, dscRNA-seq of many samples in parallel remains challenging to implement. If the genetic identity of each cell could be determined, pooling cells from different individuals in one microfluidic run would result in lower per-sample library preparation cost and eliminate confounding effects. Furthermore, if droplets containing multiple cells from different individuals could be detected, pooled cells could be loaded at higher concentrations, enabling additional reduction in per-cell library preparation cost.

Here we develop an experimental protocol for multiplexed dscRNA-seq and a computational algorithm, demuxlet, that harnesses genetic variation to determine the genetic identity of each cell (demultiplex) and identify droplets containing two cells from different individuals (Fig. 1a). While strategies to demultiplex cells from different species[1, 10, 17] or host and graft samples[17] have been reported, simultaneously demultiplexing and detecting doublets from more than two individuals has not been possible. Inspired by models and algorithms developed for detecting contamination in DNA sequencing[18], demuxlet is fast, accurate, scalable, and compatible with standard input formats[17, 19, 20].

Figure 1

Demuxlet: demultiplexing and doublet identification from single cell data

a) Pipeline for experimental multiplexing of unrelated individuals, loading onto droplet-based single-cell RNA-sequencing instrument, and computational demultiplexing (demux) and doublet removal using demuxlet. Assuming equal mixing of 8 individuals, b) 4 genetic variants can recover the sample identity of a cell, and c) 87.5% of doublets will contain cells from two different samples.

Demuxlet implements a statistical model for evaluating the likelihood of observing RNA-seq reads overlapping a set of single nucleotide polymorphisms (SNPs) from a single cell. Given a set of best-guess genotypes or genotype probabilities obtained from genotyping, imputation or sequencing, demuxlet uses maximum likelihood to determine the most likely donor for each cell using a mixture model. A small number of reads overlapping common SNPs is sufficient to accurately identify each cell. For a pool of 8 individuals and a set of uncorrelated SNPs each with 50% minor allele frequency (MAF), 4 reads overlapping SNPs are sufficient to uniquely assign a cell to the donor of origin (Fig. 1b) and 20 reads overlapping SNPs can distinguish every sample with >98% probability in simulation (Supplementary Fig. 1). We note that by multiplexing even a small number of individuals, the probability that a doublet contains cells from different individuals is very high (1 – 1/N, e.g., 87.5% for N=8 samples) (Fig. 1C). For example, if a 1,000-cell run without multiplexing results in 990 singlets with a 1% undetected doublet rate, multiplexing 1,570 cells each from 63 samples can theoretically achieve the same rate of undetected doublets, producing up to a 37-fold more singlets (36,600) if the sample identity of every droplet can be perfectly demultiplexed (Supplementary Fig. 2, see Methods for details). To minimize the effects of sequencing doublets, profiling 22,000 cells multiplexed from 26 individuals generates 23-fold more singlets at the same effective doublet rate (Supplementary Fig. 3).

We first assess the performance of multiplexed dscRNA-seq through simulation. The ability to demultiplex cells is a function of the number of individuals multiplexed, the depth of sequencing or number of read-overlapping SNPs, and relatedness of multiplexed individuals. We simulated 6,145 cells (5,837 singlets and 308 doublets) from 2 – 64 individuals from the 1000 Genomes Project[21]. We show that 50 SNPs per cell allows demultiplexing of 97% of singlets and identification of 92% of doublets in pools of up to 64 individuals (Supplementary Figs. 4–5, see Methods for details). Simulating a range of sequencing depths, we determined that 50 SNPs can be obtained with as few as 1,000 unique molecular identifiers (UMIs) per cell (Supplementary Fig. 6), and recommended sequencing depths of standard dscRNA-seq workflows would capture hundreds of SNPs. To assess dependence on the relatedness of multiplexed individuals, we simulated 6,145 cells from a set of 8 related individuals from 1000 Genomes[21]. In this simulation, 50 SNPs per cell would allow demuxlet to correctly assign over 98% of cells (Supplementary Fig. 7). These results suggest optimal multiplexed designs where cells from tens of unrelated individuals should be pooled, loaded at concentrations 2–10x higher than standard workflows, and sequenced to at least 1,000 UMIs per cell.

We evaluate the performance of demuxlet by analyzing a pool of peripheral blood mononuclear cells (PBMCs) from 8 lupus patients. By sequential pairwise pooling, three pools of equimolar concentrations of cells were generated (W1: patients S1-S4, W2: patients S5-S8 and W3: patients S1-S8) and each loaded in a well on a 10X Chromium instrument (Fig. 2a). 3,645 (W1), 4,254 (W2) and 6,205 (W3) cell-containing droplets were sequenced to an average depth of 51,000, 39,000 and 28,000 reads per droplet.

Figure 2

Performance of demuxlet

a) Experimental design for equimolar pooling of cells from 8 unrelated samples (S1-S8) into three wells (W1-W3). W1 and W2 contain cells from two disjoint sets of 4 individuals. W3 contains cells from all 8 individuals. b) Demultiplexing single cells in each well recovers the expected individuals. c) Estimates of doublet rates versus previous estimates from mixed species experiments. d) Cell type identity determined by prediction to previously annotated PBMC data. e) t-SNE plot of two individuals (S1 and S5) from different wells are qualitatively concordant.

In wells W1, W2 and W3, demuxlet identified 91% (3332/3645), 91% (3864/4254), and 86% (5348/6205) of droplets as singlets (likelihood ratio test, L(singlet)/L(doublet) > 2), of which 25% (+/− 2.6%), 25% (+/− 4.6%) and 12.5% (+/− 1.4%) mapped to each donor, consistent with equal mixing of individuals in each well. From wells W1 and W2, each containing cells from two disjoint sets of 4 individuals, we estimated a demultiplexing error rate (number of cells assigned to individuals not in the pool) of less than 1% of singlets (W1: 2/3332, W2: 0/3864) (Fig. 2b).

We next assess the ability of demuxlet to detect doublets in both simulated and real data. 466/3645 (13%) droplets from W1 were simulated as synthetic doublets by setting the cellular barcodes of 466 cells each from individuals S1 and S2 to be the same. Applied to simulated data, demuxlet identified 91% (426/466) of synthetic doublets as doublets or ambiguous, correctly recovering the sample identity of both cells in 403/426 (95%) doublets (Supplementary Fig. 8). Applied to real data from W1, W2 and W3, demuxlet identified 138/3645, 165/4254, and 384/6205 doublets corresponding to doublet rates of 5.0%, 5.2% and 7.1%, consistent with the expected doublet rates estimated from mixed species experiments (Fig. 2c).

Demultiplexing of pooled samples allows for the statistical and visual comparisons of individual-specific dscRNA-seq profiles. Singlets identified by demuxlet in all three wells cluster into known immune cell types (Fig. 2d) and are correlated with bulk RNA-sequencing of sorted cell populations (R=0.76–0.92) (Supplementary Fig. 9). For the same individuals from different wells, t-distributed stochastic neighbor embedding (t-SNE) of dscRNA-seq data are qualitatively consistent, and estimates of cell type proportions are highly correlated (R = 0.99) (Fig. 2e and Supplementary Fig. 10). Further, t-SNE projections of the pool and each individual are not confounded by well-to-well effects (Supplementary Fig. 11a). While 6 genes were differentially expressed between wells W1 and W2 (DESeq2 on pseudobulk counts, FDR < 0.05), only 2 genes were differentially expressed between W1 and W2 individuals in well W3 (FDR < 0.05) (Supplementary Fig. 11b), suggesting multiplexing reduces technical effects due to separate sample processing[22, 23].

We used multiplexed dscRNA-seq to characterize the cell type specificity and inter-individual variability of response to IFN-β, a potent cytokine that induces genome-scale changes in the transcriptional profiles of immune cells[24, 25]. From each of 8 lupus patients, PBMCs were activated with recombinant IFN-β or left untreated for 6 hours, a time point we previously found to maximize the expression of interferon-sensitive genes (ISGs) in dendritic cells (DCs) and T cells[26, 27]. Two pools, IFN-β-treated and control, were prepared with the same number of cells from each individual and loaded onto the 10X Chromium instrument.

We obtained 14,619 (control) and 14,446 (stimulated) cell-containing droplets, of which demuxlet identified 83% (12,138) and 84% (12,167) as singlets. The estimated doublet rate of 10.9% in each condition is consistent with predicted rates (Fig. 2C) and the observed and expected frequencies of doublets for each pair of individuals are highly correlated (R=0.98) (Supplementary Fig. 12). Detected doublets form distinct clusters near the periphery of other clusters defined by cell type (Supplementary Fig. 13).

Demultiplexing individuals enables the use of the 8 individuals within each pool as biological replicates to quantitatively assess cell type-specific IFN-β responses in PBMCs. Consistent with previous reports from bulk RNA-sequencing data, IFN-β stimulation induces widespread transcriptomic changes observed as a shift in the t-SNE projections of singlets[24] (Fig. 3A). As expected, IFN-β did not affect cell type proportions between control and stimulated cells (Supplementary Fig. 14), and these were consistent with flow cytometry measurements (R=0.88) (Supplementary Fig. 15). Estimates of abundances for ~2000 homologous genes in each cell type and condition correlated with similar data from mice (Supplementary Fig. 16). We identified 3,055 differentially expressed genes (logFC > 2, FDR < 0.05) in at least one cell type (Supplementary Table 1). For 709 genes, estimates of fold change in response to IFN-β stimulation in myeloid and CD4+ cells are consistent with estimates in monocyte derived dendritic cells[28] and CD4+ T cells[27], respectively (Supplementary Fig. 17) and correlated with qPCR results of sorted CD4+ T cells (Supplementary Fig. 18). Differentially expressed genes cluster into modules of cell type-specific responses enriched for distinct gene regulatory programs (Fig. 3B, Supplementary Table 2). For example, genes upregulated in all leukocytes (Cluster III: 401 genes, logFC > 2, FDR < 0.05) or only in myeloid cells (Cluster I: 767 genes, logFC > 2, FDR < 0.05) are enriched for general antiviral response (e.g. KEGG Influenza A: Cluster III P < 1.6×10−5), chemokine signaling (Cluster I P < 7.6×10−3) and pathways active in systemic lupus erythematosus (Cluster I P < 4.4×10−3). The five clusters of downregulated genes are enriched for antibacterial response (KEGG Legionellosis: Cluster II monocyte down P < 5.5×10−3) and natural killer cell mediated toxicity (Cluster IV NK/Th cell down: P < 3.6×10−2). The analysis of multiplexed dscRNA-seq data recovers cell type-specific gene regulatory programs affected by interferon stimulation consistent with published IFN-β signatures in mouse and humans[29].

Figure 3

Inter-individual variability in IFN-β response

a) t-SNE plot of unstimulated (blue) and IFN-β-stimulated (red) PBMCs and the estimated cell types. b) Cell type-specific expression in stimulated (left) and unstimulated (right) cells. Differentially expressed genes shown (FDR < 0.05, |log(FC)| > 1). Each column represents cell type-specific expression for each individual from demuxlet. c) Observed variance (y-axis) in mean expression over all PBMCs from each of the 8 individuals versus expected variance (x-axis) over synthetic replicates sampled across all cells (light blue, pink) or replicates matched for cell type proportion (blue, red). d) Cell type proportions for each individual in unstimulated and stimulated cells. e) Correlation between sample replicates in control and stimulated cells. f) Number of significantly variable genes in each cell type and condition.

Over all PBMCs, the variance of mean expression across individuals is higher than the variance across synthetic replicates whose cells were randomly sampled (Lin’s concordance = 0.022, Pearson correlation= 0.69, Fig. 3C). The variance across synthetic replicates whose cells were sampled matching for cell type proportions is more concordant with the variance across individuals (Lin’s concordance = 0.54, Pearson correlation = 0.78, Fig. 3C–D), suggesting a contribution of cell type composition on expression variability. However, for each cell type, the variance across individuals[22, 30] is also higher than the variance across synthetic replicates (Lin’s concordance = 0.007–0.20) suggesting additional inter-individual variability not due to cell type composition (Supplementary Fig. 19). In CD14+CD16− monocytes, the correlation of mean expression between pairs of synthetic replicates from the same individual (>99%) is greater than from different individuals (~97%), further indicating inter-individual variation beyond sampling (Fig. 3E). We found between 15 to 827 genes with statistically significant inter-individual variability in control cells and 7 to 613 in stimulated cells (Pearson correlation, FDR < 0.05), with most found in classical monocytes (cM) and CD4+ helper T (Th) cells. Inter-individual variable genes in stimulated cM and to a lesser extent in Th cells (P < 9.3×10−4 and 4.5×10−2, hypergeometric test, Fig. 3F) are enriched for differentially expressed genes, consistent with our previous discovery of more IFN-β response-eQTLs in monocyte-derived dendritic cells than CD4+ T cells[26, 27]. Comparing to 407 genes previously profiled in bulk monocyte-derived dendritic cells, the proportion of variance explained by inter-individual variability is more correlated in myeloid cells after stimulation (R = 0.26 – 0.3) than before (R = 0.05 – 0.19).

To map genetic variants associated with cell type proportions and cell type-specific expression using multiplexed dscRNA-seq, we sequenced an additional 15,250 (7 donors), 22,619 (8 donors) and 25,918 cells (15 donors; 8 lupus patients, 5 rheumatoid arthritis patients, and 2 healthy controls). Demuxlet identified 71% (10,766/15,250), 73% (16,618/22,619) and 60% (15,596/25,918) of droplets as singlets, correctly assigning 99% of singlets from the first two pools, W1 and W2 (10,740/10,766 and 16,616/16,618). The estimated doublet rates of 18%, 18% and 25% are consistent with the increased concentrations of loaded cells (Fig. 2C). Similar to the IFN-β stimulation experiment, we found that expression variability was determined by variability in cell type proportion (Fig. 4A) and reproducible between batches (Supplementary Fig. 20). Associating >150,000 genetic variants (MAF > 20%) with the proportion of 8 major immune cell populations, we identified a SNP (chr10:3791224) significantly associated (P = 1.03×10−5, FDR < 0.05) with the proportion of NK cells (Fig. 4B).

Figure 4

Genetic control over cell type proportion and gene expression (N=23)

a) Observed variance (y-axis) in mean expression over all PBMCs from each individual versus expected variance (x-axis) over synthetic replicates sampled across batch 1 (left, N=8) and batch 3 (right, N=15). b) Association of chr10:3791224 with NK cell type proportions. c) Genome-wide and chromosome 6 Manhattan plots across all major cell types. Horizontal lines correspond to FDR < 0.1 (blue) and FDR < 0.05 (red). d) Q-Q plots across all genes and subsets of previously published eQTLs in relevant cell types are shown for B, cM, and Th populations. e) Notable cis-eQTLs across all major immune cell types are marked with *(FDR < 0.25), ** (FDR < 0.1), and *** (FDR < 0.05). Lack of association is marked with NS (not significant).

Across 23 donors, we conducted an expression quantitative trait loci (eQTL) analysis to map genetic variants associated with expression variability in each major immune cell type. We found a total of 32 local eQTLs (+/− 100kb, FDR < 0.1), 22 of which were detected in only one cell type (Fig. 4C, Supplementary Table 3). Previously reported local eQTLs from bulk CD14+ monocytes, CD4+ T cells and lymphoblastoid cell lines are more significantly associated with gene expression in the most similar cell types (cM, Th and B cells, respectively) than other cell types (Fig. 4D). We used an inverse variance weighted meta-analysis to identify genes with pan-cell type eQTLs, including those in the major histocompatibility complex (MHC) class I antigen presentation pathway including ERAP2 (P < 3.57×10−32, meta-analysis), encoding an aminopeptidase known to cleave viral peptides[34], and HLA-C (P < 1.74×10−29, meta-analysis), which encodes the MHC class I heavy chain (Fig. 4E). HLA-DQA1 has local eQTLs only in some cell types (P <2.11×10−15, Cochran’s Q) while HLA-DQA2 has local eQTLs in all antigen presentation cells (P < 1.02e10−43, Cochran’s Q). Among other cell type-specific local eQTLs are CD52, a gene ubiquitously expressed in leukocytes that only has eQTLs in monocyte populations, and DIP2A, a gene with an eQTL only in NK cells that is associated with immune response to vaccination in peripheral blood[35]. These results demonstrate the ability of multiplexed dscRNA-seq to characterize inter-individual variation in immune response and when integrated with genetic data, reveal cell type-specific genetic control of gene expression, which would be undetectable when bulk tissues are analyzed.

The capability to demultiplex and identify doublets using natural genetic variation reduces the per-sample and per-cell library preparation cost of single-cell RNA-sequencing, does not require synthetic barcodes or split-pool strategies[36-40], and captures biological variability among individual samples while limiting unwanted technical variability. We find the optimal number of samples to multiplex is approximately 20, based on sample processing time and empirical doublet rates of current microfluidic devices and anticipate that number to increase with automated sample handling and lower doublet rates.

Compared to sorting known cell types followed by bulk RNA-seq, multiplexed dscRNA-seq is a more efficient and unbiased method for obtaining cell type-specific immune traits[41]. Demuxlet enables reliable estimation of cell type proportion, recovers cell type-specific transcriptional response to stimulation, and could facilitate further genetic and longitudinal analyses in relevant cell types and conditions across a range of sampled individuals, including between healthy controls and disease patients[42-44]. While demuxlet could in principle be applied to sequencing solid tissue, standardizing sample processing and preservation remain major challenges. Although we developed demuxlet specifically for RNA-sequencing, we anticipate that the computational framework could be easily extended to other single cell assays where synthetic barcodes or natural genetic variation are measured by sequencing.

Methods

Identifying the sample identity of each single cell

We first describe the method to infer the sample identity of each cell in the absence of doublets. Consider RNA-sequence reads from C barcoded droplets multiplexed across S different samples, where their genotypes are available across V exonic variants. Let d be the number of unique reads overlapping with the v-th variant from the c-th droplet. Let b ∈ {R, A, O}, i ∈ {1, ⋯ , d} be the variant-overlapping base call from the i-th read, representing reference (R), alternate (A), and other (O) alleles respectively. Let e ∈ {0,1} be a latent variable indicating whether the base call is correct (0) or not (1), then given e = 0, b ∈ {R = 0, A = 1} and ~ Binomial when g ∈ {0,1,2} is the true genotype of sample corresponding to c-th droplet at v-th variant. When e = 1, we assume that Pr(b|g, e) follows Supplementary Table 4. e is assumed to follow Bernoulli where q is a phred-scale quality score of the observed base call. We use the standard 10X pipeline to process the raw reads which estimates the phred-scale quality score based on the alignment of each read to the reference human transcriptome using the STAR aligner[49].

We allow uncertainty of observed genotypes at the v-th variant for the s-th sample using , the posterior probability of a possible genotype g given external DNA data Data (e.g. sequence reads, imputed genotypes, or array-based genotypes). If genotype likelihood Pr(Data|g) is provided (e.g. unphased sequence reads) instead, it can be converted to a posterior probability scale using where Pr(g) ~ Binomial(2, p) and p is the population allele frequency of the alternate allele. To allow errors ε in the posterior probability, we replace it with . The overall likelihood that the c-th droplet originated from the s-th sample is

In the absence of doublets, we use the maximum likelihood to determine the best-matching sample as argmax[L (s)].

Screening for droplets containing multiple samples

To identify doublets, we implement a mixture model to calculate the likelihood that the sequence reads originated from two individuals, and the likelihoods are compared to determine whether a droplet contains cells from one or two samples. If sequence reads from the c-th droplet originate from two different samples, s1, s2 with mixing proportions (1 − α) : α, then the likelihood in (1) can be represented as the following mixture distribution[18],

To reduce the computational cost, we consider discrete values of α ∈ {α1, ⋯ , α}, (e.g. 5 – 50% by 5%). We determine that it is a doublet between samples s1, s2 if and only if and the most likely mixing proportion is estimated to be argmaxL(s1, s2, α). We determine that the cell contains only a single individual s if , and less confident droplets are classified as ambiguous. While we consider only doublets for estimating doublet rates, we remove all doublets and ambiguous droplets to conservatively estimate singlets. Supplementary Fig. 8 illustrates the distribution of singlet, doublet likelihoods and the decision boundaries when t = 2 was used.

Theoretical expectation of deconvoluting singlets

The theoretical distribution of expected singlets with multiplexing (presented in Supplementary Fig. 2) is as follows. Let d (e.g. 0.01) be the proportion of true multiplets when x (1,000) cells are loaded when multiplexing was not used. Then the expected multiplet rates when x cells are loaded can be modeled exponentially as . Let α be the fraction of true singlets incorrectly classified as non-singlets (i.e. doublet or ambiguous), and β be the fraction of multiplets correctly classified as non-singlets. When multiplexing x cells equally from n samples, the expected multiplet rates are d(x), and are expected to be undetectable doublets mixed between the cells from the same sample. Therefore, the overall effective multiplet rate is . Similarly, the expected number of correctly identified singlets becomes . Given α, β the expected number of singlets can be calculated by fixing the multiplet rate d(x) = d0. We used d0 = 0.01, x0 = 1000 for the simulation in Supplementary Fig. 2.

Dependence of demultiplexing performance on experimental design parameters

The demuxlet ‘plp’ option was used to generate a pileup format of 6,145 cells from one well of PBMC 10x data. The reads in the pileup were then modified to reflect the genotypes of individuals sampled from the 1000 Genomes Phase 3 cohort. The pileup was downsampled to obtain different numbers of read-overlapping exonic SNPs (ranging from 5,000 to 100,000) for the whole cohort. To create simulated doublets, we randomly sampled and merged pairs of barcodes within a dataset, resulting in a 5% doublet rate in the original data. For simulations with related individuals, we simulated transcriptomes from 8 individuals in 1000 Genomes with varying degrees of relatedness, ranging from unrelated to parent-child (HG00146, HG00147, HG00500, HG00501, HG00502, HG00512, HG00514, and HG00524).

Isolation and preparation of PBMC samples

Informed consent was obtained from all patients sequenced in this study. Peripheral blood mononuclear cells were isolated from patient donors, Ficoll separated, and cryopreserved by the UCSF Core Immunologic Laboratory (CIL). PBMCs were thawed in a 37°C water bath, and subsequently washed and resuspended in EasySep buffer (STEMCELL Technologies). Cells were treated with DNAseI and incubated for 15 min at RT before filtering through a 40um column. Finally, the cells were washed in EasySep and resuspended in 1x PBMS and 0.04% bovine serum albumin. Cells from 8 donors were then re-concentrated to 1M cells per mL and then serially pooled. At each pooling stage, 1M cells per mL were combined to result in a final sample pool with cells from all donors.

IFN-β stimulation and culture

Prior to pooling, samples from 8 individuals were separated into two aliquots each. One aliquot of PBMCs was activated by 100 U/mL of recombinant IFN-β (PBL Assay Science) for 6 hrs according to the published protocol[26]. The second aliquot was left untreated. After 6 hrs, the 8 samples for each condition were pooled together in two final pools (stimulated cells and control cells) as described above.

Fluorescence-activated cell sorting and analysis

1M PBMCs from each donor were stained using standard procedure (30 min, 4 C) with the following surface antibody panel (CD3-PerCP clone SK7 (BioLegend), CD4-APC clone OKT4 (BioLegend), CD8-BV570 clone RPA-T8 (BioLegend), CD14-FITC clone 63D3 (BioLegend), CD19-BV510 clone SJ25C1 (BD), and Ghost dye A710 viability stain (Tonbo)) (Life Sciences Reporting Summary). Samples were then analyzed and sorted using a BD FACSAria Fusion instrument at the UCSF flow cytometry core. To calculate cell type proportions, the number of events in each of CD3+ CD4+ CD8− (CD4+ T cells), CD3+ CD4− CD8+ (CD8+ T cells), CD3− CD19+ (B cells), and CD3− CD14+ (monocytes) were divided by the sum of events in these gates (Supplementary Fig. 21).

Quantitative polymerase chain reaction analysis

RNA was isolated from sorted CD4+ T cells following the RNeasy micro kit protocol (QIAGEN), and cDNA was prepared using MultiScribe Reverse Transcriptase (Applied Biosystems cat #4368814). The qPCR primers were chosen from the PrimerBank reference when available [50]. Each sample was run in triplicate with the Luminaris HiGreen qPCR kit (Thermo Scientific #K0992) according to standard protocol using a Roche Light Cycler 96 instrument and fold change was calculated from ΔΔCT between control and stimulated samples with GAPDH as a reference gene.

Droplet-based capture and sequencing

Cellular suspensions were loaded onto the 10x Chromium instrument (10x Genomics) and sequenced as described in Zheng et al[17]. The cDNA libraries were sequenced using a custom program on 10 lanes of Illumina HiSeq2500 Rapid Mode, yielding 1.8B total reads and 25K reads per cell. At these depths, we recovered >90% of captured transcripts in each sequencing experiment.

Bulk isolation and sequencing

PBMCs from lupus patients were isolated and prepared as described above. Once resuspended in EasySep buffer, the EasyEights Magnet was used to sequentially isolate CD14+ (using the EasySep Human CD14 positive selection kit II, cat #17858), CD19+ (using the EasySep Human CD19 positive selection kit II, cat #17854), CD8+ (EasySep Human CD8 positive selection kitII, cat#17853), and CD4+ cells (EasySep Human CD4 T cell negative isolation kit (cat #17952) according to the kit protocol. RNA was extracted using the RNeasy Mini kit (#74104), and reverse transcription and tagmentation were conducted according to Picelli et al. using the SmartSeq2 protocol[51, 52]. After cDNA synthesis and tagmentation, the library was amplified with the Nextera XT DNA Sample Preparation Kit (#FC-131-1096) according to protocol, starting with 0.2ng of cDNA. Samples were then sequenced on one lane of the Illumina Hiseq4000 with paired end 100bp read length, yielding 350M total reads.

Alignment and initial processing of single cell sequencing data

We used the CellRanger v1.1 and v1.2 software with the default settings to process the raw FASTQ files, align the sequencing reads to the hg19 transcriptome, and generate a filtered UMI expression profile for each cell[17]. The raw UMI counts from all cells and genes with nonzero counts across the population of cells were used to generate t-SNE profiles.

Cell type classification and clustering

To identify known immune cell populations in PBMCs, we used the Seurat package to perform unbiased clustering on the 2.7k PBMCs from Zheng et al., following the publicly available Guided Clustering Tutorial[17, 53]. The FindAllMarkers function was then used to find the top 20 markers for each of the 8 identified cell types. Cluster averages were calculated by taking the average raw count across all cells of each cell type. For each cell, we calculated the Spearman correlation of the raw counts of the marker genes and the cluster averages, and assigned each cell to the cell type to which it had maximum correlation.

Differential expression analysis

Demultiplexed individuals were used as replicates for differential expression analysis. For each gene, raw counts were summed for each individual. We used the DESeq2 package to detect differentially expressed genes between control and stimulated conditions[54]. Genes with baseMean > 1 were filtered out from the DESeq2 output, and the qvalue package was used to calculate FDR < 0.05 [55].

Estimation of inter-individual variability in PBMCs

For each individual, we found the mean expression of each gene with nonzero counts. The mean was calculated from the log2 single cell UMI counts normalized to the median count for each cell. To measure inter-individual variability, we then calculated the variance of the mean expression across all individuals. Lin’s concordance correlation coefficient was used to compare the agreement of observed data and synthetic replicates. Synthetic replicates were generated by sampling without replacement either from all cells or cells matched for cell type proportion. Cell type-specific variability estimated as the correlation between synthetic replicates was compared to variability estimates from 23 biological replicates of bulk IFN-stimulated monocyte-derived dendritic cells. Protein coding genes (407/414) originally measured using Nanostring (a hybridization based PCR-free quantification method) were assessed, and variability in the bulk dataset was estimated as repeatability using a linear mixed model[56,26].

Estimation of inter-individual variability within cell types

For each cell type, we generated two bulk equivalent replicates for each individual by summing raw counts of cells sampled without replacement. We used DESeq2 to generate variance-stabilized counts across all replicates. To filter for expressed genes, we performed all subsequent analyses on genes with 5% of samples with > 0 counts. The correlation of replicates was performed on the log2 normalized counts. Pearson correlation of the two replicates from each of the 8 individuals was used to find genes with significant inter-individual variability.

Quantitative trait mapping in major immune cell types

Genotypes were imputed with EAGLE[57] and filtered for MAF > 0.2, resulting in a total of 189,322 SNPs. Cell type proportions were calculated as number of cells for each cell type divided by the number of total cells for each person. Linear regression was used to test associations between each genetic variant and cell-type proportion with the Matrix eQTL software[58]. Cis-eQTL mapping was conducted in each cell type separately. All genes with at least 50 UMI counts in 20% of the individuals in all PBMCs were tested for each cell type, resulting in a total of 4,555 genes. Variance-stabilized and log-normalized gene expression was calculated using the ‘rlog’ function of the DESeq2 package[54]. All variants within a window of 100kbp of each gene were tested with linear regression using Matrix eQTL[58]. Batch information for each sample as well as the first 3 principal components of the expression matrix were used as covariates.