Worked Paired Pseudobulk Analysis

From Cells to Paired Donor Comparisons

A paired pseudobulk analysis asks whether expression in a specified cell type differs between two conditions while accounting for the fact that both conditions were measured in the same donors. The essential unit of replication is the donor, not the individual cell. Aggregating counts within each donor-condition combination preserves that replication structure and avoids treating cells from one donor as independent biological replicates.

The workflow has four linked decisions: define the cell type and paired samples, sum raw counts within each sample, encode donor and condition in the model, and test a clearly defined gene-level contrast. The worked numbers below are hypothetical. They illustrate the calculations; they do not represent a real study.

Aggregation is performed separately for each gene and donor-condition sample. If a sample contains \(m_{dc}\) cells, then \(Y_{gdc}=\sum_{j=1}^{m_{dc}} y_{gdcj}\), where \(y_{gdcj}\) is the raw count for gene \(g\) in cell \(j\). The same cell assignments must be used consistently across genes. Counts are summed, not averaged; the resulting library size or a normalized size factor is accounted for in the count model.

Model and Paired Design Matrix

For each gene, a standard pseudobulk count model is a negative binomial generalized linear model. Write \(Y_{gdc}\) for the observed count and \(\mu_{gdc}\) for its expected count. One common parameterization has variance \(\mu_{gdc}+\phi_g\mu_{gdc}^2\), where \(\phi_g\) is the gene-specific dispersion. A log link relates the expected count to a sample-specific size factor and the design:

Here \(s_{dc}\) is a size factor, \(\alpha_{gd}\) is the gene-specific donor effect, \(T_c=0\) for control and \(T_c=1\) for treatment, and \(\beta_g\) is the common treatment-versus-control log mean contrast for gene \(g\), conditional on the size factors and donor effects. A positive \(\beta_g\) indicates a higher expected normalized count under treatment. The corresponding fold change is \(\exp(\beta_g)\), and the log base 2 fold change is \(\beta_g/\log(2)\).

This model assumes that the negative binomial mean-variance relationship is adequate, the size factors reasonably account for library-depth differences, and the modeled donor and condition effects describe the design. The paired fixed-effect version assumes a common condition effect across donors; it does not estimate a separate treatment effect for every donor. It also requires the condition contrast to be identifiable: paired samples must exist across donors, and condition must not be inseparable from another design feature such as batch.

For four complete pairs, one equivalent design-matrix parameterization uses an intercept, indicators for donors D2 through D4, and a treatment indicator. D1 is the reference donor. Rows are ordered by donor, with control then treatment within each donor:

With this coding, the coefficient of \(T\) is the treatment-versus-control contrast. In matrix notation, if the fitted coefficient vector is \(\hat{\boldsymbol\beta}_g\), the contrast is \(c^\mathsf{T}\hat{\boldsymbol\beta}_g\), where \(c=(0,0,0,0,1)^\mathsf{T}\). The null hypothesis for no condition difference is \(H_0:\beta_g=0\); a two-sided alternative is \(H_a:\beta_g\ne 0\).

Gene-Wise Contrasts and Multiplicity

A fitted model gives an estimate and uncertainty for the condition coefficient for each gene. For a Wald test, the statistic is the estimated coefficient divided by its standard error. The two-sided p-value is calculated under the null reference distribution used by the analysis. A Wald statistic using a standard normal reference is an asymptotic approximation; with few donors, its calibration and uncertainty deserve particular scrutiny. The choice of test and any small-sample calibration should be specified rather than assumed to be interchangeable.

Checks, Alternatives, and Limits

Before interpreting results, verify that each aggregate corresponds to the intended donor, condition, and cell type; that each donor’s pair is matched correctly; and that raw counts, not transformed expression values, were summed. Inspect the design matrix for duplicated or missing samples, aliased columns, and condition-batch confounding. Check whether size factors and the gene-wise mean-variance fit are plausible, and review residual, influence, and outlier diagnostics. These checks matter because a correctly calculated contrast from a poorly specified or poorly fitting model is still unreliable.

A paired model is appropriate when the same biological donors contribute both conditions and the scientific target is a condition contrast adjusted for donor baseline. An independent-group model instead treats the two condition groups as separate sets of donors; it does not use within-donor pairing and is not a replacement for a paired design when matching is present. Cell-level mixed models offer another way to represent donor clustering, but they have different assumptions and computational demands. Neither approach justifies counting cells as independent biological replicates.

The fixed donor-effect model estimates a shared condition effect. If treatment responses vary substantially among donors, that common-effect summary may not describe the variation of individual responses; the paired analysis here does not resolve that question. Incomplete pairs, very few donors, confounding, and influential samples can also limit interpretation. Statistical significance should be considered alongside effect size, uncertainty, model diagnostics, and the biological context.