From Pseudobulk Profiles to a QTL Regression
A cell-type QTL asks whether genetic variation is associated with a molecular measurement in a particular cell type. In a pseudobulk analysis, the observation used for regression is not an individual cell. It is a donor-level expression measurement formed by aggregating counts across that donor’s cells of the specified type and then normalizing the resulting profile. This makes the donor, rather than each captured cell, the independent sampling unit.
For a particular gene and cell type, the analysis has one outcome per donor with usable data. The predictor of primary interest is the donor’s genotype dosage at a genetic variant. Covariates can account for measured features such as age, sex, technical batch, or ancestry-related population structure. The regression describes an adjusted association; by itself, it does not establish that the variant causes an expression change.
The outcome is specific to a gene and cell type, so the model is fitted separately for each tested gene–variant–cell-type combination (or in an equivalent joint modeling framework). A covariate with a categorical value, such as batch, is represented by appropriate indicator variables. The intercept \(\beta_0\) is the expected outcome at genotype dosage zero and covariate reference values; it is usually less scientifically central than \(\beta_G\).
What Is Being Regressed, and What Does the Coefficient Mean?
Pseudobulk construction sums the raw counts for each gene across the cells assigned to the same donor and cell type. Normalization is then applied to the donor-by-cell-type profile. Common outcome choices include log-transformed counts per million (often with a small prior count), log-transformed size-factor-normalized counts, or a variance-stabilized expression measure. The choice should be made consistently across donors and stated when reporting results.
Normalization is not a cosmetic step: it determines the scale on which the coefficient is expressed. If \(Y_d\) is log base 2 counts per million with a prior count, then \(\beta_G\) is in log base 2 expression units per allele. When the prior count has little influence, \(2^{\beta_G}\) can be described approximately as the multiplicative change in normalized abundance per additional allele. It is not an absolute change in raw molecules, and the approximation is not exact when the prior count materially affects the transformation.
Library-size scaling aims to make profiles comparable when donors contribute different total numbers of reads or cells. Other normalization methods may also adjust for composition or mean–variance behavior. They do not all produce the same numerical outcome or estimand. An analysis should therefore report the aggregation and normalization method, any transformation, and the units or scale of the outcome. A regression on a transformed normalized outcome estimates an association on that transformed scale.
Worked Example: Normalizing Two Pseudobulk Counts
Suppose two hypothetical donors have pseudobulk counts of 80 and 120 for a gene in the same cell type. Their total pseudobulk library sizes are 400,000 and 600,000 counts, respectively. For illustration, calculate counts per million (CPM) and log base 2 CPM after adding a prior count of 1: $$ \mathrm{CPM}_d=\frac{C_d}{L_d}\times 1{,}000{,}000,\qquad Y_d=\log_2(\mathrm{CPM}_d+1), $$ where \(C_d\) is the gene count and \(L_d\) is the library size.
For donor 1, \(\mathrm{CPM}_1=(80/400{,}000)\times1{,}000{,}000=200\), so \(Y_1=\log_2(201)\approx7.651\). For donor 2, \(\mathrm{CPM}_2=(120/600{,}000)\times1{,}000{,}000=200\), so \(Y_2=\log_2(201)\approx7.651\).
The raw counts differ, but the library-size-scaled outcomes are equal. This illustrates what the scaling does; it does not establish that CPM is the best normalization for every data set. In a QTL regression using this outcome, the genotype coefficient would be measured in log base 2 CPM units per allele, with the prior count included in the outcome definition.
Estimating and Interpreting the Genotype Effect
In the model, the estimate \(\hat{\beta}_G\) summarizes the fitted change in expression associated with one additional allele, holding the included covariates fixed. Genotype coding matters: dosage 0, 1, 2 counts copies of a chosen allele, so the sign of the coefficient depends on which allele is designated. The linear dosage model also imposes an additive pattern: the expected difference from dosage 0 to 1 is the same as from 1 to 2.
With no covariates, the ordinary least squares slope is the sample covariance of genotype and outcome divided by the sample variance of genotype, in its sum-of-squares form. With covariates, the genotype coefficient is the association remaining after adjustment for those covariates. The following calculation shows the unadjusted fit and its uncertainty under a simple set of assumptions.
Worked Example: An Unadjusted Per-Allele Estimate
Consider six hypothetical independent donors with the following normalized expression values, in log base 2 CPM units, for one gene in one cell type. Assume no covariates are included.
| Genotype dosage \(G_d\) | Expression values \(Y_d\) |
|---|---|
| 0 | 1.0, 1.2 |
| 1 | 1.8, 2.0 |
| 2 | 2.5, 2.7 |
The genotype mean is \(\bar{G}=1\), and the expression mean is \(\bar{Y}=11.2/6=1.8667\). The genotype sum of squares is $$ S_{GG}=\sum_d(G_d-\bar{G})^2=2(0-1)^2+2(1-1)^2+2(2-1)^2=4. $$ The cross-product sum is $$ S_{GY}=\sum_d(G_d-\bar{G})(Y_d-\bar{Y})=3.000, $$ so the fitted slope is \(\hat{\beta}_G=S_{GY}/S_{GG}=3/4=0.7500\) log base 2 CPM units per allele. The fitted intercept is \(\hat{\beta}_0=\bar{Y}-\hat{\beta}_G\bar{G}=1.1167\).
For the usual homoskedastic linear-model standard error, the residuals are approximately \(-0.1167,\ 0.0833,\ -0.0667,\ 0.1333,\ -0.1167,\ 0.0833\). Their squared sum is \(\mathrm{SSE}=0.06333\), giving residual variance estimate \(s^2=\mathrm{SSE}/(n-2)=0.01583\). Thus $$ SE(\hat{\beta}_G)=\sqrt{\frac{s^2}{S_{GG}}} =\sqrt{\frac{0.01583}{4}}\approx0.06291. $$ Under independent, normally distributed errors with constant variance, a 95% \(t\)-interval with \(df=6-2=4\) uses \(t_{0.975,4}=2.776\). It is \(0.7500\pm2.776(0.06291)=[0.5753,0.9247]\) log base 2 CPM units per allele.
Under those model assumptions, the interval describes a range of plausible values for the population per-allele association on this outcome scale. It does not mean there is a 95% probability that the fixed coefficient lies inside this particular interval. The six-donor example is deliberately small: its residual pattern cannot establish that normality or constant variance is reasonable, so the interval is an illustration of the calculation, not strong evidence for a biological conclusion.
Covariate Adjustment and a Short Derivation
Covariates are included when they are relevant to the outcome, the genotype association, or the study design. For example, ancestry-related covariates may help address population structure, and technical batch may account for systematic differences in measurement. Including covariates changes the comparison represented by \(\hat{\beta}_G\): it is an association conditional on the modeled covariates. It does not automatically remove unmeasured confounding or correct a misspecified model.
The adjusted genotype coefficient has a useful interpretation through residualization. Let \(W\) contain the intercept and all covariates. First regress genotype \(G\) on \(W\) and retain its residuals \(r_G\). Separately regress expression \(Y\) on \(W\) and retain residuals \(r_Y\). Then the coefficient of \(G\) in the regression of \(Y\) on \(G\) and \(W\) equals the slope from regressing \(r_Y\) on \(r_G\): $$ \hat{\beta}_G=\frac{\sum_d r_{G,d}r_{Y,d}}{\sum_d r_{G,d}^2}. $$
This follows by writing the least-squares residual sum of squares for a candidate genotype coefficient \(b\) after optimizing the coefficients of \(W\). The parts of \(Y\) and \(G\) explained by \(W\) are removed, leaving the residual sum of squares \(\sum_d(r_{Y,d}-b r_{G,d})^2\). Differentiating with respect to \(b\) and setting the derivative to zero gives \(-2\sum_d r_{G,d}r_{Y,d}+2b\sum_d r_{G,d}^2=0\), which yields the stated slope when \(\sum_d r_{G,d}^2>0\). Adjustment therefore uses only genotype variation not explained by the included covariates.
Worked Example: Adjusting for a Covariate
Suppose four hypothetical donors have genotypes \(G=(0,1,1,2)\), a centered covariate \(Z=(0,0,1,1)\), and normalized expression \(Y=(2.0,3.7,4.3,5.8)\), in log base 2 CPM units. Fit a model with an intercept, genotype, and \(Z\). To calculate the adjusted genotype slope, regress \(G\) on the intercept and \(Z\). Within \(Z=0\), the genotype mean is \(0.5\); within \(Z=1\), it is \(1.5\). The fitted genotype values are therefore \((0.5,0.5,1.5,1.5)\), giving \(r_G=(-0.5,0.5,-0.5,0.5)\).
Next, regress \(Y\) on the intercept and \(Z\). The mean expression is \(2.85\) when \(Z=0\) and \(5.05\) when \(Z=1\), so the fitted values are \((2.85,2.85,5.05,5.05)\). Thus \(r_Y=(-0.85,0.85,-0.75,0.75)\). The adjusted slope is $$ \hat{\beta}_G= \frac{(-0.5)(-0.85)+(0.5)(0.85)+(-0.5)(-0.75)+(0.5)(0.75)} {(-0.5)^2+(0.5)^2+(-0.5)^2+(0.5)^2} =\frac{1.60}{1.00}=1.60. $$
The fitted association is 1.60 log base 2 CPM units per allele, conditional on \(Z\). This example demonstrates the adjustment calculation, not reliable inference: with four observations and three fitted coefficients, only one residual degree of freedom remains. It also requires that there be genotype variation after adjustment, as the nonzero denominator confirms here.
Assumptions and Limits of Interpretation
The regression coefficient can be computed without every assumption being true, but standard errors, intervals, and tests rely on assumptions about the model and sampling. For the usual ordinary least squares inference, the key conditions include:
- Independent donor-level errors. Donors contribute independent observations, or dependence is handled appropriately. Multiple cells from one donor do not create additional independent genotype observations.
- A suitable conditional mean model. The expected normalized expression is represented adequately by the included covariates and the assumed linear dosage effect. The additive coding should be scientifically and statistically defensible.
- Residual variance for the chosen inference. The conventional standard error assumes constant residual variance across fitted values. If this is not credible, a suitable robust variance estimator may be considered; it does not repair other model problems.
- For exact small-sample \(t\)-inference, normally distributed errors. With larger samples, inference may be less sensitive to departures from normality, but outliers or strong misspecification still matter.
- Usable genotype information and no exact collinearity. There must be variation in genotype dosage not fully explained by the covariates. A coefficient cannot be separated from covariate effects if the predictors are perfectly collinear.
Population structure can induce genotype–expression associations if ancestry is related to both allele frequencies and expression. Relevant ancestry covariates can help, but their adequacy depends on the data and design. Likewise, batch adjustment is useful only when batch and genotype effects can be distinguished; if genotype is completely confounded with batch, the data do not identify separate effects. A QTL association is not automatically causal, even when the regression is correctly fitted.
Why the Donor-Level Formulation Matters
The regression matches the level at which genotype varies: in a typical study, each donor has one genotype dosage for the variant, while many cells from that donor share it. Treating every cell as an independent genotype observation would artificially inflate the apparent sample size and can make uncertainty look much smaller than it is. Pseudobulk regression instead uses one expression outcome per donor, gene, and cell type, so its sample size is based on donors with usable observations for that analysis.
The result is specific to the selected cell type, gene, genotype coding, normalization, and covariates. A positive coefficient means that greater dosage is associated with higher expected expression on the modeled scale, conditional on the covariates. It does not by itself show that the effect differs from that in another cell type, nor does it distinguish changes in within-cell expression from changes in cell composition. Those are different scientific questions.
Common Mistakes
- Counting cells as independent genotype observations. This confuses the number of captured cells with the number of independent donors. Use donor-level pseudobulk observations for this regression.
- Reporting a coefficient without its expression scale. “An increase of 0.5” is incomplete. State, for example, that the fitted association is 0.5 log base 2 CPM units per additional allele.
- Calling an adjusted association causal. Covariate adjustment does not guarantee that all confounding or selection effects are absent. Describe the result as an association unless a causal design and assumptions justify stronger language.
- Assuming all normalization choices estimate the same numerical effect. Transformations and scaling change the outcome and hence the coefficient’s units and interpretation. Report the chosen method.
- Ignoring dosage coding or lack of genotype variation. The sign depends on the designated allele, and adjustment cannot identify a genotype effect when no residual genotype variation remains.
Check Your Understanding
Answer each question using the donor-level regression framework.
- In a model with log base 2 CPM as the outcome, what are the units and interpretation of a genotype coefficient of \(-0.3\)?
- Why does aggregating cells into donor-level pseudobulk profiles not make the number of cells the regression sample size?
- How can a covariate affect the genotype coefficient, and what must be true for the adjusted coefficient to be identifiable?
- What does adding a prior count before taking log base 2 change about the interpretation of normalized expression?
- Name two assumptions needed for conventional ordinary least squares inference, and explain one way population structure can threaten interpretation.