Estimating Proportions From Count References

From Reference Profiles to Spot Proportions

A spatial spot can contain transcripts from several cell types, but its observed counts combine those contributions. Estimating the proportions therefore requires both a reference describing expression by cell type and a model for how the spot’s counts arise from a mixture of those references. A regression on normalized profiles and a likelihood for counts can use similar mixture structure, yet they weight discrepancies differently and handle library size and variation differently.

Here, \(G\) genes are measured in a spot and \(K\) cell types are represented in a reference. Let \(y_g\) be the count for gene \(g\) in the spot, and let \(N=\sum_{g=1}^G y_g\) be its total count over the selected genes. For cell type \(k\), let \(R_{gk}\) be the reference count for gene \(g\), and \(M_k=\sum_g R_{gk}\) the corresponding total over those genes. A simple plug-in reference profile is \(q_{gk}=R_{gk}/M_k\).

Definition: A reference profile \(q_{\cdot k}\) is a vector of nonnegative gene proportions for cell type \(k\), with \(\sum_g q_{gk}=1\). A spot mixture proportion vector \(p=(p_1,\ldots,p_K)\) satisfies \(p_k\geq 0\) and \(\sum_k p_k=1\). Under the basic mixture model, the expected gene proportions in a spot are \(q_g(p)=\sum_k p_k q_{gk}\).

This setup treats \(p_k\) as the contribution of type \(k\) to the selected-gene transcript mixture. It need not equal the fraction of cells of type \(k\): cell types can differ in RNA content, capture, and expression. Translating transcript contributions into cell fractions requires additional assumptions or calibration.

Normalized-Profile Regression and Its Count Likelihood

The constrained regression approach considered in the earlier material on constrained regression for cell-type proportions fits the observed normalized profile \(y_g/N\) to the reference profiles. In matrix form, let \(Q\) have entries \(q_{gk}\). Ordinary least squares estimates \(p\) by minimizing the sum of squared differences between \(y/N\) and \(Qp\), subject to the simplex constraints. This is convenient, but gives equal squared-error weight to each gene-proportion residual, regardless of its expected count or sampling variability.

\widehat p_{\mathrm{LS}} = \underset{p_k\geq 0,\ \sum_k p_k=1}{\operatorname{argmin}} \sum_{g=1}^G\left(\frac{y_g}{N}-\sum_{k=1}^K p_kq_{gk}\right)^2

A count likelihood instead specifies a sampling distribution. Conditional on \(N\), the simplest model is \(y\sim\operatorname{Multinomial}(N,q(p))\). Its assumptions are that the selected-gene counts are categorical draws with probabilities \(q_g(p)\), that the total \(N\) is fixed by conditioning, and that the draws are independent given \(p\) and the reference profiles. In practice, the independence assumption can be too strong because biological and technical variation may make spots more variable than this model allows.

P(y\mid N,p,Q) = \frac{N!}{\prod_{g=1}^G y_g!} \prod_{g=1}^G q_g(p)^{y_g}, \qquad \ell(p)=C+\sum_{g=1}^G y_g\log q_g(p)

Here \(C\) does not depend on \(p\). The multinomial estimate maximizes \(\ell(p)\) over the simplex. Because each \(q_g(p)\) is linear in \(p\) and the logarithm is concave, \(\ell(p)\) is concave wherever the probabilities are positive. Thus a constrained maximum is a global maximum, although it may occur on the boundary. If a gene has \(q_g(p)=0\) but \(y_g>0\), that mixture has zero likelihood.

Worked Example: Likelihood and Least Squares Give Different Fits

Suppose two reference types have three-gene profiles \(q_{\cdot A}=(0.6,0.3,0.1)\) and \(q_{\cdot B}=(0.1,0.3,0.6)\). A hypothetical spot has counts \(y=(45,30,25)\), so \(N=100\). Write \(p\) for the proportion of type A; the type B proportion is \(1-p\).

The mixture probabilities are \(q_1(p)=0.1+0.5p\), \(q_2(p)=0.3\), and \(q_3(p)=0.6-0.5p\). For least squares, the first and third observed proportions are \(0.45\) and \(0.25\). Each gives \(p=0.70\): \(0.1+0.5p=0.45\) and \(0.6-0.5p=0.25\). The second gene has no ability to distinguish the types because both references assign it probability \(0.3\). Therefore \(\widehat p_{\mathrm{LS}}=0.70\).

For the multinomial likelihood, differentiate the terms that depend on \(p\):

\frac{d\ell}{dp} = \frac{45(0.5)}{0.1+0.5p} - \frac{25(0.5)}{0.6-0.5p}

Setting this derivative to zero gives \(45/(0.1+0.5p)=25/(0.6-0.5p)\). Thus \(45(0.6-0.5p)=25(0.1+0.5p)\), so \(27-22.5p=2.5+12.5p\), and \(\widehat p_{\mathrm{ML}}=24.5/35=0.7000\). In this particular example the estimates agree. The likelihood weights the fit through the observed counts and predicted probabilities; this agreement follows from the chosen numbers, not from general equivalence between the methods.

How Library Size Enters

Conditioning on \(N\) means the multinomial model does not try to explain why one spot has more total reads than another. For fixed observed proportions, multiplying every count by the same factor leaves the maximum-likelihood proportion estimate unchanged: it multiplies the log-likelihood’s \(p\)-dependent part by that factor. But a larger \(N\) makes the count composition more informative under the model, so the likelihood is more sharply concentrated around its maximizer.

An equivalent Poisson formulation makes the exposure role explicit when the expected total is specified. If \(L\) is a spot-level library-size exposure and the expected counts are \(\mu_g=Lq_g(p)\), independent Poisson counts have a total expected count \(L\), since \(\sum_g q_g(p)=1\). Conditioning those counts on their total produces the multinomial model. If the total is not adequately represented by this exposure, a model for total count or other offsets may be needed; simply normalizing counts does not diagnose that issue.

Worked Example: Same Composition, Different Library Sizes

Consider two hypothetical spots with two-gene counts \(y^{(1)}=(70,30)\) and \(y^{(2)}=(700,300)\). Both have observed proportions \((0.7,0.3)\). Let type A have profile \((0.8,0.2)\), type B have profile \((0.2,0.8)\), and let \(p\) denote the type A contribution. Then \(q_1(p)=0.2+0.6p\) and \(q_2(p)=0.8-0.6p\).

For either spot, maximizing the multinomial likelihood matches the observed first-gene proportion: \(0.2+0.6p=0.7\). Thus \(p=0.5/0.6=0.8333\). Increasing all counts tenfold does not change this estimate. It does increase the log-likelihood’s \(p\)-dependent terms tenfold, so under the multinomial assumptions it provides stronger evidence about the composition than the smaller library. This does not mean a tenfold larger library removes systematic reference bias.

Allowing Extra-Multinomial Variation

If spots vary more than a multinomial model predicts, a Dirichlet-multinomial model allows the gene probabilities themselves to vary around \(q(p)\). Let \(\phi>0\) be a concentration parameter and set the Dirichlet parameters to \(\phi q_g(p)\). Conditional on those probabilities, counts are multinomial; after integrating over the Dirichlet distribution, the count vector has a Dirichlet-multinomial distribution. Smaller \(\phi\) means greater overdispersion, while the limit \(\phi\to\infty\) recovers the multinomial model.

P(y\mid N,p,Q,\phi) = \frac{N!}{\prod_g y_g!} \frac{\Gamma(\phi)}{\Gamma(N+\phi)} \prod_{g=1}^G \frac{\Gamma(y_g+\phi q_g(p))}{\Gamma(\phi q_g(p))}

For a single gene count, the Dirichlet-multinomial variance is \(Nq_g(p)(1-q_g(p))(N+\phi)/(1+\phi)\), larger than the multinomial variance \(Nq_g(p)(1-q_g(p))\) when \(N>1\) and \(\phi\) is finite. The count likelihood now depends on \(\phi\) as well as the mixture proportions. Estimating proportions with this model requires specifying or estimating \(\phi\); it is not merely a change to the standard error of the multinomial fit.

Worked Example: Quantifying Extra-Multinomial Variation

Suppose a hypothetical spot has \(N=100\) counts across two categories, each with fitted probability \(q=0.5\). Under a multinomial model, the variance of the first category’s count is \(Nq(1-q)=100(0.5)(0.5)=25\).

Now use a Dirichlet-multinomial model with \(\phi=20\). The variance is \(100(0.5)(0.5)(100+20)/(1+20)=25(120/21)=142.9\) counts squared. Its standard deviation is \(\sqrt{142.9}=11.95\) counts, compared with \(\sqrt{25}=5.000\) under the multinomial model. The larger variance reflects the possibility that the spot-level probabilities vary around their mean. It does not change the fitted mean probability \(0.5\) in this calculation; it changes the distribution of counts around that mean.

Choosing and Checking the Fit

Use normalized-profile least squares when a simple constrained approximation is adequate and the scale of residuals across genes is not the main concern. Use a multinomial likelihood when the selected-gene counts can reasonably be treated as draws from the mixture profile and conditioning on the total is appropriate. Use an overdispersed likelihood when replicate spots or other diagnostics show more variability than the multinomial model predicts. Independent negative-binomial gene models are another alternative when gene-specific count means and dispersion are central, but they require a coherent mean and exposure specification; they do not automatically enforce a single multinomial composition.

Check whether fitted profiles reproduce the observed gene proportions, whether estimates repeatedly land on simplex boundaries, and whether residual discrepancies are concentrated in particular genes or spots. A boundary estimate can be a genuine consequence of the data, but it can also indicate poor separation among reference profiles or an inadequate reference. Very low reference counts can produce unstable \(q_{gk}\); zero reference entries paired with positive spot counts can make some mixtures impossible under an unmodified likelihood. Filtering, smoothing, or adding pseudocounts changes the model and should be reported rather than treated as a neutral computational fix.

A common mistake is to conclude that a larger library must produce a different proportion estimate, or that normalization makes library size irrelevant. Under the conditional multinomial model, scaling all counts equally leaves the point estimate unchanged, but changes the information in the likelihood. Another mistake is to interpret overdispersion as a change in the mean mixture: in the Dirichlet-multinomial model it describes additional variation around the mean probabilities. Neither model repairs a systematically mismatched reference; a likelihood can produce a precise fit to the wrong profiles.

Key takeaway: Count-based deconvolution combines a simplex-constrained mixture of reference profiles with an explicit sampling model. Conditioning on library size yields the multinomial fit; allowing spot-level probability variation yields an overdispersed alternative. These choices affect the weighting and variability of the fit, not the need for suitable references.

Check Your Understanding

Answer each question, then compare with the brief solutions.

If every count in a spot is multiplied by five, what happens to the multinomial maximum-likelihood proportion estimate under the same reference profiles?
Two reference profiles assign the same probability to a gene. What can that gene contribute to distinguishing those two types?
For \(N=50\), \(q=0.4\), and \(\phi=24\), calculate the Dirichlet-multinomial variance of that gene’s count.
What does an observed count in a gene with zero fitted mixture probability imply for the ordinary multinomial likelihood?

Solutions. (1) It is unchanged, because scaling all counts multiplies the \(p\)-dependent log-likelihood by a constant; the larger library does increase its curvature. (2) Nothing, because that gene’s mixture probability is constant as the proportions of those types vary. (3) \(50(0.4)(0.6)(50+24)/(1+24)=12(74/25)=35.52\) counts squared. (4) That mixture has zero likelihood, since it assigns probability zero to an observed count.