Worked Independent-Group Pseudobulk Analysis

Why the Independent-Group Model Is Different

In a two-condition study with independent donors, each donor contributes a pseudobulk count profile for the cell type being analyzed. The scientific comparison is between donor groups, not between individual cells: cells from the same donor share biological and technical features, so treating them as independent replicates would overstate the information in the data. The independent-group analysis therefore models one count per donor, gene, and condition.

For a gene, the response \(Y_i\) is the summed count for donor \(i\). The library size, or an effective library size that includes a normalization factor, accounts for differences in sequencing exposure. A condition coefficient then describes the difference in expected normalized expression between the two donor groups. It is an association between condition group and expression; by itself it does not establish a causal effect.

Model, Assumptions, and the Contrast

A standard gene-wise model is negative binomial (NB). Use the parameterization in which the conditional mean is \(\mu_i\) and the variance is \(\mu_i+\phi\mu_i^2\), where \(\phi\) is the dispersion. This allows the variance to exceed the mean, as is common for biological count data. Let \(x_i=0\) for the reference condition and \(x_i=1\) for the other condition.

The term \(\log(E_i)\) is an offset: its coefficient is fixed at one rather than estimated. The condition coefficient \(\beta_1\) is the log ratio of expected expression per unit of effective library size in the second condition relative to the reference, under this model. Thus \(\exp(\beta_1)\) is the normalized expression ratio, and \(100\{\exp(\beta_1)-1\}\%\) is the corresponding percent difference.

The comparison is usually framed as \(H_0:\beta_1=0\) against \(H_a:\beta_1\ne 0\). A prespecified one-sided alternative is possible when scientifically justified, but the direction must not be selected after inspecting the result. The contrast is \(c^\mathsf{T}\beta\), where \(\beta=(\beta_0,\beta_1)^\mathsf{T}\) and, for this comparison, \(c=(0,1)^\mathsf{T}\). Its estimated standard error is \(\sqrt{c^\mathsf{T}\widehat{\operatorname{Cov}}(\hat\beta)c}\). A Wald statistic divides the estimated contrast by that standard error and is compared with a standard normal reference distribution in a large-sample analysis.

The model relies on independent donors, a suitable log-linear mean structure, correctly specified exposures, and an adequate NB variance model. Counts are conditionally independent across donors given the model covariates; within-donor cells have already been aggregated and are not counted as independent replicates. The condition groups must not be confounded with unmodeled factors that explain the expression difference. Dispersion must be estimable or stabilized appropriately for the available donor count. With a small number of donors, asymptotic Wald p-values can be poorly calibrated; the earlier material on pseudobulk inference with few donors and calibration of gene-level test statistics is relevant to that limitation.

Estimation and a Short Justification of the Standard Error

The mean model and dispersion are fitted for each gene. In a simple method-of-moments illustration, a group-specific dispersion estimate is \((s^2-\bar{Y})/\bar{Y}^2\), where \(\bar{Y}\) and \(s^2\) are the group sample mean and sample variance. A negative value is incompatible with a nonnegative NB dispersion and is truncated at zero in this simple estimator. Production analyses generally estimate dispersion using a likelihood-based procedure, often stabilizing gene-wise estimates across genes; the moment calculation is useful for understanding what dispersion represents.

A useful large-sample approximation for the variance of a log group mean under the NB model follows directly from the variance function. For \(n_g\) independent donors in group \(g\), the variance of the sample mean is approximately \((\mu_g+\phi\mu_g^2)/n_g\). Applying the delta method to the logarithm gives:

For equal exposures and independent groups, the approximate variance of the log mean ratio is the sum of the group variances, so its standard error is the square root of that sum. With unequal exposures, the GLM uses the offsets when fitting the coefficient and its covariance; the calculation should not be replaced by a raw count ratio. This approximation is an explanatory calculation, not a substitute for the fitted GLM covariance in a general design.

Defining and Reporting the Contrast

The contrast must match the scientific comparison and coding of the design. With reference coding, \(x_i=0\) for the reference and \(x_i=1\) for the second group, the contrast vector \((0,1)^\mathsf{T}\) selects the condition coefficient. Reversing the reference level changes the coefficient's sign and reciprocates its expression ratio, but it does not change the underlying comparison or a two-sided test.

A design may include prespecified nuisance covariates, such as a technical batch variable, when those covariates are estimable and not perfectly confounded with condition. In that case, \(\beta\) contains additional coefficients, but the condition comparison remains a specific linear contrast \(c^\mathsf{T}\beta\). Its variance must use the fitted coefficient covariance, including any covariance terms. The coefficient should not be interpreted in isolation if the contrast is not simply one coefficient.

Diagnostics, Alternatives, and Limitations

A defensible result requires checking whether the model's ingredients fit the data. Review donor-level sample identity, condition labels, cell-type definitions, and effective library sizes before interpreting coefficients. Examine fitted means and residual behavior across donors; a single influential donor can have a large effect when replication is limited. Check whether the observed mean-variance relationship is compatible with the fitted dispersions and whether there are extreme counts or excess zeros that the model does not capture. These checks complement, rather than replace, a valid design.

Warning signs include an estimated dispersion at a boundary for many genes, a few donors dominating the condition difference, systematic residual patterns, unstable estimates under small changes in normalization, or condition being completely aligned with batch. A model cannot recover a condition effect separately from a perfectly confounded batch effect. Likewise, offset adjustment cannot correct a misspecified cell type, poor sample annotation, or a lack of independent donors. When these failures occur, the correct response is to investigate the data and design, not simply to report a smaller p-value from another fitting option.

The paired pseudobulk approach is appropriate when the same donors are measured under both conditions; it uses within-donor comparisons and is not interchangeable with an independent-group model. A cell-level mixed model can represent cell-level observations and donor dependence directly, but has different assumptions and computational requirements. Treating all cells as independent is not a valid alternative. A Gaussian analysis of transformed or normalized expression may be useful under its own assumptions, but it does not retain the NB count mean-variance model and offset structure described here.

Finally, a small donor count limits precision and can make dispersion and Wald approximations unreliable even when there are many cells per donor. The effective replication for the condition comparison is the number of donors, not the number of cells. For genome-wide testing, per-gene p-values also require appropriate multiplicity control; a small unadjusted p-value alone is not a claim that a gene is differentially expressed with that probability.