Why Library Size Enters the Mean Model
A pseudobulk sample with more sequenced reads will usually contain more observed counts for a gene, even when the gene has the same relative expression rate as in a shallower sample. A count model that compares raw counts without accounting for this difference can mistake sequencing depth for a biological difference. A library-size offset accounts for the expected scaling of counts with sequencing exposure while keeping the response as an integer count.
For gene \(g\) in pseudobulk sample \(j\), let \(Y_{gj}\) be the observed count and \(L_j\) the sample's library size, measured as the total relevant pseudobulk reads or counts. An effective library size may also incorporate a sample-specific normalization factor. The model then compares expression rates relative to that effective exposure, rather than treating every sample's raw count as equally exposed.
A common count-model choice is the negative binomial generalized linear model (NB GLM), with $$ Y_{gj}\sim \operatorname{NB}(\mu_{gj},\phi_g),\qquad \operatorname{Var}(Y_{gj}\mid \mathbf{x}_j)=\mu_{gj}+\phi_g\mu_{gj}^{2}, $$ where \(\phi_g\) is the dispersion for gene \(g\). The offset specifies the mean's exposure scaling; it does not change the fact that the response is modeled as a count. The variance assumption and the design matrix are part of the model too. The offset alone does not make samples comparable if the exposure measure is inappropriate or if other model assumptions fail.
Deriving the Offset Form
Suppose a gene has a rate \(r_{gj}\) per unit of effective library size in sample \(j\). The expected count is proportional to exposure: $$ \mu_{gj}=E_j r_{gj}. $$ If the log rate is modeled as \(\log(r_{gj})=\mathbf{x}_j^\mathsf{T}\boldsymbol{\beta}_g\), taking logarithms gives $$ \log(\mu_{gj})=\log(E_j)+\mathbf{x}_j^\mathsf{T}\boldsymbol{\beta}_g. $$ Thus the log-library-size term enters additively, and its coefficient is one because doubling exposure doubles the expected count when the rate is unchanged. Fixing that coefficient is the defining feature of the offset; estimating a separate coefficient for library size would instead allow the data to determine how strongly exposure scales the mean.
This derivation assumes that, conditional on modeled covariates, expected counts scale proportionally with effective library size. It also assumes that the library size and any normalization factor are available as fixed quantities for the analysis, rather than outcomes being modeled jointly with the gene count. If the linear predictor includes a condition indicator \(C_j\), for example, $$ \log(\mu_{gj})=\log(E_j)+\beta_{0g}+\beta_{1g}C_j, $$ then the modeled rate in the \(C_j=1\) group divided by the modeled rate in the \(C_j=0\) group is \(\exp(\beta_{1g})\), holding other covariates fixed. The offset handles exposure; it does not decide which covariates belong in the design.
Worked Example: Equal Rates at Different Sequencing Depths
Consider two hypothetical pseudobulk samples for one gene. Sample A has \(Y_A=100\) counts and \(L_A=10{,}000{,}000\) reads; sample B has \(Y_B=300\) counts and \(L_B=30{,}000{,}000\) reads. Suppose \(f_A=f_B=1\), so effective library size equals library size. Compare the observed counts per million reads.
For sample A, the observed rate per million is $$ \frac{10^6Y_A}{E_A} =\frac{10^6(100)}{10{,}000{,}000} =10.0\ \text{counts per million}. $$ For sample B, it is $$ \frac{10^6Y_B}{E_B} =\frac{10^6(300)}{30{,}000{,}000} =10.0\ \text{counts per million}. $$ The raw count in B is three times the count in A, but its library size is also three times as large. The observed rate ratio is \(10.0/10.0=1\), with log ratio \(\log(1)=0\). Under a model with a condition indicator separating these samples, this pattern is consistent with a condition log-rate difference of zero, not a threefold increase. With only one sample per group, this calculation is descriptive and does not provide a reliable estimate of between-sample variability or a condition test.
Using Effective Library Sizes
A raw total library size may not always be the exposure used by the model. Some workflows multiply \(L_j\) by a sample-specific factor \(f_j\), producing \(E_j=L_jf_j\). The factor can adjust the exposure scale when raw totals alone are not considered adequate. The precise procedure for choosing normalization factors is separate from the offset's role: once \(E_j\) is specified, the offset is \(\log(E_j)\). A factor above one increases effective exposure and therefore increases the expected count at a fixed rate.
For example, suppose sample A has \(L_A=8{,}000{,}000\) and \(f_A=1.25\), while sample B has \(L_B=12{,}000{,}000\) and \(f_B=0.90\). Then \(E_A=10{,}000{,}000\) and \(E_B=10{,}800{,}000\), not 8 million and 12 million. On the log scale, the offsets are \(\log(10{,}000{,}000)\) and \(\log(10{,}800{,}000)\). A common scaling of every effective library size changes the intercept's units but not a comparison of condition rates.
Worked Example: Applying an Effective-Library-Size Factor
Suppose two hypothetical samples have \(Y_A=100\) and \(Y_B=216\) counts for a gene. Their raw library sizes are 8 million and 12 million reads, and their specified factors are \(f_A=1.25\) and \(f_B=0.90\). Calculate effective library sizes and compare observed rates per million effective reads.
First calculate exposure: $$ E_A=(8{,}000{,}000)(1.25)=10{,}000{,}000,\qquad E_B=(12{,}000{,}000)(0.90)=10{,}800{,}000. $$ The rates are $$ \frac{10^6(100)}{10{,}000{,}000}=10.0,\qquad \frac{10^6(216)}{10{,}800{,}000}=20.0 $$ counts per million effective reads. The observed rate ratio is \(20.0/10.0=2.00\), and its log is \(\log(2.00)\approx0.693\). This is a descriptive comparison of these two samples. If the rate ratio were represented by a condition coefficient in a log-link model, a coefficient of about \(0.693\) would correspond to a twofold rate. It does not follow that an NB GLM fit to noisy data must return exactly this value.
Interpreting the Model and Checking Its Fit
Under the offset model, the expected count is the product of effective library size and modeled rate: $$ \mu_{gj}=E_j\exp(\mathbf{x}_j^\mathsf{T}\boldsymbol{\beta}_g). $$ For a simple two-condition model, \(\exp(\beta_{1g})\) is the modeled rate ratio for condition 1 versus condition 0, at the same values of other covariates. It is not a raw-count ratio unless the compared samples have equal effective library sizes. In particular, the offset does not itself establish a condition effect, supply a standard error, or account for donor structure; those depend on the design and the rest of the model.
For a specified rate, predicted counts scale directly with effective exposure. If a gene's modeled rate is 12 counts per million effective reads, a sample with 5 million effective reads has expected count \(12(5)=60\), while one with 15 million has expected count \(12(15)=180\). These are expected counts, not promises about observed counts. Count variability means individual observations can lie above or below their means.
Worked Example: Predicting Counts From a Rate and Offset
For a hypothetical gene, suppose the fitted rate for a particular covariate profile is 8 counts per million effective reads. Two pseudobulk samples with that profile have effective library sizes of 6 million and 9 million reads. Use the offset mean model to calculate their expected counts.
Writing exposure in millions, the rate is 8 counts per million. For the first sample, $$ \mu_1=(6\ \text{million reads}) \left(\frac{8\ \text{counts}}{1\ \text{million reads}}\right) =48\ \text{counts}. $$ For the second, $$ \mu_2=(9\ \text{million reads}) \left(\frac{8\ \text{counts}}{1\ \text{million reads}}\right) =72\ \text{counts}. $$ The ratio \(72/48=1.50\) matches the exposure ratio \(9/6=1.50\), because the modeled rates are equal. The calculation checks the offset's intended scaling; it does not claim that future observed counts will equal 48 and 72.
Diagnostics should ask whether the specified offset is behaving plausibly alongside the count model. Compare observed counts or model residuals with effective library size. Larger counts overall in deeper samples are expected; a systematic depth pattern in residuals after accounting for the offset can indicate that exposure scaling or other parts of the model are inadequate. Also inspect samples with unusually small or large effective libraries, factors, or gene counts. A single extreme sample can have substantial influence, especially when biological replication is limited.
Check that each offset corresponds to the pseudobulk sample represented by its count row, and that no sample is inadvertently assigned another sample's library size. A log is defined only for positive exposure, so zero or missing library sizes require investigation rather than inserting an arbitrary value. Do not use the target gene's observed count as its own library size; exposure should represent the sample's overall relevant library, not a quantity that mechanically adjusts away the gene-specific response.
Alternatives, Failure Modes, and Limits
| Approach | What it does | Important distinction |
|---|---|---|
| Log-library-size offset in an NB GLM | Models raw counts with a fixed exposure term in the mean. | Retains the count response and imposes proportional scaling with effective library size. |
| Compare raw counts | Compares observed counts without exposure adjustment. | Confounds rate differences with unequal library sizes when exposure differs. |
| Normalize counts first and use a different response model | Transforms or scales the response before modeling. | Not automatically equivalent to fitting the count likelihood with an offset; its error assumptions and uncertainty treatment differ. |
| Estimate a coefficient for log library size | Includes library size as an ordinary predictor with an estimated slope. | Does not impose the known proportional-exposure relationship and may estimate a slope affected by the biological design. |
A common mistake is to reason, “Sample B has three times as many counts, so the gene is expressed three times as strongly,” without comparing library sizes. In the first worked example, that conclusion is wrong: both samples have 10 counts per million, and the threefold count difference is exactly matched by a threefold exposure difference. The correction is to model or calculate rates using a justified effective library size.
An offset also has important limits. Library-size adjustment addresses sequencing exposure, not every source of between-sample variation. In particular, it does not turn relative counts into absolute transcript counts per cell, correct for changes in cell composition, or guarantee that a global shift in RNA abundance is identified. If nearly all genes change together, a relative normalization scale may not reveal that absolute shift. Nor does a library-size offset correct omitted covariates, dependence among repeated samples, or an unsuitable variance model. These issues require appropriate design and model checks rather than changing the offset arbitrarily.
The central assumption is that, at a fixed modeled expression rate, doubling effective library size doubles the expected count. If residual checks or technical knowledge make that assumption implausible, the response is not simply to add an estimated depth slope and treat the issue as solved. First check the count aggregation, the exposure calculation, the normalization factors, and sample identity; then assess whether the model's exposure interpretation is appropriate for the scientific estimand.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Two hypothetical pseudobulk samples have counts 60 and 180 and library sizes 6 million and 18 million. What are their rates per million reads, and what does the raw-count ratio alone suggest incorrectly?
- Write the mean model for a gene with effective library size \(E_j\) and a condition indicator \(C_j\). State the rate-ratio interpretation of the condition coefficient.
- A sample has \(L_j=5\) million reads and normalization factor \(f_j=1.2\). What is its effective library size, and what expected count follows from a modeled rate of 10 counts per million effective reads?
- Name one assumption addressed by the offset and two limitations it does not address.
Solutions. (1) The rates are \(10^6(60)/(6{,}000{,}000)=10\) and \(10^6(180)/(18{,}000{,}000)=10\) counts per million. The raw-count ratio is \(180/60=3\), which would misleadingly suggest a threefold rate difference; exposure also differs threefold. (2) The model is \(\log(\mu_{gj})=\log(E_j)+\beta_{0g}+\beta_{1g}C_j\). The modeled condition rate ratio is \(\exp(\beta_{1g})\), holding other covariates fixed. (3) \(E_j=(5\ \text{million})(1.2)=6\) million reads, so the expected count is \(6(10)=60\). (4) The offset assumes expected counts scale proportionally with effective library size at a fixed rate. It does not establish absolute transcript abundance per cell or correct for omitted covariates, cell-composition changes, or dependence among repeated samples.