Single cell transcriptional analysis reveals novel innate immune cell types

Primary cell extraction and single-cell sorting

Neutrophils and T lymphocytes were extracted from 5 ml whole blood from 6 healthy donors and isolated based on phenotype by negative selection using magnetic beads (EasySep neutrophil extraction kit, Stem Cell Technologies, or Dynabeads untouched T cells, Life Technologies). One donor’s neutrophil count was too low for further processing, therefore all results presented for neutrophils consist of n = 5. Negative selection was chosen so as to avoid cellular activation due to receptor cross-linking. For each purified cell type, flow cytometry sorting with a BD FACS Aria II gated by forward- and side scatter was utilized to deposit single cells into a 96-well PCR plate preloaded with 5 µl of lysis buffer with 0.05U Superase RNase inhibitor (Life Technologies) per well. The plates were centrifuged for 1 min at 4 °C and immediately frozen and stored at −80 °C. All donors were individuals enrolled in The Center for Health Discovery and Well-Being at Emory Midtown Hospital and provided written consent for participation in the study. The protocol for blood collection was approved by the Georgia Tech Institutional Review Board (approval #H09364).

Single cell qRT-PCR

The cellular lysates were converted to cDNA and 96 target genes per cell type were pre-amplified with a pool of 96 primer pairs targeting genes representing pattern recognition, cell-type markers, intracellular signaling, transcription, and immune response. For each donor, amplified cDNA samples from 48 cells of each type were then randomized and re-plated across 5 Fluidigm 96 × 96 microfluidic arrays, in order to avoid any plate effects confounding the analysis of single donors. Gene-specific quantitative real-time PCR reactions were performed using the Fluidigm BioMark I nano-scale platform. Negative controls (without cDNA) and samples of 10 and 100 cells were used as controls for single-cell loading. The mean difference in Ct value between 1 and 10 cells and between 10 and 100 cells per sample was determined in independent assays, providing a measurable control for single cell loading of each sample.

To enable reproducible comparison of gene expression between qRT-PCR samples, data is usually normalized with respect to data obtained for an internal or endogenous reference gene. Housekeeping genes such as β-actin and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) are most often used because their expression levels are expected to remain constant. Unfortunately, single cells exhibit large heterogeneity in housekeeping gene expression levels, and this method cannot be used as control for reproducible comparison of gene expression between single cell samples (Liss et al., 2001; Suzuki, Higgins & Crawford, 2000).

Quality control and data exclusion

In order to control for single-cell sensitivity and consistency in sample loading, single cell readings were compared to multi-cell controls. Mean expression levels from 10 randomly chosen single cell samples were calculated for each gene. The values obtained were compared by regression to the mean expression levels for the corresponding genes from the 10-cell samples, and R² values above 0.65 were observed in all cases, indicating good concordance between single cell and multicell measurements.

Raw data for gene expression were obtained as Ct values between 1 and 40, with lower Ct value indicating higher abundance of gene-specific product. Missing data points were coded as Ct values of 999; such values can either be due to null or very low expression of the target gene in question or due to a failed reaction (truly missing data). Single missing measures may indicate technical failures, but consistent absence of a similar set of lowly expressed transcripts is more likely to imply coordinated loss of expression. Downstream methods differ largely with respect to how the missing data is handled. Three different sets of criteria were used for data exclusion for each of the two (neutrophil and T lymphocyte) data sets.

(A) Supervised Data Exclusion. For the neutrophil data set, an empirical cutoff was set to transcripts present in at least 70% of cells, and subsequently to cell samples expressing at least 70% of these most uniformly expressed genes. We reasoned that the absence of the same set of genes in a common set of cells would imply true absence of expression, and used hierarchical clustering to provide a preliminary indication of such clusters of non-expressed genes. 23 such co-regulated low-abundance genes were identified, for which missing values were re-assigned a Ct value of 40 (the maximum number of cycles). Subsequently, for 36 genes, sporadic missing data was assumed to represent technical error and these values were reassigned to the average Ct for the gene in question in the data set. 34 genes and 18 cells were excluded in their entirety. Expression was evaluated for 59 genes in 202 cells. Because the T lymphocyte data set did not contain a natural cutoff for transcript presence, this method of analysis was not implemented for the T lymphocyte data.

(B) Data exclusion based on median standard deviation cutoff. All missing data values were initially set to Ct 40, and the mean Ct and number of missing data points were calculated for all genes. The second and third highest expressed genes in the data set were selected and their mean Ct and standard deviation calculated. (The highest expressed gene in both the neutrophil and T lymphocyte data set were treated as outliers and ignored for the purposes of calculating mean Ct, due to expression levels far higher than all other genes). Any gene whose average expression was within a cutoff of three standard deviations of the mean Ct value for the two chosen genes was included. All cells expressing less than half of these genes were then excluded. A plot of the maximum Ct across all cells for all 96 genes in the neutrophil data set showed a bimodal distribution of maximum Ct values, with a second peak starting at Ct 37 that corresponds to cells deemed not to express the target gene. The limit of detection (LOD) was thus set to Ct 37 for neutrophils, and the LOD Ct was set to 38 for T lymphocytes by the same methodology. All data values above LOD Ct, including Ct 999, were replaced with 37, and the LOD Ct value was then subtracted from all other Ct values, according to the Log₂EX method (Log₂EX = LOD Ct–Ct (Abdel-Rahman et al., 2008)). Consequently, the adjusted expression measure for this method is inverted and ranges from 0 to LOD Ct, with more highly expressed transcripts having higher values, more in line with intuition and with microarray or RNA-Seq data analysis. For the T cell data set, Ct values above LOD were interpreted and analyzed in two alternate ways; either as representing no expression of the target gene, with Ct values set to 0 or, alternatively, as missing data points due to technical error, with missing values replaced with average Ct for the gene (analogous to the Supervised data analysis method for neutrophils). Subsequently, entire cells were excluded, if the two gene targets with highest expression in our data set were more than three standard deviation units lower than the median. Additionally, any genes that were not expressed in any cell sample were excluded from the data set. For the T lymphocyte data set, 2 genes that were only expressed in one cell were also excluded from analysis. This resulted in the exclusion of 12 neutrophils and 31 genes in the neutrophil data set, and 7 T lymphocytes and 63 genes in the T lymphocyte data set. Expression was evaluated for 62 genes in 208 cells in the neutrophil data set and just 29 genes in 244 cells in the T lymphocyte data set.

(C) Inclusion of all data points. All data points were initially included in analysis, with the exception of genes not expressed in any of the control samples (cDNA, tRNA, 10-cell samples). This excluded 12 genes from analysis in the neutrophil data set and 13 in the lymphocyte data set. In addition, any transcripts missing from all samples in an array were excluded. This excluded 3 genes in the neutrophil set, for a total of 15 excluded genes. No cells were excluded. LOD Ct was subtracted according to the Log₂EX method as described above. Expression was evaluated for 81 genes in 220 cells for the neutrophil data set and 85 genes in 247 cells for the T lymphocyte data set.

Data normalization

Three different sets of criteria were used for data normalization for each of the data sets generated from the three methods for data exclusion.

(1) Mean centering. The mean Ct value for each cell was calculated and subtracted from each data point for the same cell. This approach removes the dependence of magnitudes, allowing for easier visualization and comparison of relative differences in expression levels.

(2) Quantile normalization. Gene expression data for each cell was re-ordered by raw Ct value, and mean Ct values for each cell were calculated. The original data was then replaced by the average quantile, such that the highest value was replaced by the mean of the highest values, the second highest value by the mean of the second highest values, and so on. This method of rank-order analysis eliminates cell-to-cell differences in data density, by making the data distributions identical.

(3) Standardization of the genes. Gene expression data were mean-centered for each cell, and then the values for each gene were standardized (converted to z-scores) by mean-centering and dividing by the standard deviation. Residuals from an ANOVA with Plate as the main effect were extracted. This method adjusts the distribution only of targets whose expression differs among plates. A further centering of residual expression values to a mean of zero for each cell ensures that no cells have artificially low or high expression of all genes.

Analysis of gene expression patterns

The single cell transcript abundance distribution for each gene was determined using SAS JMP Pro 10 (Cary, NC). For each gene, several models, including Normal, Gamma, Johnson Su, Johnson SI, Lognormal, and Weibull, were tested in order to find the model best fitting the data. Modality was assessed for the two best fitting models by Akaike information criteria (AIC) score and was further verified by calculating deltaAIC, comparing scores of bimodality, trimodality and unimodality, as well as visual observation. Genes exhibiting bimodality were tracked and cluster membership was determined in the raw data set as well as after data exclusion and normalization methods deemed most suitable, using the criteria above, namely exclusion by missing data cutoff and normalization by standardization of the genes. In addition, the number of cells included in each cluster was determined for each donor. Known gene product functionality was obtained from three data bases: ToppFun (Chen et al., 2009), DAVID v6.7 (Huang da, Sherman & Lempicki, 2009a; Huang da, Sherman & Lempicki, 2009b), and KEGG Pathway (Nakao et al., 1999).

Analysis of donor-to-donor variability

For each donor, the cell count was determined for each of the cell clusters defined within the population across all donors. Following this, the observed frequencies were compared to expected frequency by Chi-square test comparison of the number of cells of each class in each of the five individuals relative to the expectation assuming equivalent proportions.

Comparison of primary analysis methods by concordance of cell clusters

Combining the methods for data exclusion and normalization generated nine alternate sets of processed data for each of the two cell types. Each data set was organized by hierarchical clustering as well as k-means clustering by cell, resulting in cell clusters based on shared gene expression patterns. Concordance, defined as the percentage of cells ascribed to the same cluster, was compared between all combinations of analysis methods for both methods of clustering. For hierarchical clustering, data was clustered using Ward’s minimum variance method (Ward, 1963), which minimizes the total within-cluster variance using the total within-cluster sum of squares, under the assumption that distances between individual objects are proportional to Euclidean distance. The k-means method of clustering aims to sort data into a pre-defined number of clusters, k, with each data point belonging to the cluster with the nearest mean (MacQueen, 1967). K-means clustering was performed on all data sets with k values of 2 or 3 for both neutrophils and T lymphocytes. The k values were evaluated using Cubic Clustering Criteria (CCC) with external cluster validation. All computations were performed in SAS JMP-Genomics v5.0 (Cary, NC).

Gene expression pattern analysis

Gene expression analysis of the raw neutrophil data revealed the existence of different expression patterns for genes, such as unimodal distribution of expression (Fig. 2A), bimodal distribution with or without the existence of low expressors (Ct₃₅–Ct₃₉) (Figs. 2B and 2C), and trimodal distribution (Fig. 2D). The existence of non-expressing cells poses the problem of how to define these data points. One approach is to assign all such values the maximum Ct of 40, but this assumes that these data represent true missing expression; they could also result from technical errors due to failed PCR reactions. If the latter is the case, apparent bimodality with on/off expression patterns would in reality represent unimodal distribution with missing data points being technical artifacts instead of biologically relevant information. An alternative approach for addressing this issue is to look at patterns of missing data within the sets. If missing data points from the same genes tend to correlate within the cells, the cause is likely to be biological, suggesting that the populations contain cellular subgroups. In order to determine whether the existence of genes with bimodal expression patterns signaled the existence of cellular subclasses, the data was clustered based on shared gene expression patterns. Clustering showed that for neutrophils, bimodal genes exhibiting on/off pattern tended to be off in the same cells, although they clustered together with unimodal genes implying that the differential expression between cell types is not restricted to bimodality (Fig. 3A). Another potential cause for missing data points is low initial concentration of RNA in the sample, owing to inefficient RNA extraction, leading to complete loss of signal for the lowest-abundance genes that share the technical inefficiency. In order to address this, we controlled for overall abundance of RNA by normalizing our data sets.

Similarly to the neutrophils, the T lymphocyte data set contained genes exhibiting bimodal gene expression. As seen in neutrophils, T lymphocyte genes with bimodal on/off expression patterns also tended to be interspersed with unimodally expressed transcripts (Fig. 3B).

Detection of neutrophil subtypes

Hierarchical clustering was applied to the datasets using Ward’s method, which has been shown to discriminate clusters efficiently on gene expression datasets (Ferreira & Hitchcock, 2009; Ma & Zhang, 2012). Figure 4 shows the results of hierarchical clustering with nine different methods combining the three methods for data exclusion and three methods for data normalization. The color coding (purple, green and orange) shows the degree of concordance of clustering relative to the method based on supervised data removal with mean centering (top left). Employing exclusion with any of the three methods, followed by either mean centering or quantile normalization, three clusters of neutrophils were observed consistently, with notable separation of the orange, and most of the green, clusters from the purple one. Concordance, defined as the percentage of cells assigned to the same cluster, ranged from 75% to 100%, prima facie supporting the presence of three cell types in our samples.

However, when a Fluidigm array plate effect was fit to the standardized gene expression z-scores, only two major clusters were observed regardless of the data exclusion method (Fig. 4G), and concordance of the two-way classification of orange versus green/purple cells was perfect. This analysis implies that a plate effect caused the splitting of the large purple/green clusters observed with the mean-centering and quantile normalization methods. That is to say, very low abundance gene expression led to loss of signal on one of the plates, generating an artificial signature of co-regulation of some cells. However, the orange cluster remains robustly detected by all methods. We conclude that there are two main clusters of cell types in neutrophils. There is also a hint of a sub-type within the orange cells defined by differential expression of a half-dozen genes, but a larger sample will be required to validate this inference.

Hierarchical clustering verified the existence of cellular subgroups

Having compared methods for data exclusion and normalization, we opted to focus on the analysis method using a two standard deviation cutoff for exclusion with normalization by standardization of the genes (Fig. 4H). Hierarchical clustering revealed 2 major subclasses in both neutrophils (Fig. 5A) and T lymphocytes (Fig. 5B). The more clear definition of neutrophil subgroups, as compared to T lymphocytes, could be due to different levels of bimodality in the gene sets, such that more bimodality in the neutrophil data set gives rise to more distinct cellular subclasses. Alternatively, the two data sets could incorporate the same level of overall bimodality but differ in the level of co-variation of bimodally expressed genes. Since the expression of many genes on the T-cell array was too low to detect consistently, the analysis is based on fewer genes which also reduces the power to detect clusters.

More refined clustering of the T-cell data was also heavily impacted by the decision as to how to handle missing data. Including genes in the analysis according to the 2 standard deviation cutoff, setting missing data to null expression resulted in 6 clusters of cells irrespective of the data normalization procedure. In contrast, when missing data was assumed to be due to technical error and thus assigned the mean value for that transcript, the number of cellular subgroups observed after clustering differed: mean centering resulted in 2 large and 6 small clusters, quantile normalization in 7 clusters, and standardization of the genes in 6 clusters. The all inclusive method of data exclusion also resulted in differing numbers of cell clusters depending on the normalization method, with mean centering indicating 4 cellular subgroups, while quantile normalization showed five groups, and standardization of the gene resulted in six groups after hierarchical clustering by visual observation. Concordance of cell clustering across methods ranges between 70% and 80%, arguing that there are multiple cell states despite the high degree of heterogeneity. Concordance of the two-state clustering indicated in Fig. 5B was 95%.

Individual differences in donor representation in cellular subgroups

We next turned to analysis of differences in cellular abundance among donors, and asked whether cells from all donors were equally distributed among the observed clusters. The results show that the frequency of cells in each neutrophil cluster differed between donors (Tables 1A and 1B) with donor 3 having a significantly lower than expected proportion of cells in cluster A, whereas donor 4 has the inverse profile. Setting the number of clusters to 2 in the analysis of the standardized data following supervised normalization, the χ² value for differences in cell type abundance is 24.5 (p = 6 × 10⁻⁵, 4 degrees of freedom) (Table 1A).

As the intrinsic variability of T cell populations is greater than that of neutrophil populations, it is perhaps not surprising that we found donor-to-donor variability to be larger for T cells than for neutrophils. Compared to neutrophils, T cells had considerable variability in cell distribution between subgroups. The counts associated with the smallest subgroup were not large enough to establish whether the donors differ, but they do suggest divergence for the other clusters. Setting the number of clusters to 2 following data normalization, the χ² value for differences in cell type abundance is 36.8 (p = 7 × 10⁻⁷, 5 degrees of freedom) (Table 1C). Sampling of more cells in more donors will be required to establish whether these differences correlate with physiological and immunological attributes of the individuals.