PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of bioinfoLink to Publisher's site
 
Bioinformatics. 2013 February 15; 29(4): 461–467.
Published online 2012 December 24. doi:  10.1093/bioinformatics/bts714
PMCID: PMC3570210

Data exploration, quality control and testing in single-cell qPCR-based gene expression experiments

Abstract

Motivation: Cell populations are never truly homogeneous; individual cells exist in biochemical states that define functional differences between them. New technology based on microfluidic arrays combined with multiplexed quantitative polymerase chain reactions now enables high-throughput single-cell gene expression measurement, allowing assessment of cellular heterogeneity. However, few analytic tools have been developed specifically for the statistical and analytical challenges of single-cell quantitative polymerase chain reactions data.

Results: We present a statistical framework for the exploration, quality control and analysis of single-cell gene expression data from microfluidic arrays. We assess accuracy and within-sample heterogeneity of single-cell expression and develop quality control criteria to filter unreliable cell measurements. We propose a statistical model accounting for the fact that genes at the single-cell level can be on (and a continuous expression measure is recorded) or dichotomously off (and the recorded expression is zero). Based on this model, we derive a combined likelihood ratio test for differential expression that incorporates both the discrete and continuous components. Using an experiment that examines treatment-specific changes in expression, we show that this combined test is more powerful than either the continuous or dichotomous component in isolation, or a t-test on the zero-inflated data. Although developed for measurements from a specific platform (Fluidigm), these tools are generalizable to other multi-parametric measures over large numbers of events.

Availability: All results presented here were obtained using the SingleCellAssay R package available on GitHub (http://github.com/RGLab/SingleCellAssay).

Contact: rgottard/at/fhcrc.org

Supplementary information: Supplementary data are available at Bioinformatics online.

1 INTRODUCTION

The development of fluorescence-based flow cytometry (FCM) revolutionized single-cell analysis. Although populations of cells sorted by FCM using surface markers may appear monolithic, mRNA expression of specific genes within these cells can be heterogeneous (Dalerba et al., 2011) and could further discriminate cell subsets. On the other hand, classical gene expression experiments [microarrays, RNA-seq, quantitative polymerase chain reactions (qPCR)] richly characterize a cellular population but at the cost of reporting a summation of expression from many individual cells. Recent advances in microfluidic technology now permit performing thousands of PCRs in a single device, enabling gene expression measurements at the single-cell level across hundreds of cells and genes (Kalisky and Quake, 2011). This provides a technology that probes the stochastic nature of biochemical processes, resulting in relatively large cell-to-cell expression variability.

This heterogeneity may carry important information; thus, single-cell expression data should not be analysed in the same fashion as population-level data. At the scale of a single cell, biological variability (the object of interest) and technical variability (a nuisance factor) are often of the same magnitude, making it difficult to distinguish between the two. No expression (i.e. the gene is off) may be detected in individual cells owing to real biological effects, resulting in zero-inflation of otherwise continuous measures. These features of single-cell data require special attention during analysis.

Here, we focus on the reverse-transcriptase qPCR (rt-qPCR)-based Fluidigm (San Francisco, CA) single-cell gene expression assay, which provides simultaneous measurements of up to 96 genes on mRNA sources as minute as a single cell. In traditional rt-qPCR, despite careful measurement of input cDNA concentrations, differences in starting material below the limit of detection require correction for reliable results (Vandesompele et al., 2002). Subtraction of internal control genes, or averages thereof is typically used (e.g. the An external file that holds a picture, illustration, etc.
Object name is bts714i1.jpg-Ct method), and results are often reported in fold increase per cell (Schmittgen and Livak, 2008). In array-based gene expression, differences in hybridization and washing of non-specific DNA between chips require additional correction.

Such normalization schemes are not directly applicable in single-cell gene expression experiments, nor is it obvious that they are needed. For single cells, the individual cell is the atomic unit of normalization and the amount of starting material naturally measured in number of cells per reaction. Even if one attempted direct application of traditional normalization approaches, the dichotomous nature of single-cell expression hinders their use.

Nonetheless, it is important to test for and address any technical biases. We present a filtering approach for removing outlying measurements at the single-cell level that accounts for the dichotomous nature of the data. Using concordance measures derived from three datasets where gene expression was measured at the single-cell and 100-cell levels, we show that classical rt-qPCR type normalization is not necessary with single-cell multiplexed PCR data, and that our filtering step removes technical artifacts that most severely impact quantitation.

A typical goal of gene expression experiments is to search for differential expression across groups. The zero-inflation of expression in Fluidigm introduces problems for testing differential representation of cell subsets characterized by expression patterns, as well. Traditional tests of differential expression such as the t-test or other approaches based on normality are likely inappropriate for zero-inflated data (Gottardo et al., 2006; Smyth, 2004). Approaches to this problem have varied. Powell et al. (2012) used a winsorized z-transformation of the expression values and then treated them as continuous. Glotzbach et al. (2011) used the non-parametric, Kolmorgov–Smirnov test for differences in distribution to find differentially expressed genes after winsorizing. Flatz et al. (2011) dichotomized the expression and worked with the binary trait. Of these authors, only Flatz et al. (2011) and Glotzbach et al. (2011) made use of formal tests of differential expression. However, as we will see later, both the continuous and discrete parts of the measurements are informative for differential expression and should be used. A parametric test allows directions of difference to be assessed.

Here, we propose a discrete/continuous model for single-cell expression data based on a mixture of a point mass at zero and a log-normal distribution. Using this model, we derive a likelihood ratio test (LRT) that can simultaneously test for changes in mean expression (conditional on the gene being expressed) and in the percentage of expressed cells.

2 METHODS

2.1 Datasets and notations

We use three Fluidigm single-cell gene expression datasets described later in the text. We offer a brief overview of the assay technology used for our data. Desired cells (e.g. antigen-specific CD8+ T cells) are selected and lysed, and a cDNA library is generated through rt-qPCR. A short (c. 15 cycle), multiplexed pre-amplification selects and enriches for the desired genes. These products are loaded onto the Fluidigm chip, and gene-specific primers are added for single-cell gene expression quantitation. For the data presented here, we used a An external file that holds a picture, illustration, etc.
Object name is bts714i2.jpg format plate, i.e. 96 genes across 96 cells. The design of the chip generates each combination of the 96 genes and 96 enriched cDNA libraries producing 9216 separate PCR reactions. After each cycle, the fluorescence is read. The cycle (or interpolated fraction thereof) at which the fluorescence crosses a pre-determined threshold is recorded, defined as the ‘ct’ value. For all datasets considered here, primers were chosen to have An external file that holds a picture, illustration, etc.
Object name is bts714i3.jpg amplification efficiency.

Dataset A: Twenty-eight An external file that holds a picture, illustration, etc.
Object name is bts714i4.jpg format plates of CMV- or HIV-specific CD8+ single cell T cells were isolated from 16 individuals. The donors’ cells were stimulated with one of four tetramers. Cells were sorted immediately after tetramer incubation (‘unstimulated’) or after 3 hours of exposure (‘stimulated’). Approximately 90 individual cells were measured for each patient-stimulation combination (‘unit’).

Dataset B: Ten subjects were considered, and ~180 activated CD4+ memory T cells were sorted per subject, with each subject crossed between two arrays.

Dataset C: Two subjects were considered. Fluorescent staining of CD4+ T cells allowed cytometric sorting into CD154+/− sub-populations. Approximately 40 cells were sorted per sub-population per subject across three arrays.

Additionally, for each individual and treatment within each dataset, aggregates of 100 cells (i.e. 100 cells per well on the array) were isolated and assayed by Fluidigm technology. The expression measured in these 100-cell aggregates, after dividing by 100, provides a ‘biological’ average of expression per cell and can be compared with an in silico average of the single-cell measurements. The concordance between these two averages serves as a measure of experimental fidelity (Lin, 1989).

Notations: The standard assumptions of qPCR-based assays apply to the Fluidigm technology, namely that the cycle threshold (ct) is inversely proportional to the log of fluorescence. The fluorescence is directly proportional to the starting concentration of mRNA (Higuchi et al., 1992; Karlen et al., 2007). The Fluidigm instrument returns the cycle threshold (ct); however, we find it more useful to work with the complement of ct, which we define as the expression threshold (et)

equation image

where An external file that holds a picture, illustration, etc.
Object name is bts714i5.jpg is the maximum number of cycles used, 40 in our case. Assuming all reactions are in the exponential amplification phase, this quantity should be directly proportional to the log-abundance of mRNA, plus an intercept term corresponding to the number of cycles it takes for the minimally detectable quantity of mRNA to cross threshold. If the fluorescence does not cross the threshold after 40 cycles, then the Fluidigm instrument records a value of N/A, and we say that the gene is not detected. As we will see in the Section 3, detected genes typically have a value of ct much less than An external file that holds a picture, illustration, etc.
Object name is bts714i6.jpg suggesting that undetected genes might be regarded as unexpressed genes. This assumption is supported by the idea that transcription of mRNA is thought to occur in bursts of activity (Kaufmann and van Oudenaarden, 2007; Levsky et al., 2002), followed by quiescence. Other authors have noted this feature in single-cell expression studies as well (Glotzbach et al., 2011). When looking at the concordance of the single-cell and 100-cell experiments, this assumption is reasonable and leads to better concordance than omitting the N/A values. As a consequence, we treat the undetected genes as unexpressed genes, and we set the corresponding et value to An external file that holds a picture, illustration, etc.
Object name is bts714i7.jpg so that the mRNA abundance is zero (i.e. An external file that holds a picture, illustration, etc.
Object name is bts714i8.jpg).

For a fixed sample or experimental unit, let us denote by An external file that holds a picture, illustration, etc.
Object name is bts714i9.jpg the expression threshold of well An external file that holds a picture, illustration, etc.
Object name is bts714i10.jpg and gene An external file that holds a picture, illustration, etc.
Object name is bts714i11.jpg, for An external file that holds a picture, illustration, etc.
Object name is bts714i12.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i13.jpg. This results in a matrix of An external file that holds a picture, illustration, etc.
Object name is bts714i14.jpgbased expression values, An external file that holds a picture, illustration, etc.
Object name is bts714i15.jpg, just as in array-based gene expression. Similarly, we will denote by An external file that holds a picture, illustration, etc.
Object name is bts714i16.jpg the matrix of untransformed expression values, where An external file that holds a picture, illustration, etc.
Object name is bts714i17.jpg. Usually, a well contains one cell, but the Fluidigm technology can be used with multiple cells per well to quantify the gene expression of a mixture of cells. As a consequence, we prefer to use the term ‘well’ instead of ‘cell’. In the three datasets used here, wells will contain either 1 or 100 cells. Finally, several biological units are typically measured in an experiment, and in this case, we will use an extra index An external file that holds a picture, illustration, etc.
Object name is bts714i18.jpg to refer to the biological units.

2.2 A model for single-cell expression

As described previously, for a given cell, a gene can be defined as on (i.e. a positive et value is recorded) or as off (i.e. the gene is undetected and An external file that holds a picture, illustration, etc.
Object name is bts714i19.jpg). To simplify our model, we will denote by An external file that holds a picture, illustration, etc.
Object name is bts714i20.jpg the indicator variable equal to one if the gene An external file that holds a picture, illustration, etc.
Object name is bts714i21.jpg is expressed in well An external file that holds a picture, illustration, etc.
Object name is bts714i22.jpg and zero otherwise. Following classical statistical conventions, we use upper cases to denote the random variables and lower cases to denote the values taken by these random variables. Using these notations, we introduce the following model of single-cell expression

equation image
(1)

equation image
(2)

equation image
(3)

where An external file that holds a picture, illustration, etc.
Object name is bts714i23.jpg denotes a point mass at zero, An external file that holds a picture, illustration, etc.
Object name is bts714i24.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i25.jpg are the An external file that holds a picture, illustration, etc.
Object name is bts714i26.jpg-based mean and variance expression-level parameters conditional on the gene being expressed (i.e. An external file that holds a picture, illustration, etc.
Object name is bts714i27.jpg), and An external file that holds a picture, illustration, etc.
Object name is bts714i28.jpg is the frequency of expression of gene An external file that holds a picture, illustration, etc.
Object name is bts714i29.jpg across all cells. In the datasets considered here, the frequency of expression greatly varies across genes from 0 to 0.99 with a median value of An external file that holds a picture, illustration, etc.
Object name is bts714i30.jpg ~0.1 (see Supplementary Fig. S1). Assuming a log-Normal model for An external file that holds a picture, illustration, etc.
Object name is bts714i31.jpg is equivalent to modeling An external file that holds a picture, illustration, etc.
Object name is bts714i32.jpg as normally distributed. The empirical distribution of the data (Fig. 1 and Supplementary Figs S8–S10) motivates our selection of a log-normal distribution and follows observations of previous authors (Bengtsson et al., 2005).

Fig. 1.
Histogram and theoretical (normal) distribution of An external file that holds a picture, illustration, etc.
Object name is bts714i33.jpg for single-cell (left, light gray) and 100-cell experiments (right, dark gray). Genes FASLG, IFN- An external file that holds a picture, illustration, etc.
Object name is bts714i34.jpg, BIRC3 and CD69 are depicted. The frequency expression of each gene in the single-cell experiments ...

Thus, in a particular gene, three parameters characterize the expression distribution: An external file that holds a picture, illustration, etc.
Object name is bts714i36.jpg, the mean and standard deviation of the An external file that holds a picture, illustration, etc.
Object name is bts714i37.jpg, and An external file that holds a picture, illustration, etc.
Object name is bts714i38.jpg, the Bernoulli probability of expression.

2.3 Quality control and filtering

The Fluidigm assay is sensitive, and owing to the exponential amplification of starting mRNA, even minute contamination can render a measurement unreliable. Similarly, variation in cell preparation can have significant impact on the resulting experiment and data, such as unintentional empty wells, which would distort estimates of An external file that holds a picture, illustration, etc.
Object name is bts714i39.jpg. This suggests identifying, and possibly removing outliers before conducting further analysis. We examine both the discrete component An external file that holds a picture, illustration, etc.
Object name is bts714i40.jpg and the continuous component An external file that holds a picture, illustration, etc.
Object name is bts714i41.jpg in screening for outliers. We define the robust z-transformed positive expression value as

equation image

where the median and median absolute deviation (MAD) are calculated, for a given gene, over expressed cells (i.e. An external file that holds a picture, illustration, etc.
Object name is bts714i42.jpg), and An external file that holds a picture, illustration, etc.
Object name is bts714i43.jpg is a scaling constant that gives the standard deviation in terms of the MAD for the normal distribution. Next, let An external file that holds a picture, illustration, etc.
Object name is bts714i44.jpg be the Bernoulli variance-stabilizing transformation of the proportion of genes expressed in well An external file that holds a picture, illustration, etc.
Object name is bts714i45.jpg. Then, we define the robust z-transformed fraction as

equation image

where the median, MAD and An external file that holds a picture, illustration, etc.
Object name is bts714i46.jpg are as defined previously. This leads to the following steps for filtering:

  1. Remove null cells with no detected genes, i.e. An external file that holds a picture, illustration, etc.
Object name is bts714i47.jpg, for all An external file that holds a picture, illustration, etc.
Object name is bts714i48.jpg.
  2. Pick threshold for z filtering (An external file that holds a picture, illustration, etc.
Object name is bts714i49.jpg); threshold for An external file that holds a picture, illustration, etc.
Object name is bts714i50.jpg filtering (An external file that holds a picture, illustration, etc.
Object name is bts714i51.jpg).
  3. Calculate An external file that holds a picture, illustration, etc.
Object name is bts714i52.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i53.jpg
  4. Remove wells in which genes have An external file that holds a picture, illustration, etc.
Object name is bts714i54.jpg OR An external file that holds a picture, illustration, etc.
Object name is bts714i55.jpg.

Step 1 removes wells where no cells were loaded, and thus all measured expression values are null. It is important to perform this step first to prevent break-down in the median and MAD estimates for the zeta values in experiments with many amplification or FCM failures. Finally, step 4 removes unreliable wells that either have an extreme proportion of expression or extreme cell An external file that holds a picture, illustration, etc.
Object name is bts714i57.jpg gene expression values. The thresholds An external file that holds a picture, illustration, etc.
Object name is bts714i58.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i59.jpg control the tolerance to outliers; therefore, typical advice for outlier thresholding applies. Biological replicates, such as the 100-cell replicates described in Section 2.1, permit the assessment of intra-class deviance, and then the thresholds can be selected to minimize this deviance. We present such a calculation in the Supplementary Material. Using this approach, we find that picking An external file that holds a picture, illustration, etc.
Object name is bts714i60.jpg works well for the datasets we consider here, see Section 3.

2.4 Testing for ET differences between experimental groups

One typical goal of gene expression analysis is to test for difference in expression patterns between experimental units. Here, we focus on testing differential gene expression between two paired-biological units, e.g. before and after stimulation. Our framework should be generalizable to other types of situations, see Section 4. The classical test for changes in mean for samples with continuous measurements is the An external file that holds a picture, illustration, etc.
Object name is bts714i61.jpg-test. Conversely, if only a change in An external file that holds a picture, illustration, etc.
Object name is bts714i62.jpg were of interest, then a contingency table test (Chi-square, Fisher’s exact or Bernoulli likelihood ratio) is appropriate. However, in our case, we would like to test for a change in An external file that holds a picture, illustration, etc.
Object name is bts714i63.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i64.jpg simultaneously, as both could be biologically relevant. Formally, we wish to test

equation image

versus the alternative

equation image

This can be accomplished using an LRT that would simultaneously test for differences in means or proportions of expression.

Suppose that An external file that holds a picture, illustration, etc.
Object name is bts714i65.jpg wells are assayed in each unit, though the unbalanced case (An external file that holds a picture, illustration, etc.
Object name is bts714i66.jpg) would be treated similarly with obvious changes of notation. Based on (1), the likelihood function for one gene across two biological units, omitting the gene index An external file that holds a picture, illustration, etc.
Object name is bts714i67.jpg for clarity, is given by

equation image
(4)

where An external file that holds a picture, illustration, etc.
Object name is bts714i68.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i69.jpg are the vectors of observations for the gene across the two groups, An external file that holds a picture, illustration, etc.
Object name is bts714i70.jpg is the vector of unknown parameters, An external file that holds a picture, illustration, etc.
Object name is bts714i71.jpg is the set of cells expressing the gene in group An external file that holds a picture, illustration, etc.
Object name is bts714i72.jpg (i.e. An external file that holds a picture, illustration, etc.
Object name is bts714i73.jpg}), An external file that holds a picture, illustration, etc.
Object name is bts714i74.jpg is the number of cells expressing the gene in group An external file that holds a picture, illustration, etc.
Object name is bts714i75.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i76.jpg is the density function of the log-normal distribution with parameters An external file that holds a picture, illustration, etc.
Object name is bts714i77.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i78.jpg. The LRT statistic An external file that holds a picture, illustration, etc.
Object name is bts714i79.jpg is then defined as the ratio of the null and alternative likelihoods obtained by replacing the unknown parameters with their null and alternative maximum likelihood estimates. Detailed derivations of the likelihood function and the LRT statistics are described in Supplementary Material.

An interesting observation is that the likelihood function given by (4) is the product of the Bernoulli likelihood for all cells and the log-normal likelihood for the expressed cells. It follows that the log-LRT statistic decomposes as a sum of a Bernoulli log-LRT test statistic and a log-normal log-LRT test statistic, as each component can be maximized independently. It thus combines the two sources of information in a natural way, and this decomposition allows post-hoc assessment of which of the component(s) drive the detected difference by simply comparing the magnitude of the two log-LRTs. In Section 3, we will show that our combined LRT test is more powerful than the Bernoulli or log-normal tests alone.

Applying classical asymptotic results about LRTs (Wilks, 1938), An external file that holds a picture, illustration, etc.
Object name is bts714i82.jpg converges to a An external file that holds a picture, illustration, etc.
Object name is bts714i83.jpg distribution with two degrees of freedom under An external file that holds a picture, illustration, etc.
Object name is bts714i84.jpg. Some care is warranted in invoking this asymptotic result, as even for large An external file that holds a picture, illustration, etc.
Object name is bts714i85.jpg, the sample size for the log-normal LRT will be An external file that holds a picture, illustration, etc.
Object name is bts714i86.jpg. We show in Supplementary Figs S2 and S3 that the An external file that holds a picture, illustration, etc.
Object name is bts714i87.jpg convergence is adequate for An external file that holds a picture, illustration, etc.
Object name is bts714i88.jpg even under departures from normality. Below this value, it is possible to derive the null distribution of this statistic through permutation procedures as is commonly done for microarray data (Ge et al., 2003). This proviso applies similarly for purpose of power calculations; hence, one may wish to conduct these through simulation.

3 RESULTS

3.1 Distributional assumptions

In Figure 1, we observe good agreement between the empirical distributions of positive et values and their postulated normal distribution for four genes in dataset A. This confirms that a log-normal model for the positive expression level, An external file that holds a picture, illustration, etc.
Object name is bts714i89.jpg, is appropriate. Even cells in the lowest quantiles of et (and lowest quantiles of expression) still have expression far away from the bound at 0, suggesting that undetected genes represent cells with null or negligible RNA abundance. It is also noteworthy that the difference between the means (shown as a heavy, vertical line) of the 100-cell replicates and single-cell replicates is approximately An external file that holds a picture, illustration, etc.
Object name is bts714i90.jpg cycles, where An external file that holds a picture, illustration, etc.
Object name is bts714i91.jpg is the expression frequency of gene An external file that holds a picture, illustration, etc.
Object name is bts714i92.jpg in the single-cell experiments. As such, in genes with An external file that holds a picture, illustration, etc.
Object name is bts714i93.jpg, such as FASLG, this difference between means is smaller than genes with An external file that holds a picture, illustration, etc.
Object name is bts714i94.jpg. As we will see the next section, inclusion of the unexpressed cells (An external file that holds a picture, illustration, etc.
Object name is bts714i95.jpg) is important to accurately relate the expression level of the single-cell experiments to the 100-cell experiments.

3.2 Concordance between 100-cell and single-cell experiments

The 100-cell aggregates (see Section 2.1) allows us to assess the accuracy and reliability of our single-cell experiments by comparing this in vitro 100-cell expression to an in silico estimate obtained by averaging the expression of 100 single-cell measurements. The in silico average of signal in a gene An external file that holds a picture, illustration, etc.
Object name is bts714i96.jpg and unit An external file that holds a picture, illustration, etc.
Object name is bts714i97.jpg from 100 single-cell wells is An external file that holds a picture, illustration, etc.
Object name is bts714i98.jpg where An external file that holds a picture, illustration, etc.
Object name is bts714i99.jpg is the expression measurement of gene An external file that holds a picture, illustration, etc.
Object name is bts714i100.jpg in cell An external file that holds a picture, illustration, etc.
Object name is bts714i101.jpg and unit An external file that holds a picture, illustration, etc.
Object name is bts714i102.jpg. We compare this with the in vitro ‘average’ of signal from a 100-cell aggregate. In this case, we just use the expression of a gene unit and divide by the number of cells (100).

The concordance here is assessed both visually by plotting An external file that holds a picture, illustration, etc.
Object name is bts714i103.jpg versus An external file that holds a picture, illustration, etc.
Object name is bts714i104.jpg (Fig. 2) and by calculating the concordance correlation coefficient (An external file that holds a picture, illustration, etc.
Object name is bts714i112.jpg) between the two variables, which is often used to quantify reproducibility (Lin, 1989). The shifted log transformation allows visualization of both the discrete and continuous components while being on the et scale.

Fig. 2.
Concordance between 100 cell An external file that holds a picture, illustration, etc.
Object name is bts714i105.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i106.jpg, the in silico average of single-cell wells for datasets A, B and C. In the top row, wells with An external file that holds a picture, illustration, etc.
Object name is bts714i107.jpg are included and treated as exact zeroes. In the middle row, they are excluded, resulting in a clear lack of concordance. ...

We first use this concordance experiment to test whether wells that do not cross the fluorescence threshold after An external file that holds a picture, illustration, etc.
Object name is bts714i113.jpg should be treated as exact zeroes or missing values. If we suppose that An external file that holds a picture, illustration, etc.
Object name is bts714i114.jpg implies an assay failure and the measurement should be discarded, we would simply compute the single-cell average over expressed cells, i.e. An external file that holds a picture, illustration, etc.
Object name is bts714i115.jpg. Figure 2 demonstrates good concordance between the 100-cell and single-cell experiments when the undetected genes are treated as zeros. However, this is not the case when the zeros are treated as missing values.

3.3 Filtering outlying cells

In addition to the concordance measure An external file that holds a picture, illustration, etc.
Object name is bts714i116.jpg, we use another goodness-of-fit measure to optimize our filtering parameters An external file that holds a picture, illustration, etc.
Object name is bts714i117.jpg, An external file that holds a picture, illustration, etc.
Object name is bts714i118.jpg defined by,

equation image
(5)

where An external file that holds a picture, illustration, etc.
Object name is bts714i119.jpg is the number of positive wells for gene An external file that holds a picture, illustration, etc.
Object name is bts714i120.jpg in unit An external file that holds a picture, illustration, etc.
Object name is bts714i121.jpg in the single-cell experiments. For a particular gene and unit, the An external file that holds a picture, illustration, etc.
Object name is bts714i122.jpg decreases as we lower the filtering threshold and extreme values are filtered. Eventually, so many cells are removed that there is zero expression (and a large deviance) for the in silico estimate. Thus, we wish to find a set of values for the filtering parameters that would lead to the lowest An external file that holds a picture, illustration, etc.
Object name is bts714i123.jpg measure across the three datasets used here. The addition of the scaling factor An external file that holds a picture, illustration, etc.
Object name is bts714i124.jpg gives higher weight to combinations with more ex ante positive observations so that the contribution to the sum of squares would be smaller in gene An external file that holds a picture, illustration, etc.
Object name is bts714i125.jpg unit combinations that have fewer expressed cells. The factor An external file that holds a picture, illustration, etc.
Object name is bts714i126.jpg can also be interpreted as the scaling factor for the variance of the mean over positive observations. Finally, the An external file that holds a picture, illustration, etc.
Object name is bts714i127.jpg is computed on the An external file that holds a picture, illustration, etc.
Object name is bts714i128.jpg scale to reduce the effect of extreme outliers.

When 100-cell aggregates are available, one can optimize the filter parameters An external file that holds a picture, illustration, etc.
Object name is bts714i129.jpg by minimizing the An external file that holds a picture, illustration, etc.
Object name is bts714i130.jpg over possible combinations. In our case, we found that setting An external file that holds a picture, illustration, etc.
Object name is bts714i131.jpg achieves the best reduction in An external file that holds a picture, illustration, etc.
Object name is bts714i132.jpg across the three datasets explored here (Supplementary Figs S4–S6 and Supplementary Table S1). Using these values, our filtering criteria moderately improve the concordance between the single-cell and 100-cell experiments in two of the datasets but dramatically improve (decrease) the weighted sum of squares. This improvement is evident graphically, as the per unit averages of et of multiple genes move toward the diagonal.

Beside improving An external file that holds a picture, illustration, etc.
Object name is bts714i133.jpg and generally improving An external file that holds a picture, illustration, etc.
Object name is bts714i134.jpg, we explore the effect of filtering on detection of control genes in the Supplementary Material (Supplementary Table S2).

3.4 Normalization and housekeeping genes

Other authors have noted that ‘the gene transcript number is ideally standardized to the number of cells’ (Vandesompele et al., 2002), which is the case with gene expression from sorted cells. Therefore, it is not entirely a surprise that we find little evidence for housekeeping genes providing useful normalization here. For a housekeeper to have good validity, it should have high cross-correlation with other housekeeping genes. This is not the case for housekeepers GAPDH and POLR2A, which in dataset A, in linear regression, have an An external file that holds a picture, illustration, etc.
Object name is bts714i135.jpg. In Supplementary Figure S7, we observe in scatter plots of housekeepers’ et that the correlation drops even further (to an An external file that holds a picture, illustration, etc.
Object name is bts714i136.jpg) after filtering outlying cells (see previous section). As the correlation between housekeepers is present primarily in cells we suspect suffered from technical error, we find little utility in normalization schemes. In fact, the use of housekeeping genes for normalization could even result in masking cellular artifacts that should be filtered.

3.5 An efficient test of differential expression for single cells

In dataset A, ~90 cells in each of 16 subjects were measured in unstimulated and stimulated states (see Section 2.1). This permits conducting a test for each gene in each subject for differences in An external file that holds a picture, illustration, etc.
Object name is bts714i137.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i138.jpg, as described in Section 2.4. We plot the number of discoveries at various false discovery rates (FDR) in Figure 3. The combined likelihood test produces the greatest number of discoveries over a wide range of FDR. For example, at an FDR of 1%, our combined test could detect more than 20 additional gene An external file that holds a picture, illustration, etc.
Object name is bts714i141.jpg unit changes across the four stimulations.

Fig. 3.
Number of discoveries (genes An external file that holds a picture, illustration, etc.
Object name is bts714i139.jpg units) versus FDR, by treatment, dataset A. The combined LRT is compared with a Bernoulli or normal-theory only LRT, as well as a t-test of the raw expression values (An external file that holds a picture, illustration, etc.
Object name is bts714i140.jpg scale), including zero measurements

Another feature of the combined LRT is its robustness to background gene frequency An external file that holds a picture, illustration, etc.
Object name is bts714i142.jpg. Of course, if An external file that holds a picture, illustration, etc.
Object name is bts714i143.jpg on average, then any test will be underpowered to detect group differences. But using only the continuous components amounts to “throwing away” data for genes with intermediate An external file that holds a picture, illustration, etc.
Object name is bts714i144.jpg. And similarly, using only the dichotomous component results in a test insensitive to differences in An external file that holds a picture, illustration, etc.
Object name is bts714i145.jpg in frequently expressed genes. This robustness to the An external file that holds a picture, illustration, etc.
Object name is bts714i146.jpg spectrum is shown in Figure 4 in which An external file that holds a picture, illustration, etc.
Object name is bts714i150.jpg P-values are shown for the Bernoulli, normal and combined LRTs versus frequency of An external file that holds a picture, illustration, etc.
Object name is bts714i151.jpg.

Fig. 4.
An external file that holds a picture, illustration, etc.
Object name is bts714i147.jpg of tests (genes An external file that holds a picture, illustration, etc.
Object name is bts714i148.jpg units) versus frequencies of expression An external file that holds a picture, illustration, etc.
Object name is bts714i149.jpg of the genes. The Bernoulli, normal-theory and combined LRTs are plotted. Asterisk indicates test is different from the combined test at 5% significance in a Wilcoxon signed-rank test

A total of 65 genes were detected at an FDR of 1% in at least one individual. We define An external file that holds a picture, illustration, etc.
Object name is bts714i152.jpg as the negative An external file that holds a picture, illustration, etc.
Object name is bts714i153.jpg P-value times an indicator variable, which is positive when stimulated groups have greater expression, and negative otherwise. Figure 5 plots a heatmap of signed An external file that holds a picture, illustration, etc.
Object name is bts714i155.jpg P-values. The selected genes are in clustered rows; the 16 individuals are arranged in columns by stimulation. The color above each column indicates which of the four antigen stimulations the individual received. From this, it is clear that genes cluster into upregulated and downregulated modules, and that there is much individual variability in response. Some genes appear to have stronger responses to particular antigens, such as the response to CMV (red and purple columns) in FASLG and CLEC2B.

Fig. 5.
Heatmap of signed An external file that holds a picture, illustration, etc.
Object name is bts714i154.jpg for selected genes (rows, see main text) and all 16 individuals (columns). The color above each column indicates the antigen stimulation applied to the cells; thus, individuals are randomly arranged in each antigen block. Red and purple ...

4 CONCLUSIONS

Current approaches for analysis of single-cell assays have incompletely used the salient features of the experiment, and the resulting inference can be suboptimal. In this article, we have presented a framework for data exploration, quality control and testing for differential expression using single-cell data. Our comparison of 100-cell and single-cell measurements shows that undetected genes in an assay should be treated as effective ‘zeroes’. Both the discrete, zero-inflated portion and continuous portion of single-cell expression data are meaningful for detecting outliers. Moreover, differences in either could be of biological interest; therefore, it is desirable to combine evidence from both to detect changes in expression. Our LRT allows just that.

Although we have suggested default parameters for the filtering of outliers, informed from several datasets, our defaults are likely conservative. They are 3–4 times larger than the most substantial difference in expression between experimental groups we observed. Acquiring forms of ground-truth besides ‘bulk’ experiments (in our case, 100-cell aggregates) could allow forming tighter bounds. As in any outlier-based filtering procedure, it is desirable to tune for the problem at hand. The thresholds An external file that holds a picture, illustration, etc.
Object name is bts714i156.jpg and An external file that holds a picture, illustration, etc.
Object name is bts714i157.jpg should be different when eliminating potential technical error is of greatest concern than when one is most interested in searching for biological heterogeneity.

In this article, we have used the An external file that holds a picture, illustration, etc.
Object name is bts714i158.jpg asymptotic distribution of the LRT to compute P-values and assess significance. This approximation is relatively accurate and robust to the distributional form of An external file that holds a picture, illustration, etc.
Object name is bts714i159.jpg when the expected number of expressed cells is greater than 8 (see Supplementary Material). Otherwise, approximating the null distribution using permutations as in Ge et al., 2003 is more appropriate.

Further work, incorporating a mixed-effects model to our LRT, could extend its applicability. The test outlined in this article may not be appropriate in cases where traits of interest are not blocked within individuals (e.g. comparing between phenotypes like HIV+ versus HIV−). In this case, one wishes to identify gene expression changes across groups, despite high individual-to-individual heterogeneity. By modeling the mean and proportion of expression as common across groups and adding specific random effects for between-individual variability, our model could be extended to address such experimental questions as well.

Single-cell gene expressions assays have already been shown to be useful in multiple studies and will become even more routine once sequencing at the single-cell level becomes practical (Ramskold et al., 2012; Varadarajan et al., 2011). As a consequence, the development of effective statistical methods to analyse such data is becoming increasingly important. This article offers a coherent framework for researchers using this nascent technology.

Funding: Intramural Research Program of the National Institute of Allergy and Infectious Diseases, National Institute of Health and the Collaboration for AIDS Vaccine Discovery [#38650]; National Institute of Health [U19 AI089986-01, R01 EB008400 to R.G., G.F. and A.M.]; and the Bill & Melinda Gates Foundation [#OPP1032325].

Conflict of Interest: none declared.

Supplementary Material

Supplementary Data:

REFERENCES

  • Bengtsson M, et al. Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels. Genome Res. 2005;15:1388–1392. [PubMed]
  • Dalerba P, et al. Single-cell dissection of transcriptional heterogeneity in human colon tumors. Nat. Biotechnol. 2011;29:1120–1127. [PMC free article] [PubMed]
  • Flatz L, et al. Single-cell gene-expression profiling reveals qualitatively distinct CD8 T cells elicited by different gene-based vaccines. Proc. Natl Acad. Sci. USA. 2011;108:5724–5729. [PubMed]
  • Ge Y, et al. Resampling-based multiple testing for microarray data analysis. TEST. 2003;12:1–77.
  • Glotzbach JP, et al. An information theoretic, microfluidic-based single cell analysis permits identification of subpopulations among putatively homogeneous stem cells. PLoS One. 2011;6:e21211. [PMC free article] [PubMed]
  • Gottardo R, et al. Bayesian robust inference for differential gene expression in microarrays with multiple samples. Biometrics. 2006;62:10–18. [PubMed]
  • Higuchi R, et al. Simultaneous amplification and detection of specific DNA sequences. Biotechnology. 1992;10:413–417. [PubMed]
  • Kalisky T, Quake SR. Single-cell genomics. Nat. Meth. 2011;8:311–314. [PubMed]
  • Karlen Y, et al. Statistical significance of quantitative PCR. BMC Bioinformatics. 2007;8:131. [PMC free article] [PubMed]
  • Kaufmann BB, van Oudenaarden A. Stochastic gene expression: from single molecules to the proteome. Curr. Opin. Genet. Dev. 2007;17:107–112. [PubMed]
  • Levsky JM, et al. Single-cell gene expression profiling. Science. 2002;297:836–840. [PubMed]
  • Lin LI. A concordance correlation coefficient to evaluate reproducibility. Biometrics. 1989;45:255–268. [PubMed]
  • Powell AA, et al. Single cell profiling of circulating tumor cells: transcriptional heterogeneity and diversity from breast cancer cell lines. PLoS One. 2012;7:e33788. [PMC free article] [PubMed]
  • Ramskold D, et al. Full-length mRNA-Seq from single-cell levels of RNA and individual circulating tumor cells. Nat. Biotech. 2012;30:777–782. [PMC free article] [PubMed]
  • Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative CT method. Nat. Protocols. 2008;3:1101–1108. [PubMed]
  • Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 2004;3 article 3. [PubMed]
  • Vandesompele J, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3 RESEARCH0034. [PMC free article] [PubMed]
  • Varadarajan N, et al. A high-throughput single-cell analysis of human CD8 T cell functions reveals discordance for cytokine secretion and cytolysis. J. Clin. Invest. 2011;121:4322–4331. [PMC free article] [PubMed]
  • Wilks SS. The large-sample distribution of the likelihood ratio for testing composite hypotheses. Ann. Math. Stat. 1938;9:60–62.

Articles from Bioinformatics are provided here courtesy of Oxford University Press