Why Regress Expression on Genotype Dosage?
An expression quantitative trait locus (eQTL) analysis asks whether genetic variation is associated with variation in the expression phenotype of a molecular feature. For a particular variant and feature, a simple starting point is to represent each donor’s genotype by a dosage and regress the measured expression phenotype on that dosage. This turns the question into a familiar one: how much does the expected expression value change for a one-allele increase in dosage?
The additive model is compact and interpretable, but it makes a particular assumption: the expected expression values at dosages 0, 1, and 2 lie on a straight line. The estimated slope is therefore a per-allele association under that model. It is not automatically a causal effect, and the quality of its uncertainty estimate depends on the regression assumptions.
Dosage is the number of copies of the chosen allele: the usual genotype categories have values 0, 1, and 2. A dosage can also be fractional when it represents an estimated expected allele count rather than a directly observed category. In either case, the allele being counted must be stated, because changing the counted allele reverses the direction of the coefficient. The response \(Y_i\) is the expression phenotype as prepared for analysis; its units may be transformed or normalized units rather than raw molecule counts.
The Model and Its Assumptions
In the unadjusted model, the fitted value for donor \(i\) is \(\hat{Y}_i=\hat{\beta}_0+\hat{\beta}_GG_i\). The intercept is the fitted expression at dosage zero. The slope is the fitted change in expression for each additional effect allele. Consequently, the fitted difference between dosage 2 and dosage 0 is \(2\hat{\beta}_G\), provided the additive relationship is appropriate.
The model can be estimated as a least-squares line. For the usual standard error and exact small-sample \(t\)-test to be justified, the model needs more than a numerical fit. In particular, the conditional mean should be linear in dosage, the errors should have mean zero, donors should contribute independent errors, and the error variance should be constant across dosage values. For exact \(t\)-based inference in a small sample, errors are additionally assumed to be normally distributed. These assumptions concern the model errors conditional on the observed dosages; they do not assert that the genotype itself is normally distributed.
A single-variant regression is an association analysis. Even if its assumptions hold, the slope alone does not establish that the variant causes the expression difference. An association can reflect other sources of dependence between genotype and expression, and the simple model does not by itself rule them out.
Estimating the Per-Allele Coefficient
Write \(\bar{G}\) and \(\bar{x}\) for the sample means of the observed dosages and expression values, respectively. The least-squares slope is the sum of cross-products of centered dosage and expression divided by the sum of squared centered dosages. The intercept then makes the fitted line pass through the two sample means.
The denominator measures how much genotype dosage varies in the sample. If all donors have the same dosage, it is zero, so a slope cannot be estimated. Even when it is positive, little dosage variation gives less information about the per-allele association and typically a larger standard error.
Worked Example: Estimating a Slope From Six Donors
Suppose a hypothetical expression phenotype is recorded in normalized units for six unrelated donors. The counted allele is specified, and the model uses dosages 0, 1, and 2. The observed data are:
| Dosage \(G_i\) | Expression \(Y_i\) |
|---|---|
| 0 | 10.1 |
| 0 | 9.9 |
| 1 | 11.2 |
| 1 | 10.8 |
| 2 | 12.1 |
| 2 | 11.9 |
The dosage mean is \(\bar{G}=(0+0+1+1+2+2)/6=1\), and the expression mean is \(\bar{x}=66/6=11\) normalized units. The centered-dosage sum of squares is
The centered cross-product sum is \(S_{GY}=(-1)(-0.9)+(-1)(-1.1)+(0)(0.2)+(0)(-0.2)+(1)(1.1)+(1)(0.9)=4\). Therefore,
The intercept is \(\hat{\beta}_0=11-(1)(1)=10.0\) normalized units. The fitted means at dosages 0, 1, and 2 are thus 10.0, 11.0, and 12.0. In this hypothetical sample, the estimated expression is one normalized unit higher per additional counted allele under the additive model. This is an estimated association, not evidence by itself that adding an allele would cause expression to rise.
Residual Variation and the Standard Error
The residual for donor \(i\) is \(e_i=Y_i-\hat{Y}_i\). In a simple regression with an intercept and one slope, two parameters have been estimated, so the residual degrees of freedom are \(df=n-2\). The residual variance estimate is the residual sum of squares divided by those degrees of freedom. Its square root, \(s\), estimates the standard deviation of the model errors. The standard error of the slope then scales this estimated error spread by the amount of genotype variation.
A short justification explains the denominator. With the observed dosages treated as fixed, the slope estimator is a weighted sum of the expression values, with weights proportional to \(G_i-\bar{G}\). Under independent errors with common variance \(\sigma^2\), the variance of that weighted sum is
Replacing the unknown \(\sigma\) by the residual estimate \(s\) gives the stated standard error. This calculation shows why the error spread and genotype spread both matter: greater residual variation increases uncertainty, while greater dosage variation decreases it. The formula relies on independent, equal-variance errors; the exact small-sample \(t\) distribution additionally relies on normal errors.
Worked Example: Calculating the Slope Standard Error
Continue with the six-donor data above. The fitted expression values are 10.0, 10.0, 11.0, 11.0, 12.0, and 12.0. The residuals are 0.1, \(-0.1\), 0.2, \(-0.2\), 0.1, and \(-0.1\) normalized units, giving
Here \(n=6\), so \(df=6-2=4\), and \(s^2=0.12/4=0.0300\) squared normalized units. Since \(S_{GG}=4\),
This standard error describes the estimated sampling variability of the slope under the specified regression assumptions. It is not the standard deviation of individual expression measurements; that is estimated here by \(s=\sqrt{0.0300}\) normalized units.
Testing for an Expression Association
The usual two-sided association test asks whether the population slope in the model is zero. The null and alternative hypotheses are \(H_0:\beta_G=0\) and \(H_a:\beta_G\ne 0\). Under the model assumptions, the test statistic is the estimated slope minus its null value, divided by its standard error. Under \(H_0\), it follows a \(t\) distribution with \(n-2\) degrees of freedom when the errors are normally distributed.
The two-sided p-value is the probability, assuming \(H_0\) and the test assumptions, of observing a \(t\) statistic at least as far from zero as the one calculated. A small p-value is evidence against a zero slope in this model; it is not the probability that \(H_0\) is true. The sign and size of \(\hat{\beta}_G\) still matter for describing the estimated direction and magnitude.
Worked Example: Testing a Dosage Association
Consider a second hypothetical six-donor dataset with dosages \(0,0,1,1,2,2\). Let \(c=\sqrt{11/24}\), approximately 0.6770. The expression values, in normalized units, are \(5+c,\ 5-c,\ 6.5,\ 6.5,\ 8+c,\ 8-c\). The exact values are specified by \(c\); their decimal approximations are 5.677, 4.323, 6.500, 6.500, 8.677, and 7.323. These data are constructed for illustration, not taken from a study.
The dosage mean is 1, the expression mean is 6.5, and \(S_{GG}=4\). The centered cross-product sum is \(S_{GY}=(-1)(-1.5+c)+(-1)(-1.5-c)+(1)(1.5+c)+(1)(1.5-c)=6\). Thus \(\hat{\beta}_G=6/4=1.50\) normalized expression units per effect allele, and the intercept is \(6.5-(1.5)(1)=5.00\).
The fitted values are 5.0, 5.0, 6.5, 6.5, 8.0, and 8.0. The residuals are \(c,-c,0,0,c,-c\), so \(RSS=4c^2=4(11/24)=11/6\). With \(df=6-2=4\), the residual variance estimate is \((11/6)/4=11/24\). The slope standard error is therefore
The test statistic for \(H_0:\beta_G=0\) is \(t=1.50/0.3385\approx 4.431\), with 4 degrees of freedom. The two-sided p-value is approximately 0.01141. At significance level \(\alpha=0.05\), this is evidence against a zero dosage slope under the model assumptions. The estimated association is positive and equals 1.50 normalized expression units per effect allele; the p-value does not establish causation or measure the probability that the null is true.
Uncertainty in the Per-Allele Estimate
A confidence interval supplements a test by showing a range of slope values compatible with the data and the model. For a two-sided 95% interval under the normal-error regression model, use the appropriate \(t\) critical value with \(n-2\) degrees of freedom:
Worked Example: A 95% Interval for the Slope
For the second hypothetical dataset, \(\hat{\beta}_G=1.50\), \(SE(\hat{\beta}_G)\approx0.3385\), and \(df=4\). The 95% critical value is approximately 2.776, so the margin of error is \(2.776(0.3385)\approx0.9397\). The interval is
Under repeated sampling and the model assumptions, this interval procedure would cover the true slope in 95% of samples. For these data, the interval gives a range of plausible per-allele slopes under the model. It does not mean there is a 95% probability that the fixed slope lies within this particular calculated interval.
Interpretation and Common Pitfalls
A regression coefficient must be interpreted in the units of the phenotype and with the counted allele made explicit. If the expression phenotype was transformed before analysis, the coefficient is on that transformed scale. A two-allele contrast is twice the fitted per-allele difference only because the additive model imposes a linear dosage relationship. The regression does not demonstrate that the observed means at the three genotype categories truly follow that relationship.
- Reversing the allele without changing the sign. If dosage is recoded to count the other allele, the direction of the slope reverses. State which allele the dosage counts.
- Calling a slope causal. A nonzero slope describes an association under the model. It does not, by itself, establish that changing the genotype would change expression.
- Reading the p-value as an effect size. A small p-value does not say that the expression difference is large or biologically important. Report the coefficient and its units as well.
- Treating a standard error as individual variability. The slope standard error quantifies uncertainty in the estimated coefficient; residual spread and coefficient uncertainty are different quantities.
- Ignoring model assumptions. The usual standard error and \(t\)-based inference can be unreliable if the mean relationship is not adequately linear or errors are dependent or heteroskedastic.
In an eQTL analysis, the result is most informative when the coefficient, allele coding, phenotype scale, standard error or interval, and test are reported together. A regression p-value for one variant-feature pair is not a complete summary of the broader evidence, and it should be interpreted alongside study design and the assumptions used to obtain it.
Check Your Understanding
Answer these questions using the model and calculations above.
- In the additive dosage model, what does \(\beta_G\) measure, and what are its units?
- Why does the least-squares slope require variation in the observed genotype dosages?
- For a simple regression with an intercept and one dosage predictor, why are the residual degrees of freedom \(n-2\)?
- What does the two-sided p-value quantify under \(H_0\), and what does it not quantify?
- If the dosage is recoded to count the alternate allele instead of the originally counted allele, what happens to the slope’s direction?