Why Adjust for Cell State?
A peak may be accessible more often in one cell state, while its nearby gene is also expressed more highly in that state. A pooled association between accessibility and expression can therefore arise even if the two measurements are not associated among cells of the same state. Regression makes the comparison explicit: it estimates the accessibility–expression relationship while holding specified cell-state variables fixed.
The resulting coefficient is a conditional association. It describes how expression differs with accessibility among observations that agree on the included covariates, under the chosen model. Covariate adjustment does not, by itself, establish that accessibility causes expression. Unmeasured factors, measurement error, selection, and the direction of biological influence can all complicate a regulatory interpretation.
Covariates might include annotated cell type, a continuous cell-state score, cell-cycle score, donor, batch, or a library-size measure. Which are appropriate depends on the question and data-generating process. In particular, adjustment is not automatically beneficial: a covariate can be measured poorly, overlap strongly with accessibility, or lie on a biological pathway one wants to describe. State what the model conditions on and why.
Regression Models and Their Assumptions
For a continuous expression summary—for example, a suitably transformed expression value—a linear model is
Here, \(\beta_A\) is the expected difference in \(Y\) for a one-unit increase in accessibility, comparing cells with the same modeled covariates. Its units are expression units per accessibility unit. For a binary accessibility indicator, it is the adjusted difference in expected expression between accessible and inaccessible observations. The linear model assumes a correctly specified conditional mean; conventional standard errors additionally rely on independent errors with constant variance. Approximate normality of errors is relevant to small-sample \(t\)-based inference, not to the definition of the coefficient.
For nonnegative RNA counts, a count model can be more suitable. One common choice is negative-binomial regression with a log link and a library-size offset:
In this expression, \(L_i\) is a specified library-size factor and \(\theta\) is a dispersion parameter; one common parameterization has \(\operatorname{Var}(Y_i\mid A_i,Z_i)=\mu_i+\mu_i^2/\theta\). The offset adjusts the expected count for the chosen library-size factor. A one-unit increase in \(A_i\) is associated with a conditional expected-count ratio of \(\exp(\beta_A)\), for equal covariates and offset. The model requires an appropriate conditional mean and dispersion structure. If counts are especially sparse or accessibility is binary, the measurement and model choice should be made deliberately rather than treating zeros as ordinary continuous observations.
Cells from the same donor are generally not independent biological replicates. If multiple donors contribute cells, a model that accounts for donor clustering—such as a suitable mixed model or cluster-robust inference—may be needed. Treating thousands of cells from a few donors as thousands of independent replicates can produce standard errors that are too small. The effective replication for a biological claim is constrained by the independent biological units.
What Adjustment Does: A Short Derivation
The role of covariate adjustment in linear regression can be seen by residualizing. Let \(X\) be the design matrix containing an intercept and the covariates \(Z\), but not accessibility. Regress \(A\) on \(X\) and call the residuals \(r_A\); regress \(Y\) on \(X\) and call the residuals \(r_Y\). The coefficient of \(A\) in the full linear model is
To justify this expression, write the full model as \(Y=X\boldsymbol{\delta}+A\beta_A+\varepsilon\). For any fixed \(\beta_A\), least squares chooses \(\boldsymbol{\delta}\) to minimize the squared residuals after subtracting \(A\beta_A\). Projecting both \(Y\) and \(A\) onto the space orthogonal to the columns of \(X\) leaves residuals \(r_Y\) and \(r_A\). The remaining least-squares problem minimizes \(\sum_i(r_{Y,i}-\beta_A r_{A,i})^2\). Differentiating with respect to \(\beta_A\), setting the derivative to zero, and solving gives the displayed ratio. Thus the adjusted slope relates the parts of expression and accessibility not linearly explained by the included covariates. If all residualized accessibility values are zero, the slope is not identifiable.
Worked Example: An Adjusted Linear Association
Suppose six hypothetical cells have a continuous expression summary \(Y\), accessibility \(A\), and a binary state indicator \(S\). Consider the linear model \(Y_i=\beta_0+\beta_A A_i+\beta_S S_i+\varepsilon_i\).
| State \(S\) | Accessibility \(A\) | Expression \(Y\) |
|---|---|---|
| 0 | 0 | 2.0 |
| 0 | 1 | 2.5 |
| 0 | 2 | 3.0 |
| 1 | 0 | 5.0 |
| 1 | 1 | 5.5 |
| 1 | 2 | 6.0 |
Within each state, the mean accessibility is \(1\), and the mean expression is \(2.5\) in state 0 and \(5.5\) in state 1. After subtracting each state’s mean, the accessibility residuals are \((-1,0,1)\) in each group, and the expression residuals are \((-0.5,0,0.5)\) in each group. Therefore,
The adjusted estimate is \(0.5\) expression units per accessibility unit, conditional on state. The state coefficient is \(3.0\) expression units, since the two state means differ by \(3.0\) at the same accessibility. These data lie exactly on the fitted lines, leaving no residual variation from which to estimate ordinary sampling uncertainty. The coefficient calculation is valid, but this tiny constructed example does not support a useful uncertainty assessment.
Worked Example: A Pooled Association Can Reflect State
Consider six hypothetical cells. In state 0, accessibility is \(0,1,2\) and expression is \(1,1,1\). In state 1, accessibility is \(3,4,5\) and expression is \(5,5,5\). Within each state, expression does not change with accessibility. Yet the pooled slope is positive.
The pooled means are \(\bar{A}=2.5\) and \(\bar{Y}=3\). The unadjusted least-squares slope is
For example, the numerator is \(9\) from state 0 and \(9\) from state 1: the state-0 expression deviations are all \(-2\), while its accessibility deviations from the pooled mean are \(-2.5,-1.5,-0.5\); the corresponding products sum to \(9\). The products for state 1 also sum to \(9\). Once state is included, the within-state expression residuals are all zero, so the adjusted accessibility slope is \(0\). This example shows how a pooled association can reflect differences between states rather than a within-state relationship. Its exact values also leave no residual variation for inferential uncertainty estimates.
Worked Example: Interpreting a Count-Model Coefficient
Suppose a hypothetical negative-binomial model for a gene’s counts uses a log library-size offset, binary peak accessibility, cell state, and donor. Its fitted accessibility coefficient is \(\hat{\beta}_A=0.70\), with standard error \(SE=0.20\). Assume the model’s mean, dispersion, and dependence assumptions are adequate for this illustration, and that the reported standard error is appropriate for the donor structure.
For accessible versus inaccessible cells with the same modeled covariates and library-size offset, the estimated conditional expected-count ratio is
Thus, under this fitted model, accessible observations have an estimated expected count about \(2.01\) times that of inaccessible observations, conditional on the included variables. A Wald 95% interval on the coefficient scale is \(0.70\pm1.96(0.20)=[0.308,1.092]\). Exponentiating its endpoints gives an approximate 95% interval for the conditional expected-count ratio of \([\exp(0.308),\exp(1.092)]\approx[1.361,2.980]\). This interval concerns repeated-sampling uncertainty under the model; it is not a 95% probability that the fixed coefficient lies inside these particular endpoints. Neither the estimate nor its interval establishes a causal regulatory effect.
Diagnostics, Alternatives, and Limits
A fitted coefficient is informative only to the extent that the model and data support it. Useful checks include:
- Overlap and collinearity: inspect whether accessible and inaccessible cells occur across the covariate strata of interest. If a peak is accessible only in one state, its separate association from state may be weakly identified or impossible to estimate. Large standard errors, unstable coefficients, or strong sensitivity to covariate choices can signal this problem.
- Functional form: a linear accessibility term assumes a constant change in the model’s linear predictor per accessibility unit. Plot or otherwise examine residual patterns and consider whether a nonlinear term is scientifically and statistically justified. For a log-link count model, linearity is on the log-mean scale.
- Residual and variance behavior: for a linear model, examine residuals against fitted values and covariates for curvature, unequal spread, and influential observations. For a count model, assess whether the mean and dispersion assumptions are plausible. An excess of zeros or unmodeled heterogeneity can make a chosen count model inadequate.
- Replication and dependence: check how many independent donors or other biological units are represented, not only how many cells. Refit or reassess inference using a method that respects clustering when the design requires it.
- Measurement and selection: sparse accessibility calls and noisy expression summaries can attenuate or distort associations. If many peak–gene pairs are tested, account for multiple testing when making collections of claims; a small nominal p-value alone does not resolve model or measurement problems.
A simple correlation is useful as a descriptive summary but does not adjust for cell state. Stratifying by state allows separate within-state summaries, but can become unstable with many sparse strata and does not automatically combine information efficiently. Regression provides a compact conditional estimate and can include continuous covariates, but depends on its functional form and inference assumptions. A mixed model or cluster-robust approach addresses dependence structures that ordinary cell-level regression ignores; it does not cure confounding or poor measurement. Each alternative answers a somewhat different question, so report the model, covariates, scale, and biological units used.
A further limitation is interpretation. Cell state may be a confounder, a consequence of a regulatory process, or a proxy for multiple unmeasured features. Conditioning on it changes the association being described and does not necessarily recover a direct or total regulatory effect. A regression coefficient is evidence of a conditional pattern in the analyzed data, not proof that a particular peak regulates a particular gene.
Common Mistake
Correction: The coefficient estimates a conditional association under the specified model. Adjustment only controls the measured covariates included in that model; it does not rule out unmeasured confounding, reverse direction, or measurement artifacts. Describe the result as an adjusted accessibility–expression association unless a separate causal design and its assumptions justify stronger language.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- In a linear model with a binary accessibility indicator, what does \(\beta_A\) represent when cell type is included?
- Why might an unadjusted accessibility–expression slope be positive even when there is no within-state association?
- A log-link count model estimates \(\hat{\beta}_A=0.40\). What is the estimated conditional expected-count ratio for a one-unit accessibility increase?
- What is the identification problem if accessibility is completely determined by cell state?
- Why can treating every cell from a small number of donors as an independent replicate give misleading inference?
Solutions. (1) The adjusted expected expression difference per accessibility unit, comparing observations with the same modeled cell type and other covariates; for binary accessibility it is the accessible-versus-inaccessible difference. (2) State can shift both average accessibility and average expression, creating a between-state association that is not present within states. (3) \(\exp(0.40)\approx1.492\), or about 1.49 times the conditional expected count. (4) There is no accessibility variation within state to separate its coefficient from the state effect, so the accessibility association is not identifiable from those data. (5) Cells from one donor can share biological and technical influences; ignoring that dependence can understate uncertainty and overstate the amount of independent replication.