Worked Donor-Level Analysis of a Cell-Type Shift

From Cell Counts to a Donor-Level Comparison

A cell-type shift is a comparison of proportions across biological samples, not a comparison of individual cells as if each were an independent replicate. For each donor, the observed proportion for a cell type is its cell count divided by the total number of cells profiled for that donor. Variation among donors can be substantially larger than binomial sampling variation alone predicts. A beta-binomial regression provides a way to model that extra variation while adjusting for a donor-level covariate.

The hypothetical data below include six donors in each condition. The four cell-type counts sum to the same total of 1,000 cells per donor. The covariate is age, measured in years. These data are illustrative, not from a real study.

Donor	Condition	Age (years)	Target type	Type B	Type C	Type D	Total cells
1	Control	42	94	306	250	350	1,000
2	Control	47	111	289	250	350	1,000
3	Control	51	98	302	250	350	1,000
4	Control	54	105	295	250	350	1,000
5	Control	58	92	308	250	350	1,000
6	Control	62	100	300	250	350	1,000
7	Treatment	43	153	287	240	320	1,000
8	Treatment	48	169	271	240	320	1,000
9	Treatment	52	158	282	240	320	1,000
10	Treatment	55	172	268	240	320	1,000
11	Treatment	59	151	289	240	320	1,000
12	Treatment	63	164	276	240	320	1,000

The Beta-Binomial Regression Model

Definition: For donor \(i\), let \(Y_i\) be the number of cells assigned to the target type and \(N_i\) the total number of profiled cells. A beta-binomial model treats \(Y_i\) as beta-binomial with total \(N_i\), mean proportion \(\mu_i\), and overdispersion parameter \(\rho\). In a regression model, the mean is linked to donor-level predictors by

\operatorname{logit}(\mu_i)=\log\left(\frac{\mu_i}{1-\mu_i}\right) =\beta_0+\beta_1 T_i+\beta_2(A_i-52.5),

Here \(T_i=1\) for treatment and \(T_i=0\) for control, and \(A_i\) is age in years. Centering age at 52.5 makes \(\beta_0\) the control log odds at age 52.5 and makes \(\beta_1\) the adjusted treatment-versus-control log odds ratio at that same age. With this model, the condition effect is assumed constant across ages; an interaction would be needed to let it vary.

The beta-binomial assumption allows donors to have different underlying proportions, beyond the binomial sampling variation conditional on a single proportion. One representation is \(P_i\sim\operatorname{Beta}(a_i,b_i)\) and \(Y_i\mid P_i\sim\operatorname{Binomial}(N_i,P_i)\), where \(a_i=\mu_i(1-\rho)/\rho\) and \(b_i=(1-\mu_i)(1-\rho)/\rho\) for \(0<\rho<1\). Here \(\mu_i\) is the predictor-dependent mean and \(\rho\) is common across donors; the binomial case is obtained as the limit \(\rho\to0\). In regression, the mean varies with predictors while the overdispersion parameter describes residual donor-to-donor variation.

The variance makes the role of overdispersion explicit. By the conditional-variance identity, \(\operatorname{Var}(Y_i)=E[\operatorname{Var}(Y_i\mid P_i)]+\operatorname{Var}(E[Y_i\mid P_i])\). Since \(\operatorname{Var}(Y_i\mid P_i)=N_iP_i(1-P_i)\) and \(E[Y_i\mid P_i]=N_iP_i\), this gives

\begin{aligned} \operatorname{Var}(Y_i) &=N_iE[P_i(1-P_i)]+N_i^2\operatorname{Var}(P_i)\\ &=N_i\mu_i(1-\mu_i)\{1+(N_i-1)\rho\}. \end{aligned}

Thus \(\rho=0\) recovers binomial variance, while positive \(\rho\) inflates the variance. Even modest overdispersion can matter when each donor contributes many cells: the effective information is not equivalent to treating every cell as an independent biological replicate.

Assumptions Before Fitting

The condition test and its uncertainty are meaningful only if the model and design support them. The principal assumptions are:

Donors are the independent biological units, and their observations are independent under the sampling design. Cells from one donor are not independent replicates.
The numerator is a count of target-type cells among the stated total, and donor-level counts are correctly recorded. The beta-binomial model accounts for extra-binomial variation, but not arbitrary dependence or data errors.
The logit mean is adequately described by the stated predictors. In particular, the age effect is linear on the log odds scale, and the condition effect is constant across age in this model.
There is enough overlap in age between conditions to distinguish an adjusted condition effect from an age effect. A covariate cannot repair a design in which condition and age are perfectly confounded.
The observed proportions are interpreted as relative proportions among profiled cells. Differences in capture or annotation can affect those proportions, and this model alone does not estimate absolute cell numbers.

The model is fitted by maximizing the beta-binomial likelihood, estimating the regression coefficients and overdispersion. A Wald interval and test use an approximately normal sampling distribution for the estimated coefficient. With few donors, that approximation can be fragile; profile-likelihood intervals or a donor-level bootstrap may be preferable when supported by the analysis plan and software.

Worked Example: Estimate the Adjusted Condition Effect

Suppose fitting the target-type model to the hypothetical donor data gives \(\hat{\beta}_1=0.620\) with \(SE=0.240\). These are stipulated illustrative model output; the arithmetic below shows how to interpret and report the fitted result.

The coefficient is an adjusted log odds ratio. Exponentiating gives the odds ratio for treatment relative to control at the same age:

\widehat{OR}=\exp(0.620)=1.859.

A Wald 95% confidence interval for the log odds ratio is

0.620\pm1.96(0.240)=0.620\pm0.4704=(0.1496,\ 1.0904).

Exponentiating the endpoints gives an odds-ratio interval of \((\exp(0.1496),\exp(1.0904))=(1.161,\ 2.975)\). Under the fitted model, at a fixed age the estimated odds of belonging to the target type among profiled cells are about 1.86 times as high in treatment as in control; the 95% interval is approximately 1.16 to 2.98. The interval describes uncertainty in the model-based odds ratio, not a 95% probability that a fixed parameter lies in this particular interval.

The Wald statistic and two-sided p-value are

z=\frac{0.620}{0.240}=2.583,\qquad p=2\{1-\Phi(2.583)\}\approx0.0098.

This p-value is the probability, under a zero condition coefficient and the test’s model assumptions, of obtaining a statistic at least this extreme in absolute value. It is not the probability that the null hypothesis is true.

Worked Example: Translate the Coefficient to Fitted Proportions

Suppose the fitted intercept is \(\hat{\beta}_0=-2.20\), with age centered at 52.5 years as above. At age 52.5, the fitted control proportion is

\hat{\mu}_{C}=\frac{1}{1+\exp(-(-2.20))}=0.09975.

For a treatment donor of the same age, the fitted log odds are \(-2.20+0.620=-1.58\), so

\hat{\mu}_{T}=\frac{1}{1+\exp(-(-1.58))}=0.1708.

The fitted difference at that reference age is \(0.1708-0.09975=0.0711\), or about 7.11 percentage points among profiled cells. This is a model-based contrast at age 52.5, not a universal absolute change: the percentage-point difference depends on the baseline proportion. It also does not imply that treatment caused the shift.

At age 62.5, ten years above the reference, the control linear predictor is \(-2.20+0.015(10)=-2.05\), using the illustrative age coefficient \(\hat{\beta}_2=0.015\). The treatment predictor is \(-2.05+0.620=-1.43\). Applying the inverse logit gives approximately 0.114 for control and 0.193 for treatment. The odds ratio remains \(\exp(0.620)=1.859\), but the fitted proportion difference is now about 7.9 percentage points. This illustrates why a log odds ratio and a proportion difference are not interchangeable effect measures.

Multiplicity Across Cell Types

A cell-type analysis commonly fits one condition test per cell type. If several null hypotheses are tested, reporting each unadjusted p-value as though it were the only test inflates the chance of at least one false discovery. The Benjamini–Hochberg procedure controls the false discovery rate under independence and certain forms of positive dependence among tests.

For \(m\) p-values, sort them from smallest to largest, \(p_{(1)}\le\cdots\le p_{(m)}\). At false discovery rate level \(q\), find the largest rank \(k\) satisfying \(p_{(k)}\le kq/m\). Reject all hypotheses with ranks 1 through \(k\); if no rank qualifies, reject none. The adjusted p-value for a test is its Benjamini–Hochberg value, computed by taking the running minimum of \(mp_{(j)}/j\) over ranks \(j\) at or above its rank, capped at 1.

Worked Example: Adjust Four Cell-Type Tests

Suppose the four condition tests from the same hypothetical analysis have p-values 0.0098 for the target type, 0.12 for Type B, 0.40 for Type C, and 0.08 for Type D. For \(m=4\) tests and \(q=0.05\), the ordered p-values are 0.0098, 0.08, 0.12, and 0.40. The rank thresholds \(kq/m\) are 0.0125, 0.0250, 0.0375, and 0.0500.

At rank 1, \(0.0098\le0.0125\), so rank 1 qualifies. Ranks 2, 3, and 4 do not qualify because \(0.08>0.0250\), \(0.12>0.0375\), and \(0.40>0.0500\). The largest qualifying rank is therefore 1, so only the first-ranked test is selected at this threshold. In general, the procedure selects all tests up to the largest qualifying rank, not just a qualifying rank considered in isolation.

For the target test, the adjusted p-value is the minimum of the rank-based values \(4(0.0098)/1=0.0392\), \(4(0.08)/2=0.160\), \(4(0.12)/3=0.160\), and \(4(0.40)/4=0.400\), which is 0.0392. The target-type result therefore remains below 0.05 after this correction. Report the unadjusted p-value and the adjusted p-value, and state the number of tests in the family.

Diagnostics, Alternatives, and Limits

A fitted model should be checked rather than accepted solely because it converged. Examine donor-level proportions against fitted values and covariates, residual patterns by condition, and whether a small number of donors drives the estimated effect. Inspect the estimated overdispersion and compare observed residual variability with model expectations. Warning signs include systematic age patterns left in residuals, pronounced condition-specific spread, fitted means near zero or one when observations are not, and highly influential donors. A strong pattern may indicate a missing nonlinear covariate effect, an unsuitable mean model, or data problems. These checks assess model adequacy; they do not prove the assumptions.

The beta-binomial model is one option rather than a universal choice. A binomial regression is simpler but can understate uncertainty when donor proportions are overdispersed. A quasi-binomial model allows variance inflation without specifying the full beta-binomial likelihood, but does not provide the same likelihood-based model comparison. Proportion-based methods such as propeller use donor-level proportions and can be useful alternatives, with their own distributional assumptions and diagnostics. The appropriate choice depends on the design, number of donors, count structure, and analysis goal; large cell counts do not compensate for very few independent donors.

The result concerns the fraction of captured and classified cells assigned to a type. Because proportions are compositional, an increase in one fraction can arise from a decrease in other components even if the target type’s absolute number is unchanged. Capture efficiency, tissue processing, and annotation can also influence the observed fractions. Consequently, the adjusted coefficient supports a statement about relative abundance in the profiled samples under the model, not an absolute cell-number claim or a causal treatment effect without a design that justifies causal interpretation.

Common pitfall: Treating thousands of cells as thousands of independent replicates can make standard errors far too small. The condition comparison must be supported by independent donors; the model uses each donor’s count and total while allowing donor-level overdispersion. Also avoid translating an odds ratio directly into a percentage-point change: calculate fitted proportions at a stated covariate value for that interpretation.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

If \(\hat{\beta}_1=0.40\) in the stated model, what is the adjusted odds ratio? What comparison does it represent?
For four tests at \(q=0.05\), the sorted p-values are 0.010, 0.030, 0.040, and 0.20. Which tests are selected by the Benjamini–Hochberg rule?
Why can a binomial model give overly narrow intervals for donor-level proportions?
Does a fitted target-type proportion increase establish that the absolute number of target cells increased? Explain.

Solutions. (1) \(\exp(0.40)=1.492\). It is the treatment-versus-control odds ratio at the same covariate value, under the model’s constant condition-effect assumption. (2) The thresholds are 0.0125, 0.0250, 0.0375, and 0.0500. Rank 1 qualifies, while ranks 2, 3, and 4 do not; the largest qualifying rank is 1, so only the test at rank 1 is selected. (3) A binomial model assumes no extra donor-to-donor variation beyond binomial sampling. When donor proportions are more variable, it understates the variance and can understate uncertainty. (4) No. A proportion is relative to the profiled total; changes in other cell types, capture, or annotation can alter it without an increase in the target type’s absolute number.