Design Matrices for Pseudobulk Differential Expression

Why the Design Matrix Matters

A pseudobulk analysis has a sample-level model: for a given cell type, each column of the data represents a donor-by-condition sample, and each row of the design matrix describes that sample's condition, donor, batch, and measured covariates. The design matrix is where the comparison becomes explicit. If the condition indicator is reversed, a reference level is misunderstood, or condition is inseparable from batch, the fitted model cannot answer the intended question merely by using a more sophisticated count distribution.

This tutorial concerns the structure and identifiability of that matrix. It does not specify a count distribution or estimation procedure. In a generalized linear model, the matrix enters the systematic component through \(X\beta_g\), where \(g\) indexes a gene; the distribution and link determine how that component relates to the observed counts. The algebraic checks here apply before selecting that distribution.

Definition: For \(n\) pseudobulk samples and \(p\) model terms, the design matrix \(X\) is an \(n \times p\) matrix. Row \(j\) contains the predictor values for sample \(j\), and column \(k\) represents one model coefficient. A model's systematic component for gene \(g\) is \(X\beta_g\), where \(\beta_g\) is the vector of gene-specific coefficients.

The rows must correspond to the actual samples analyzed, not to individual cells. If one donor contributes multiple condition samples, those are multiple rows, but they are not automatically independent observations. The design matrix records predictors; it does not by itself specify the covariance structure or account for dependence. Those modeling assumptions must be handled separately. Covariates must also be defined at the sample level and recorded for the samples to which their rows refer.

Reference Coding and the Meaning of Columns

For a categorical predictor with \(K\) levels, reference coding uses \(K-1\) indicator columns and selects one level as the reference. With an intercept, the reference level is represented by zeros in all of that factor's indicator columns. Each other level receives a 1 in its corresponding column. The intercept and these indicators together represent the factor without including a redundant indicator for every level.

For example, if control is the reference level for condition, a treatment indicator \(T_j\) is 0 for control samples and 1 for treated samples. In a model containing only an intercept and this indicator, the intercept corresponds to the reference condition's systematic component, and the treatment coefficient describes its difference from that reference on the model's coefficient scale. The precise interpretation of coefficients in a count model depends on its link; that is distinct from how the design matrix is encoded.

The same principle applies to donor and batch indicators. State the reference levels and column order explicitly, especially when interpreting software output. A different valid coding, such as sum-to-zero coding, changes the coefficient parameterization but need not change the comparisons the model can represent. Coding does not create information absent from the sample layout.

Worked Example: A Paired Design With Donor Fixed Effects

Suppose three hypothetical donors each contribute one control and one treated pseudobulk sample. Use control as the condition reference and donor 1 as the donor reference. The intended columns are an intercept, donor 2 indicator, donor 3 indicator, and treatment indicator \(T\). The six rows are:

X = \begin{pmatrix} 1&0&0&0\\ 1&0&0&1\\ 1&1&0&0\\ 1&1&0&1\\ 1&0&1&0\\ 1&0&1&1 \end{pmatrix}.

The row order is donor 1 control, donor 1 treated, donor 2 control, donor 2 treated, donor 3 control, donor 3 treated. The fourth column is 1 only for treated samples; the second and third columns identify donors 2 and 3 relative to donor 1.

The matrix has full column rank. To see this, suppose a linear combination of its columns is zero. Subtracting the control row from the treated row within donor 1 forces the coefficient of \(T\) to be zero. The donor 1 control row then forces the intercept coefficient to be zero. The donor 2 and donor 3 control rows force the coefficients of the donor indicators to be zero. Thus only the zero coefficient vector maps to zero. The paired layout therefore separates a common treatment term from donor indicators in this model. It does not, by itself, assert that treatment responses are identical across donors or specify how within-donor dependence is modeled.

Rank, Confounding, and Estimability

The columns of \(X\) are linearly independent when no nonzero vector \(v\) satisfies \(Xv=0\). In that case, \(X\) has full column rank, and its coefficient vector is uniquely determined by the systematic component. If a column is a linear combination of other columns, the matrix is rank deficient: different coefficient vectors can give exactly the same fitted systematic component.

Definition: A coefficient or linear combination \(c^\mathsf{T}\beta_g\) is estimable if it has the same value for every coefficient vector producing the same systematic component \(X\beta_g\). Equivalently, \(c\) must lie in the row space of \(X\), or be orthogonal to every vector \(v\) in the null space of \(X\).

The equivalence follows directly. If \(Xv=0\), then \(X(\beta_g+v)=X\beta_g\), so the data's systematic component cannot distinguish \(\beta_g\) from \(\beta_g+v\). The target remains unchanged across those coefficient vectors exactly when \(c^\mathsf{T}v=0\) for every such \(v\). The vectors orthogonal to the null space of \(X\) form its row space. Thus a rank-deficient matrix may still identify some combinations of coefficients, even though it cannot identify all coefficients separately.

Worked Example: Condition Completely Confounded With Batch

Consider four hypothetical samples: two controls, both processed in batch A, and two treated samples, both processed in batch B. With an intercept, treatment indicator \(T\), and batch B indicator \(B\), the design matrix is:

X = \begin{pmatrix} 1&0&0\\ 1&0&0\\ 1&1&1\\ 1&1&1 \end{pmatrix}.

Here the treatment and batch columns are identical: \(T=B\). Consequently, \(X(0,1,-1)^\mathsf{T}=0\). Increasing the treatment coefficient by any amount and decreasing the batch coefficient by the same amount leaves \(X\beta_g\) unchanged. The treatment coefficient and batch coefficient are not separately estimable; only their sum is identified by this design.

Adding more samples with the same pattern does not resolve the confounding. Samples that break the alignment—for example, controls and treated samples represented in both batches—are needed to separate condition from batch. If the study did not collect such samples, no choice of reference level or analysis software can recover the separate effects from these data alone.

A Direct Rank Check With Several Predictors

In small examples, full rank can be checked by identifying a square submatrix with a nonzero determinant. For larger matrices, numerical rank-revealing methods, such as a QR decomposition with pivoting, are more practical. The determinant itself is not a stable general-purpose diagnostic for large matrices; the useful principle is whether the columns contain redundant information.

Worked Example: Verifying Full Rank for Five Columns

For five hypothetical samples, use an intercept, treatment indicator, donor B indicator, batch B indicator, and a binary covariate \(Z\). The rows below follow those column definitions:

X = \begin{pmatrix} 1&0&0&0&0\\ 1&1&0&0&0\\ 1&0&1&1&1\\ 1&0&1&0&1\\ 1&0&1&1&0 \end{pmatrix}.

To check the determinant of this square matrix, subtract the first row from each of rows 2 through 5. Expansion along the first column then leaves the determinant of:

A = \begin{pmatrix} 1&0&0&0\\ 0&1&1&1\\ 0&1&0&1\\ 0&1&1&0 \end{pmatrix}.

Expanding \(A\) along its first column leaves the \(3 \times 3\) determinant:

\det \begin{pmatrix} 1&1&1\\ 1&0&1\\ 1&1&0 \end{pmatrix} = 1(0-1)-1(0-1)+1(1-0) =-1+1+1=1.

The displayed matrix is the appropriate minor. Its determinant is 1; alternatively, subtract its first row from each of the other rows to obtain a triangular matrix with determinant 1. Therefore \(\det(X)=1\), which is nonzero, and the five columns are linearly independent. This example also illustrates why each determinant step should be checked rather than inferred from a visual pattern.

Design Checks Before Fitting a Model

A useful review starts with the sample table, not just a software-generated matrix. Confirm that each row maps to the intended donor, condition, batch, and covariate values. Then check the factor levels and columns: a treatment indicator should mark the intended condition, and a donor indicator should use the stated reference. Empty or unexpected categories, mislabeled factors, and accidental numeric coding can produce a matrix that is technically computable but does not represent the scientific comparison.

Check redundant columns. With an intercept, including an indicator for every level of a factor creates an exact linear dependence. One level must be omitted or another nonredundant coding used.
Check confounding. Compare condition against batch and other categorical predictors in a cross-tabulation. A condition represented in only one batch may be inseparable from batch.
Check sample-level covariates. A donor-level covariate is constant within a donor. If donor fixed indicators are also included, that covariate may be exactly redundant with them or their linear combinations.
Check numerical dependence. Continuous covariates on very different scales can make numerical calculations unstable even without exact rank deficiency. Centering or rescaling can help computation; it does not repair confounding or add information.

Full rank is a necessary algebraic condition for estimating every coefficient in the specified parameterization, not a guarantee of precise estimates. A nearly redundant column can lead to large uncertainty even when rank is technically full. Likewise, a design with as many samples as columns can be full rank yet leave no residual degrees of freedom for some modeling tasks. Adequate biological replication and the information limits of donor-level covariates remain essential.

How This Check Differs From Other Modeling Choices

Reference coding and sum-to-zero coding are alternative parameterizations of categorical predictors. They can give different coefficient values and coefficient interpretations while representing the same set of fitted systematic components, provided each coding is full rank and the same model is intended. Choosing a reference level is therefore a matter of transparent interpretation, not a method for resolving a missing comparison.

Similarly, treating donor as a fixed effect or representing donor variation with a random effect makes different modeling assumptions. The design matrix for fixed donor indicators explicitly includes donor columns. A random-effects formulation has additional distributional and covariance structure; its design matrices and assumptions are not interchangeable with fixed-effect coding. Neither approach can create condition-versus-batch information when the observed sample layout perfectly confounds those factors.

A practical failure signal is a software message about aliased, non-estimable, or dropped coefficients. It should prompt inspection of the sample layout and design columns, not an automatic acceptance of the remaining output. A coefficient silently dropped by a fitting routine does not become estimable; the reported parameterization or target may have changed. Before proceeding, state which comparison is identifiable and whether it matches the scientific question.

Key takeaway: A design matrix is an explicit map from pseudobulk samples to model terms. Verify its row identities, reference levels, rank, and estimable targets before fitting; otherwise the model may be unable to distinguish the biological condition from donor, batch, or covariate structure.

Check Your Understanding

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

With an intercept and a two-level condition factor, why is it redundant to include indicators for both conditions?
If treatment and batch B indicators are identical in every sample, which quantities can the model distinguish?
What does full column rank guarantee, and what does it not guarantee?
Why might a donor-level covariate be redundant when donor fixed indicators are in the same model?

Solutions. (1) The two indicators sum to the intercept column, so the columns are linearly dependent; use one reference level or another full-rank coding. (2) The model can identify their combined contribution, but not the treatment and batch coefficients separately. (3) Full column rank means every coefficient in that parameterization is algebraically identifiable from the systematic component. It does not guarantee adequate precision, correct model assumptions, or sufficient biological replication. (4) The covariate is constant within each donor and can therefore be represented by the donor indicators, or by a linear combination of them, leaving no independent column information for its separate coefficient.