Why Pair Pseudobulk Samples?
In a paired experiment, each donor contributes a pseudobulk sample in both conditions for the same cell type. The two observations from a donor are related: they share biological features that can affect expression, even when the condition changes. A paired analysis uses that structure to estimate a condition contrast within donors, rather than comparing two groups as if their observations came from unrelated people.
The unit of replication is still the donor, not the individual cell. For a particular gene and cell type, the data consist of donor-condition pseudobulk counts. Pairing does not turn the two samples from each donor into independent biological replicates, nor does it increase the number of donors. The introductory material on pseudobulk estimands establishes why these donor-level observations are the basis for inference.
Organizing the Matched Data
Index donors by \(d\), conditions by \(c\), and genes by \(g\). Let \(Y_{dgc}\) be the summed count for gene \(g\) in donor \(d\), condition \(c\), and the cell type being analyzed. Let \(L_{dc}\) be the library-size measure used as the offset. A paired design has both condition records for each included donor, with a reliable donor identifier linking them.
The table below is a small hypothetical example for one gene and one cell type. Library sizes are shown in millions of counts for readability; in an analysis, the offset uses the corresponding positive library-size value on the model's count scale.
| Donor | Condition | Gene count | Library size (millions) |
|---|---|---|---|
| D1 | Control | 120 | 1.0 |
| D1 | Treated | 180 | 1.5 |
| D2 | Control | 80 | 0.8 |
| D2 | Treated | 100 | 1.0 |
| D3 | Control | 150 | 1.2 |
| D3 | Treated | 210 | 1.4 |
Each donor contributes one row per condition. Raw counts should not be compared without regard to library size: for example, D1's treated count is higher, but its library is also larger. The library-size offset adjusts the expected count for this exposure difference. The pairing variable is the donor ID; a condition label alone cannot establish which observations belong together.
Paired Negative-Binomial Model and Its Estimand
For a given gene and cell type, a standard count model assumes conditionally independent pseudobulk counts with negative-binomial variation. A common parameterization has mean \(\mu_{dgc}\) and variance \(\mu_{dgc}+\phi_g\mu_{dgc}^2\), where \(\phi_g\) is the gene's dispersion. With control coded \(x_{d,\mathrm{control}}=0\) and treated coded \(x_{d,\mathrm{treated}}=1\), a paired log-link model is
Here \(\log(L_{dc})\) is an offset with coefficient fixed at one, \(\alpha_{dg}\) represents the donor-specific baseline for this gene, and \(\beta_g\) is the common treated-versus-control effect on the log scale. The model's donor terms encode the pairing. The later treatment of donor fixed effects develops their estimation and design-matrix representation; for the present purpose, the key point is what cancels in the within-donor contrast.
The model states that the expected count per library-size unit is \(\exp(\alpha_{dg})\) in control and \(\exp(\alpha_{dg}+\beta_g)\) in treated. Dividing expected counts by their library sizes removes the offset contribution. The within-donor ratio of these expected normalized counts is therefore
Taking logarithms gives \(\beta_g\). Thus the coefficient is a within-donor log rate ratio for treated versus control, conditional on the model. It is not a raw difference in counts, a percentage of cells expressing the gene, or evidence that treatment caused the change. If donor-specific condition effects vary, one common coefficient summarizes them under the model's fitting and weighting; it need not equal every donor's effect.
Worked Example: Reading the Paired Design
Suppose the hypothetical table above represents one cell type, and the scientific contrast is treated versus control. Identify the paired units, explain the role of the offset, and state the condition coefficient's interpretation.
Pairing. D1's control observation is paired with D1's treated observation, and likewise for D2 and D3. There are three donor pairs, not six independent donors. The six rows are sample records, while the three donors provide the independent biological replication for the paired comparison.
Offset. The D1 treated library has 1.5 million counts compared with 1.0 million in control. The model uses \(\log(L_{dc})\) as an offset so this difference in library size is not mistaken for a condition effect. In particular, the raw gene-count ratio \(180/120=1.5\) is not by itself the normalized expression ratio; the library sizes must also be accounted for.
Effect. With control coded zero and treated coded one, \(\beta_g\) is the treated-versus-control log rate ratio within donors. Its exponent, \(\exp(\beta_g)\), is the corresponding rate ratio. Reversing the condition coding changes the sign of the log contrast and reciprocates the rate ratio.
Estimation, Inference, and Interpretation
After fitting the paired count model, the estimated condition coefficient \(\hat{\beta}_g\) quantifies the contrast; its standard error measures uncertainty under the fitted model. A Wald statistic is the estimate divided by its standard error, and a two-sided normal-approximation p-value compares the absolute statistic with a standard normal reference distribution under the null hypothesis \(H_0:\beta_g=0\). This tests no condition-associated log rate difference, conditional on the model. The validity of the approximation depends on adequate model fitting and an appropriate inferential calibration, especially the number of independent donors.
Worked Example: From a Paired Coefficient to a Fold Change
Consider a hypothetical paired pseudobulk fit for a gene in one cell type. The fitted treated-versus-control coefficient is \(\hat{\beta}_g=0.4055\), with model-based \(SE=0.1800\). Assume for this illustration that the large-sample Wald normal approximation is an appropriate calibration.
Fold change. The estimated normalized expression ratio is
The fitted model therefore estimates about 1.50 times the normalized expression rate in treated relative to control, or approximately 50.0% higher. This is the estimated within-donor contrast under the model, not a claim that every donor has exactly that ratio.
Wald test. The statistic is
The two-sided normal-approximation p-value is \(2P(Z\geq |2.2528|)\approx 0.02427\), where \(Z\sim N(0,1)\). At \(\alpha=0.05\), this is evidence against \(H_0:\beta_g=0\) under the stated approximation. The p-value is not the probability that the null is true.
Interval. A Wald 95% confidence interval on the log scale is
Exponentiating its endpoints gives an approximate 95% interval of \([\exp(0.0527),\exp(0.7583)]\approx[1.054,\;2.135]\) for the treated-to-control rate ratio. Under repeated use of the same valid interval procedure, 95% of such intervals would cover the true model parameter. This example illustrates the arithmetic; a small donor count or poor model fit can make a normal approximation unreliable.
Checks, Alternatives, and Failure Modes
Before interpreting a paired result, verify that every sample is assigned to the right donor and condition, and that the same cell type and comparable processing rules were used for both conditions. Examine library sizes, count distributions, fitted values, and influential donors. A striking coefficient driven by one donor, strong model misfit, or an effect that changes substantially when one pair is removed deserves investigation rather than automatic acceptance. These checks diagnose potential problems; they do not prove the model assumptions.
Pairing is not possible merely because the data contain the same number of samples in two conditions. It requires genuine matching: both observations must come from the same donor. If control and treated samples come from different donors, labeling rows as pairs creates a false design. Also check whether condition is entangled with another factor: if all control samples are processed in one batch and all treated samples in another, the data may not distinguish a condition effect from a batch effect.
The main alternative is an unpaired model that treats the two conditions as independent donor groups. When measurements truly are matched, that analysis discards within-donor information and may fail to represent the sampling structure. A simple paired-difference analysis is another option when a suitable continuous, approximately normal donor-level measurement is available: compute one treated-minus-control difference per donor and analyze those differences. For overdispersed pseudobulk counts, a negative-binomial model with a library-size offset retains the count structure and estimates a log rate ratio directly. The paired-difference approach is not interchangeable without checking its measurement scale and assumptions.
Pairing has limits. It controls donor-level baseline differences represented by the model, but it cannot remove time-varying confounding, measurement errors in the matching, or condition effects that are inseparable from batch. It also cannot recover a within-donor contrast for a donor missing one condition without additional modeling assumptions. If effects vary greatly across donors, a single common coefficient may be an incomplete summary. Finally, many genes are tested in a typical analysis, so a gene-level p-value should not be mistaken for a conclusion that accounts for the full multiple-testing problem.
Common Mistake
A second interpretive error is to call \(\beta_g\) a difference in raw counts. The coefficient is on the log scale and compares expected counts after accounting for library size. Report the contrast direction, the gene and cell type, and the rate ratio or log rate ratio with its units understood as a relative expression scale rather than an absolute count.
Check that donor IDs link the correct two condition samples and identify which donors have complete pairs.
Choose a reference condition and coding so the sign of the condition coefficient has an unambiguous direction.
Use the appropriate pseudobulk count and offset, represent donor pairing, and examine fit and influential observations.
Interpret the coefficient as a log rate ratio, exponentiate it for a rate ratio, and report uncertainty and the inferential method.
Check Your Understanding
Answer the questions, then compare your reasoning with the solutions.
- Eight donors each contribute a control and treated pseudobulk sample. How many donor pairs are represented, and why are there not 16 independent biological replicates?
- In a model with control coded zero and treated coded one, what does \(\exp(\beta_g)\) represent?
- A fitted coefficient for control versus treated is \(-0.30\). What is the estimated control-to-treated rate ratio, approximately?
- All control samples were processed in batch A and all treated samples in batch B. What concern does this raise for interpreting the condition coefficient?
Solutions. (1) There are eight pairs. The two observations from each donor share donor-level features and are linked; they are not independent donors. (2) It is the treated-to-control ratio of expected normalized expression rates within donors, under the common-effect model. (3) \(\exp(-0.30)\approx0.741\), so the estimated control rate is about 0.741 times the treated rate. (4) Condition and batch are perfectly confounded in this design, so their separate effects cannot be identified from these observations alone.