stochprofML: Stochastic Profiling Using Maximum Likelihood Estimation in \proglangR

Suppose there are $k$ (tissue) samples, indexed by $i\in\{1,\ldots,k\}$. From each tissue sample $i$, we collect a pool of a known number of cells. The cells are either indexed by $j\in\{1,\ldots,n\}$ if the cell pool size is the same in all measurements, or, as possible in the latest implementation, by $j_{i}\in\{1,\ldots,n_{i}\}$ in case cell pool sizes vary between measurements. In the latter, more general case, the cell numbers are variable over the $k$ cell pools and summarized by $\vec{n}=(n_{1},\ldots,n_{k})$. From each sample, the gene expression of $m$ genes is measured, indexed by $g\in\{1,\ldots,m\}$. We assume that each cell stems from one out of $T$ cell populations, indexed by $h\in\{1,\ldots,T\}$. If $T>1$ in the set of all cells of interest, the tissue is called heterogeneous. The notation is illustrated in Figure 1.

Biologically, the different cell populations correspond to different regulatory states or — especially in the context of cancer — to different (sub-)clones. For example, there might be two populations within a considered tissue: one occupying a basal regulatory state, where the expression of genes is at a low level, and one from a second regulatory state, where genes are expressed at a higher level.

As described in Section 1, there are various technologies to measure gene expression. In case of microarray techniques as done in previous applications of stochastic profiling, see Janes et al., 2010 and Bajikar et al., 2014, the measured gene expression is typically described in terms of continuous probability distributions. Conditioned on the cell population, \pkgstochprofML provides two choices for the distribution of the expression of one gene:

stochprofML is tailored to analyze gene expression measurements of small pools of cells, beyond the analysis of standard single-cell gene expression data. In other words, the single-cell gene expression $X_{ij_{i}}^{(g)}$ described in Section 2.1 is assumed latent. Instead, we consider observations

for $i=1,\ldots,k$, which represent the overall gene expression of the $i$th cell pool for gene $g$. In the first version of \pkgstochprofML, pools had to be of equal size $n$, i.e. for each measurement $Y_{i}^{(g)}$ one had to extract the same number of cells from each tissue sample. This was a restrictive assumption from the experimental point of view. The recent extension of \pkgstochprofML allows each cell pool $i$ to contain a different number $n_{i}$ of cells (see also Figures 1 and 2).

The algorithm aims to estimate the single-cell population parameters despite the fact that measurements are available only in convoluted form. To that end, we derive (in Section 2.3) the likelihood function of the parameters in the convolution model (2), where we assume the gene expression of the single cells to be independent within a tissue sample. For better readability, we suppress for now the superscript $(g)$ and introduce it again in Section 2.3.

The derivation of the distribution of $Y_{i}$ is described in Appendix A. The corresponding PDF $f_{n_{i}}(y_{i}|\boldsymbol{\theta},\boldsymbol{p})$ of an observation $y_{i}$ which represents the overall gene expression from sample $i$ (consisting of $n_{i}$ cells) is given by

where $\ell_{T}=n_{i}-\sum_{h=1}^{T-1}\ell_{h}$ and $p_{T}=1-\sum_{h=1}^{T-1}p_{h}$. Here, $f_{(\ell_{1},\ell_{2},\ldots,\ell_{T})}$ describes the PDF of a pool of $n_{i}$ cells with known composition of the single populations, i.e. it is known that there are $\ell_{1}$ cells from population $1$, $\ell_{2}$ cells from population $2$ etc. ${n_{i}\choose{\ell_{1},\ell_{2},\ldots,\ell_{T}}}p_{1}^{\ell_{1}}p_{2}^{\ell_{2}}\cdots p_{T}^{\ell_{T}}$ represents the multinomial probability of obtaining exactly this composition $(\ell_{1},\ldots,\ell_{T})$ using the multinomial coefficient ${n_{i}\choose{\ell_{1},\ell_{2},\ldots,\ell_{T}}}=n_{i}!/(\ell_{1}!\ldots\ell_{T}!)$. Equation (3) sums up over all possible compositions $(\ell_{1},\ldots,\ell_{T})$ with $\ell_{1},\ldots,\ell_{T}\in\mathbb{N}_{0}$ and $\ell_{1}+\ldots+\ell_{T}=n_{i}$. Taken together, $f_{n_{i}}(y_{i}|\boldsymbol{\theta},\boldsymbol{p})$ determines the PDF of $y_{i}$ with respect to each possible combination of $n_{i}$ cells of $T$ populations.

Thus, the calculation of $f_{n_{i}}(y_{i}|\boldsymbol{\theta},\boldsymbol{p})$ requires knowledge of $f_{(\ell_{1},\ell_{2},\ldots,\ell_{T})}(y_{i}|\boldsymbol{\theta})$ . The derivation of this PDF depends on the choice of the single-cell model (LN-LN, rLN-LN, or EXP-LN; see Section 2.1) that was made for $X_{ij_{i}}$ .

1. 1.

   LN-LN:

   $$f_{(\ell_{1},\ldots,\ell_{h},\ldots,\ell_{T})}(y_{i}|\boldsymbol{\theta})=f^{\text{LN-LN}}_{(\ell_{1},\ldots,\ell_{h},\ldots,\ell_{T})}(y_{i}|\mu_{1},\ldots,\mu_{h},\ldots,\mu_{T},\sigma^{2})$$

   is the density of a sum $Y_{i}=X_{i1}+\ldots+X_{in_{i}}$ of ${n_{i}}$ independent random variables with

   $$X_{ij_{i}}\sim\begin{cases}\mathcal{L}\mathcal{N}(\mu_{1},\sigma^{2})&\text{if }1\leq j_{i}\leq J_{1}\\ \vdots&\\ \mathcal{L}\mathcal{N}(\mu_{h},\sigma^{2})&\text{if }J_{h-1}<j_{i}\leq J_{h}\\ \vdots&\\ \mathcal{L}\mathcal{N}(\mu_{T},\sigma^{2})&\text{if }J_{T-1}<j_{i}\leq J_{T}=n_{i},\end{cases}$$

   with $J_{1}=\ell_{1},\ldots,J_{h}=\ell_{1}+\ell_{2}+\ldots+\ell_{h},\ldots,J_{T}=\ell_{1}+\ell_{2}+\ldots+\ell_{T}=n_{i}$. $Y_{i}$ is the convolution of random variables $X_{i1},\ldots,X_{in_{i}}$, which is here the convolution of $T$ sub-convolutions: a convolution of $\ell_{1}$ times $\mathcal{L}\mathcal{N}(\mu_{1},\sigma^{2})$, plus a convolution of $\ell_{2}$ times $\mathcal{L}\mathcal{N}(\mu_{2},\sigma^{2})$, and so on, up to a convolution of $\ell_{T}$ times $\mathcal{L}\mathcal{N}(\mu_{T},\sigma^{2})$.

   There is no analytically explicit form for the convolution of lognormal random variables. Hence, $f^{\text{LN-LN}}_{(\ell_{1},\ldots,\ell_{h},\ldots,\ell_{T})}$ is approximated using the method by Fenton (1960). That is, the distribution of the sum $A_{1}+\ldots+A_{m}$ of independent random variables$A_{i}\sim\mathcal{L}\mathcal{N}(\mu_{A_{i}},\sigma_{A_{i}}^{2})$ is approximated by the distribution of a random variable $B\sim\mathcal{L}\mathcal{N}(\mu_{B},\sigma_{B}^{2})$ such that

   $$\E(B)=\E(A_{1}+\ldots+A_{m})\quad\text{and}\quad\text{Var}(B)=\text{Var}(A_{1}+\ldots+A_{m}).$$

   According to Equation (1), that means that $\mu_{B}$ and $\sigma_{B}$ are chosen such that the following equations are fulfilled:

   $$\exp\left(\mu_{B}+\frac{\sigma_{B}^{2}}{2}\right)=\exp\left(\mu_{A_{1}}+\frac{\sigma_{A_{1}}^{2}}{2}\right)+\ldots+\exp\left(\mu_{A_{m}}+\frac{\sigma_{A_{m}}^{2}}{2}\right)=:\Gamma$$

   and

   	$\displaystyle\exp\left(2\mu_{B}+\sigma_{B}^{2}\right)\left(\exp\left(\sigma_{B}^{2}\right)-1\right)=$
   	$\displaystyle\quad\exp\left(2\mu_{A_{1}}+\sigma_{A_{1}}^{2}\right)\left(\exp\left(\sigma_{A_{1}}^{2}\right)-1\right)+\ldots+\exp\left(2\mu_{A_{m}}+\sigma_{A_{m}}^{2}\right)\left(\exp\left(\sigma_{A_{m}}^{2}\right)-1\right)=:\Delta.$

   That is achieved by choosing

   $$\mu_{B}=\log(\Gamma)-\frac{1}{2}\sigma_{B}^{2}\quad\text{and}\quad\sigma_{B}^{2}=\log\left(\frac{\Delta}{\Gamma^{2}}+1\right).$$

   This approximation is implemented in the function \coded.sum.of.lognormals(). The overall PDF is computed through \coded.sum.of.mixtures.LNLN().

2. 2.

   rLN-LN:

   $$f_{(\ell_{1},\ldots,\ell_{h},\ldots,\ell_{T})}(y_{i}|\boldsymbol{\theta})=f^{\text{rLN-LN}}_{(\ell_{1},\ldots,\ell_{h},\ldots,\ell_{T})}(y_{i}|\mu_{1},\ldots,\mu_{h},\ldots,\mu_{T},\sigma_{1}^{2},\ldots,\sigma_{h}^{2},\ldots,\sigma_{T}^{2})$$

   is the PDF of a sum $Y_{i}=X_{i1}+\ldots+X_{in_{i}}$ of ${n_{i}}$ independent random variables with

   $$X_{ij_{i}}\sim\begin{cases}\mathcal{L}\mathcal{N}(\mu_{1},\sigma_{1}^{2})&\text{if }1\leq j_{i}\leq J_{1}\\ \vdots&\\ \mathcal{L}\mathcal{N}(\mu_{h},\sigma_{h}^{2})&\text{if }J_{h-1}<j_{i}\leq J_{h}\\ \vdots&\\ \mathcal{L}\mathcal{N}(\mu_{T},\sigma_{T}^{2})&\text{if }J_{T-1}<j_{i}\leq J_{T}=n_{i},\end{cases}$$

   with $J_{1}=\ell_{1},\ldots,J_{h}=\ell_{1}+\ell_{2}+\ldots+\ell_{h},\ldots,J_{T}=\ell_{1}+\ldots+\ell_{T}=n_{i}$. Again, $f^{\text{rLN-LN}}_{(\ell_{1},\ldots,\ell_{h},\ldots,\ell_{T})}$ is approximated using the method by Fenton (1960), analogously to the LN-LN model. It is implemented in \coded.sum.of.mixtures.rLNLN().

3. 3.

   EXP-LN:

   $$f_{(\ell_{1},\ell_{2},\ldots,\ell_{T})}(y_{i}|\boldsymbol{\theta})=f^{\text{EXP-LN}}_{(\ell_{1},\ell_{2},\ldots,\ell_{T})}(y_{i}|\lambda,\mu_{1},\ldots,\mu_{T-1},\sigma^{2})$$

   is the density of a sum $Y_{i}=X_{i1}+\ldots+X_{in_{i}}$ of $n_{i}$ independent random variables with

   $$X_{ij_{i}}\sim\begin{cases}\mathcal{L}\mathcal{N}(\mu_{1},\sigma^{2})&\text{if }1\leq j_{i}\leq J_{1}\\ \vdots&\\ \mathcal{L}\mathcal{N}(\mu_{h},\sigma^{2})&\text{if }J_{h-1}<j_{i}\leq J_{h}\\ \vdots&\\ \mathcal{L}\mathcal{N}(\mu_{T-1},\sigma^{2})&\text{if }J_{T-2}<j_{i}\leq J_{T-1}\\ \mathcal{E}\!\mathcal{X}\!\mathcal{P}(\lambda)&\text{if }J_{T-1}<j_{i}\leq J_{T}=n_{i},\end{cases}$$

   with $J_{1}=\ell_{1},\ldots,J_{h}=\ell_{1}+\ell_{2}+\ldots+\ell_{h},\ldots,J_{T}=\ell_{1}+\ldots+\ell_{T}=n_{i}$. The sum of independent exponentially distributed random variables with equal intensity parameter follows an Erlang distribution (Feldman and Valdez-Flores, 2010), which is a gamma distribution with integer-valued shape parameter that represents the number of exponentially distributed summands. Thus, the PDF for the EXP-LN mixture model is the convolution of one Erlang (or gamma) distribution (namely the sum of all exponentially distributed summands) and one lognormal distribution (namely the sum of all lognormally distributed summands, again using the approximation method by Fenton, 1960). The PDF for this convolution is not known in analytically explicit form but expressed in terms of an integral that is solved numerically through the function \codelognormal.exp.convolution(). The overall PDF of the EXP-LN model is implemented in \coded.sum.of.mixtures.EXPLN().

Example: Mixture of two populations - Part 2. In this example of the two-population model, let each observation consist of the same number of $n=10$ cells. Then $Y_{i}$ is a $10$-fold convolution, and the PDF (3) simplifies to

where $f_{(\ell,10-\ell)}$ is the PDF of the sum $Y_{i}$ of ten independent random variables, that is$Y_{i}=X_{i1}+\ldots+X_{i10}$. This PDF depends on the particular chosen model:

1. 1.

   LN-LN

   $$f_{(\ell,10-\ell)}(y_{i}|\boldsymbol{\theta})=f^{\text{LN-LN}}_{(\ell,10-\ell)}(y_{i}|\mu_{1},\mu_{2},\sigma^{2})$$

   is the PDF of a sum $Y_{i}=X_{i1}+\ldots+X_{i10}$ of ten independent random variables with

   $$X_{ij}\sim\begin{cases}\mathcal{L}\mathcal{N}(\mu_{1},\sigma^{2})&\text{if }1\leq j\leq\ell\\ \mathcal{L}\mathcal{N}(\mu_{2},\sigma^{2})&\text{if }\ell<j\leq 10.\end{cases}$$

2. 2.

   rLN-LN

   $$f_{(\ell,10-\ell)}(y_{i}|\boldsymbol{\theta})=f^{\text{rLN-LN}}_{(\ell,10-\ell)}(y_{i}|\mu_{1},\mu_{2},\sigma_{1}^{2},\sigma_{2}^{2})$$

   is the PDF of a sum $Y_{i}=X_{i1}+\ldots+X_{i10}$ of ten independent random variables with

   $$X_{ij}\sim\begin{cases}\mathcal{L}\mathcal{N}(\mu_{1},\sigma_{1}^{2})&\text{if }1\leq j\leq\ell\\ \mathcal{L}\mathcal{N}(\mu_{2},\sigma_{2}^{2})&\text{if }\ell<j\leq 10.\end{cases}$$

3. 3.

   EXP-LN

   $$f_{(\ell,10-\ell)}(y_{i}|\boldsymbol{\theta})=f^{\text{EXP-LN}}_{(\ell,10-\ell)}(y_{i}|\lambda,\mu,\sigma^{2})$$

   is the PDF of a sum $Y_{i}=X_{i1}+\ldots+X_{i10}$ of ten independent random variables with

   $$X_{ij}\sim\begin{cases}\mathcal{L}\mathcal{N}(\mu,\sigma^{2})&\text{if }1\leq j\leq\ell\\ \mathcal{E}\!\mathcal{X}\!\mathcal{P}(\lambda)&\text{if }\ell<j\leq 10.\\ \end{cases}$$

Overall, after re-introducing the superscript $(g)$ for measurements of genes $g=1,\ldots,m$, we obtain the PDF

with model-specific choice of $f_{(\ell_{1},\ell_{2},\ldots,\ell_{T})}$. While $\boldsymbol{n}=(n_{1},\ldots,n_{k})$ is considered known, we aim to infer the unknown model parameters $\boldsymbol{\theta}=\{\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(m)},\boldsymbol{p}\}$ by maximum likelihood estimation. Assuming independent observations ${\boldsymbol{y}}=\{y_{i}^{(g)}|i=1,\ldots,k;g=1,\ldots,m\}$ of $Y_{i}^{(g)}$ for $m$ genes and $k$ tissue samples, where sample $i$ contains $n_{i}$ cells, the likelihood function is given by

Consequently, the log-likelihood function of the model parameters reads

Example: Mixture of two populations - Part 3. Returning to the two-population example with 10-cell pools, the log-likelihood for $k=100$ tissue samples and $m=5$ genes is given by

where $f_{10}\left(y_{i}^{(g)}|\boldsymbol{\theta}^{(g)},\boldsymbol{p}\right)$ is given by Equation (4).

The \pkgstochprofML algorithm aims to infer the unknown model parameters using maximum likelihood estimation. As input, we expect an $m\times k$ data matrix of pooled gene expression, known cell numbers $\vec{n}$, the assumed number of populations $T$ and the choice of single-cell distribution (LN-LN, rLN-LN, EXP-LN). Based on this input, the algorithm aims to find parameter values of $\boldsymbol{\theta}=\{\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(m)},\boldsymbol{p}\}$ that maximize $\ell(\boldsymbol{\theta}|\boldsymbol{y})$ as given by Equation (6). This section describes practical aspects of the optimization procedure.

Example: Mixture of two populations - Part 4. Several challenges occur during parameter estimation. We explain these on the two-population LN-LN example: First, we aim to ensure parameter identifiability. This is achieved for the two-population LN-LN model by constraining the parameters to fulfil either $p\leq 0.5$ or $\mu_{1}>\mu_{2}$. Otherwise, the two combinations $(p,{\boldsymbol{\mu}}_{1},{\boldsymbol{\mu}}_{2},\sigma)$ and $(1-p,{\boldsymbol{\mu}}_{2},{\boldsymbol{\mu}}_{1},\sigma)$ would yield identical values of the likelihood function and could cause computational problems. For our implementation, we preferred the second possibility, i. e.  $\mu_{1}>\mu_{2}$. The alternative, i. e.  requiring $p\leq 0.5$, led to switchings between $\mu_{1}$ and $\mu_{2}$ in case of $p\approx 0.5$. As a second measure, we implement unconstrained rather than constrained optimization: Instead of estimating $(p,{\boldsymbol{\mu}}_{1},{\boldsymbol{\mu}}_{2},\sigma)$ under the constraints $p\in[0,1]$, $\mu_{1}>\mu_{2}$ and $\sigma>0$, the parameters are transformed to $(\text{logit}(p),{\boldsymbol{\mu}}_{1},{\boldsymbol{\mu}}_{2},\log(\sigma))$, and an unconstrained optimization method is used. This is substantially faster.

The aforementioned transformations are likewise employed for all other models (rLN-LN and EXP-LN) and population numbers. In particular, $\sigma$ and $\lambda$ are log-transformed, and the lognormal populations are ordered according to the log-means $\mu_{h}^{(1)}$ of the first gene in the gene list. The population probabilities are transformed to $\mathbb{R}$ such that they still sum up to one after back-transformation. For details, see Appendix C.

The log-likelihood function is multimodal. Thus, a single application of some gradient-based optimization method does not suffice to find a global maximum. Instead, two approaches are combined which are alternately executed: First, a grid search is performed, where the log-likelihood function is computed at randomly drawn parameter values. This way, high likelihood regions can be identified with low computational cost. In the second step, the Nelder-Mead algorithm (Nelder and Mead, 1965) is repeatedly executed. Its starting values are randomly drawn from the high likelihood regions found before. The following grid search then again explores the regions around the obtained local maxima, and so on.

If a dataset contains gene expressions for $m$ genes, and if we assume $T$ populations, there are at minimum $T(m+1)$ parameters which one seeks to estimate depending on the model framework. This is computationally difficult, because the number of modes of the log-likelihood function increases with the number of parameters. The performance of the numerical optimization crucially depends on the quality of the starting values, and a large number of restarts is required. When analyzing a large gene cluster, it is advantageous to start by considering small clusters and use the derived estimates as initial guesses for larger clusters. This is implemented in the function \codestochprof.loop() (parameter \codesubgroups) and demonstrated in \codeanalyze.toycluster().

Approximate marginal 95% confidence intervals for the parameter estimates are obtained as follows: We numerically compute the Hessian matrix of the negative log-likelihood function on the unrestricted parameter space and evaluate it at the (transformed) maximum likelihood estimator. Denote by $d_{i}$ the $i$th diagonal element of the inverse of this matrix. Then the confidence bounds for the $i$th transformed parameter $\theta_{i}$ are

We obtain respective marginal confidence intervals for the original true parameters by back-transformation of the above bounds. This approximation is especially appropriate in the two-population example for the parameters $p$ and $\sigma$ when conditioning on ${\boldsymbol{\mu}}_{1}$ and ${\boldsymbol{\mu}}_{2}$. In this case, in practice, the profile likelihood is seemingly unimodal.

By increasing the number $T$ of populations, we can describe the observed data more precisely, but this comes at the cost of potential overfitting. For example, a three-population LN-LN model may lead to a larger likelihood at the maximum likelihood estimator than a two-population LN-LN model on the same dataset. However, the difference may be small, and the additional third population may not lead to a gain of knowledge. For example, the estimated population probability $\hat{p}_{3}$ may be tiny, or the log-means of the second and third population, $\hat{\mu}_{2}$ and $\hat{\mu}_{3}$ might hardly be distinguished from each other.

To objectively find a trade-off between necessary complexity and sufficient interpretability, we employ the Bayesian information criterion (BIC, Schwarz, 1978):

where $\hat{{\boldsymbol{\theta}}}$ is the maximum likelihood estimate of the respective model, $\dim(\hat{{\boldsymbol{\theta}}})$ the number of parameters and $k$ the size of the dataset. From the statistics perspective, the model with smallest BIC is considered most appropriate among all considered models.

In practice, it is required to estimate all models of interest separately with the \pkgstochprofML algorithm, e. g. the LN-LN model with one, two and three populations, and/or the respective rLN-LN and EXP-LN models. The BIC values are returned by the function \codestochprof.loop().

We first generate a dataset of $k=1000$ sample observations, where each sample consists of $n=10$ cells. We choose a single-cell model with two populations, both of lognormal type, i.e. we use the LN-LN model. Let us assume that the overall population of interest is a mixture of $62\%$ of population $1$ and $38\%$ of population $2$, i.e. $p_{1}=0.62$. As population parameters we choose $\mu_{1}=0.47$, $\mu_{2}=-0.87$ and $\sigma=0.03$. Synthetic gene expression data for one gene is generated as follows: {Schunk} {Sinput} R> library("stochprofML") R> set.seed(10) R> k <- 1000 R> n <- 10 R> TY <- 2 R> p <- c(0.62, 0.38) R> mu <- c(0.47, -0.87) R> sigma <- 0.03 R> gene_LNLN <- r.sum.of.mixtures.LNLN(k = k, n = n, p.vector = p, + mu.vector = mu, sigma.vector = rep(sigma, TY))

Figure 3 shows a histogram of the simulated data as well as the theoretical PDF of the 10-cell mixture. The following code produces this figure: {Schunk} {Sinput} R> x <- seq(from = min(gene_LNLN), to = max(gene_LNLN), length = 500) R> stochprofML:::set.model.functions("LN-LN") R> y <- d.sum.of.mixtures(x, n, p, mu,rep(sigma,TY), logdens = FALSE) R> hist(gene_LNLN, main = paste("Simulated Gene"), breaks = 50, + xlab = "Sum of mixtures of lognormals", ylab = "Density", + freq = FALSE, col = "lightgrey") R> lines(x, y, col="blue", lwd = 2) R> legend("topright", legend = "data generating pdf", col = "blue", + lwd = 2, bty = "n")

Next, we show how the parameters used above can be back-inferred from the generated dataset using maximum likelihood estimation. {Schunk} {Sinput} R> set.seed(20) R> result <- stochprof.loop(model = "LN-LN", + dataset = matrix(gene_LNLN, ncol = 1), n = n, TY = TY, + genenames = "SimGene", fix.mu = FALSE, loops = 10, + until.convergence = FALSE, print.output = FALSE, show.plots = TRUE, + plot.title = "Simulated Gene", use.constraints = FALSE) When the fitting is done, pressing <enter> causes \proglangR to show plots of the estimation process, see Figure 4, and displays the results in the following form.

Maximum likelihood estimate (MLE): p_1 mu_1_gene_SimGene mu_2_gene_SimGene sigma 0.6146 0.4710 -0.8720 0.0310

Value of negative log-likelihood function at MLE: 1204.371

Violation of constraints: none

BIC: 2436.373

Approx. 95lower upper p_1 0.60501813 0.6240938 mu_1_gene_SimGene 0.46972264 0.4722774 mu_2_gene_SimGene -0.87827704 -0.8657230 sigma 0.02967451 0.0323847

Top parameter combinations: p_1 mu_1_ge_SimGene mu_2_gene_SimGene sigma target p_1 mu_1_gene_SimGene mu_2_gene_SimGene sigma target [1,] 0.6146 0.471 -0.872 0.031 1204.371 [2,] 0.6146 0.470 -0.872 0.031 1204.371 [3,] 0.6146 0.471 -0.872 0.031 1204.371 [4,] 0.6146 0.470 -0.872 0.031 1204.371 [5,] 0.6145 0.471 -0.872 0.031 1204.371 [6,] 0.6146 0.471 -0.872 0.031 1204.371

Hence, the marginal confidence intervals cover the true parameter values.

Stochastic profiling, i.e. the analysis of small-pool gene expression measurements, is a compromise between the analysis of single cells and the consideration of large bulks: Single-cell information is most immediate, but a fixed number $k$ of samples will only cover $k$ cells. In pools of cells, on the other hand, information is convoluted, but $k$ pools of size $n$ cover $n$ times as much material. An obvious question is the optimal pool size $n$. The answer is not available in analytically closed form. We hence study this question empirically.

For this simulation study, first, we generate synthetic data for different pool sizes with identical parameter values and settings. Then, we re-infer the model parameters using the \pkgstochprofML algorithm. This is repeated 1,000 times for each choice of pool size, enabling us to study the algorithm’s performance by simple summary statistics of the replicates.

The fixed settings are as follows: We use the two-population LN-LN model to generate data for one gene with $p_{1}=0.2$, $\mu_{1}=2$, $\mu_{2}=0$ and $\sigma=0.2$. For each pool size we simulate $k=50$ observations. The pool sizes are chosen in nine different ways: In seven cases, pool sizes are identical for each sample, namely $n\in\{1,2,5,10,15,20,50\}$. In two additional cases, pool sizes are mixed, i.e.  each of the $k$ samples within one dataset represents a pool of different size $n_{i}\in\{1,2,5,10\}$ or $n_{i}\in\{10,15,20,50\}$.

Figure 5 summarizes the point estimates of the 1,000 datasets for each of the nine pool size settings. It seems that (for this particular choice of model parameter values) parameter estimation works reliably for pool sizes up to ten cells, with smaller variance from single-cells to 5-cells. This applies also for the mixture of pool sizes for the small cell numbers. For cell numbers larger than ten, the range of estimated values becomes considerably larger, but without obvious bias, which also applies to the mixture of the larger pool sizes. Appendix F shows repetitions of this study for different choices of population parameters. The results there confirm the observations made here.

Figure 5 suggests $n=5$ or varying small pool sizes as ideal choices since its estimates show smaller variance than the other pool sizes. This simulation study, however, has been performed in an idealized in silico setting: We did not include any measurement noise. In practice, however, it is well known that single-cells suffer more from such noise than samples with many cells. The ideal choice of pool size may hence be larger in practice.

The underlying data-generating model obviously influences the ability of the maximum likelihood estimator to re-infer the true parameter values: Values of $p_{1}$ close to $0.5$, small differences between $\mu_{1}$ and $\mu_{2}$ and large $\sigma$ blur the data and complicate parameter inference in practice. In the simulation study of this section, we investigate the sensitivity of parameter inference and which scenarios could be realistically identified.

We use the same datasets as in the previous simulation study: The parameter choices from Section 4.1 are considered as the standard and compared to those from Appendix F. In detail, $p_{1}$ is reduced from $0.2$ to $0.1$ in one setting and increased to $0.4$ in the next. $\mu_{2}$ is increased from $0$ to $1$, and $\sigma$ increases from $0.2$ to $0.5$. $\mu_{1}$ is kept fixed to $2$ in all settings. As before, we consider 1,000 data sets for every parameter setting and compare the resulting estimates to the true values. This was done for all pool sizes considered in Section 4.1, but here we only comment on the results of the 10-cell pools and refer to Appendix F for all other pool size settings.

Figure 6 shows the results of the study. In each row of the plot, we compare the estimates of the datasets that were simulated with the standard parameters to the estimates of the datasets that were simulated with one of the parameters changed. Even if only one parameter is changed all parameters are estimated. Each violin accumulates the estimates of 1,000 datasets. For easier comparison, each of the twelve tiles shows the standard setting as turquoise violin, which means those are repeated in each row.

When changing the parameter values, they can still be derived without obvious additional bias, but accuracy decreases for increasing $p$, decreasing $\mu_{2}-\mu_{1}$ and increasing $\sigma$ (with few exceptions). Result for other pool sizes (see Appendix F) show that these observations can be transferred to other pool sizes with some additions: Larger pool sizes infer parameters more accurately if $p$ is smaller. In an increased first population setting ($p=40\%$), $\mu_{1}$ can be better inferred if the data set consists of smaller pools. For larger pools, the estimation of $\mu_{1}$ and $\mu_{2}$ works comparably well after increasing $\mu_{2}$. In general, the estimation of $\sigma$ is the most difficult one; already in pools of ten cells with increased $\mu_{2}$, $\sigma$ is underestimated. This worsens in larger pools.

One key assumption of the \pkgstochprofML algorithm is that the exact number of cells in each cell pool is known. In Janes et al. (2010), accordingly, ten cells were randomly taken from each sample by experimental design. However, different experimental protocols may not reveal the exact cell number: In Tirier et al. (2019), for example, tissue samples were taken as whole cancer spheroids. Here, the cell numbers were experimentally unknown but estimated using light sheet microscopy and 3D image analysis. Since the \pkgstochprofML algorithm requires the pool sizes as input parameter, some estimate has to be passed to it. It is intuitively obvious that the better the prior knowledge about the cell pool sizes, the better the final model parameter estimate. In this simulation study, we investigate the consequences of misspecification.

In a first simulation study, we reuse from Section 4.1 the 1,000 synthetic 10-cell datasets. Each of these contains 50 10-cell samples, simulated with underlying model parameters $p=0.2$, $\mu_{1}=2$, $\mu_{2}=0$ and $\sigma=0.2$. As before, we re-infer the population parameters using the \pkgstochprofML algorithm. This time, however, we use varying pool sizes from $5$ to $15$ as input parameters of the algorithm. This is a misspecification except for the true value $10$. The resulting parameter estimates (empirical median and 2.5%-/97.5%-quantiles across the 1,000 datasets) are depicted in Figure 7. Estimates are optimal or at least among the best in terms of empirical bias and variance when using the correct pool size. With increasing assumed cell number, the estimates of $p$ decrease, i. e.  the fraction of cells from the higher expressed population is assumed to be smaller. This is a reasonable consequence of overestimating $n$, because in this case the surplus cells are assigned to the second population with lower (or even close-to-zero) expression. Consequently, at the same time the estimates of $\mu_{2}$ decrease to be even smaller.

In a second simulation study, we use the two settings with mixed cell pool sizes as introduced in Section 4.1. One setting embraces cell pools with rather small cell numbers (single-, 2-, 5- and 10-cell samples), the other one pools with larger cell numbers (10-, 15-, 20- and 50-cell samples). For each of the two scenarios, we generate one dataset with 50 samples. We denote the true 50-dimensional pool size vectors by $\vec{n}_{\text{small}}$ and $\vec{n}_{\text{large}}$ and employ these vectors for re-estimating the model parameters $p$, $\mu_{1}$, $\mu_{2}$ and $\sigma$. Then, we estimate the parameters again for the same two datasets for 1,000 times, but this time using perturbed pool size vectors as input to the algorithm, introducing artificial misspecification. These 50-dimensional pool size vectors are generated as follows: For each component, we draw a Poisson-distributed random variable with intensity parameter equal to the respective component of the true vectors $\vec{n}_{\text{small}}$ or $\vec{n}_{\text{large}}$. Zeros are set to one, the minimum pool size. Figure 8 shows these $2\times 1,000$ parameter estimates as compared to the true parameter values and those for which the true size vectors $\vec{n}_{\text{small}}$ and $\vec{n}_{\text{large}}$ were used as input.

The violins of the estimates for the smaller cell pools (based on $\vec{n}_{\text{small}}$) indicate that the estimates of $p$ and $\mu_{1}$ are fairly accurate, but the estimates of $\mu_{2}$ have large variance, and $\sigma$ is overestimated in all 1,000 runs. This is plausible as population 1 (the one with higher log-mean gene expression) is only present on average in 20% of the cells; even when misspecifying the pool sizes, the cells of population 1 are still detectable since this is the population responsible for most gene expression. Consequently, all remaining cells are assigned to population 2, which has lower or even almost no expression. If the pool size is assumed too low, this second population will be estimated to have on average a higher expression; if it is assumed too large, the second population will be estimated to have a lower expression. This leads to a broader distribution and thus an overestimation of $\sigma$.

The results for the larger cell pools (based on $\vec{n}_{\text{large}}$) show a similar pattern. In this case, however, the impact of misspecification is less visible, as also confirmed by additional simulations in Appendix F. For large cell pools, the averaging effect across cells is strong anyway and in that sense more robust. In the study here, the $\sigma$ parameter is often even better estimated when using a misspecified pool size vector than when using the true one.

Taken together, \pkgstochprofML can be used even if exact pool sizes are unknown. In that case, the numbers should be approximated as well as possible.

The \pkgstochprofML algorithm estimates the assumed parameterized single-cell distributions underlying the samples and; as described in Section 2.5, we can select the most appropriate number of cell populations using the BIC. Assume we have performed this estimation for samples from two different groups, cases and controls. One may in practice then want to know whether the inferred single-cell populations are substantially different between the two groups, e.g.  in case the estimated log-means $\hat{\mu}_{\text{cases}}$ and $\hat{\mu}_{\text{controls}}$ are close to each other. A related question is whether the difference is biologically relevant.

We hence seek a method that can judge statistical significance and potentially reject the null hypothesis that two single-cell populations are the same; and at the same time allow the interpretation of similarity. Direct application of Kolmogorov-Smirnov or likelihood-ratio tests to the observed data is impossible here since the single-cell data is unobserved: We only measure the overall gene expression of pools of cells. Calculation of the Kullback-Leibler divergence of the two distributions would be possible; however, it is not target-oriented for our application where we seek an interpretable measure of similarity rather than a comparison between more than two population densities.

For our purposes, we use a simple intuitive measure of similarity — the overlap of two PDFs, that is the intersection of the areas under both PDF curves:

for two continuous one-dimensional PDFs $f$ and $g$ (see also Pastore and Calcagnì, 2019). The overlap lies between zero and one, with zero indicating maximum dissimilarity and one implying (almost sure) equality. In our case, we are particularly interested in the overlap of two lognormal PDFs: {Code} OVL_LN_LN <- function(mu_1, mu_2, sigma_1, sigma_2) f1 <- function(x)dlnorm(x, meanlog = mu_1, sdlog = sigma_1) f2 <- function(x)dlnorm(x, meanlog = mu_2, sdlog = sigma_2) f3 <- function(x)pmin(f1(x), f2(x)) integrate(f3, lower = 0, upper = Inf, abs.tol = 0)$value$

Figure 9 shows examples of such overlaps. Here, the overlap ranges from 12% for two quite different distributions to 86% for two seemingly similar distributions. The question is where to draw a cutoff, that is, at what point we decide to label two distributions as different. Current literature considers two cases: Either the parametric case (e.g. Inman and Bradley, 1989), where both distributions are given by their distribution families and parameter values; or the nonparametric case (e.g. Pastore and Calcagnì, 2019), where observations (but no theoretical distributions) are available for the two populations. Our application builds a third case: On the one hand, we want to compare two parametric distributions, but the model parameters are just given as estimates based on (potentially small) datasets, thus they are uncertain; on the other hand, we do not directly observe the single-cell gene expression but just the pooled one.

To address this issue, we suggest to again take into account the original data that led to the estimated parametric PDFs. As an example, assume that we consider two sets of pooled gene expression, one for a group of cases and one for a group of controls. In both groups, pooled gene expression is available as 10-cell measurements, but the two groups differ in sample size. Let’s say the cases contain 50 samples and the controls 100. We assume the LN-LN model with two populations and estimate the mixture and population parameters using the \pkgstochprofML algorithm separately for each group, leading to estimates $\hat{p}_{\text{cases}}$, $\hat{\mu}_{1,\text{cases}}$, $\hat{\mu}_{2,\text{cases}}$, $\hat{\sigma}_{\text{cases}}$ and $\hat{p}_{\text{controls}}$, $\hat{\mu}_{1,\text{controls}}$, $\hat{\mu}_{2,\text{controls}}$, $\hat{\sigma}_{\text{controls}}$. We now aim to assess whether the first populations in both groups have identical characteristics, i. e.  whether $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{cases}},\hat{\sigma}^{2}_{\text{cases}})$ and $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{controls}},\hat{\sigma}^{2}_{\text{controls}})$ are the same.

Figure 9 displays the single-cell PDFs of the first population and their overlaps for various values of the estimates. For example, in Figure 9D, the orange curve shows the single-cell PDF of population 1 inferred from the cases, yielding $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{cases}}=2.10,\hat{\sigma}^{2}_{\text{cases}}=0.19^{2})$, and the blue one shows the inferred single-cell PDF of population 1 from the controls, $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{controls}}=2.03,\hat{\sigma}^{2}_{\text{controls}}=0.20^{2})$. The overlap of these two inferred PDFs equals 86%.

We now aim to test the null hypothesis that these populations are the same versus the experimental hypothesis that they are different. We perform a sampling-based test: Taking into account the inferred population probabilities $\hat{p}_{\text{cases}}$ and $\hat{p}_{\text{controls}}$ and the number of samples and cells in the data, we can estimate the number of cells which the estimates $\hat{\boldsymbol{\theta}}_{\text{cases}}$ and $\hat{\boldsymbol{\theta}}_{\text{controls}}$ relied on. The larger this cell number, the less expected uncertainty about the estimated population distributions $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{cases}},\hat{\sigma}^{2}_{\text{cases}})$ and $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{controls}},\hat{\sigma}^{2}_{\text{controls}})$ (neglecting the impact of pool sizes).

In our example, let $\hat{p}_{\text{cases}}=12\%$. Then, approximately 12% of the 500 cells from the cases group ($50$ $\times$ 10-cell samples) belonged to population 1, that is 60 cells. For $\hat{p}_{\text{controls}}=20\%$, 200 cells were expected to be from the first population (that is 20% of 1,000 cells, coming from the 100 $\times$ 10-cell measurements for the controls). In our procedure, we compare parameter estimates that are based on the respective numbers of single cells, i. e.  60 cells for cases and 200 cells for controls. We perform the following steps:

• •

  Calculate $\text{OVL}_{\text{original}}$, the overlap of the PDFs of $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{cases}}=2.10,\hat{\sigma}^{2}_{\text{cases}}=0.19^{2})$ and $\mathcal{L}\mathcal{N}(\hat{\mu}_{1,\text{controls}}=2.03,\hat{\sigma}_{\text{controls}}^{2}=0.20^{2})$.

• •

  Under the null hypothesis, the two distributions are identical. We approximate the parameters of this identical distribution as $\tilde{\mu}_{1,\text{mean}}=(\hat{\mu}_{1,\text{cases}}+\hat{\mu}_{1,\text{controls}})/2$ and $\tilde{\sigma}_{\text{mean}}=(\hat{\sigma}_{\text{cases}}+\hat{\sigma}_{\text{controls}})/2$.

• •

  Repeat $N=1,000$ times:

  • –

    Draw dataset $A$ of size $60$ from $\mathcal{L}\mathcal{N}(\tilde{\mu}_{1,\text{mean}},\tilde{\sigma}_{\text{mean}}^{2})$.

  • –

    Draw dataset $B$ of size $200$ from $\mathcal{L}\mathcal{N}(\tilde{\mu}_{1,\text{mean}},\tilde{\sigma}_{\text{mean}}^{2})$.

  • –

    Estimate the log-mean and log-sd for these two datasets using the method of maximum likelihood, yielding $\hat{\mu}_{A}$, $\hat{\sigma}_{A}$, $\hat{\mu}_{B}$ and $\hat{\sigma}_{B}$.

  • –

    Calculate $\text{OVL}\left(f_{\mathcal{L}\mathcal{N}(\hat{\mu}_{A},\hat{\sigma}_{A}^{2})},f_{\mathcal{L}\mathcal{N}(\hat{\mu}_{B},\hat{\sigma}_{B}^{2})}\right)$.

• •

  Sort the $N$ overlap values and select the empirical 5% quantile $\text{OVL}_{0.05}$.

• •

  Compare the overlap from the original data to this quantile:

  • –

    If $\text{OVL}_{original}\leq\text{OVL}_{0.05}$, the null hypothesis that both populations are the same can be rejected.

  • –

    If $\text{OVL}_{original}>\text{OVL}_{0.05}$, the null hypothesis cannot be rejected.

This proceeding is related to the idea of parametric bootstrap with the difference that our original data is on the $n$-cell level and the parametrically simulated data is on the single-cell level.

The left panel of Figure 10 shows one outcome of the above-described procedure (i. e.  the stochastic, sampling-based algorithm was run once) with the above-specified values of the parameter estimates. Here, $\text{OVL}_{original}$ lies in the critical range such that we reject the null hypothesis that the gene expression of the populations in question stem from the same lognormal distribution. We thus assume a difference here. The right panel of Figure 10 demonstrates the importance of taking into account the number of cells which the original estimates were based on: Here, we show one outcome of the above described steps, but this time we assume that for the control group there were only 30 10-cell samples (i.e. 300 cells in total). With the same population fraction as before ($\hat{p}_{\text{controls}}=20\%$), the datasets $B$ now contain only 60 cells. Here, the value $\text{OVL}_{original}$ does not fall into the critical range, and therefore we would not reject the null hypothesis that the two populations of interest are the same.

When testing for heterogeneity for several genes simultaneously, multiple testing issues should be taken into account. However, genes will in general not be independent from each other.

The \pkgstochprofML algorithm estimates the parameters of the mixture model, i. e. — in case of at least two populations — the probability for each cell within a pool to fall into the specific populations. It does not reveal the individual pool compositions. In some applications, however, exactly this information is of particular interest. Here, we present how one can infer likely population compositions of a particular cell pool.

This is done in a two-step approach via conditional prediction: First, one estimates the model parameters from the observed pooled gene expression, i. e.  one obtains an estimate $\hat{\boldsymbol{\theta}}$ of ${\boldsymbol{\theta}}$. Then, one assumes that $\boldsymbol{\theta}$ equals $\hat{\boldsymbol{\theta}}$ and derives the most probable population composition via maximizing the conditional probability of a specific composition given the pooled gene expression (for calculations, see Appendix D).

We evaluate this procedure via a simulation study. As before, we simulate data using the \pkgstochprofML package. In particular, we use the LN-LN model with two populations with parameters $\boldsymbol{p}=(0.2,0.8)$, $\boldsymbol{\mu}=(2,0)$ and $\sigma=0.2$. Each simulated measurement shall contain the pooled expression of $n=5$ cells, and we sample $k=100$ such measurements. We store the original true cell pool compositions from the data simulation step in order to later compare the composition predictions to the ground truth.

Having generated the synthetic data, we apply \pkgstochprofML to estimate the model parameters $\boldsymbol{p}$, $\boldsymbol{\mu}$ and $\sigma$. Figure 11 shows a histogram of one simulated data set along with the PDF of the true population mixture and the PDF of the estimated population mixture (that is the LN-LN model with parameters $\boldsymbol{\hat{p}}=(0.14,0.86)$, $\boldsymbol{\hat{\mu}}=(2.04,0)$ and $\hat{\sigma}=0.20$).

Next, we calculate the conditional probability mass function (PMF; see Appendix D for details) for each possible population composition conditioned on the particular pooled gene expression measurement. Figure 12A and Table 1 show results for the first six (out of 100) pooled measurements.

In particular, Figure 12A displays the conditional PMF of all possible compositions (i. e.  $k$ times population 1 and $5-k$ times population 2 for $k\in\{0,1,\ldots,5\}$). Blue bars stand for these probabilities when $\hat{\boldsymbol{\theta}}$ is used as model parameter value. Orange stands for the hypothetical case where the true value $\boldsymbol{\theta}$ is known and used. These two scenarios are in good agreement with each other.