QTLs Across Continuous Cell States

Why Model Genetic Effects Across a Continuous State?

Cells of the same broad type can occupy a range of activation, differentiation, or other molecular states. A genetic variant may have a different association with gene expression at different points along that range. Dividing cells into a few categories can conceal a gradual change, make results depend on arbitrary cut points, or leave some categories with little data. A continuous-state model instead asks whether the expression association changes with a measured state score.

The key is to distinguish two questions. A state-dependent expression effect asks whether genotype is associated with different expression levels at different values of the state score. A state-abundance effect asks whether genotype is associated with having more cells at some state values than others. These are different outcomes and different estimands: one concerns expression conditional on state, while the other concerns the distribution of states.

A Genotype-by-State Regression Model

Let \(i\) index donors and \(j\) index cells from donor \(i\). Let \(Y_{ij}\) be the expression outcome for a gene in cell \(j\), \(G_i\) the donor's genotype dosage coded as 0, 1, or 2 copies of the effect allele, and \(S_{ij}\) the cell's continuous state score. For interpretation, the score is often centered so that \(S=0\) represents a meaningful reference, such as the mean score in the analyzed cells. It may also be scaled, for example to have standard deviation 1.

Definition: A genotype-by-state interaction model includes the product of genotype dosage and the continuous state score. In a linear model with donor-level covariates \(C_i\), one form is

Y_{ij}=\beta_0+\beta_G G_i+\beta_S S_{ij}+\beta_{GS}(G_iS_{ij})+C_i^\top\gamma+b_i+\epsilon_{ij}.

Here, \(\beta_G\) is the genotype-associated expression difference per additional effect allele when \(S=0\); \(\beta_S\) is the state-expression slope at \(G=0\); and \(\beta_{GS}\) measures how much the genotype slope changes for a one-unit increase in state score. The term \(b_i\) is a donor-level random intercept, and \(\epsilon_{ij}\) is the remaining cell-level error. A random intercept or another appropriate donor-level method accounts for dependence among cells from the same donor; cells are not independent genetic replicates.

This linear specification assumes that the conditional mean expression changes linearly with state and that the genotype effect changes linearly with state over the observed range. For a continuous expression outcome, a common working model also assumes suitably behaved residual variation for the chosen inferential procedure. Covariates should be included when justified by the study design and analysis plan. If expression is modeled on a log scale, coefficients are differences on that scale; they are not raw-unit expression changes.

Key result: At a fixed state value \(s\), the model-implied expression difference for one additional effect allele is

\Delta_G(s)=\beta_G+\beta_{GS}s.

Thus, \(\beta_G\) alone is not the genetic effect at every state. It is the effect at the reference score \(S=0\). The interaction coefficient \(\beta_{GS}\) determines whether and how that effect changes along the score.

Worked Example: Calculating Genetic Effects at Two State Values

Suppose a hypothetical study fits the model above to a log-expression outcome, using a centered, standardized state score. Its estimated genotype coefficient is \(\hat{\beta}_G=-0.20\) log-expression units per allele, and its estimated interaction coefficient is \(\hat{\beta}_{GS}=0.12\) log-expression units per allele per one-unit increase in state score. Calculate the estimated per-allele expression association at \(S=-1\), \(S=0\), and \(S=1\).

Use the conditional slope \(\widehat{\Delta}_G(s)=\hat{\beta}_G+\hat{\beta}_{GS}s\). At \(S=-1\),

\widehat{\Delta}_G(-1)=-0.20+(0.12)(-1)=-0.32.

At the reference state \(S=0\), the slope is \(-0.20+(0.12)(0)=-0.20\). At \(S=1\),

\widehat{\Delta}_G(1)=-0.20+(0.12)(1)=-0.08.

The estimated association is therefore negative at each of these three score values, but it becomes less negative as state score increases. The difference between the slopes at \(S=1\) and \(S=-1\) is \(-0.08-(-0.32)=0.24\) log-expression units per allele. These are model-based conditional associations, not evidence that genotype causes a cell to move along the state axis.

Testing Whether the Genetic Effect Varies With State

The direct test of state dependence concerns the interaction coefficient. The null hypothesis is that the per-allele expression association is constant over the modeled state range. For a single interaction term, the hypotheses are

H_0:\beta_{GS}=0 \qquad\text{versus}\qquad H_a:\beta_{GS}\ne 0.

A Wald test divides the estimated interaction coefficient by its standard error. With an appropriate \(t\) reference distribution and degrees of freedom for the fitted model, the statistic is \(t=\hat{\beta}_{GS}/SE(\hat{\beta}_{GS})\). The reference distribution and degrees of freedom must reflect the analysis method and donor-level dependence; a cell count alone does not provide the number of independent genetic observations.

Worked Example: Testing the Genotype-by-State Interaction

In a hypothetical donor-aware analysis, the estimated interaction is \(0.10\) log-expression units per allele per one-unit increase in state score, with \(SE=0.04\) in the same units. Suppose the fitted test uses \(df=40\). Test for a nonzero interaction at \(\alpha=0.05\), using a two-sided \(t\)-test.

The statistic is

t=\frac{0.10}{0.04}=2.50.

For a \(t\) distribution with 40 degrees of freedom, the two-sided p-value is \(2P(T_{40}\ge 2.50)\approx 0.01662\). Since \(0.01662<0.05\), reject \(H_0\) at the 5% significance level. The data provide evidence that the per-allele expression association varies with state under this model. This p-value is not the probability that \(H_0\) is true.

For reference, the 95% confidence interval is approximately

0.10\pm t_{0.975,40}(0.04) =0.10\pm(2.021)(0.04) =[0.0192,\ 0.181].

Under the model and its assumptions, this interval gives a range of plausible values for the change in the per-allele expression association for each one-unit increase in state score. Its units are log-expression units per allele per state-score unit. It does not describe the genetic effect at a particular state; that requires \(\beta_G+\beta_{GS}s\).

Why the Interaction Formula Holds

The interaction interpretation follows by comparing the model-predicted expression at two genotypes while holding state, covariates, and donor effect fixed. Increasing genotype dosage by one changes the genotype term by \(\beta_G\) and the product term by \(\beta_{GS}s\); all terms that do not involve genotype cancel.

\begin{aligned} E[Y\mid G+1,S=s,C,b]-E[Y\mid G,S=s,C,b] &=\beta_G+\beta_{GS}s. \end{aligned}

This also shows why the interaction tests a change in the slope rather than a difference in average expression between state values. If \(\beta_{GS}=0\), the per-allele association is the same at all state values in the model, even if expression itself changes substantially with state through \(\beta_S\).

Worked Example: Expression Association Versus State Abundance

Consider a hypothetical dataset with three genotype groups and two representative regions of a continuous state axis, described here only to make the distinction easy to see. Suppose the mean log expression is 2 units in the low-state region and 6 units in the high-state region for every genotype. Within either region, genotype has no expression association. However, the proportions of cells in the high-state region are 20% for \(G=0\), 50% for \(G=1\), and 80% for \(G=2\).

The mean expression across all cells in each genotype group is a weighted average. For \(G=0\), it is

$$ (0.80)(2)+(0.20)(6)=1.60+1.20=2.80. $$

For \(G=1\), the mean is \((0.50)(2)+(0.50)(6)=1.00+3.00=4.00\). For \(G=2\), it is \((0.20)(2)+(0.80)(6)=0.40+4.80=5.20\). The overall means rise by \(1.20\) log-expression units per genotype step, even though the within-region expression association is zero in both regions.

The pattern arises because genotype groups have different distributions of cells across states, not because genotype is associated with different expression at a fixed state. In a continuous-state analysis, the abundance question concerns how the distribution of \(S\) varies with genotype; the expression interaction question concerns how \(Y\) varies with genotype conditional on \(S\). An overall expression comparison can combine both patterns and should not, by itself, be described as evidence for a state-dependent expression effect.

Assumptions, Interpretation, and Practical Checks

A genotype-by-state coefficient is interpretable only relative to the chosen state score, its scale, and the range of states represented in the data. Centering makes the main genotype coefficient refer to \(S=0\); scaling makes the interaction coefficient refer to a one-unit change on the scaled score. If the state score is not centered, the genotype main effect refers to the arbitrary score value zero, which may lie outside the observed range.

The linear interaction also imposes a shape: the genetic association changes at a constant rate per score unit. Before interpreting that rate, inspect the observed score range and the genotype representation across it. A fitted slope can be misleading when some genotypes are observed only at low scores and others only at high scores, because the model may rely on extrapolation rather than comparisons supported by data. A nonlinear relationship may require a more flexible state function and corresponding genotype-by-function terms; the simple linear interaction does not establish that all state dependence is linear.

The response model must suit the expression measurement. A linear mixed model can be a reasonable working model for a suitably transformed expression outcome, while sparse counts may call for a count model. The repeated-cell structure must also be represented appropriately, as in the earlier material on cell-level mixed models and repeated cell types within donors. In either case, genotype is measured at the donor level, so many cells from a small number of donors do not amount to the same evidence as many independent donors.

Finally, a cell-state score is a measured covariate, not automatically an externally fixed cause. If genotype is related to the state itself, conditioning on state estimates a conditional expression association and need not equal a total genetic association across all cells. The model alone does not establish a causal mechanism or separate all pathways linking genotype, state, and expression.

Common Mistakes

Interpreting \(\beta_G\) as the effect everywhere. It is the per-allele association only at \(S=0\). At another score \(s\), use \(\beta_G+\beta_{GS}s\).
Calling an abundance shift an expression interaction. A genotype-related change in the fraction of cells at high scores is about the state distribution. It does not show a different genotype-expression association among cells at the same score.
Treating cells as independent genetic replicates. Genotype is shared within donor. Account for donor clustering or use an appropriate donor-level summary and inferential procedure.
Interpreting an interaction without checking scale and support. State what one score unit means, report the observed range, and avoid extrapolating a linear slope beyond well-supported values.
Turning a conditional association into a causal claim. The interaction describes how a modeled association varies with measured state; it does not show that genotype causes state changes or that state causes the expression pattern.

Key takeaway: A genotype-by-continuous-state term tests whether the conditional per-allele expression association changes along the state score. Its coefficient is a change in that association per score unit. A genotype-related shift in the distribution of cells across states is a separate abundance pattern and must not be substituted for an expression interaction.

Check Your Understanding

Try these questions using the model and distinctions introduced above.

If \(\hat{\beta}_G=0.05\) and \(\hat{\beta}_{GS}=-0.03\), what is the estimated per-allele association at \(S=2\), and what are its units if expression is measured in log units?
What does \(H_0:\beta_{GS}=0\) say about the genetic expression association across state values?
Why can a difference in overall mean expression between genotype groups occur even when expression has no genotype association at a fixed state?
What is the interpretation of \(\beta_G\) when the state score has not been centered?
Why does observing thousands of cells from a few donors not provide thousands of independent genotype observations?