Parameterizing cell-to-cell regulatory heterogeneities via stochastic transcriptional profiles

Probability Models for Heterogeneous Transcriptional Regulation.

To make reliable single-cell inferences, it was critical to start with simple probabilistic models of gene expression that were biologically accurate. Our method considers genes that exhibit two distinct regulatory states in a population of cells (19, 20, 27). Within each state, the cell-to-cell variation of expression was originally described by a lognormal distribution according to measurements of high-copy transcripts in single mammalian cells (28, 29). We tested whether there was a mechanistic foundation for using two lognormal subpopulations by examining a standard model of regulated gene expression (30, 31) (SI Appendix, Fig. S1). In this model, transcript levels per cell are determined by the kinetics of polymerase binding–unbinding, transcriptional elongation, and mRNA degradation. The relative magnitudes of the kinetic rate parameters together govern the steady-state distribution of transcripts in the population (32), allowing different regulatory states to be simulated.

For parameter sets where the probability of observing zero transcripts per cell was near zero, we found that the lognormal distribution was a suitable approximation of basal expression (Fig. 1A, blue). Parameter sets yielding median expression levels as low as 20 copies per cell showed only minor skewness in quantile–quantile comparisons with a lognormal distribution (Fig. 1A, blue inset). Starting with this basal distribution, we simulated a second cellular regulatory state by increasing the rate of polymerase binding, decreasing the rate of mRNA degradation, or both (Fig. 1A, orange). The apparent rate of polymerase binding increases upon recruitment by transcription factors that are expressed or activated heterogeneously within a population of cells (4, 5). Conversely, mRNA stabilization occurs posttranscriptionally through dedicated signal-transduction pathways activated by environmental stresses and proinflammatory stimuli (33). We found that either mechanism of gene up-regulation led to right-shifted distributions that were lognormal (Fig. 1A, orange inset). These simulations indicated that lognormal random variables were appropriate for the regulated expression of mid- to high-abundance transcripts.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Simple probability models capture regulated changes in gene expression. (A and B) Probability densities for the number of transcripts per cell were calculated using a kinetic model (30, 31) whose parameters led to basal regulatory states (blue) with either nonzero copies per cell in A or with near-zero copies per cell in B. The basal-regulatory states were compared with a lognormal distribution in A or an exponential distribution in B through a quantile–quantile (QQ) plot (blue insets). A second, induced regulatory state (orange) was created by increasing the polymerase binding rate (Lower Left), decreasing the transcript degradation rate (Upper Right), or both (Lower Right) in the model (SI Appendix, Fig. S1). All induced regulatory states were compared with a lognormal distribution through a QQ plot (orange insets).

One drawback of the lognormal distribution is that it has no support at zero copies (34), making it poor for capturing low-abundance genes that are completely silenced in some cells. To identify an alternative in this circumstance, we reconfigured the parameters of the model and defined a steady-state population where most cells would contain close to zero transcripts (Fig. 1B, blue). As noted before (32), this regulatory state was best captured by an exponential distribution (Fig. 1B, blue inset). Importantly, we found that when the kinetic parameters of a basal exponential state were modified to create a second right-shifted state (Fig. 1B, orange), the resulting distributions were lognormal (Fig. 1B, orange inset). Together, we conclude that the basic mechanisms of gene expression lead to steady-state distributions described by probability models that are relatively simple.

Deconvolution of Random 10-Cell Averages by Maximum-Likelihood Inference.

Our results from the gene-expression model suggested that single-cell regulatory heterogeneities could be depicted as a mixture of two lognormal states or as a mixture of an exponential state and a lognormal state (Fig. 1). Either mixture gives rise to a probability distribution characterized by four key parameters. The lognormal–lognormal (LN–LN) mixture requires the log-mean expression of the two regulatory states (µ₁ and µ₂), the log-SD for biological noise (σ), and the expression frequency (F) describing the probability that cells will occupy the higher regulatory state (step 1, Fig. 2A). (For simulations, the two lognormal states are assumed to share a common σ, but in practice we test whether inferences are improved when each lognormal state is allowed its own noise parameter; see below.) Thus, an LN–LN gene that is expressed at an approximately eightfold higher level in 20% of the population with a coefficient of variation (CV) of ∼50% is captured by µ₁ – µ₂ = 2, F = 20%, and σ = 0.48.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Inferring cellular subpopulations by maximum-likelihood inference of stochastic 10-cell samples from an LN–LN mixture of regulatory states. (A) The maximum-likelihood approach involves four steps. (1) A model of heterogeneous gene regulation is posed, where single cells are assumed to express genes at a low or high level with a common coefficient of variation for both subpopulations. The weight of each subpopulation is defined by the integrated single-cell expression distribution of the subpopulation (ϕ₁ and ϕ₂). The four parameters of the model are the log-mean expression for each subpopulation (μ₁ and μ₂), the proportion of cells in the high subpopulation (F), and the common log-SD of expression (σ). (2) Random 10-cell samples are collected to build a distribution of measurements for inference by the model. (3) Based on the model in step 1, a likelihood function is derived (Methods). (4) The likelihood function is then maximized by searching through the four parameters of the model to identify those that are most likely given the experimental observations. Additionally, we obtain measures of confidence for each estimated parameter (gray). (B–E) Accurate prediction of single-cell parameters from simulated 10-cell samples. Ten-cell expression data were simulated using different values of (B) μ₁, (C) μ₂, (D) F, and (E) σ (solid gray line) and then estimated by maximum likelihood. For each group of simulations, the remaining three model parameters were kept fixed at (C–E) μ₁ = 0.5, (B, D, and E) μ₂ = –2.5, (B–C and E) F = 22.5%, and (B–D) σ = 0.225 (dashed gray line). Solid gray line shows the one-to-one mapping of inferred-to-known parameter value. Off-diagonal plots are categorical plots of the fixed parameter estimates for a given value of the perturbed parameter. Graphs show the parameter estimates together with 95% maximum-likelihood CIs from three independent sets of 50 10-cell samples. (F–H) Prediction and validation of expression frequency for the heterogeneous transcript, SOD2, during breast-epithelial acinar morphogenesis. (F) Distribution of 81 10-cell qPCR measurements of SOD2 in outer ECM-attached epithelial cells and estimated subpopulation distribution (red line). Maximum-likelihood parameters (red box) are shown with 95% CI in brackets. (G) Representative RNA FISH image of endogenous SOD2 expression. A pseudocolored image (Left) is shown alongside a two-color image with DRAQ5 counterstain to visualize nuclei (Right). Arrows indicate ECM-attached cells with high SOD2 expression. (Scale bar, 20 μm.) (H) Percentage of cells showing high expression of SOD2 by RNA FISH (gray bar) compared with the maximum-likelihood estimate of F (white dashed line). RNA FISH data are shown as the mean percentage with 95% CI of ECM-attached cells showing high expression of SOD2. Maximum-likelihood predictions are shown as the parameter point estimate (white) with 95% CI (red).

The exponential–lognormal (EXP–LN) mixture also requires σ and F, along with a single log-mean for the high lognormal state (µ) and a rate parameter for the low exponential state (λ) (step 1, SI Appendix, Fig. S2). The rate parameter relates to how quickly the lower distribution decays above zero copies per cell. For example, a rate parameter of λ = 1 creates a distribution that has ∼37% overlap with that of a high lognormal state of µ = 0.5 and σ = 0.225 whereas λ = 3 causes only a ∼6.3% overlap. We modeled two distinct regulatory states by restricting the simulations to rate parameters that caused negligible overlap with the high lognormal state (λ > 3). Together, the different mixture models enabled us to simulate stochastic-profiling data by summing the expression of 10 cells randomly sampled from the appropriate two-state distribution (step 2, Fig. 2A and SI Appendix, Fig. S2).

To infer the most-likely parameters from a collection of random 10-cell samples, we derived maximum-likelihood estimators for the LN–LN and EXP–LN mixtures (Methods and SI Appendix, Methods). Maximum-likelihood estimation requires a defined probability density function (pdf). The stochastic-sampling pdf is the convolution of 10 binomial choices drawn from the two underlying distributions in the mixture (step 3, Fig. 2A and SI Appendix, Fig. S2). The pdf has a ≤11-modal shape where each mode corresponds to choosing 0–10 cells from the high regulatory state. The most-likely parameter combination was calculated by maximizing the likelihood function (Methods), yielding parameters with interval estimates that best explained the data (step 4, Fig. 2A and SI Appendix, Fig. S2). By performing this maximum-likelihood estimation, we could “invert” stochastic profiling data to infer single-cell characteristics from 10-cell samples.

Theoretical and Experimental Validation of Maximum-Likelihood Inference.

We evaluated the performance of our approach by using computational simulations of 10-cell samples with known distribution parameters. First, it was important to identify the minimum number of random samples needed to ensure accurate parameter estimation. Given hundreds to thousands of samples, we found that robust and accurate estimates were obtained for all model parameters irrespective of the mixture type (SI Appendix, Fig. S3 A and B). Conversely, with very few samples (∼20 or fewer), the convolved distributions were incompletely populated and our resulting estimates were highly uncertain and sometimes inaccurate for the LN–LN and EXP–LN mixtures. The transition between the two regimes occurred at 50–100 samples, which we defined as the approximate number of data points required for effective maximum-likelihood inference of single transcripts.

We next used simulations to identify the parameter ranges where maximum-likelihood inference makes accurate estimates of each regulatory state. Starting with the LN–LN mixture, we perturbed µ₁, µ₂, σ, or F individually while keeping the other three parameters fixed and simulated 50 random 10-cell samples. For a wide range of subpopulation log-means (µ₁ and µ₂), maximum-likelihood inference accurately inferred model parameters with negligible bias (Fig. 2 B and C). We also observed good performance when altering the expression frequency (F). Accuracy declined near F = 50%, when the two subpopulations offset one another and disguise as a distribution with large σ (Fig. 2D). Nevertheless, the estimation procedure still accurately and confidently captured ∼70% of the total parameter space (F = 0–35% over the range of 0–50%). For the log-SD (σ), performance declined only when this parameter was extremely high (Fig. 2E). Parameter estimates were accurate until σ reached ∼0.8, corresponding to a ∼95% CV that is higher than nearly all genes examined thus far (35, 36). None of the mixture parameters could be reliably inferred from higher-order moments of the 10-cell distributions, although low F or high σ correlated with a slight increase in skewness (SI Appendix, Fig. S4). These results indicated that maximum-likelihood inference could extract parameters that were otherwise inaccessible by descriptive statistics.

We repeated the simulations for the EXP–LN mixture and arrived at very similar conclusions. As long as λ and μ were large enough to prevent overlap of the two regulatory states, we found that parameter estimates were accurate, although the variance of inferred σ was somewhat higher than in the LN–LN mixture (SI Appendix, Fig. S5). Together, these simulations suggested that maximum-likelihood inference is able to deconvolve a wide range of regulatory heterogeneities from 10-cell samples.

To examine the accuracy of maximum-likelihood inference with real 10-cell samples, we focused on expression of the detoxifying enzyme superoxide dismutase 2 (SOD2) during breast-epithelial acinar morphogenesis. We used a culture model in which immortalized human breast epithelial cells are grown as single-cell clones in reconstituted basement-membrane ECM to form 3D organotypic spheroids (37, 38). Earlier stochastic-profiling studies of developing spheroids had suggested that there were two SOD2 regulatory states among the ECM-attached cells (27, 39). To apply maximum-likelihood inference, we deeply sampled SOD2 expression by quantitative PCR (qPCR) in 81 random samples of 10 ECM-attached cells (Fig. 2F, Left). Using these data, we maximized the likelihood of the LN–LN and EXP–LN models, as well as that of a relaxed LN–LN model, which allowed each regulatory state to have its own log-SD (σ₁ and σ₂). The three estimates were compared by using the Bayesian information criterion (BIC) score to calculate the quality of the fit relative to the number of inferred parameters (SI Appendix, Table S1). The best overall estimate was the mixture model that parameterized two distinct regulatory states with the lowest BIC score.

For the 10-cell measurements of SOD2, we found that the LN–LN mixture was slightly preferred over the EXP–LN mixture (Fig. 2F, Right and SI Appendix, Table S1). The inability to discriminate clearly between these two models was likely caused by the basal regulatory state, which could be described as an exponential distribution (λ = 46) or a lognormal distribution with a very small log-mean (µ₂ = –4.1) given the sampling data. Regardless, the two models predicted similar SOD2 expression frequencies among ECM-attached cells: 23% (13–33%) for the LN–LN mixture vs. 19% (12–27%) for the EXP–LN mixture. To determine the accuracy of this shared prediction, we directly measured F in 3D spheroids by RNA FISH (Fig. 2G). Scoring individual cells with high SOD2 fluorescence intensity, we calculated an expression frequency of ∼26%. This measurement closely agreed with the inferred parameter of the LN–LN mixture (the better-scoring model; Fig. 2H and SI Appendix, Table S1) and lay within the estimated confidence interval (CI) of the EXP–LN mixture. By resampling the 10-cell SOD2 data, we found that at least 50 observations were required to arrive at an accurate result (SI Appendix, Fig. S3C), confirming our earlier estimates using simulated data (SI Appendix, Fig. S3 A and B). The SOD2 parameterization suggested that maximum-likelihood inference could correctly extract single-cell information from 10-cell sampling data.

Maximum-Likelihood Inference of Coordinated Stochastic Transcriptional Profiles.

Programs of gene expression are often controlled by common upstream factors that enforce the regulatory state. We reasoned that coordinated single-cell gene programs would be the product of an overarching regulatory heterogeneity characterized by a shared F. If true, then it should be possible to estimate the expression frequency more confidently and with fewer samples by extending maximum-likelihood inference to gene clusters with coordinated 10-cell fluctuations.

We extended the approach as follows (Fig. 3A). First, each gene within the cluster was assigned its own µ₁ and µ₂ (or µ and λ for the EXP–LN mixture) to account for gene-to-gene differences in expression level and detection sensitivity. Next, we assumed that the genes within a cluster share a common F and σ (or F, σ₁, and σ₂ in the relaxed LN–LN mixture), implying a shared mechanism of regulation (39, 40). Therefore, each mixture model of a cluster of g genes involved 2g + 2 or 2g + 3 parameters. Even for small gene programs (g ≤ 10), this parameter search space was too large for nonconvex optimization methods to maximize the global likelihood function quickly (Methods). To increase the speed and efficiency of estimation, the cluster was broken down into smaller four-gene subgroups to infer log-means for each gene in the subgroup together with local estimates of F, σ , and λ (steps 1 and 2, Fig. 3A). After log-means were locally estimated, the remaining parameters were globally inferred by remaximizing the likelihood function for the entire gene cluster while retaining the local gene-specific estimates of µ₁ and µ₂ (LN–LN mixture) or µ (EXP–LN mixture) (steps 3 and 4, Fig. 3A). As before, selection of the LN–LN, relaxed LN–LN, and EXP–LN mixture model was made according to the lowest BIC score (SI Appendix, Table S1). This revised formulation of maximum-likelihood inference enabled accurate and confident estimates of the expression frequency while requiring only approximately one-third of the sample size (SI Appendix, Fig. S6).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Maximum-likelihood inference accurately estimates subpopulation frequencies from 10-cell gene-expression clusters. (A) The maximum-likelihood approach was modified for gene clusters with coordinated 10-cell sampling fluctuations as follows. (1) Global gene measurements are grouped and assumed to share a common F and σ. (2) An expression cluster of interest is divided into four-gene subsets for the first round of parameter estimation of μ₁ and μ₂ for each gene in the subset. (3) A maximum-likelihood estimator is derived based on an expanded version of the model in Fig. 2A, where each gene k in a group of genes, {1, …, g}, has its own and . The likelihood function is maximized to infer and locally. (4) The likelihood function is then remaximized for the entire dataset keeping the log-mean estimates ( and ) fixed to provide clusterwide estimates of F and σ. Note that each gene has a different range of gene expression to reflect differences in overall expression levels, which are captured in the model predictions as well. (B–G) Prediction and validation of expression frequency for heterogeneously expressed gene programs during breast-epithelial acinar morphogenesis. (B and C) Heat map of clustered 10-cell transcriptional profiles (27). Gray labels indicate the 10-cell sample numbers from Fig. 3B. Maximum-likelihood estimate of expression frequency (red box) is shown with 95% CI in brackets for each cluster. Note that the two gene clusters are predicted to have substantially different frequencies of high expression based on their 10-cell sampling fluctuations. (D and E) Representative RNA FISH images of transcripts from (D) the infrequent cluster and (E) the rare cluster. Images are shown with DRAQ5 counterstain to visualize nuclei. Arrows show ECM-attached cells with high expression. (Scale bar, 10 μm.) (F and G) Percentage of cells showing high expression by RNA FISH (gray bar) of a subset of genes in each cluster compared with the maximum-likelihood estimate of F (white dashed line). RNA FISH data are shown as the mean percentage with 95% CI of ECM-attached cells showing high expression. Maximum-likelihood predictions are shown as the parameter point estimate (white) with 95% CI (red).

We tested our extension of maximum-likelihood inference by extracting from an earlier study two coexpression clusters that were completely uncharacterized (27) (SI Appendix, Fig. S7). These clusters contained one to two dozen genes with strongly coordinated expression fluctuations across 16 samples of 10 ECM-attached cells, but the patterns of fluctuation were markedly different (Fig. 3 B and C). Accordingly, when we inferred the parameters for the two clusters, the model predicted two very different expression frequencies. The first “infrequent” gene cluster was predicted to be up-regulated in ∼25% of the ECM-attached population (Fig. 3B). The LN–LN mixture model was preferred over the EXP–LN or relaxed LN–LN mixtures (SI Appendix, Table S1), although all three models converged upon similar values for F. By contrast, the expression frequency of the “rare” second cluster was predicted to be ∼10% by the LN–LN mixture (Fig. 3C), which was the best-scoring model of the three (SI Appendix, Table S1). Our parameterization of the two clusters emphasizes the mosaicked regulatory states that evolve even in a very simple model of tissue architecture (27, 38, 41).

To test whether the predicted values of F were accurate within the coexpressed clusters, we designed and validated riboprobes for four or five genes in each cluster and quantified their frequency of high expression by RNA FISH (SI Appendix, Fig. S8 A and B). We found that transcripts in the infrequent expression cluster were strongly expressed in three to five ECM-attached cells per acinus cross-section (Fig. 3 D and F and SI Appendix, Fig. S9A), yielding an average expression frequency of ∼25%. Conversely, genes in the rare expression cluster (Fig. 3 E and G and SI Appendix, Fig. S9B) were strongly expressed in one or two ECM-attached cells per acinus cross-section, consistent with an expression frequency of ∼10%. The expression frequencies of both clusters closely agreed with the inferred F parameters, suggesting that our extended inference approach was effective and accurate.

We evaluated the estimates of expression frequency more broadly by selecting four additional clusters from the same dataset for parameterization (SI Appendix, Fig. S7) (27). The clusters showed distinct fluctuation patterns and consequently led to F estimates that ranged from less than 5% to greater than 25% (Fig. 4 A–D, Upper). We validated riboprobes for a representative gene in each cluster and scored the expression frequency (Fig. 4 A–D, Lower and SI Appendix, Fig. S8C). Together with the earlier clusters, we observed a strong correlation between the expression frequency inferred computationally and the manual counts obtained by RNA FISH (R = 0.89, Fig. 4E). The accuracy of the manual counts was further confirmed by correlation with an expression-frequency index derived from digital image analysis of segmented acini (Methods and SI Appendix, Fig. S10). Taken together, these data indicate that maximum-likelihood inference accurately infers single-cell expression frequencies from cluster-wide patterns of 10-cell fluctuations.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Widespread parameterization of single-cell expression frequency by maximum-likelihood inference. (A–D, Upper) Clusters of 10-cell expression fluctuations among ECM-attached cells (27). A complete list of transcripts in each cluster is shown in SI Appendix, Fig. S7. Maximum-likelihood estimate of expression frequency (red box) is shown with 95% CI in brackets for each cluster. (A–D, Lower) RNA FISH images of a representative transcript in each cluster. Images are shown with DRAQ5 counterstain to visualize nuclei. Arrows show ECM-attached cells with high expression. (Scale bars, 10 μm.) (E) Percentage of cells scored for high expression by RNA FISH compared with the maximum-likelihood estimate of F. RNA FISH data are shown as the mean percentage with 95% CI of ECM-attached cells showing high expression. Maximum-likelihood predictions are shown as the parameter point estimate with 95% CI. The gray bar shows a one-to-one correspondence with 5% measurement tolerance. Estimates for SOD2 are reprinted from Fig. 2H. Estimates for Fig. 3 D and E were calculated by pooling all scored transcripts within each cluster. Pearson correlation (R) between measured and inferred expression frequencies is shown.

Identification of a Peculiar, Very Rare Transcriptional Regulatory State.

Maximum-likelihood inference provides critical information about the state distribution and expression frequency of any gene cluster identified by stochastic profiling to be heterogeneously regulated. As a proof-of-concept application, we screened gene clusters from the 3D profiling data (27) to identify unusual regulatory states that warranted follow-up study. One cluster was notable among those surveyed because the predicted expression frequency of the high regulatory state was very rare (F = 2.3%). The very rare cluster was also distinguished by its strong concordance with the relaxed LN–LN mixture compared with the alternative mixture models (SI Appendix, Table S1). Moreover, the log-mean of the low regulatory state was extremely low (μ₂ ∼–3.3), suggesting that the cluster was at or below detection in the population. Within this coexpression cluster, we recognized several genes that were strongly associated with breast cancer, including the breast cancer susceptibility gene BRIP1 [alternatively called FANCJ or BACH1 (42)], the breast cancer-associated gene IRF2 (43), and the zinc-finger gene HIVEP2, which is frequently down-regulated or mutated in breast cancer (44, 45) (Fig. 5A). We speculated that genes within the cluster were tightly regulated so that they could be activated in a restricted cellular context where their expression was critical.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

A unique, very rare regulatory state is marked by PIK3CD, which is important for normal suppression of proliferation during breast-epithelial acinar morphogenesis. (A) Heat map of clustered 10-cell transcriptional profiles (27). Gray labels indicate the 10-cell sample numbers from Fig. 3B. Maximum-likelihood estimate of expression frequency (red box) is shown with 95% CI in brackets. (B) PIK3CD expression is up-regulated during 3D morphogenesis. Relative PIK3CA (red), PIK3CB (blue), and PIK3CD (black) expression was measured by qPCR at various time points during 3D morphogenesis. Data are shown as mean expression ± SEM normalized to the day-4 expression of PIK3CD of three independent experiments. PIK3CG was not expressed in MCF10A-5E cells (SI Appendix, Fig. S11). (C) Representative RNA FISH image of PIK3CD expression is shown with DRAQ5 counterstain to visualize nuclei. Arrow shows ECM-attached cells with high expression of PIK3CD. (Scale bar, 20 µm.) (D) Knockdown of p110δ by shRNA. MCF10A-5E cells were infected with either shGFP (lane 1) or with one of two shRNA sequences targeting p110δ (lanes 2 and 3). Lysates were analyzed by immunoblotting with tubulin used as the loading control. Densitometry of p110δ abundance is shown relative to the shGFP control. (E and F) Disruption of normal PIK3CD regulation elicits a hyperproliferative phenotype in 3D culture. shGFP, shPIK3CD #1, and shPIK3CD #2 cells or shGFP cells + 20 μM p110δ inhibitor IC87114 were fixed at day 15 of 3D morphogenesis, stained for pRb (red), and analyzed by confocal immunofluorescence. Cells were counterstained with DRAQ5 (blue) to label nuclei. Arrows in F highlight pRb-positive cells. (Scale bars, 20 μm.) Quantification of proliferating acini in each condition is shown in E as the mean ± SEM of eight independent experiments.

Among the genes in the very rare cluster, we were most intrigued by the phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit δ isoform (PIK3CD alternatively called p110δ). Three-dimensional breast epithelial cultures abundantly express two other PI3K isoforms, PIK3CA and PIK3CB (Fig. 5B and SI Appendix, Fig. S11), and it is generally thought that any PI3K isoform can support proliferation and survival (46). Nevertheless, we found that the low-copy expression of PIK3CD was transcriptionally up-regulated with delayed kinetics compared with the other PI3K isoforms (Fig. 5B), suggesting a unique regulatory mechanism. When PIK3CD abundance was visualized in single cells by RNA FISH, we observed a striking pattern. Most cells lacked PIK3CD or expressed it at very low levels; however, we consistently identified a sporadic subpopulation of cells (roughly one or two cells every other acinus cross-section) with high PIK3CD expression (Fig. 5C and SI Appendix, Fig. S8D). The overall frequency of cells in the PIK3CD^hi state was somewhat higher than the estimates of F for the cluster, but the inferred frequency agreed with the very rare expression of two other members of the cluster, FEM1A and IRF2 (SI Appendix, Fig. S12). Together, these observations pointed to an acute (and likely transient) regulatory event triggering the selective induction of cluster genes in single ECM-attached cells.

We next asked whether PIK3CD was specifically important for normal acinar morphogenesis. To eliminate the very rare PIK3CD^hi subpopulation, we perturbed p110δ by two independent methods: RNA interference and the p110δ-specific small-molecule inhibitor, IC87114 (ref. 47, Fig. 5D, and SI Appendix, Fig. S13). When shPIK3CD cells were placed in 3D culture, we found that acini were larger and distorted, suggesting a defect in proliferation arrest. Using phosphorylated Rb (pRb) as a proliferative marker, we observed that shPIK3CD acini were still cycling after 15 d of 3D culture when shGFP control acini had quiesced (Fig. 5 E and F). Furthermore, when control cells were cultured with IC87114, we observed sustained proliferation that phenocopied PIK3CD knockdown (Fig. 5 E and F). These data together indicate that p110δ activity stemming from the very rare PIK3CD^hi regulatory state is critical for normal proliferation arrest of breast epithelia in 3D culture. More broadly, our results with the very rare cluster illustrate how maximum-likelihood inference can be used to hone in on gene programs with an expression frequency or regulatory pattern of interest.

Comparison with Alternative Deconvolution Methods.

We compared the performance of our method to other computational approaches for deconvolving mixed populations (48⇓–50). The alternative methods invoked different mathematical formalisms—Bayesian statistics (48), nonnegative matrix factorization (49), and principal component analysis (50)—and none had previously been applied to transcriptional profiles of small samples. Using the sampling fluctuations within the infrequent, rare, and very rare clusters, we attempted estimates of expression frequency and found that all were substantially less accurate than our approach (Table 1). The comparison illustrates that our method is uniquely effective at parameterizing transcriptional regulatory states within cell populations.

Direct Comparison of Single-Cell and 10-Cell Sampling Strategies.

Maximum-likelihood inference reconstructs the single-cell expression distribution without the need to measure single cells. Ignoring the technical challenges of global single-cell methods (17, 20, 21, 27), it should also be theoretically possible to recreate the complete expression distribution by measuring many individual cells. However, it was not clear whether single-cell profiling would be as effective as stochastic profiling when reconstructing from a limited number of 1- or 10-cell samples. We anticipated that low expression frequencies would be particularly difficult for single-cell methods to characterize because of uncertainty in reliably capturing the rare regulatory state.

To compare single-cell profiling with stochastic profiling, we repeatedly simulated 1- or 10-cell measurements of gene clusters with similar characteristics to those previously examined (Figs. 3 F and G and 5A and Methods). The three 12-gene clusters varied in their expression fraction—infrequent (F = 25%), rare (F = 10%), and very rare (F = 5%)—and the very rare cluster was simulated as an LN–LN mixture or an EXP–LN mixture. When the number of observations was limited to 16 (as in the actual data), we found that maximum-likelihood inference provided superior estimates of F when using 10-cell groups (Table 2). Estimates from simulated observations of 16 single cells showed substantially higher mean squared error (MSE) for all gene clusters compared with 16 10-cell observations. The larger MSE of one-cell estimates was caused by increases in both the bias and variance of the estimate, whose magnitudes depended on the cluster characteristics and mixture model. These computational simulations provide an upper bound on performance, because experimental error from actual single-cell experiments (17, 19) should blur the data much more. By collecting a greater total number of cells when observations are limited, maximum-likelihood inference of stochastic 10-cell profiles provides a more accurate picture of the single-cell distribution than single-cell profiles.