From Cells to a Donor-Level cis-eQTL
A single-cell data set may contain thousands of cells, but a genetic association analysis must respect the unit at which genotype varies. Ordinarily, every cell from the same donor has the same inherited genotype at a given variant. Treating those cells as independent genetic observations would therefore overstate the amount of independent information. A pseudobulk analysis instead summarizes the cells of a chosen type within each donor and relates those donor-level summaries to donor genotypes.
This tutorial works through that calculation for one gene and one variant in a specified cell type. The example uses a simple regression so that the arithmetic is visible. In an applied analysis, researchers would generally consider relevant donor-level covariates, prespecified quality criteria, and the study’s multiple-testing plan. The association estimate alone does not establish that the variant causes a change in expression.
Worked Example: Aggregating Counts Within Donors
Worked Example: Aggregating Counts Within Donors
Suppose a hypothetical study selects one cell type and one gene, Gene G. For each of three donors, the table gives the Gene G counts in the assigned cells. It also gives the number of assigned cells, to make clear what is being summed.
| Donor | Gene G counts across assigned cells | Number of cells | Pseudobulk count \(K_{gdc}\) |
|---|---|---|---|
| A | 4, 0, 3, 1 | 4 | 8 |
| B | 2, 5, 1 | 3 | 8 |
| C | 0, 2, 2, 1, 3 | 5 | 8 |
For donor A, the pseudobulk count is \(4+0+3+1=8\). For donor B it is \(2+5+1=8\), and for donor C it is \(0+2+2+1+3=8\). Thus, all three donor-level counts are 8, even though the donors contributed different numbers of cells and the counts were distributed differently across cells.
The units for the genetic association are the three donors, not the 12 individual cells. The sum retains the total observed gene count for the selected donor and cell type. It is not an average across cells, and the cell count itself is not the number of independent genotype observations.
Normalize Before Comparing Expression
Pseudobulk libraries can differ in size because donors contribute different numbers of cells, cells have different sequencing depths, or both. A raw count of 100 is not directly comparable across libraries if one library contains many more total reads than another. One simple illustration of library-size normalization is counts per million (CPM).
Worked Example: Calculating a Transformed Expression Value
For a hypothetical donor and cell type, suppose Gene G has pseudobulk count \(K_{gdc}=8\), and the total pseudobulk library size is \(L_{dc}=400\) counts. Calculate its CPM and transformed expression value.
First, the CPM is
Then the log-transformed value is
The result is in log base 2 units of CPM plus 1. These small library totals are chosen only to keep the arithmetic simple; real pseudobulk libraries are typically much larger. CPM and this transformation are an illustration of a possible expression scale, not a universal requirement for every count-based analysis.
Specify the Genotype Regression
Code genotype as the number of copies of a chosen effect allele: \(G_d=0\), \(1\), or \(2\). This is an allele-dosage coding. For a particular gene and cell type, a basic additive cis-eQTL model is
Here, \(y_{gdc}\) is the donor-level expression value, \(\beta_0\) is the intercept, and \(\beta_G\) is the expected change in expression for one additional copy of the coded allele. The residual \(\varepsilon_d\) represents donor-level variation not explained by genotype in this model. If the expression scale is \(\log_2(\mathrm{CPM}+1)\), then \(\beta_G\) is measured in log base 2 units of CPM plus 1 per allele copy.
A practical model may include donor-level covariates, such as ancestry principal components or technical variables, when appropriate to the study design. In that case the regression includes those covariates, and the genotype coefficient describes the association conditional on them. Covariates should be chosen and justified rather than added solely because they change the estimate.
Worked Example: Estimating the Allele Effect
Worked Example: Estimating the Allele Effect
Consider eight hypothetical donors. For this hand calculation, the table reports their cell-type pseudobulk expression values on the log2(CPM + 1) scale and their genotype dosages at one prespecified cis variant. No covariates are included in this deliberately simple example.
| Donor | Genotype dosage \(G_d\) | Expression \(y_{gdc}\) |
|---|---|---|
| 1 | 0 | 5.2 |
| 2 | 0 | 4.8 |
| 3 | 0 | 5.0 |
| 4 | 1 | 5.4 |
| 5 | 1 | 5.8 |
| 6 | 1 | 5.6 |
| 7 | 2 | 6.2 |
| 8 | 2 | 6.2 |
The genotype mean is \(\bar{G}=7/8=0.875\). The expression mean is \(\bar{y}=44.2/8=5.525\). The genotype sum of squares is
The cross-product is
Therefore, the least-squares genotype slope and intercept are
The fitted model is \(\hat{y}=5.000+0.600G\). The fitted values for genotypes 0, 1, and 2 are 5.0, 5.6, and 6.2, respectively. The residuals from the table are \(0.2,-0.2,0,-0.2,0.2,0,0,0\), so the residual sum of squares is \(0.16\). With two fitted parameters, the residual degrees of freedom are \(df=n-2=6\), giving residual variance estimate \(s^2=0.16/6=0.02667\).
For this simple regression, the standard error of the slope is
The test statistic for \(H_0:\beta_G=0\) is \(t=\hat{\beta}_G/SE(\hat{\beta}_G)=0.600/0.07396\approx 8.11\), with 6 degrees of freedom under the model assumptions. The two-sided \(p\)-value is approximately 0.0002. This is evidence against a zero genotype slope under this model; it is not the probability that the null hypothesis is true.
Using the 97.5th percentile of a \(t\) distribution with 6 degrees of freedom, approximately 2.447, the nominal 95% confidence interval is
Under the model, the estimated association is an increase of 0.600 log2(CPM + 1) units per additional copy of the coded allele, with a nominal 95% confidence interval from 0.419 to 0.781 such units per copy. Exponentiating the slope gives \(2^{0.600}\approx1.52\), a ratio on the CPM + 1 scale per allele copy. It should not be described as proof that the variant itself causes a 1.52-fold change in biological expression.
Why the Slope and Its Uncertainty Have These Forms
The simple-regression formulas follow directly from least squares. For a proposed slope \(b\), the intercept that minimizes the sum of squared residuals is \(\bar{y}-b\bar{G}\). Substituting this intercept into the residual sum of squares gives \(\sum_d[(y_{gdc}-\bar{y})-b(G_d-\bar{G})]^2\). Differentiating with respect to \(b\), setting the derivative to zero, and solving gives the slope \(S_{Gy}/S_{GG}\).
If the residual variance is \(\sigma^2\), the variance of this slope under the regression assumptions is \(\sigma^2/S_{GG}\): greater genotype spread across donors provides more information about the slope. Replacing \(\sigma^2\) by the residual mean square \(s^2\) gives \(SE(\hat{\beta}_G)=\sqrt{s^2/S_{GG}}\). With normally distributed, independent, constant-variance residuals, the standardized difference between the estimated slope and its true value follows a \(t\) distribution with \(n-2\) degrees of freedom. That result justifies the \(t\)-based confidence interval and test used above.
For a regression with covariates, the genotype coefficient is estimated while accounting for those covariates; its standard error depends on the full design matrix, not just \(S_{GG}\). The simple formula shown here applies to the intercept-and-genotype model in the worked example.
Interpretation and Common Pitfalls
A cis-eQTL result is an association between genotype and expression for a particular gene, cell type, and expression scale. A positive coefficient means that the coded allele dosage is associated with higher expected expression on that scale, conditional on the fitted model. The interpretation depends on which allele was coded as the effect allele; reversing the coding reverses the sign of the slope.
- Counting cells as independent donors. Thousands of cells from a small number of donors do not provide thousands of independent genotype observations. Aggregate within donor and cell type, and base genetic inference on independent donors.
- Interpreting a raw count as normalized expression. Raw pseudobulk counts can reflect library size as well as expression. State the normalization and expression scale before interpreting a coefficient.
- Reading an interval as a probability statement about the parameter. A nominal 95% confidence procedure has 95% coverage over repeated samples when its assumptions hold; it does not mean there is a 95% probability that this particular fixed parameter lies in the observed interval.
- Calling the association causal. The regression estimates an association, not the biological mechanism. Linkage disequilibrium, unmodeled structure, or other sources of confounding can contribute to an observed association.
- Treating a single-variant \(p\)-value as a study-wide conclusion. The example’s \(p\)-value is for one test. When many variants, genes, or cell types are examined, inference must account for the testing plan and multiplicity.
Check Your Understanding
Use the model and definitions above to answer these questions.
- Why is the donor, rather than each individual cell, the genetic observation in a donor-level pseudobulk cis-eQTL analysis?
- A donor has a gene count of 12 and a total pseudobulk library size of 600. Calculate the CPM and the log2(CPM + 1) value.
- In the worked regression, what does the estimated slope of 0.600 mean, including its units and allele-dosage interpretation?
- What model assumptions support the \(t\)-based confidence interval in the eight-donor example?
- Why does the small \(p\)-value in the worked example not demonstrate that the variant causes altered expression?