From Analysis Decisions to an Auditable Report
A pseudobulk result is useful only if a reader can tell what biological samples entered the analysis, how cell-level counts became sample-level counts, which model was fitted, and what comparison its coefficients represent. Reporting these steps is not administrative detail: changing the aggregation unit, filtering rule, reference level, or contrast can change the analysis and its interpretation.
The goal here is a compact, reproducible account of one gene-level differential-expression analysis. The examples are hypothetical. They use a paired design and a negative-binomial count model with a donor random intercept; the reporting principles also apply when a different design or model is justified.
Report the Analysis Set Before the Model
Start by naming the biological unit, the grouping used for aggregation, and the inclusion rules. A reader should be able to reconstruct the mapping from cells to pseudobulk libraries: for example, aggregate raw gene counts separately for each donor, condition, and cell type. State how cell identities and sample labels were assigned, and whether any cells or samples were excluded before aggregation.
Filtering should be stated as an operational rule, not merely described as “low-expression genes were removed.” Give the threshold, the level at which it was applied, and the number of genes before and after filtering. Report the number of biological samples in each condition, the number of donors, and the number of pseudobulk observations actually analyzed. If the rule requires a gene to pass a threshold in a specified number of libraries, make clear that the libraries are pseudobulk samples, not individual cells.
Worked Example: Documenting Aggregation and Filtering
Suppose a hypothetical paired experiment has 8 donors, each measured in control and treated conditions, for one cell type. The cell-level matrix contains raw UMI counts. After applying a pre-specified cell-quality rule, each donor-condition group has at least 20 cells. For gene A in donor 1, control-cell counts are 2, 0, 3, and 1; treated-cell counts are 4, 1, 2, and 3.
Aggregation. Sum counts within each donor-condition-cell-type group. The two gene A pseudobulk counts for donor 1 are \(2+0+3+1=6\) in control and \(4+1+2+3=10\) in treated. The resulting analysis has 16 pseudobulk libraries, not 16 independent donors and not one independent observation per cell. There are 8 biological replicates with paired measurements.
Filtering. Suppose the pre-specified rule retains a gene if its raw pseudobulk count is at least 10 in at least 6 of the 16 libraries. Gene B has counts of at least 10 in 7 libraries, so it passes; gene C reaches 10 in only 5 libraries, so it does not. The report should state the rule and these outcomes, as well as the total number of genes tested after filtering. This rule is an example, not a universal threshold.
A reproducible report would also give the cell type definition, aggregation software or code version, and a sample-level table or equivalent mapping. It would state whether a pseudobulk library with too few cells was excluded and whether that exclusion left any donor without a usable paired observation.
State the Model and the Estimand
A model report must identify the response, distribution, link, offset, fixed effects, random effects, and reference levels. For a paired count analysis, one possible gene-specific model is
Here \(Y_{gdj}\) is the summed count for gene \(g\), donor \(d\), and condition-specific sample \(j\); \(s_{dj}\) is the effective library-size factor used as an offset; and \(I(\text{treated})\) equals 1 for treated and 0 for control. The model uses the convention \(\operatorname{Var}(Y\mid b)=\mu+\phi_g\mu^2\). The donor-specific random intercept satisfies \(b_{gd}\sim N(0,\tau_g^2)\). Conditional on the random effects and model covariates, the count observations are assumed independent.
This model assumes that the chosen offset adequately accounts for library-size differences, that the negative-binomial mean-variance form is adequate, and that donor random intercepts are appropriately represented by the specified distribution. It also requires an identifiable design: the condition effect cannot be separated from a perfectly confounded covariate. State any additional fixed effects, such as batch, only when they are part of the planned and estimable design. Describe the normalization and fitting procedure sufficiently for another analyst to reproduce it.
With control as the reference and treatment coded as 1, \(\beta_{1g}\) is the treated-versus-control log fold-change on the natural-log scale, conditional on the model’s other terms. Its exponential, \(\exp(\beta_{1g})\), is the fitted treated-to-control fold change. State the scale explicitly; a natural-log coefficient is not itself a log base 2 fold-change.
Worked Example: Reporting a Contrast and Its Uncertainty
Suppose a hypothetical fitted model for gene D estimates \(\hat{\beta}_{1D}=0.405\) with \(SE=0.120\), using the treatment-minus-control contrast and a large-sample Wald normal approximation. The contrast is the coefficient of the treated indicator, so the model-based log fold-change estimate is \(0.405\).
The Wald statistic for \(H_0:\beta_{1D}=0\) is
The two-sided normal-reference p-value is approximately \(2P(Z\geq 3.375)=0.00074\), where \(Z\sim N(0,1)\). A 95% Wald interval on the natural-log scale is
Exponentiating gives a fold-change estimate of \(\exp(0.405)=1.499\), with interval \([\exp(0.1698),\exp(0.6402)]=[1.185,\ 1.897]\). Thus the fitted expression is estimated to be 1.50 times as high in treated as in control, with a model-based 95% interval from 1.185 to 1.897 for that ratio. This is an estimate of association under the specified model, not evidence by itself that treatment caused the change. The report must identify the inference method; a Wald interval is not interchangeable with an interval produced by a different small-sample procedure.
For a more general contrast, write the target as \(c^\mathsf{T}\beta_g\), where \(c\) is a vector specifying the coefficient combination. Its estimate is \(c^\mathsf{T}\hat{\beta}_g\), and its model-based standard error is \(\sqrt{c^\mathsf{T}\widehat{\operatorname{Cov}}(\hat{\beta}_g)c}\). This makes the reported effect traceable to the model coefficients. For interactions or multiple conditions, state the actual coefficient combination and the condition levels to which it applies; a coefficient label alone may not describe the scientific comparison.
Diagnostics, Adjusted Findings, and a Reporting Checklist
A report should summarize checks relevant to the fitted model rather than claim that diagnostics “passed” without evidence. Depending on the analysis, report residual patterns, dispersion and mean-variance checks, fitted means and excess-zero checks, influential donor samples, outlier-count handling, convergence, singularity, and calibration of test statistics. These checks should be interpreted in relation to the data and the model. A concern may motivate sensitivity analysis, a revised model, or a qualified conclusion; it should not disappear from the report because it complicates the result.
For gene-level inference across many tested genes, report the adjustment method and the testing family. Under the Benjamini–Hochberg (BH) procedure, order the \(m\) raw p-values as \(p_{(1)}\leq \cdots \leq p_{(m)}\). The adjusted value for rank \(i\) is the smaller of 1 and the minimum of \(mp_{(j)}/j\) over ranks \(j\geq i\). State the chosen decision threshold and whether the rule is “at or below” or strictly “below” that threshold.
Worked Example: Applying BH and Reporting the Cutoff
Suppose five hypothetical genes have ordered raw p-values \(0.001,\ 0.010,\ 0.025,\ 0.040,\ 0.200\). For \(m=5\), the rank-specific BH quantities are \(5p_{(i)}/i\): \(0.005,\ 0.025,\ 0.04167,\ 0.05000,\ 0.200\). These are already nondecreasing, so the cumulative minimum-from-the-largest-rank adjustment leaves them unchanged.
| Gene | Raw p-value | BH-adjusted p-value | At or below 0.05? |
|---|---|---|---|
| A | 0.001 | 0.00500 | Yes |
| B | 0.010 | 0.02500 | Yes |
| C | 0.025 | 0.04167 | Yes |
| D | 0.040 | 0.05000 | Yes |
| E | 0.200 | 0.2000 | No |
At the stated adjusted-p-value rule of \(\leq 0.05\), genes A through D meet the threshold, including D, whose adjusted value is exactly 0.05000. If the rule were strictly \(<0.05\), D would not meet it. The report should give the number of genes tested, the adjustment method, the threshold, and both effect estimates and adjusted p-values; an adjusted p-value alone does not communicate the size or direction of an effect.
A concise report can present, for each reported gene, the contrast and direction, log fold-change scale, standard error or interval, raw p-value, and adjusted p-value. Include the filtering and testing counts, model formula and reference levels, normalization and software details, and a summary of diagnostics or sensitivity analyses. The earlier material on residual diagnostics, influence diagnostics, and calibration of gene-level test statistics gives the specialized checks; the report’s task is to record what was checked and what consequences followed.
Different approaches answer related but not identical questions. A paired model accounts for donor-level pairing; treating the two conditions as independent discards that structure. A donor fixed-effect model can control for donor-specific baselines in a paired design, while a random-intercept model represents donor variation through a distribution and can support estimation across donors under its assumptions. A simpler independent-group model is appropriate only when the biological samples are independent and the design supports that assumption. In every case, the report should explain why the chosen structure matches the sampling design rather than selecting a model solely because software fits it.
Pseudobulk reporting also has limits. Aggregation discards cell-level variation within each sample, and a model cannot repair a design with too few independent donors, confounded conditions, or poorly measured covariates. Results can depend on the cell-type definition, aggregation, filtering, normalization, and model assumptions. A transparent report makes those dependencies visible; it does not turn observational comparisons into causal claims or imply that a non-significant result proves no effect.
State the sample and cell-type definitions, cell exclusions, aggregation rule, numbers of donors and pseudobulk libraries, and gene-filtering rule.
Give the response distribution, offset, formula, reference levels, random-effects structure, assumptions, and exact coefficient combination that represents the target comparison.
Summarize diagnostics and any sensitivity analysis, then report effect estimates with uncertainty, raw p-values, adjustment method, testing family, threshold, and adjusted findings.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- A study has 10 donors, each with control and treated samples. How many independent biological replicates are there, and how many pseudobulk libraries are there if every sample is retained?
- A coefficient is \(0.693\) on the natural-log scale. What is the estimated fold change, to three significant figures?
- Five genes have ordered raw p-values \(0.002,\ 0.012,\ 0.030,\ 0.060,\ 0.400\). Calculate the BH-adjusted p-value for the second-ranked gene, assuming the adjusted sequence is obtained by the usual cumulative minimum rule.
- Why is “adjusted p-value = 0.05; significant if below 0.05” an insufficient basis for calling that result significant?
Solutions. (1) There are 10 independent biological replicates (donors) and 20 pseudobulk libraries, with two repeated condition measurements per donor. (2) \(\exp(0.693)\approx 2.00\), approximately a 2.00-fold increase. (3) At rank 2, \(mp_{(2)}/2=5(0.012)/2=0.030\). The rank-3, rank-4, and rank-5 quantities are \(0.050,\ 0.075,\ 0.400\), so the cumulative minimum from rank 2 onward is \(0.030\). (4) “Below” means strictly less than 0.05, so a value equal to 0.05 does not qualify; a rule using “at or below 0.05” would include it.