Count Likelihoods for RNA and Other Sequencing Blocks

Why Sequencing Counts Need Their Own Likelihood

RNA-sequencing and other sequencing assays record nonnegative integer counts: for example, the number of reads assigned to a gene in a sample. Such counts are not continuous measurements, so a Gaussian likelihood for an omics profile is generally not a natural starting point. In a clustering model, the count likelihood describes how plausible a sample's observed counts are under each cluster's profile, after accounting for sequencing depth.

Write $Y_{igb}$ for the count for sample $i$, feature $g$, and sequencing block $b$. A block might be an RNA-sequencing assay or another assay whose observations are counts. Let $z_i=k$ indicate that sample $i$ belongs to cluster $k$. Use $E_{ib}$ for a known, positive sample-specific exposure or library-size scaling factor, and $\lambda_{gbk}>0$ for the expected count rate per unit exposure for feature $g$ in block $b$ and cluster $k$. The expected count is then

\mu_{igb}=E_{ib}\lambda_{gbk}\quad\text{when }z_i=k.

The exposure is treated as known here: it adjusts the expected count for differences in sequencing depth. This setup assumes that the exposure is a suitable scaling factor and that the cluster-specific rate captures the remaining systematic differences in expected abundance. It does not, by itself, resolve every compositional issue that can arise when sequencing produces a fixed total number of reads.

Poisson Likelihood Within a Cluster

Definition: Under a Poisson count model, conditional on the cluster assignment and model parameters, $Y_{igb}$ is Poisson with mean $\mu_{igb}$: $$ P(Y_{igb}=y\mid z_i=k,E_{ib},\lambda_{gbk}) =\frac{\exp(-\mu_{igb})\mu_{igb}^{y}}{y!}, \qquad y=0,1,2,\ldots. $$ The model has conditional mean and variance both equal to $\mu_{igb}$.

For a cluster's sequencing block, a common working assumption is conditional independence of counts across samples and features given the cluster assignments, exposures, and rates. Under that assumption, the block likelihood is the product of the relevant Poisson probabilities. For one sample, its log likelihood contribution is a sum over features:

\ell_{ib}(k) =\sum_g\left[ Y_{igb}\log(E_{ib}\lambda_{gbk}) -E_{ib}\lambda_{gbk} -\log(Y_{igb}!) \right].

Terms that do not depend on $k$, such as $\log(Y_{igb}!)$, may be omitted when comparing cluster scores for that same sample. They must be retained if an absolute likelihood is needed. Conditional independence is a modeling assumption, not a consequence of using a Poisson distribution: correlated features or remaining sample-level effects can make the product likelihood too concentrated.

Worked Example: Estimating a Cluster-Specific Poisson Rate

Suppose three hypothetical samples are assigned to the same cluster for one RNA feature. Their counts are $12,18,10$, and their known exposures are $1,1.5,0.75$. Assume independent Poisson counts with means $E_i\lambda$, sharing one unknown rate $\lambda$. The exposures are positive, the values are nonnegative integers, and the shared-rate and conditional-independence assumptions are stipulated for this calculation.

The log likelihood, omitting no terms that depend on $\lambda$, is

\ell(\lambda)=\sum_{i=1}^{3} \left[Y_i\log(E_i\lambda)-E_i\lambda-\log(Y_i!)\right].

Differentiating with respect to $\lambda$ gives

\frac{d\ell}{d\lambda} =\frac{Y_1+Y_2+Y_3}{\lambda}-(E_1+E_2+E_3).

Setting this derivative to zero yields $\hat{\lambda}=\sum_iY_i/\sum_iE_i$. Substituting the hypothetical observations,

\hat{\lambda}=\frac{12+18+10}{1+1.5+0.75} =\frac{40}{3.25}=12.31

expected counts per unit exposure. The fitted mean counts are $1(12.31)=12.31$, $1.5(12.31)=18.46$, and $0.75(12.31)=9.23$ reads, respectively. This is a maximum-likelihood estimate under the stated Poisson model; with only three observations, the example does not establish that the Poisson variance assumption is adequate.

Overdispersion and the Negative-Binomial Model

The Poisson assumption ties conditional variance to conditional mean. Sequencing counts can vary more than that, even after accounting for exposure and cluster-specific mean profiles. This is called overdispersion relative to the Poisson model. A negative-binomial likelihood allows the conditional variance to exceed the mean.

Definition: Parameterize the negative-binomial distribution by its conditional mean $\mu>0$ and dispersion $\phi\geq 0$. For $\phi>0$, $$ P(Y=y\mid\mu,\phi) =\frac{\Gamma(y+1/\phi)}{\Gamma(1/\phi)y!} \left(\frac{1}{1+\phi\mu}\right)^{1/\phi} \left(\frac{\phi\mu}{1+\phi\mu}\right)^y, \qquad y=0,1,2,\ldots. $$ Its conditional mean is $\mu$, and its conditional variance is $\mu+\phi\mu^2$. The Poisson distribution is recovered as the limit when $\phi$ tends to zero.

The parameter $\phi$ controls extra variability: for any $\phi>0$, the variance exceeds the mean. In a cluster model, one may use $\mu_{igb}=E_{ib}\lambda_{gbk}$ in this likelihood and specify a dispersion parameter for the feature, block, or another modeling level. The appropriate structure depends on the data and model; estimating many separate dispersion parameters can be unstable when observations are sparse.

A useful justification for this form is a Poisson rate that varies across observations. Let $U$ be a positive random multiplier with a gamma distribution of shape $1/\phi$ and rate $1/\phi$, so $E(U)=1$ and $\operatorname{Var}(U)=\phi$. Suppose $Y\mid U\sim\operatorname{Poisson}(\mu U)$. Integrating over $U$ gives

\begin{aligned} P(Y=y) &=\int_0^\infty \frac{\exp(-\mu u)(\mu u)^y}{y!} \frac{(1/\phi)^{1/\phi}}{\Gamma(1/\phi)} u^{1/\phi-1}\exp(-u/\phi)\,du\\ &=\frac{\Gamma(y+1/\phi)}{\Gamma(1/\phi)y!} \left(\frac{1}{1+\phi\mu}\right)^{1/\phi} \left(\frac{\phi\mu}{1+\phi\mu}\right)^y. \end{aligned}

For the variance, the law of total variance gives a direct check:

\begin{aligned} E(Y)&=E\{E(Y\mid U)\}=E(\mu U)=\mu,\\ \operatorname{Var}(Y) &=E\{\operatorname{Var}(Y\mid U)\} +\operatorname{Var}\{E(Y\mid U)\}\\ &=E(\mu U)+\operatorname{Var}(\mu U) =\mu+\phi\mu^2. \end{aligned}

Thus the negative-binomial model can represent heterogeneity in the underlying rate. The derivation does not prove that gamma rate variation is the true biological or technical mechanism; it explains the distribution and variance form used by the likelihood.

Worked Example: Extra-Poisson Variation and a Negative-Binomial Probability

Suppose repeated hypothetical observations for a feature have estimated conditional mean $\mu=4$ reads and variance $12$ reads squared. These summaries are assumed to describe comparable observations after exposure adjustment. The Poisson model would predict variance $4$, so the variance is larger than its Poisson value. Matching the negative-binomial variance gives

12=4+\phi(4)^2 \quad\Longrightarrow\quad \hat{\phi}=\frac{12-4}{16}=0.5.

For a count of $y=8$, $1/\phi=2$. The negative-binomial probability is

\begin{aligned} P(Y=8\mid\mu=4,\phi=0.5) &=\frac{\Gamma(10)}{\Gamma(2)8!} \left(\frac{1}{3}\right)^2 \left(\frac{2}{3}\right)^8\\ &=9\cdot\frac{1}{9}\cdot\left(\frac{2}{3}\right)^8 \approx 0.03902. \end{aligned}

The mean-variance calculation is a moment-based illustration, not a claim that two summary statistics identify a reliable dispersion estimate in a small data set. In practice, dispersion estimation and checking require adequate replication and attention to how the model pools information across features.

Using Count Likelihoods to Compare Cluster Profiles

For a given sample, each candidate cluster supplies a mean for each feature through its rate profile and the sample's exposure. Under conditional feature independence, multiply the feature probabilities, or equivalently add their log probabilities, to obtain a block score. Cluster scores are likelihood contributions; a posterior cluster probability also depends on cluster prevalence and any other model components. The calculation below uses equal prior probabilities solely to make that distinction concrete.

Worked Example: Comparing Two Negative-Binomial Cluster Scores

Consider one hypothetical sample with exposure $1$ and two counts, $Y_1=8$ and $Y_2=1$. Suppose two candidate clusters have mean profiles $(4,2)$ and $(8,1)$, respectively, and both features use $\phi=0.5$. Assume counts are conditionally independent across these features and the exposure and profile means are correctly specified. For the first feature, the probabilities under means 4 and 8 are

\begin{aligned} P(Y_1=8\mid\mu=4,\phi=0.5) &=9(1/3)^2(2/3)^8\approx 0.03902,\\ P(Y_1=8\mid\mu=8,\phi=0.5) &=9(1/5)^2(4/5)^8\approx 0.06040. \end{aligned}

For the second feature, the probability of $Y_2=1$ under mean 2 is $2(1/2)^2(1/2)=0.2500$. Under mean 1 it is $2(2/3)^2(1/3)=8/27\approx0.2963$. Therefore the two block likelihoods are

\begin{aligned} L_A&=(0.03902)(0.2500)\approx0.009755,\\ L_B&=(0.06040)(0.2963)\approx0.01790. \end{aligned}

The likelihood ratio is $L_B/L_A\approx1.834$, so these two counts are about 1.83 times as likely under cluster B's specified profile as under cluster A's. If the two clusters have equal prior probabilities and these are the only modeled data contributions, the corresponding posterior probability for B is $L_B/(L_A+L_B)\approx0.647$. This is a model-based comparison, not evidence that the sample truly belongs to B with certainty.

Assumptions, Diagnostics, and Alternatives

Choosing between Poisson and negative-binomial likelihoods is not just a choice of formula. It requires checking whether the mean structure, exposure adjustment, and dependence assumptions are plausible. Useful diagnostics compare observed and model-predicted variability within groups that are comparable in exposure and fitted mean. A simple standardized residual is

r_{igb}=\frac{Y_{igb}-\hat{\mu}_{igb}} {\sqrt{\hat{\mu}_{igb}+\hat{\phi}\hat{\mu}_{igb}^{\,2}}},

using $\hat{\phi}=0$ for a Poisson fit. Residual patterns across samples, features, fitted means, and sequencing depth can reveal systematic lack of fit. In particular, a Poisson model may show more large residuals than expected if it understates variability. A negative-binomial model may still fail when there are excess zeros, unusually heavy tails, unmodeled sample effects, or correlated features. A residual diagnostic is a prompt to investigate, not a formal proof that the model is correct or incorrect.

Model choice	Conditional variance	When it can be useful	Main caution
Poisson	$\mu$	Counts are compatible with equidispersion after exposure and profile effects.	Underestimated variation can make cluster likelihoods too sharp.
Negative binomial	$\mu+\phi\mu^2$	Counts show extra variation that grows with the mean.	Dispersion assumptions and estimates need checking; it does not model every dependence pattern.
Gaussian model for transformed measurements	Specified by a continuous-data model	May suit appropriately transformed or normalized continuous measurements.	It is not the same as a raw-count likelihood and can handle zeros or discreteness poorly.

The Gaussian likelihood for continuous omics profiles is therefore not automatically interchangeable with a count likelihood. Likewise, Poisson and negative-binomial models treat counts as conditionally independent in their basic product forms; an alternative that explicitly models dependence or zero inflation may be needed if diagnostics indicate those features. Such alternatives change the model rather than simply correcting the variance formula.

A common mistake is to see a large count and assume the Poisson likelihood must be suitable because counts are integers. The correction is that count-valued support does not establish equidispersion: under the Poisson model, the conditional variance must equal the mean. Another mistake is to treat a negative-binomial fit as proof of a particular biological mechanism. Its extra variance is a flexible distributional feature, and its good fit alone does not identify why that variation arose.

Key takeaway: Use the Poisson likelihood when its conditional mean and variance assumptions are credible; use a negative-binomial likelihood when extra-Poisson variation is supported and its mean-variance form is reasonable. In either case, state the exposure and conditional-independence assumptions, and evaluate fit rather than relying on the count scale alone.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For a Poisson count with mean $6$, what conditional variance does the model imply? What variance would a negative-binomial model with $\phi=0.25$ imply at the same mean?
A sample has twice the exposure of another, while both belong to the same cluster and have the same feature rate. How do their expected counts compare?
For a negative-binomial count with mean $\mu=3$ and dispersion $\phi=1/3$, calculate the conditional variance.
Why can a negative-binomial likelihood be preferable to a Poisson likelihood, and what does that preference not establish?

Solutions. (1) Poisson variance is $6$. Negative-binomial variance is $6+0.25(6^2)=15$. (2) Since $\mu=E\lambda$, doubling $E$ doubles the expected count when $\lambda$ is unchanged. (3) The variance is $3+(1/3)(3^2)=6$. (4) It can represent variance greater than the mean; that fit does not identify the cause of extra variation or guarantee that all dependence and zero patterns are modeled correctly.

Model choice	Conditional variance	When it can be useful	Main caution
Poisson	\(\mu\)	Counts are compatible with equidispersion after exposure and profile effects.	Underestimated variation can make cluster likelihoods too sharp.
Negative binomial	\(\mu+\phi\mu^2\)	Counts show extra variation that grows with the mean.	Dispersion assumptions and estimates need checking; it does not model every dependence pattern.
Gaussian model for transformed measurements	Specified by a continuous-data model	May suit appropriately transformed or normalized continuous measurements.	It is not the same as a raw-count likelihood and can handle zeros or discreteness poorly.