Massively parallel digital transcriptional

Full text

Turn on search term navigation

ARTICLE

Received 20 Sep 2016 | Accepted 23 Nov 2016 | Published 16 Jan 2017

Grace X.Y. Zheng1, Jessica M. Terry1, Phillip Belgrader1, Paul Ryvkin1, Zachary W. Bent1, Ryan Wilson1, Solongo B. Ziraldo1, Tobias D. Wheeler1, Geoff P. McDermott1, Junjie Zhu1, Mark T. Gregory2, Joe Shuga1, Luz Montesclaros1, Jason G. Underwood1,3, Donald A. Masquelier1, Stefanie Y. Nishimura1,Michael Schnall-Levin1, Paul W. Wyatt1, Christopher M. Hindson1, Rajiv Bharadwaj1, Alexander Wong1, Kevin D. Ness1, Lan W. Beppu4, H. Joachim Deeg4, Christopher McFarland5, Keith R. Loeb4,6,William J. Valente2,7,8, Nolan G. Ericson2, Emily A. Stevens4, Jerald P. Radich4, Tarjei S. Mikkelsen1, Benjamin J. Hindson1 & Jason H. Bielas2,6,8,9

Characterizing the transcriptome of individual cells is fundamental to understanding complex biological systems. We describe a droplet-based system that enables 30 mRNA counting of tens of thousands of single cells per sample. Cell encapsulation, of up to 8 samples at a time, takes place in B6 min, with B50% cell capture efciency. To demonstrate the systems technical performance, we collected transcriptome data from B250k single cells across 29 samples. We validated the sensitivity of the system and its ability to detect rare populations using cell lines and synthetic RNAs. We proled 68k peripheral blood mononuclear cells to demonstrate the systems ability to characterize large immune populations. Finally, we used sequence variation in the transcriptome data to determine host and donor chimerism at single-cell resolution from bone marrow mononuclear cells isolated from transplant patients.

1 10x Genomics Inc., Pleasanton, California, 94566, USA. 2 Translational Research Program, Public Health Sciences Division, Fred Hutchinson Cancer Research Center, Seattle, Washington 98109, USA. 3 Department of Genome Sciences, University of Washington, Seattle, Washington 98195, USA. 4 Clinical Research Division, Fred Hutchinson Cancer Research Center, Seattle, Washington 98109, USA. 5 Seattle Cancer Care Alliance Clinical Immunogenetics Laboratory, Seattle, Washington 98109, USA. 6 Department of Pathology, University of Washington, Seattle, Washington 98195, USA. 7 Medical Scientist Training Program, University of Washington School of Medicine, Seattle, Washington 98195, USA. 8 Molecular and Cellular Biology Graduate Program, University of Washington, Seattle, Washington 98195, USA. 9 Human Biology Division, Fred Hutchinson Cancer Research Center, Seattle, Washington 98109, USA. Correspondence and requests for materials should be addressed to B.J.H. (email: mailto:[email protected]

Web End [email protected] ) or to J.H.B. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

DOI: 10.1038/ncomms14049 OPEN

Massively parallel digital transcriptional proling of single cells

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049

Understanding of biological systems requires the knowledge of their individual components. Single-cell RNA-sequencing (scRNA-seq) can be used to dissect transcriptomic

heterogeneity that is masked in population-averaged measurements1,2. scRNA-seq studies have led to the discovery of novel cell types and provided insights into regulatory networks during development3. However, previously described scRNA-seq methods face practical challenges when scaling to tens of thousands of cells or when it is necessary to capture as many cells as possible from a limited sample49. Commercially available, microuidic-based approaches have limited throughput5,6. Plate-based methods often require time-consuming uorescence-activated cell sorting (FACS) into many plates that must be processed separately4,9. Droplet-based techniques have enabled processing of tens of thousands of cells in a single experiment7,8, but current approaches require generation of custom microuidic devices and reagents.

To overcome these challenges, we developed a droplet-based system that enables 30 messenger RNA (mRNA) digital counting of thousands of single cells. Approximately 50% of cells loaded into the system can be captured, and up to eight samples can be processed in parallel per run. Reverse transcription takes place inside each droplet, and barcoded complementary DNAs (cDNAs) are amplied in bulk. The resulting libraries then undergo Illumina short-read sequencing. An analysis pipeline, Cell Ranger, processes the sequencing data and enables automated cell clustering. Here we rst demonstrated comparable sensitivity of the system to existing droplet-based methods by performing scRNA-seq on cell lines and synthetic RNAs. Next, we proled 68k fresh peripheral blood mononuclear cells (PBMCs) and demonstrated the scRNA-seq platforms ability to dissect large immune populations. Last, we developed a computational method to distinguish donor from host cells in bone marrow transplant samples by genotype. We combined this method with clustering analysis to compare subpopulation changes in acute myeloid leukemia (AML) patients. This analysis enables transplant monitoring of the complex interplay between donor and host cells.

ResultsDroplet-based platform enables barcoding of cells. The scRNA-seq microuidics platform builds on the GemCode technology, which has been used for genome haplotyping, structural variant analysis and de novo assembly of a human genome1012. The core of the technology is a Gel bead in EMulsion (GEM). GEM generation takes place in an 8-channel microuidic chip that encapsulates single gel beads at B80% ll rate (Fig. 1ac). Each gel bead is functionalized with barcoded oligonucleotides that consists of: (i) sequencing adapters and primers, (ii) a 14 bp barcode drawn from B750,000 designed sequences to index GEMs, (iii) a 10 bp randomer to index molecules (unique molecular identier, UMI) and (iv) an anchored 30 bp oligo-dT to prime polyadenylated RNA transcripts (Fig. 1d). Within each microuidic channel, B100,000 GEMs are formed per B6-min run, encapsulating thousands of cells in GEMs. Cells are loaded at a limiting dilution to minimize co-occurrence of multiple cells in the same GEM.

Cell lysis begins immediately after encapsulation. Gel beads dissolve and release their oligonucleotides for reverse transcription of polyadenylated RNAs. Each resulting cDNA molecule contains a UMI and shared barcode per GEM, and ends with a template switching oligo at the 30 end (Fig. 1e). Next, the emulsion is broken and barcoded cDNA is pooled for PCR amplication, using primers complementary to the

switch oligos and sequencing adapters. Finally, amplied cDNAs are sheared, and adapter and sample indices are incorporated into nished libraries, which are compatible with next-generation short-read sequencing. Read1 contains the cDNA insert while Read2 captures the UMI. Index reads, I5 and I7, contain the sample indices and cell barcodes, respectively. This streamlined approach enables parallel capture of thousands of cells in each of the 8 channels for scRNA-seq analysis.

Technical demonstration with cell lines and synthetic RNAs. To assess the technical performance of our system, we loaded a mixture of B1,200 human (293T) and B1,200 mouse (3T3)

cells and sequenced the library on the Illumina NextSeq 500 to yield B100k reads per cell. Sequencing data were processed by CellRanger (Supplementary Methods and Fig. 1f). Briey, 98 nucleotides (nt) of Read1s were aligned against the union of human (hg19) and mouse (mm10) genomes with STAR. Barcodes and UMIs were ltered and corrected (Supplementary Methods). PCR duplicates were marked using the barcode, UMI and gene ID. Only condently mapped, non-PCR duplicates with valid barcodes and UMIs were used to generate a gene-barcode matrix for further analysis. Thirty-eight per cent and 33% of reads mapped to human and mouse exonic regions, respectively, and o6% of reads mapped to intronic regions (Supplementary

Table 1). The high mapping rate is comparable to previously reported scRNA-seq systems49.

Based on the distribution of total UMI counts for each barcode (Supplementary Methods), we estimated that 1,012 GEMs contained cells, of which 482 and 538 contained reads that mapped primarily to the human and mouse transcriptome, respectively (and will be referred to as human and mouse GEMs) (Fig. 2a). Greater than eighty-three per cent of UMI counts were associated with cell barcodes, indicating low background of cell-free RNA. Eight cell-containing GEMs had a substantial fraction of human and mouse UMI counts (the UMI count is 41% of each species UMI count distribution), yielding an inferred multiplet rate (rate of GEMs containing 41 cell) of1.6% (Supplementary Methods and Fig. 2a). A cell titration experiment across six different cell loads showed a linear relationship between the multiplet rate and the number of recovered cells ranging from 1,200 to 9,500 (Supplementary Fig. 1a). The multiplet rate and trend are consistent with Poisson loading of cells, and have been validated by independent imaging experiments (Supplementary Methods and Supplementary Fig. 1b). In addition, we observed B50% cell capture rate, which is the ratio of the number of cells detected by sequencing and the number of cells loaded. The capture rate is consistent across four types of cells with cell loads ranging from B1,000 to B23,000 (Supplementary Table 2), a key improvement over some existing scRNA-seq methods57. Last, the mean fraction of UMI counts from the other species was 0.9% in both human and mouse GEMs, indicating a low level of cross-talk between cell barcodes (Online Supplementary Methods). This, coupled with the low multiplet rate and high cell capture rate, is particularly important when processing samples that are in extreme limited supply.

The sensitivity of scRNA-seq methods is critical to many applications. At 100k reads per cell, we detected a median of B4,500 genes and B27,000 transcripts (UMI counts) in each human and mouse cell (Fig. 2b,c). UMI counts showed a standard deviation of B43% of the mean (coefcient of variation (CV)) in human cells, and B33% of the mean in mouse cells, where the trend is consistent in four independent human and mouse mixture experiments (Supplementary Fig. 1c,d). Genes of different guanine-cytosine (GC) composition and length show

2 NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049 ARTICLE

Eight-channel microfluidics chip

Barcoded cDNA

Oil

Cells + reagents

Beads RT

Construct Library

Break Emulsion

AmplifycDNA Sequence

GEM outlet

Single-cell GEMs

f Extract cell-barcode, UMI, RNA read Correct cell-barcode sequences

Align reads using STAR

Tag reads with gene, transcript hits

Correct UMI sequences

Count UMIs by (cell, gene)

Select cell-associated barcodes

Gene-barcode matrix

1 Hamming distance

Barcoded primer gel beads

Cells reagents

Oil

Uniquely mapped reads only

Barcodes (T)30VN

poly(A)

cDNA

UMI

P7 R2

% of GEMs

RNA

P7 R2 P5

10 Barcodes

UMI (T)30VN

cDNA

Insert

Sample Index

0 N=0 N=1 N>1

No. of gel beads per GEM

Figure 1 | GemCode single-cell technology enables 30 proling of RNAs from thousands of single cells simultaneously. (a) scRNA-seq workow on

GemCode technology platform. Cells were combined with reagents in one channel of a microuidic chip, and gel beads from another channel to form

GEMs. RT takes place inside each GEM, after which cDNAs are pooled for amplication and library construction in bulk. (b) Gel beads loaded with primers and barcoded oligonucleotides are rst mixed with cells and reagents, and subsequently mixed with oil-surfactant solution at a microuidic junction. Single-cell GEMs are collected in the GEM outlet. (c) Percentage of GEMs containing 0 gel bead (N 0), 1 gel bead (N 1) and 41 gel bead (N41). Data

include ve independent runs from multiple chip and gel bead lots over 470k GEMs for each run, n 5, means.e.m. (d) Gel beads contain barcoded

oligonucleotides consisting of Illumina adapters, 10x barcodes, UMIs and oligo dTs, which prime RT of polyadenylated RNAs. (e) Finished library molecules consist of Illumina adapters and sample indices, allowing pooling and sequencing of multiple libraries on a next-generation short read sequencer.(f) CellRanger pipeline workow. Gene-barcode matrix (highlighted in green) is an output of the pipeline.

similar UMI count distributions, suggesting low transcript bias (Supplementary Fig. 1eh).

We also directly measured cDNA conversion rate by loading External RNA Controls Consortium (ERCC) synthetic RNAs into GEMs in place of cells. We found that mean UMI counts from sequencing was highly correlated (r 0.96)

with molecule counts calculated from the loading concentration of ERCC (Fig. 2d and Supplementary Fig. 2a). Furthermore, we inferred 6.78.1% efciency from both ERCC RNA Spike-in Mix1 and Mix2 in a 1:50 dilution (Supplementary Fig. 2b), with minimal evidence of GC bias, and limited bias for transcripts longer than 500 nt (Supplementary Fig. 2c,d). Additionally, we estimated the conversion rate of cell transcripts in Jurkat cells by droplet digital PCR (ddPCR)13. The amount of cDNA of eight genes obtained from single cells after reverse transcription in GEMs was compared with the expected RNA inferred from bulk proling. The conversion rates among eight genes were between2.5 and 25.5%, which is consistent with the ERCC data (Supplementary Methods and Supplementary Fig. 2e).

The ERCC experiments also allowed us to estimate the relative proportion of biological and technical variation. Since ERCCs are in solution, they do not introduce biological variation related to differences in cell size, RNA content or transcriptional activity. Thus, technical variation is the only source of variation. When the ERCCs are dilute (UMI counts are small), sampling noise dominates; when the UMI counts increase, technical variations become dominant14 (Supplementary Fig. 2f). These

variations include variation in droplet size, variation in concentration of reverse transcription reagents in the droplets, variation in the concentration of sample in the droplets and variation in RT and/or PCR efciency of the distinct gel bead barcode sequences. The squared CV (CV2) is B7% among all the ERCC experiments. In comparison, CV2 in samples of mouse and human cells is B1119% (Supplementary Fig. 1d), suggesting that technical variance accounts for B50% of total variance, consistent with the observations from Klein et al.8

Detection of individual populations in mixed samples. We tested the ability of the system to accurately detect heterogeneous populations by mixing two cell lines, 293T and Jurkat cells, at different ratios (Supplementary Table 1). We performed principal component analysis (PCA) on UMI counts from all detected genes after pooling all the samples (Supplementary Fig. 3a). In the sample where an equal number of 293T and Jurkat cells was mixed, PC 1 separated cells into two clusters of equal size (Fig. 2e, Supplementary Fig. 4a and Supplementary Data 1). Based on the expression of cell type-specic markers, we infer that one cluster corresponds to Jurkat cells (preferentially expressing CD3D), and the other corresponds to 293T cells (preferentially expressing XIST, as 293T is a female cell line, and Jurkat is a male cell line) (Fig. 2e and Supplementary Fig. 4b). Points located between the two clusters are likely multiplets, as they expressed both CD3D and XIST (Fig. 2e and

NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049

Mouse-only

Human-only

Human:mouse

Mouse (3T3)

Mouse UMI counts (k)

Median genes per cell (k)

Median UMIs per cell (k)

Mean UMI counts (log 10)

1,012 cells recovered

0 20 40 60 80

Human UMI counts (k)

50 75

Raw reads per cell (k)

100

50 75

100

0 1 2 3 4

Raw reads per cell (k)

ERCC molecules per GEM (log10)

293T-only Jurkat-only

293T:Jurkat (50:50) 293T:Jurkat (1:99)

293T-only Jurkat-only

150

293T:Jurkat

293T-only

100 Jurkat-only

3,396 cells recovered Multiplet rate=3.1%

293T SNV counts293T SNV counts

100

Principal component 3

CD3D

0 50 100 150

150

293T:Jurkat (50:50)

150

293T:Jurkat (1:99)

100

Principal component 1

0 50 100 150 Jurkat SNV counts

0 50 100 150

Jurkat SNV counts

Figure 2 | Demonstration of technical performance of GemCode single-cell technology with cell lines and ERCC. (a) Scatter plot of human and mouse UMI counts detected in a mixture of 293T and 3T3 cells. Cell barcodes containing primarily mouse reads are colored in cyan and termed Mouse only; cell barcodes with primarily human reads are colored in red and termed Human only; and cell barcodes with signicant mouse and human reads are coloured in grey and termed Human:Mouse. A multiplet rate of 1.6% was inferred. Median number of genes (b) and UMI counts (c) detected per cell in a mixture of 293T (red) and 3T3 (cyan) cells at different raw reads per cell. Data from three independent experiments were included, means.e.m. (d) Mean observed UMI counts for each ERCC molecule is compared with expected number of ERCC molecules per GEM. A straight line was tted to summarize the relationship. (e) Principal component analysis was performed on normalized scRNA-seq data of Jurkat and 293T cells mixed at four different ratios (100% 293T, 100% Jurkat, 50:50 293T:Jurkat and 1:99 293Tand Jurkat). PC1 and PC3 are plotted, and each cell is colored by the normalized expression of CD3D. (f) SNV analysis was performed, and 293T- and Jurkat-enriched SNVs were plotted for each sample. A 3.1% multiplet rate was inferred from the 50:50 293T: Jurkat sample.

Supplementary Fig. 4b). In contrast, PC1 did not separate cells into two clusters in the 293T- and the Jurkat-only samples (Fig. 2e). Furthermore, in the sample with 1% 293T and 99% Jurkat cells, the number of cells in each of the two clusters was observed at the expected ratio (Fig. 2e, Supplementary Fig. 4a and Supplementary Fig. 4b). A similar trend was observed for 12 independent samples where 293T and Jurkat cells were mixed at ve different proportions, demonstrating the systems ability to perform unbiased detection of rare single cells (Supplementary Fig. 4a).

Our scRNA-seq data not only provides a digital transcript count but it also provides B250 nt sequence for each cDNA that could be used for single nucleotide variant (SNV) detection. On average, there are B350 SNVs detected in each 293T or

Jurkat cell (Supplementary Fig. 4c and Supplementary Table 3), and we investigated whether the SNVs could be used independently to distinguish cells in the mixture. We selected a set of high-quality SNVs that were only observed in 293T or Jurkat cells, but not both (Supplementary Methods).

We then scored cells in the mixed samples based on the number of 293T or Jurkat-enriched SNVs (Supplementary Methods). In the 1:1 mixed sample, 45% 293T cells primarily (96%) harbored 293T-enriched SNVs, whereas 50% Jurkat cells primarily (94%) harbored Jurkat-enriched SNVs (Fig. 2f and Supplementary Data 2). Jurkat and 293T cells inferred from marker-based analysis are 99% consistent with SNV-based assignment. We observed a multiplet rate of B3%, accounting for multiplets from Jurkats:293Ts as well as Jurkats:Jurkats and 293Ts:293Ts. The multiplet rate is consistent with that predicted from the human and mouse mixing experiment, when B3,000 cells were recovered (Supplementary Fig. 1a). Our result demonstrates that SNVs detected from scRNA-seq data can be used to classify individual cells.

Subpopulation discovery from a large immune population. The GemCode single-cell technology can also be used for scRNA-seq

4 NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049 ARTICLE

a b j

10,000

10 (0.5%)

1(9.3%)

3(17.0%)

4 (11.2%)

7 (5.7%)

9 (6.6%)

6 (12.5%)

2(31.2%)

5 (5.7%)

8(0.3%)

CD19+

B cells

CD4+/CD25+ Reg T CD4+/CD45RO+ memory T

CD8+/CD45 RA+ Naive cytotoxic

Dendritic

CD4+/CD25+ Reg T

CD56+ NK

CD8+ cytotoxic T

Megakaryocytes

CD14+ Monocytes

CD4+/CD45 RA+/ CD25- Naive T CD8+/CD45 RA+ Naive cytotoxic

Counts

1,000

tSNE2

Genes/Cell

UMIs Counts/Cell

tSNE1

d e

CD3D

CD8A

NKG7

Clusters

1 2

1: Activated CD8+

8: Megakaryocytes

6: CD8+

10 4 3 GZMK

CD8A

PTCRA

TNFRSF18

tSNE2

SIGLEC7

GNLY

5: NK

9: Monocytes and dendritic

10: B, dendritic, T

LGALS3

CD79A

hFCER1A

CD16

S100A8

3: Memory and Reg T

CCR10 CD4

CLEC4C

CD8B ID3

1 (37%) 3(16%)

2 (47%)

1 (37%) 3(16%)

2 (47%)

1 (37%) 3(16%)

2 (47%)

PF4

tSNE2

2: Naive CD8+

4: Naive CD4+

7: B

Normalized expression

tSNE1

Figure 3 | Distinct populations can be detected in fresh 68k PBMCs. (a) Distribution of number of genes (left) and UMI counts (right) detected per 68k PBMCs. (b) tSNE projection of 68k PBMCs, where each cell is grouped into one of the 10 clusters (distinguished by their colours). Cluster number is indicated, with the percentage of cells in each cluster noted within parentheses. (c) Normalized expression (centred) of the top variable genes (rows) from each of 10 clusters (columns) is shown in a heatmap. Numbers at the top indicate cluster number in (b), with connecting lines indicating the hierarchical relationship between clusters. Representative markers from each cluster are shown on the right, and an inferred cluster assignment is shown on the left. (di) tSNE projection of 68k PBMCs, with each cell coloured based on their normalized expression of CD3D, CD8A, NKG7, FCER1A, CD16 and S100A8. UMI normalization was performed by rst dividing UMI counts by the total UMI counts in each cell, followed by multiplication with the median of the total UMI counts across cells. Then, we took the natural log of the UMI counts. Finally, each gene was normalized such that the mean signal for each gene is 0, and standard deviation is 1. (j) tSNE projection of 68k PBMCs, with each cell coloured based on their correlation-based assignment to a puried subpopulation of PBMCs. Subclusters within T cells are marked by dashed polygons. NK, natural killer cells; reg T, regulatory T cells.

of primary cells. To study immune populations within PBMCs, we obtained fresh PBMCs from a healthy donor (Donor A). 89k cells were captured from each of 8 channels and pooled to obtain B68k cells. Data from multiple sequencing runs were merged using the Cell Ranger pipeline. At B20k reads per cell, the median number of genes and UMI counts detected per cell was B525 and 1,300, respectively (Fig. 3a and

Supplementary Fig. 5a). The UMI count is roughly 10% of that from 293T and 3T3 samples at B20k reads per cell, likely reecting the differences in cells RNA content (B1 pg RNA per cell in PBMCs versus B15 pg RNA per cell in 293T and 3T3 cells) (Supplementary Fig. 5a,b).

We performed clustering analysis to examine cellular heterogeneity among PBMCs. We applied PCA on the top 1,000 variable genes ranked by their normalized dispersion, following a similar approach to Macosko et al.7 (Supplementary Figs 3b and 5c and Supplementary Methods). K-means15 clustering on the rst 50 PCs identied 10 distinct cell clusters, which were visualized in two-dimensional projection of t-distributed stochastic neighbour embedding (tSNE)16 (Supplementary Methods, Fig. 3b and Supplementary Fig. 5d).

To identify cluster-specic genes, we calculated the expression difference of each gene between that cluster and average of the rest of clusters. Examination of the top cluster-specic genes revealed major subtypes of PBMCs at expected ratios17: 480% T cells (enrichment of CD3D, part of the T-cell receptor complex, in clusters 13 and 6), B6% NK cells (enrichment of NKG7 (ref. 18) in cluster 5), B6% B cells (enrichment of CD79A (ref. 19) in cluster 7) and B7% myeloid cells (enrichment of S100A8 and S100A9 (ref. 20) in cluster 9 (Supplementary Methods, Fig. 3bf, Supplementary Fig. 5e and Supplementary Data 3). Finer substructures were detected within the T-cell cluster; clusters 1, 4 and 6 are CD8 cytotoxic T cells, whereas

clusters 2 and 3 are CD4 T cells (Fig. 3e and Supplementary

Fig. 5f). The enrichment of NKG7 on cluster 1 cells implies a cluster of activated cytotoxic T cells21 (Fig. 3f). Cells in cluster 3 showed high expression of CCR10 and TNFRSF18, markers for memory T cells22 and regulatory T cells23 respectively, and likely consisted of a mixture of memory and regulatory T cells (Fig. 3c and Supplementary Fig. 5g). The presence of ID3, which is important in maintaining a naive T-cell state24, suggests that cluster 2 represents naive CD8 T cells, whereas

NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049

cluster 4 represents naive CD4 T cells (Fig. 3c). To identify subpopulations within the myeloid population, we further applied k-means clustering on the rst 50 PCs of cluster 9 cells. At least three populations were evident: dendritic cells (characterized by the presence of FCER1A25), CD16 mono

cytes and CD16 /low monocytes26 (Fig. 3gi and Suppleme

ntary Data 3). Overall, these results demonstrate that our scRNA-seq method can detect all major subpopulations expected to be present a PBMC sample.

Our analysis also revealed some minor cell clusters, such as cluster 8 (0.3%) and cluster 10 (0.5%) (Fig. 3b). Cluster 8 showed preferential expression of megakaryocyte markers, such as PF4, suggesting that it represents a cluster of megakaryocytes (Fig. 3b,c and Supplementary Fig. 5h). Cells in cluster 10 express markers of B, T and dendritic cells, suggesting a likely cluster of multiplets (Fig. 3b,c). The size of the cluster suggests the multiplets comprises mostly B:dendritic and B:T:dendritic cells (Supplementary Methods). With B9k cells recovered per channel, we expect a B9% multiplet rate and that the majority of multiplets would only contain T cells. More sophisticated methods will be required to detect multiplets from identical or highly similar cell types.

To further characterize the heterogeneity among 68k PBMCs, we generated reference transcriptome proles through scRNA-seq of 10 bead-enriched subpopulations of PBMCs from Donor A (Supplementary Figs 6 and 7 and Supplementary Table 4). Clustering analysis revealed a lack of substructure in most samples, consistent with the samples being homogenous populations, and in agreement with FACS analysis (Supplementary Methods and Supplementary Figs 6 and 7). However, substructures were observed in CD34 and CD14 monocyte samples

(Supplementary Methods and Supplementary Fig. 7b,j). In the CD34 sample, B70% cell clusters show expression of

CD34 (Supplementary Fig. 7j). In the CD14 sample, the minor

population showed marker expression for dendritic cells (for example, CLEC9A), providing another reference transcriptome to classify the 68k PBMCs (Supplementary Fig. 7b). This result also demonstrates the power of scRNA-seq in selecting appropriate cells for further analysis.

We classied 68k PBMCs based on their best match to the average expression prole of 11 reference transcriptomes (Supplementary Methods and Fig. 3j). Cell classication was largely consistent with previously described marker-based classication, although the boundaries among some of the T-cell subpopulations were blurred. Namely, part of the inferred CD4 naive T population was classied as CD8 T cells. We

have also tried to cluster the 68k PBMC data with Seurat27 (Supplementary Methods). While it was able to distinguish inferred CD4 naive from inferred CD8 naive T cells, it

was not able to cleanly separate out inferred activated cytotoxic T cells from inferred NK cells (Supplementary Fig. 5i). Such populations have overlapping functions, making separation at the transcriptome level particularly difcult and even unexpected. However, the complementary results of Seurats and our analysis suggest that more sophisticated clustering and classication methods can help address these problems.

Single-cell RNA proling of cryopreserved PBMCs. To determine the effect that a freeze-thaw might have on gene expression and thus on the ability of our scRNA-seq pipeline to classify cell type in frozen repository specimens, we froze the remaining fresh PBMCs from Donor A, and made a scRNA-seq library from gently thawed cells 3 weeks later where B3k cells were recovered (Supplementary Methods). The two data sets (fresh and frozen) showed a high similarity between their average

gene expression (r 0.96; Supplementary Methods and

Supplementary Fig. 8a). Fifty-seven genes showed twofold upregulation in the frozen sample, with B50% being ribosomal protein genes, and the rest not enriched in any pathways (Supplementary Table 5). In addition, the number of genes and UMI counts detected from fresh and frozen PBMCs was very similar (P 0.8 and 0.1, respectively), suggesting that the

conversion efciency of the system is not compromised when proling frozen cells (Supplementary Fig. 8b). Furthermore, subpopulations were detected from frozen PBMCs at a similar proportion to that of fresh PBMCs, demonstrating the applicability of our method on frozen samples (Supplementary Methods and Supplementary Fig. 8c).

Genotype-based method to detect individual cell populations. Next, we applied the GemCode technology to study host and donor cell chimerism in an allogeneic hematopoietic stem cell transplant (HSCT) setting. Following a stem cell transplant, it is important to monitor the proportion of donor and host cells in major cell lineages to ensure complete engraftment and as a sensitive means of detecting impending relapse. Currently, the amount of host and donor chimerism is often measured from ow-sorted cell populations using PCR assays with a panel of SNV-specic primers. Current clinical chimerism tests have a number of limitations28, namely (1) the ow-sorted cell populations are limited by cell surface markers, (2) only populations with sufcient cell counts can be used for PCR assays and (3) they are not intended for the detection of minimal residual disease. Here we present a simple method that addresses these limitations, resolves host and donor chimerism at single-cell resolution and enables extensive characterization of cell subtypes by integrating scRNA-seq with de novo SNV calling.

While previous studies have used existing SNVs from DNA sequencing or large-scale copy number changes in the transcriptome data to distinguish cells by genotype2932, these methods cannot be applied to transplant samples where donor and host genotype is not known a priori, and when donor and host are closely matched in genotype. To address these limitations, we rst developed a method to infer the relative presence of host and donor genotypes in a mixed population based on SNVs directly predicted from the transcriptome data. The method identies SNVs and infers a genotype at each SNV. It then classies cells based on their genotypes across all SNVs (Supplementary Methods).

To evaluate the technical performance of this method, we generated scRNA-seq libraries from PBMCs of two healthy Donors B and C, with B8k cells captured for each sample (Supplementary Table 1). We rst performed in silico mixing of PBMCs B and C at 12 mixing ratios ranging from 0 to 50%. Only condently mapped reads from samples B and C were used, and a total of 6,000 cells were selected (Supplementary Methods). There were B15k reads per cell, with B50 ltered SNVs per cell (Supplementary Methods, Supplementary Fig. 9a,b, and Supplementary Tables 1 and 3). We then classied the cells based on variants detected from the mixed transcriptome. Sensitivity and positive predictive values (PPV) were calculated by comparing the predicted call of each cell against its true labelling. Our method was able to identify minor genotypes as low as 3% at 495% sensitivity and PPV (Fig. 4a,b). A minor population could not be detected when the mixed ratio was below 3% (Fig. 4c). The accuracy of this method is affected by the number of observed SNVs per cell, which is dependent on cell types, diversity between subjects and variant calling sensitivity. Nevertheless, the base error rate and

6 NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049 ARTICLE

a b

1.0

Positive predictive value (PPV)

0.9

Called mix fraction (%)

0.9

0.8

Sensitivity

0.7

0.6

0.5

0 1 2 3 4 5 6 7 8 9 10 50

Minor population %

0 1 2 3 4 5 6 7 8 9 10 50

Minor population %

Not confident call Confident call

0 1 2 3 4 5 6 7 8 9 10 Actual mix fraction (%)

Figure 4 | Genotype analysis of in silico and in vitro mixing of PBMCs. (a) Sensitivity versus percentage of minor population, where sensitivity is evaluated against the true labelling of in silico mixed PBMCs from Donors B and C. Red line indicates that the major population comes from Donor B PBMCs. Blue line indicates that the major population comes from Donor C PBMCs. (b) Positive predictive value (PPV) versus percentage of minor population, where PPV is evaluated against the true labelling of in silico mixed PBMCs from Donors B and C. Red line indicates that the major population comes from Donor B cells. Blue line indicates that the major population comes from Donor C cells. (c) Called mix fraction versus actual mix fraction in in silico mixing of PBMCs from Donors B and C. Fifty per cent actual mix fraction is correctly called, but omitted from the plot so that the rest of the ratios can be clearly displayed.

Table 1 | Genotype comparison of predicted genotype groups to puried populations.

Sample Observed % of minor population

Expected % of minor population

% of Genotype overlap with pre-transplant sample (host)

Likely identity

Genotype group

Genotype overlap with Donor B

PBMCs

Genotype overlap with Donor C

PBMCs

AML027 post-transplant

B only 0 0 1 100 77

C only 0 0 1 77 100

B:C 50:50 43 50

1 63 94

2 96 58

B:C 90:10 12 10

1 13.8 52 Donor

2 86.2 97 Host

AML035post-transplant 1 100 78 Donor

AML, acute myelod leukaemia.

1 47 97

2 82 74

B:C 99:1 Not

detected

1 1 97 77

Table 2 | Predicted genotype groups and their genotype overlap with pre-transplant samples.

Sample Genotype group

% of Genotype group

variant calling errors have a limited effect on the accuracy of the method, as the method uses all instead of a small subset of SNVs (Supplementary Fig. 9c).

We further validated the performance of the method in experiments where PBMCs from Donors B and C were mixed at three ratios, 50:50, 90:10 and 99:1, before scRNA-seq. In the 50:50 mixture sample, cells from Donors B and C were almost indistinguishable by RNA expression (Supplementary Fig. 9d,e). However, they can be separated by their genotype at the correct proportion (Table 1). In addition, the genotype overlap between genotype group 1 and Donor C is 94%, whereas the overlap between genotype group 1 and Donor B is only 63%, both within the range of positive and negative controls, suggesting that group 1 comes from Donor C (Supplementary Methods and Table 1). Similarly, genotype group 2 was inferred to be from Donor B (Supplementary Methods and Table 1). The proportions of the minor genotype were accurately predicted at the 90:10 mixing ratio. Consistent with the in silico mixing results, the minor population could not be detected when B and C were mixed at 99:1 ratio (Table 1).

Single-cell analysis of transplant bone marrow samples. Single-cell RNA-seq libraries were generated from cryopreserved bone marrow mononuclear cell (BMMC) samples obtained from

two patients before and after undergoing HSCT for AML (AML027 and AML035) (Supplementary Table 2). Since HSCT samples are fragile, cells were carefully washed in PBS with 20% fetal bovine serum (FBS) before loading them into chips. Relative to BMMCs from two healthy controls, we found the median number of UMI counts per cell to be 35 times higher in AML samples at B15k reads per cell, suggesting their vastly abnormal transcriptional programs (Supplementary Fig. 10a). Approximately 35 and 60 SNVs per cell were detected from AML027 and AML035 pre-transplant samples, respectively (Supplementary Table 3 and Supplementary Fig. 10b,c). Our SNV analysis detected the presence of two genotypes in the post-transplant sample of AML027: one at 13.8% and one at86.2% (Table 2). As expected, there was no evidence of multiple genotype groups in the pre-transplant host sample. We compared the major and minor inferred genotypes present in the post-transplant sample to the genotype found in the host cells. The major inferred genotype in the post-transplant sample was 97% similar to that inferred from the host sample, while the minor inferred genotype was only 52% similar to that of the host sample (Table 2). The observed range of genotype overlap between the same individuals is B98% (Supplementary

Methods), indicating errors in the genotypes inferred from individual SNVs. Ninty-seven per cent is within the observed range, and this result suggests that the post-transplant sample consists mainly (86.2%) of host cells. This observation is

NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049

Healthy control 1

AML027 pre-transplant (host)

AML035 pre-transplant (host)

Blasts and

Immature Ery 1

Immature

Granulocytes

Monocytes

tSNE2

Immature Ery 2

Immature Ery 3

AML027 post-transplant (donor)

AML027 post-transplant (host)

AML035 post-transplant (donor)

Mature Ery

tSNE2

tSNE1

tSNE1 tSNE1

Healthy control 1

Healthy control 2

AML027 pre-transplant (host)

AML027 post-transplant (host)

AML027 post-transplant (donor)

AML035 pre-transplant (host)

AML035 post-transplant (donor)

0 25 50 75 100

Percent of cells

Figure 5 | Genotype and single-cell expression analysis of transplant BMMCs. (a) tSNE projection of scRNA-seq data from a healthy control, AML027 pre- and post-transplant samples (post-transplant sample is separated into host and donor) and AML035 pre- and post-transplant samples. tSNE projection was also performed on a second healthy control, but the plot is not included here as it is very similar to that of the rst healthy control. Each cell is coloured by their classication, which is labelled next to the cell clusters. (b) Proportion of subpopulations in each sample.

consistent with the clinical chimerism assay, which demonstrated only 12% donor in the post-transplant sample. In contrast, SNV analysis on the post-HSCT sample from AML035 did not detect the presence of two genotype groups. The sample only shared 78% similarity with AML035 host cells, suggesting that the post-HSCT sample was all donor-derived (Table 2). This nding was validated by the independent clinical chimerism assay.

SNV and scRNA-seq analyses enable subpopulation comparison between individuals within and across multiple samples. We applied these analyses on BMMC scRNA-seq data from healthy controls and AML patients (Supplementary Methods), and observed subpopulation differences in AML patients after HSCT. First, while T cells dominate the healthy BMMCs and donor cells of the AML027 post-transplant sample as expected, erythroids constitute the largest population among AML samples (Fig. 5a). Different sets of progenitor and differentiation markers (for example, CD34, GATA1, CD71 and HBA1) were detected among the erythroids, indicating populations at various stages of erythroid development (Supplementary Methods and Supplementary Fig. 10df). AML027 showed the highest level of erythroid cells (480%, consist of mostly mature erythroids) before transplant, consistent with

the erythroleukaemia diagnosis of AML027 (Fig. 5b). In contrast, after transplant, AML027 showed the highest level of blast cells and immature erythroids (CD34 , GATA1 ), consistent

with the relapse diagnosis and return of the malignant host AML (Fig. 5b). These observations would have been difcult to make with FACS analysis, with limited number of markers for early erythroid lineages. Second, B20% cells in the AML027 post-transplant sample show markers of immature granulocytes (AZU1, IL8; Fig. 5b and Supplementary Fig. 10df), which are absent in AML035 post-transplant sample, and generally low among AML patients31. These cells lack marker expression for mature cells, suggesting the presence of residual precursor cells that may be part of the leukemic clone. Third, monocytes are abundant in both AML patients before transplant (10% and 25% in AML027 and AML035 respectively), but are not detectable after transplant (Fig. 5b). Monocytes have been previously identied in post-transplant samples, and the unexpected monocytopenia needs to be followed up with additional studies. Taken together, the analysis provided insights into the cellular composition and possible presence of residual disease in the bone marrows of HSCT recipients that was not available from routine clinical assays.

8 NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049 ARTICLE

DiscussionHere we present a droplet-based scRNA-seq technology that enables encapsulation of tens of thousands of single cells within minutes. We demonstrated the scalability and robustness of the system through transcriptome analysis of B250k single cells across 29 samples. scRNA-seq of cell lines and synthetic RNAs showed the systems comparable sensitivity to other droplet-based methods7,8.

The GemCode technology platform enables high-throughput scRNA-seq with rapid cell encapsulation and a high cell capture rate that addresses the challenges associated with existing scRNA-seq platforms (Supplementary Table 6). Single gel beads are encapsulated into GEMS at B80% ll rate. This ll rate combined with Poisson loading of cells results in B50% cell capture rate, enabling the processing of samples with limited cell input material. We demonstrate the ability to load from 1,000 to 23,000 cells per well, from four different cell lines and two primary cell types (PBMCs and BMMCs), illustrating the applicability of the GemCode platform to a wide variety of cell types. The GEM-based encapsulation of single cells within the microuidics platform reduces the need for expensive sorting equipment and complicated workows involving large numbers of plates (Supplementary Table 6). The scalability and high-throughput nature of the GemCode platform is achieved in two ways: hundreds to thousands of cells can be encapsulated per channel, and each chip has eight channels. Therefore, a large number of cells can be processed within a very short period of time, minimizing the perturbation of the cellular transcriptome. In addition, multiple samples can be processed simultaneously, a key advantage for experimental setups that involve a time course or multiple treatments.

The speed, reproducibility and high cell capture rate of GemCode technology allowed us to prole the transcriptome of B17,000 single BMMCs, which are notoriously fragile and difcult to work with, as these cells are isolated from patients that have undergone multiple rounds of high-intensity chemotherapy, transplant conditioning and immunosuppression. Other existing single-cell RNA-seq methods would take signicantly more time to capture the comparable number of cells (Supplementary Table 6), which could lead to RNA degradation and vast cell lysis, compromising the ability to analyse the cellular transcriptome of leukemic BMMCs.

We also developed a novel method to infer cell origin in transplant bone marrow samples by using SNVs identied from scRNA-seq data, without prior knowledge of either individuals genotype. To our knowledge, this is the rst report of a high-throughput method to determine chimerism of immune cell populations. This allowed us to discover new insights into the disease state of the host before and after transplant that were not readily achievable with traditional PCR, FACS-based analysis or any other methodology described to date. For example, in one patient, AML027, we detected a population of atypical blast and granulocyte precursors that were likely refractory to chemotherapy before transplantation and thus responsible for disease relapse. As information on erythrocyte precursors is not easy to obtain because of limited number of surface markers for these populations, the current clinical means to assess the presence of atypical myeloid blasts by ow and morphology failed to detect the expanded erythroid precursor population. While this proof-of-principle study only contained BMMCs from two transplant patients, it highlights the potential clinical impact of the technology, and lays the foundation for more extensive studies involving larger numbers of patient samples. It is our belief that the GemCode single-cell RNA-seq technology coupled with de novo chimerism testing will, in the near future, greatly expand research possibilities for clinicians and

basic scientists, and ultimately lead to improved patient care and survival.

Methods

High-speed imaging of gel beads and cells in GEMs. A microscope (Nikon Ti-E, 10 objective) and a high-speed video camera (Photron SA5,

frame rate 4,000 s 1) were used to image every GEM as they were generated

in the microuidic chip. A custom analysis software was used to count the number of GEMs generated and the number of beads present in each GEM, based on edge detection and the contrast between bead edges and GEM edges and the adjacent liquid. The results of the analysis are summarized in Fig. 1c. To estimate the distribution of cells in GEMs, manual counting was used for B28k frames of one video on a subset of GEMs. The results indicate an approximate adherence to a Poisson distribution. However, the percentage of multiple cell encapsulations was 16% higher than the expected value, possibly due to subsampling error or to cellcell interactions (some two-cell clumps were observed during the manual count) (Supplementary Fig. 1b).

Cell lines and transplant patient samples. Jurkat (ATCC TIB-152), 293T (ATCC CRL-11268) and 3T3 (ATCC CRL-1658) cells were acquired from ATCC and cultured according to ATCC guidelines. Fresh PBMCs, frozen PBMCs and BMMCs were purchased from ALLCELLS. Frozen PBMCs from Donor A were made from fresh PBMCs from Donor A by mixing 1e6 cells in freezing medium (15% dimethylsulphoxide (DMSO) in Iscoves modied Dulbeccos media containing 20% FBS) gently, and chilled in CoolCell FTS30 (BioCision) in

80 C for at least 4 h before transferring to liquid nitrogen for storage for 3 weeks.

The Institutional Review Board at the Fred Hutchinson Cancer Research Center approved the study on transplant samples. The procedures followed were in accordance with the Declaration of Helsinki of 1975 and the Common Rule. Samples were obtained after patients had provided written informed consent on molecular analyses. We identied patients with AML undergoing allogeneic hematopoietic stem cell transplant at the Fred Hutchinson Cancer Research Center. The diagnosis of AML was established according to the revised criteria of the World Health Organization33.

Bone marrow aspirates were obtained for standard clinical testing2030 days before transplant and serially post-transplanted according to the treatment protocol. Bone marrow aspirate aliquots were processed within2 h of the draw. The BMMCs were isolated using centrifugation through a Ficoll gradient (Histopaque-1077; Sigma Life Science, St Louis, MO, USA).

The BMMCs were collected from the serum-Ficoll interface with a disposable Pasteur pipette and transferred to the 50 ml conical tube with 2% patient serum in 1 PBS. The BMMCs were counted using a haemacytometer and viability

was assessed using Trypan blue. The BMMCs were resuspended in 90% FBS, 10% DMSO freezing media and frozen using a Thermo Scientic Nalgene Mr Frosty (Thermo Scientic) in a 80 C freezer for 24 h before being

transferred to liquid nitrogen for long-term storage.

Estimation of RNA content per cell. The amount of RNA per cell type was determined by quantifying (Qubit; Invitrogen) RNA extracted (Maxwell RSC simplyRNA Cells Kit) from several different known numbers of cells.

Cell preparation. Fresh cells were harvested, washed with 1 PBS and

resuspended at 1 106 cells per ml in 1 PBS and 0.04% bovine serum albumin.

Fresh PBMCs were frozen at 10 by resuspending PBMCs in DMEM 40%

FBS 10% DMSO, freezing to by C in a CoolCells FTS30 (BioCision) and

then placed in liquid nitrogen for storage.

Frozen cell vials from ALLCELLS and transplant studies were rapidly thawed in a 37 C water bath for B2 min. Vials were removed when a tiny ice crystal was left. Thawed PBMCs were washed twice in the medium and then resuspendedin 1 PBS and 0.04% bovine serum albumin at room temperature. Cells were

centrifuged at 300 r.c.f. for 5 min each time. Thawed BMMCs were washed and resuspended in 1 PBS and 20% FBS. The nal concentration of thawed cells was

1 106 cells per ml.

Sequencing library construction using the GemCode platform. Cellular suspensions were loaded on a GemCode Single-Cell Instrument (10x Genomics, Pleasanton, CA, USA) to generate single-cell GEMs. Single-cell RNA-Seq libraries were prepared using GemCode Single-Cell 30 Gel Bead and Library

Kit (now sold as P/N 120230, 120231, 120232, 10x Genomics). GEM-RT was performed in a C1000 Touch Thermal cycler with 96-Deep Well Reaction Module (Bio-Rad; P/N 1851197): 55 C for 2 h, 85 C for 5 min; held at 4 C.

After RT, GEMs were broken and the single-strand cDNA was cleaned up with DynaBeads MyOne Silane Beads (Thermo Fisher Scientic; P/N 37002D)and SPRIselect Reagent Kit (0.6 SPRI; Beckman Coulter; P/N B23318).

cDNA was amplied using the C1000 Touch Thermal cycler with 96-DeepWell Reaction Module: 98 C for 3 min; cycled 14 : 98 C for 15 s, 67 C for 20 s,

and 72 C for 1 min; 72 C for 1 min; held at 4 C. Amplied cDNA product

NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049

was cleaned up with the SPRIselect Reagent Kit (0.6 SPRI). The cDNA

was subsequently sheared to B200 bp using a Covaris M220 system (Covaris; P/N 500295). Indexed sequencing libraries were constructed using the reagents inthe GemCode Single-Cell 30 Library Kit, following these steps: (1) end repair and

A-tailing; (2) adapter ligation; (3) postligation cleanup with SPRIselect; (4) sample index PCR and cleanup. The barcode sequencing libraries were quantied by quantitative PCR (KAPA Biosystems Library Quantication Kit for Illumina platforms P/N KK4824). Sequencing libraries were loaded at 2.1 pM on an Illumina NextSeq500 with 2 75 paired-end kits using the following read length: 98 bp

Read1, 14 bp I7 Index, 8 bp I5 Index and 10 bp Read2. Some earlier libraries were made with 5 nt UMI, and 5 bp Read2 was obtained instead. These libraries have been documented in Supplementary Table 1.

ERCC assay. ERCC synthetic spike-in RNAs (Thermo Fisher Scientic;

P/N 4456740) were diluted (1:10 or 1:50) and loaded into a GemCode Single-Cell Instrument, replacing cells normally used to generate GEMs. Spike-in Mix1 and Mix2 were both tested. A slightly modied protocol was used as only a small fraction of GEMs were collected for RT and cDNA amplication. After the completion of GEM-RT, 1.25 ml of the emulsion was removed and addedto a biphasic mixture of Recovery Agent (125 ml) (P/N 220016) and 25 mM additive 1 (30 ml) (P/N 220074, 10x Genomics). The recovery agent was then removed and the remaining aqueous solution was cleaned up with the SPRISelect Reagent Kit (0.8 SPRI). cDNA was amplied using the C1000 Touch

Thermal cycler with 96-Deep Well Reaction Module: 98 C for 3 min; cycled14 : 98 C for 15 s, 67 C for 20 s, and 72 C for 1 min; 72 C for 1 min; held at

4 C. Amplied cDNA product was cleaned up with the SPRIselect Reagent Kit(0.8 ) cDNA was subsequently sheared to B200 bp using a Covaris M220

system to construct sample-indexed libraries with 10x Genomics adapters. Expected ERCC molecule counts were calculated based on the amount ofERCC molecules used and sample dilution factors. The counts were compared to detected molecule counts (UMI counts) to calculate conversion efciency.

ddPCR assay. Jurkat cells were used in ddPCR assays to estimate conversion efciency as follows: (1) the amount of RNA per Jurkat cell was determined by quantifying (Qubit, Invitrogen) RNA extracted (Maxwell RNA Purication Kits) from several different known number of Jurkat cells. (2) Bulk RT-ddPCR (Bio-Rad One-Step RT-ddPCR Advanced Kit for Probes 1864021) was performed on the extracted RNA to determine the copy number per cell of eight selected genes. (3) Approximately 5,000 Jurkat cells were processed using the GemCode Single-Cell 30 platform, and single-stranded cDNA was collected after RT in GEMs following the protocols listed in the section Sequencing library construction using the GemCode platform. cDNA copies of the eight genes were determined using ddPCR (Bio-Rad ddPCR Supermix for Probes (no dUTP) P/N 1863024). The actual Jurkat cell count was found by sequencing a subset of the GEM-RT reactions on a MiSeq. The conversion efciency is the ratio between cDNA copies percell (step 3) and RNA copies per cell from bulk RT-ddPCR (step 2), assuming a 50% efciency in RT-ddPCR34.

The probe sequences for the ddPCR assay are as follows: SERAC1_f, 50-CACGAGCCGCCAGC-30 and SERAC1_r, 50-TCTGCAACAGATGAC

GCAATAAG-30; SERAC1_p: /56-FAM/CGCCTGCCG/ZEN/GCAGAATGTC/ 3IABkFQ/. AP1S3_f, 50-GAAGCAGCCATGGTCTAAGC-30 and AP1S3_r, 50-CCTTGTCGACTGAAGAGCAATATG-30; AP1S3_p: /56-FAM/CGGCCC AGC/ZEN/CACGATGATACAT/3IABkFQ/OR. AOV1_f, 50-CCGGAAGTG GGTCTCGTOR-30 and AOV1_r, 50-TTCTTCATAGCCTTCCCGATACCOR-30; AOV1_p: /56-FAM/TCGTGATGG/ZEN/CGGATGAGAGGTTTCA/3IABkFQ/. DOLPP1_f, 50-ATGGCAGCGGACGGA-30 and DOLPP1_r, 50-GGCTCAG

GTAGGCAAGGA-30; DOLPP1_p: /56-FAM/CCACGTCGA/ZEN/ATATCCT GCAGGTGATCT/3IABkFQ/. KPNA6_f, 50-TGAAAGCTGCCGCTGAAG-30and KPNA6_r, 50-CCCTGGGCTCGCCAT-30; KPNA6_p: /56-FAM/CGGACCCG C/ZEN/GATGGAGACC/3IABkFQ/. ITSN2_f, 50-GTGACAGGCTACGCAA CAG-30 and ITSN2_r, 50-TCCTGAGTTTTCCTTGCTAGCT-30;

ITSN2_p: /56-FAM/AGGGCGCCA/ZEN/GATGGCTGA/3IABkFQ/. LCMT1_f, 50-GTCGACCCCGCTTCCA-30 and LCMT1_r, 50-GGTCATGCCAGTA

GCCAATG-30; LCMT1_p: /56-FAM/ATGCTTCCC/ZEN/TGTGCAAGAGG TTTGC/3IABkFQ/. AP2M1_f, 50-GCAGCGGGCAGACG-30 and AP2M1_r, 50-ATGGCGGCAGATCAGTCT-30; AP2M1_p: /56-FAM/CATCGCTCT/ZEN/ GAGAACAGACCTGGTG/3IABkFQ/.

Cell capture efciency calculation. The efciency is calculated by taking the ratio of the number of cells detected by sequencing versus the number of cells loaded into the chip. The latter is determined from (volume added input concentration

of cells). The input concentration of cells was determined using a Countess II Automated Cell Counter (Thermo Fisher Scientic). It is worth noting that there is a 1520% error in cell counts, which could account for at least some of the variability in the calculated efciencies.

Chimerism assay. PowerPlex 16 System (Promega) was used in conjunction with an Applied Biosystems (Life Technologies) 3130xl Genetic Analyzer. Donor BMMCs were used as the reference baseline.

Alignment, barcode assignment and UMI counting. The Cell Ranger Single-Cell Software Suite was used to perform sample demultiplexing, barcode processing and single-cell 30 gene counting (http://software.10xgenomics.com/single-cell/overview/welcome

Web End =http://software.10xgenomics.com/single-cell/overview/

http://software.10xgenomics.com/single-cell/overview/welcome

Web End =welcome ). First, sample demultiplexing was performed based on the 8 bp sample index read to generate FASTQs for the Read1 and Read2 paired-end reads, as well as the 14 bp GemCode barcode. Ten basepair UMI tags were extracted from Read2 (14 libraries were made with 5 bp UMI tags, as noted in Supplementary Table 1, due to an earlier iteration of the methods. For these samples, 5 bp UMI tags were extracted from Read2.). Then, Read1, which contains the cDNA insert, was aligned to an appropriate reference genome using STAR35. For mouse cells, mm10 was used and for human cells, hg19 was used. For samples with mouse and human cell mixtures, the union of hg19 and mm10 were used. ForERCC samples, ERCC reference (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt

Web End =https://tools.thermosher.com/content/sfs/ https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt

Web End =manuals/cms_095047.txt ) was used.

Next, GemCode barcodes and UMIs were ltered. All of the known listed of barcodes that are 1-Hamming-distance away from an observed barcode are considered. Then, the posterior probability that the observed barcode was produced by a sequencing error is computed, given the base qualities of the observed barcode and the prior probability of observing the candidate barcode (taken from the overall barcode count distribution). If the posterior probability for any candidate barcode is at least 0.975, then the barcode is corrected to the candidate barcode with the highest posterior probability. If all candidate sequences are equally probable, then the one appearing rst by lexicalorder is picked.

UMIs with sequencing quality score 410 were considered valid if they were not homopolymers. Qual 10 implies 90% base call accuracy. A UMI that is

1-Hamming-distance away from another UMI (with more reads) for the same cell barcode and gene is corrected to the UMI with more reads. This approach is nearly identical to that in Jaitin et al.4, and is similar to that in Klein et al.8 (although Klein et al.8 also used UMIs to resolve multimapped reads, which was not implemented here).

Last, PCR duplicates were marked if two sets of read pairs shared a barcode sequence, a UMI tag, and a gene ID (Ensembl GTFs GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/ Homo_sapiens.GRCh37.82.gtf.gz and GRCm38.84, ftp://ftp.ensembl.org/pub/ release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz, were used). Only condently mapped (MAPQ 255), non-PCR duplicates with valid

barcodes and UMIs were used to generate gene-barcode matrix.

Cell barcodes were determined based on distribution of UMI counts. All top barcodes within the same order of magnitude (410% of the top nth barcode, where n is 1% of the expected recovered cell count) were considered cell barcodes. Number of reads that provide meaningful information is calculated as the product of four metrics: (1) valid barcodes; (2) valid UMI; (3) associated with a cell barcode; and (4) condently mapped to exons.

In the mouse and human mixing experiments, multiplet rate was dened as twice the rate of cell barcodes with signicant UMI counts from both mouse and human, where top 1% of UMI counts was considered signicant. The extent of barcode crosstalk was assessed by the fraction of mouse reads in human barcodes, or vice versa.

Samples processed from multiple channels can be combined by concatenating gene-cell-barcode matrices. This functionality is provided in the Cell RangerR Kit (http://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit

Web End =http://support.10xgenomics.com/single-cell/software/pipelines/latest/rkit). Sequencing data from multiple sequencing runs of a library can be combinedby counting non-duplicated reads. This functionality is provided in the Cell Ranger pipeline. In addition, sequencing data can be subsampled to obtain a given number of UMI counts per cell. This functionality is also provided in the Cell Ranger R Kit, and is useful when combining data from multiple samples for comparison.

PCA analysis of mixing of Jurkat and 293T cells. Gene-cell-barcode matrix from each of the four samples was concatenated. Only genes with at least one UMI count detected in at least one cell are used. UMI normalization was performed by rst dividing UMI counts by the total UMI counts in each cell, followed by multiplication with the median of the total UMI counts across cells. Then, we took the natural log of the UMI counts. Finally, each gene was normalized such that the mean signal for each gene is 0, and standard deviation is 1. PCA was run on the normalized gene-barcode matrix. The normalized UMI counts of each gene are used to show expression of amarker in a tSNE plot.

SNV analysis of Jurkat and 293T scRNA-seq data. SNVs were called by running Freebayes 1.0.2 (ref. 36) on the genome BAM produced by Cell Ranger. High-quality SNVs (SNV calling Qual4 100 with at least 10 UMI counts from at

least two cells; indels are ignored) that were only observed in Jurkat or 293T cells (but not both) were selected. Cells were labelled as Jurkat or 293T basedon Jurkat- and 293T-specic SNV counts, where the fraction of countsfrom the other species is o0.2. Cells with a fraction of SNV from eitherspecies between 0.2 and 0.8 are considered multiplets. The inferred multiplet rate is 2* observed multiplet rate (to account for Jurkat:Jurkat and293T:293T multiplets).

10 NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049 ARTICLE

PCA and tSNE analysis of PBMCs. Genes with at least one UMI count detected in at least one cell are used. The top 1,000 most variable genes were identied based on their mean and dispersion (variance/mean), which is similar to the approach used by Macoscko et al.7 Genes were placed into 20 bins based on their mean expression. Normalized dispersion is calculated as the absolute difference between dispersion and median dispersion of the expression mean, normalized by median absolute deviation within each bin.

PCA was run on the normalized gene-barcode matrix of the top1,000 most variable genes to reduce the number of feature (gene) dimensions. UMI normalization was performed by rst dividing UMI counts by the total UMI counts in each cell, followed by multiplication with the median of the total UMI counts across cells. Then, we took the natural log of the UMI counts. Finally, each gene was normalized such that the mean signal for each gene is 0, and standard deviation is 1. PCA was run on the normalized gene-barcode matrix. After running PCA, Barnes-hut37 approximation to t-SNE16 was performed on the rst 50 PCs to visualize cells in a two-dimensional space. Fifty PCs were used because: (1) using all PCs would take a very long time with tSNE analysis;(2) they explained B25% of total variance. K-means15 clustering was run to group cells for the clustering analysis. k 10 was selected based on the sum of

squared error scree plot (Supplementary Fig. 5d).

Identication of cluster-specic genes and marker-based classication. To identify genes that are enriched in a specic cluster, the mean expression of each gene was calculated across all cells in the cluster. Then each gene from the cluster was compared to the median expression of the same gene from cells in all other clusters. Genes were ranked based on their expression difference, and the top10 enriched genes from each cluster were selected. For hierarchical clustering, pair-wise correlation between each cluster was calculated, and centred expression of each gene was used for visualization by heatmap.

Classication of PBMCs was inferred from the annotation of cluster-specic genes. In the case of cluster 10, marker expression of multiple cell types(for example, B, dendritic and T) was detected. Since the relative clustersize of B, dendritic and T is 5.7%, 6.6% and 81%, respectively, we would expect cluster 10 (which is only 0.5%) to contain multiplets consisting mostly from B:dendritic (0.36%) and B:dendritic:T (0.3%).

Selection of puried subpopulations of PBMCs. Each population of puried PBMCs was downsampled to B16k reads per cell. PCA, tSNE and k-means clustering were performed for each downsampled matrix, following the same steps outlined in PCA and t-SNE analysis of PBMCs. Only one cluster was detected in most samples, consistent with the FACS analyses (Supplementary Fig. 6).

For samples with more than one cluster, only clusters that displayed the expected marker gene expression were selected for downstream analysis. For CD14 monocytes, two clusters were observed and identied as CD14

monocytes and dendritic cells based on expression of marker genes FTL and CLEC9A, respectively.

Cell classication analysis using puried PBMCs. Each population of puried PBMCs was downsampled to B16k condently mapped reads per cell. Then, an average (mean) gene expression prole across all cells was calculated. Next, gene expression from every cell of the complex population was compared with the gene expression proles of puried populations of PBMCs by Spearmans correlation. The cell was assigned the ID of the puried population if it had the highest correlation with that population. Note that the difference between the highest and second highest correlation was small for some cells (for example,the difference between cytotoxic T and NK cells), suggesting that the cell assignment was not as condent for these cells. A few of the puried PBMC populations overlapped with each other. For example, CD4 T-helper cells include all

CD4 cells. This means that cells from this sample will overlap with cells

from samples that contain CD4 cells, including CD4 /CD25 T reg,

CD4 /CD45RO T memory, CD4 /CD45RA /CD25 naive T. Thus, when

a cell was assigned the ID of CD4 T-helper cell based on the correlation score,

the next highest correlation was checked to see if it was one of the CD4 samples.

If it was, the cells ID was updated to the cell type with the next highest correlation. The same procedure was performed for CD8 cytotoxic T and

CD8 /CD45RA naive cytotoxic T (which is a subset of CD8 cytotoxic T).

The R code used to analyse 68k PBMCs and puried PBMCs can be found here: https://github.com/10XGenomics/single-cell-3prime-paper

Web End =https://github.com/10XGenomics/single-cell-3prime-paper .

Cell clustering and classication with Seurat. The gene-cell-barcode matrix of 68k PBMCs was log-transformed as an input to Seurat. The top 469 most variable genes selected by Seurat were used to compute the PCs. The rst 22 PCs were signicant (Po0.01) based on the built-in jackstraw analysis, and usedfor tSNE visualization. Cell classication was taken from Cell classication analysis using puried PBMCs.

Comparison between fresh and frozen PBMCs. The sequencing data of 68k fresh PBMCs and 3k frozen PBMCs were down-sampled such that each sample has

B14k condently mapped reads per cell. Only genes that are detected in atleast one cell were included for the comparison, which uses the mean of each gene across all cells.

For cell classication comparison between puried and frozen PBMCs, we pooled all the cells labelled as T or natural killer cells together. This is because the subpopulations within T and between T and natural killer cells are sometimes difcult to cluster separately. We did not want the comparison between fresh and frozen cells to be affected by the clustering methods used.

SNV-based genotype assignment. SNVs were called by running Freebayes 1.0.2 (ref. 36) on the genome BAM produced by Cell Ranger. Only SNVs with support from at least two cell barcodes, with a minimal SNV Qual score 4 30, minimal

SNV base Qual4 1 were included. Reference (R) and alternate (A) allele counts

were computed at each SNV, producing a matrix of cell-reference UMI counts and cell-alternate-allele UMI counts. These matrices were modeled as a mixture of two genomes where the likelihood of any of the three genotypes (R/R, R/A or A/A) at a site was taken to be binomially distributed with a xed error rate of 0.1%. For each sample, two models were inferred in parallel, one where only one genome is present (K 1) and another where two genomes are present (K 2). Inference of the

model parameters (cell-to-genome assignments and the K sets of genotypes) was performed by using a Gibbs sampler to approximate their posterior distributions. To ameliorate the label-switching problem in Monte Carlo inference of mixture models, relabelling of the sampled cell-to-genome assignments was performed as per Stephens et al.38

In in silico cell mixing experiments, when the K 2 model failed to adequately

separate the two genomes, it reported a distribution of posterior probabilities near0.5 for the cell-genome calls, indicating a lack of condence in those calls. We applied a requirement that 90% of the cells have a posterior probability 475% to select the K 2 model over the K 1 model. Selecting K 1 indicates that the

mixture fraction is below the level of detection of the method, which in in silico mixing experiments was determined to be 4% of 6,000 cells.

Genotype comparison with the pure sample. To ascertain the assignment of genotypes to individuals, only shared SNVs between the genotype group and the pure sample were considered. Then, the average genotype of all the cells was compared to that of the pure sample. To obtain some baseline for the % genotype overlap among different individuals, we performed pairwise comparison of genotypes called from the same individuals (11 pairwise comparisons) or from different individuals (15 pairwise comparisons). The per cent genotype overlap between the same individuals averages B980.3%, whereas the per cent genotype overlap between the different individuals averages B732%.

PCA and tSNE analysis of BMMCs. Data from six samples were used: two healthy controls, AML027 pre- and post-transplant, and AML035 pre- and post-transplant. Each sample was downsampled to B10k condently mapped reads per cell. Then the gene-cell barcode matrix from each sample was concatenated. PCA, tSNE and k-means clustering were performed on the pooled matrix, following the same steps outlined in PCA and tSNE analysis of PBMCs. For k-means clustering, K 10 was

used based on the bend in the sum of squared error scree plot.

Cluster-specic genes were identied following the steps outlined in Identication of cluster-specic genes and marker-based classication. Classication was assigned based on cluster-specic genes, and based on expression of some well-known markers of immune cell types. Blasts and Immature Ery 1 refers to cluster 4, which expresses CD34, a marker of hematopoietic progenitors39, and Gata2, a marker for early erythroids40. Immature Ery 2 refers toclusters 5 and 8, which show expression of Gata1, a transcription factor essential for erythropoiesis41, but not CD71, which are often found in more committed erythroid cells39. Immature Ery 3 refers to cluster 1, which show expressionof CD71. Mature Ery refers to cluster 2. HBA1, a marker of mature erythroid cells, is preferentially detected in cluster 2. Cluster 3 was assigned as Immature Granulocytes because of the expression of early granulocyte markers such as AZU1 and IL8 (ref. 42), and the lack of expression of CD16. Cluster 7 was assigned as Monocytes because of the expression of CD14 and FCN1, for example. B refers clusters 6 and 9 because of markers such as CD19 and CD79A. T refers to cluster 10, because of markers such as CD3D and CD8A.

Data availability. All relevant data are available from the authors. Single-cell RNA-seq data have been deposited in the Short Read Archive under accession number SRP073767. Data are also available at http://support.10xgenomics.com/single-cell/datasets

Web End =http://support.10xgenomics.com/ http://support.10xgenomics.com/single-cell/datasets

Web End =single-cell/datasets . The analysis code for the 68k PBMC analysis is available at https://github.com/10XGenomics/single-cell-3prime-paper

Web End =https://github.com/10XGenomics/single-cell-3prime-paper .

References

1. Shalek, A. K. et al. Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells. Nature 498, 236240 (2013).

2. Wills, Q. F. et al. Single-cell gene expression analysis reveals genetic associations masked in whole-tissue experiments. Nat. Biotechnol. 31, 748752 (2013).

NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14049

3. Liu, S. & Trapnell, C. Single-Cell Transcriptome Sequencing: Recent Advances and Remaining Challenges Vol. 5 F1000 Research (2016).

4. Jaitin, D. A. et al. Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types. Science 343, 776779 (2014).

5. Pollen, A. A. et al. Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex. Nat. Biotechnol. 32, 10531058 (2014).

6. Fluidigm. Single-cell whole genome sequencing on the C1 System: a performance evaluation https://www.fluidigm.com/binaries/content/documents/fluidigm/marketing/single-cell-whole-genome-sequencing/single-cell-whole-genome-sequencing/fluidigm&percnt;3Afile

Web End =https://www.uidigm.com/binaries/content/ https://www.fluidigm.com/binaries/content/documents/fluidigm/marketing/single-cell-whole-genome-sequencing/single-cell-whole-genome-sequencing/fluidigm&percnt;3Afile

Web End =documents/uidigm/marketing/single-cell-whole-genome-sequencing/single- https://www.fluidigm.com/binaries/content/documents/fluidigm/marketing/single-cell-whole-genome-sequencing/single-cell-whole-genome-sequencing/fluidigm&percnt;3Afile

Web End =cell-whole-genome-sequencing/uidigm%3Ale (2016).

7. Macosko, E. Z. et al. Highly parallel genome-wide expression proling of individual cells using nanoliter droplets. Cell 161, 12021214 (2015).

8. Klein, A. M. et al. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell 161, 11871201 (2015).

9. Soumillon, M., Cacchiarelli, D., Semrau, S., van Oudenaarden, A. & Mikkelsen, T. S. Characterization of directed differentiation by high-throughput single-cell RNA-seq. Preprint at http://biorxiv.org/content/early/2014/03/05/003236

Web End =http://biorxiv.org/content/early/2014/03/05/003236 (2016).

10. Zheng, G. X. et al. Haplotyping germline and cancer genomes with high-throughput linked-read sequencing. Nat. Biotechnol. 34, 303311 (2016).

11. Narasimhan, V. M. et al. Health and population effects of rare gene knockouts in adult humans with related parents. Science 352, 474477 (2016).

12. Mostovoy, Y. et al. A hybrid approach for de novo human genome sequence assembly and phasing. Nat. Methods 13, 587590 (2016).

13. Hindson, B. J. et al. High-throughput droplet digital PCR system for absolute quantitation of DNA copy number. Anal. Chem. 83, 86048610 (2011).

14. Brennecke, P. et al. Accounting for technical noise in single-cell RNA-seq experiments. Nat. Methods 10, 10931095 (2013).

15. Sherlock, G. Analysis of large-scale gene expression data. Curr. Opin. Immunol. 12, 201205 (2000).

16. van der Maaten, L. J. P. & Hinton, G. E Visualizing high-dimensional data using t-SNE. J. Mach. Learn. Res. 9, 25792605 (2008).

17. Stem Cell Technologies. Frequencies of cell types in human peripheral blood. Available at: http://www.stemcell.com/media/files/wallchart/WA10006-Frequencies_Cell&percnt;20Types_Human_Peripheral_Blood.pdf

Web End =http://www.stemcell.com/media/les/wallchart/WA10006- http://www.stemcell.com/media/files/wallchart/WA10006-Frequencies_Cell&percnt;20Types_Human_Peripheral_Blood.pdf

Web End =Frequencies_Cell%20Types_Human_Peripheral_Blood.pdf (2016).

18. Borrego, F., Masilamani, M., Marusina, A. I., Tang, X. & Coligan, J. E The CD94/NKG2 family of receptors: from molecules and cells to clinical relevance. Immunol. Res. 35, 263278 (2006).

19. Chu, P. G. & Arber, D. A CD79: a review. Appl. Immunohistochem. Mol. Morphol. 9, 97106 (2001).

20. Schiopu, A. & Cotoi, O. S S100A8 and S100A9: DAMPs at the crossroads between innate immunity, traditional risk factors, and cardiovascular disease. Mediat. Inamm. 2013, 828354 (2013).

21. Turman, M. A., Yabe, T., McSherry, C., Bach, F. H. & Houchins, J. P Characterization of a novel gene (NKG7) on human chromosome 19 that is expressed in natural killer cells and T cells. Hum. Immunol. 36, 3440 (1993).

22. Lubberts, E The IL-23-IL-17 axis in inammatory arthritis. Nat. Rev. Rheumatol. 11, 562 (2015).

23. Ronchetti, S. et al. Glucocorticoid-induced tumour necrosis factor receptor-related protein: a key marker of functional regulatory T cells. J. Immunol. Res. 2015, 171520 (2015).

24. Lin, Y. Y. et al. Transcriptional regulator Id2 is required for the CD4 T cell immune response in the development of experimental autoimmune encephalomyelitis. J. Immunol. 189, 14001405 (2012).

25. Greer, A. M. et al. Serum IgE clearance is facilitated by human FcepsilonRI internalization. J. Clin. Invest. 124, 11871198 (2014).

26. Harman, A. N. et al. Identication of lineage relationships and novel markers of blood and skin human dendritic cells. J. Immunol. 190, 6679 (2013).

27. Satija, R. Seurat: R toolkit for single cell genomics. http://www.satijalab.org/seurat.html

Web End =http://www.satijalab.org/ http://www.satijalab.org/seurat.html

Web End =seurat.html (2016).

28. Seattle Cancer Care Alliance. Chimerism testing/engraftment analysis. http://www.seattlecca.org/healthcare-professionals/clinical-labs/clinical-immunogenetics-laboratory/chimerism-testing

Web End =http://www.seattlecca.org/healthcare-professionals/clinical-labs/clinical- http://www.seattlecca.org/healthcare-professionals/clinical-labs/clinical-immunogenetics-laboratory/chimerism-testing

Web End =immunogenetics-laboratory/chimerism-testing (2016).

29. Patel, A. P. et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344, 13961401 (2014).

30. Tirosh, I. et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189196 (2016).

31. Lee, M. C. et al. Single-cell analyses of transcriptional heterogeneity during drug tolerance transition in cancer cells by RNA sequencing. Proc. Natl. Acad. Sci. USA 111, E4726E4735 (2014).

32. Kim, K. T. et al. Single-cell mRNA sequencing identies subclonal heterogeneity in anti-cancer drug responses of lung adenocarcinoma cells. Genome Biol. 16, 127 (2015).

33. Vardiman, J. W. et al. The 2008 revision of the World Health Organization (WHO) classication of myeloid neoplasms and acute leukemia: rationale and important changes. Blood 114, 937951 (2009).

34. Zhong, J. F. et al. A microuidic processor for gene expression proling of single human embryonic stem cells. Lab Chip 8, 6874 (2008).

35. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 1521 (2013).

36. Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. Preprint at https://arxiv.org/abs/1207.3907

Web End =https://arxiv.org/abs/1207.3907 (2012).

37. van der Maaten, L. Barnes-Hut-SNE Preprint at arXiv:1301.3342 [cs.LG] (2013).38. Stephens, M. Dealing with label switching in mixture models. J. R. Stat. Soc. Ser. B 62, 795809 (2000).

39. Liu, Q. et al. Signicance of CD71 expression by ow cytometry in diagnosis of acute leukemia. Leuk. Lymphoma 55, 892898 (2014).

40. Novershtern, N. et al. Densely interconnected transcriptional circuits control cell states in human hematopoiesis. Cell 144, 296309 (2011).

41. Bonora, M. et al. Molecular mechanisms of cell death: central implication of ATP synthase in mitochondrial permeability transition. Oncogene 34, 14751486 (2015).

42. Schinke, C. et al. IL8-CXCR2 pathway inhibition as a therapeutic strategy against MDS and AML stem cells. Blood 125, 31443152 (2015).

Acknowledgements

We thank Deanna Church for critical reading of the manuscript, and members of the Bielas laboratory and 10x Genomics for helpful discussions. We thank members of the clinical immunogenetics laboratory at the Fred Hutchinson Cancer Research Center for their assistance in sample preparation and data review: David Wu, Debra Cordell, Aida Guzman, Reena Patel, Ada Ng, Chuck Kellum and Gana Balgansuren. We acknowledge support from the Canary Foundation (to J.H.B.), an Ellison Medical Foundation New Scholar Award (AG-NS-0577-09 to J.H.B.), the National Cancer Institute (R01 CA 175215-01 and U10CA180861 to J.P.R.), an Outstanding New Environmental Scientist Award (ONES) (R01) from the National Institute of Environmental Health Sciences (R01ES019319 to J.H.B.), a grant from the Congressionally Directed Medical Research Programs/US Department of Defense (W81XWH-10-1-0563 to J.H.B.) and the Pacic Ovarian Cancer Research Consortium Ovarian Cancer SPORE Award (P50 CA083636). W.J.V is supported by an Achievement Rewards for College Scientists (ARCS) Foundation Fellowship, and a Ruth L. Kirschstein National Research service F30 Award for Individual Predoctoral MD/PhD Degree Fellows (NCI F30CA200247).

Author contributions

G.X.Y.Z., J.M.T., P.B., P.R., Z.W.B., T.S.M., B.J.H., J.H.B., E.A.S. and J.P.R. designed experiments. J.M.T., P.B., Z.W.B., S.B.Z., T.D.W., G.P.M., J.S., L.M., S.Y.N., E.A.S., N.G.E., L.W.B., H.J.D., C.M., K.R.L. and W.J.V. conducted experiments. T.D.W., D.A.M., R.B., K.D.N. and B.J.H. designed the instrument. G.P.M., Z.W.B., S.Y.N., C.M.H., P.W.W. and K.D.N. designed the reagents. P.R., R.W., A.W., G.X.Y.Z., J.J.Z., T.S.M. and M.S.L. wrote the analysis software. G.X.Y.Z., P.R., T.S.M., J.Z., K.R.L. and M.T.G. analysed the data. G.X.Y.Z., E.A.S., J.P.R., T.S.M., B.J.H. and J.H.B. wrote the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: G.X.Y.Z., J.M.T., P.B., P.R., Z.W.B., R.W., S.B.Z., T.D.W., J.J.Z., G.P.M., J.S., L.M., D.A.M., S.Y.N., M.S.L., P.W.W., C.M.H., R.B., A.W., K.D.N., T.S.M. and B.J.H. are employees of 10x Genomics.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Zheng, G. X. Y. et al. Massively parallel digital transcriptional proling of single cells. Nat. Commun. 8, 14049 doi: 10.1038/ncomms14049 (2017).

Publishers note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afliations.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2017

12 NATURE COMMUNICATIONS | 8:14049 | DOI: 10.1038/ncomms14049 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

Word count: 12514

Show less

Abstract

Translate

Characterizing the transcriptome of individual cells is fundamental to understanding complex biological systems. We describe a droplet-based system that enables 3' mRNA counting of tens of thousands of single cells per sample. Cell encapsulation, of up to 8 samples at a time, takes place in ∼6 min, with ∼50% cell capture efficiency. To demonstrate the system's technical performance, we collected transcriptome data from ∼250k single cells across 29 samples. We validated the sensitivity of the system and its ability to detect rare populations using cell lines and synthetic RNAs. We profiled 68k peripheral blood mononuclear cells to demonstrate the system's ability to characterize large immune populations. Finally, we used sequence variation in the transcriptome data to determine host and donor chimerism at single-cell resolution from bone marrow mononuclear cells isolated from transplant patients.

Details

Title

Massively parallel digital transcriptional profiling of single cells

Author

Zheng, Grace X Y; Terry, Jessica M; Belgrader, Phillip; Ryvkin, Paul; Bent, Zachary W; Wilson, Ryan; Ziraldo, Solongo B; Wheeler, Tobias D; Mcdermott, Geoff P; Zhu, Junjie; Gregory, Mark T; Shuga, Joe; Montesclaros, Luz; Underwood, Jason G; Masquelier, Donald A; Nishimura, Stefanie Y; Schnall-levin, Michael; Wyatt, Paul W; Hindson, Christopher M; Bharadwaj, Rajiv; Wong, Alexander; Ness, Kevin D; Beppu, Lan W; Deeg, H Joachim; Mcfarland, Christopher; Loeb, Keith R; Valente, William J; Ericson, Nolan G; Stevens, Emily A; Radich, Jerald P; Mikkelsen, Tarjei S; Hindson, Benjamin J; Bielas, Jason H

Pages

14049

Publication year

2017

Publication date

Jan 2017

Publisher

Nature Publishing Group

e-ISSN

20411723

Source type

Scholarly Journal

Language of publication

English

DOI

https://doi.org/10.1038/ncomms14049

ProQuest document ID

1858738069

Massively parallel digital transcriptional profiling of single cells

Jump to:

Full text

Abstract

Details

Suggested sources