Count models

limma models log counts; edgeR and DESeq2 model the counts themselves, with a negative binomial distribution whose variance grows with the mean. CytePool implements edgeR 4’s quasi-likelihood pipeline and DESeq2’s analysis, each checked against the R package on the same data.

Normalization¶

Counts depend on how deeply each sample was sequenced and on its composition: a few very highly expressed genes in one sample take reads from all others. Two corrections are standard. edgeR multiplies each library size by a normalization factor, from the trimmed mean of M values (tmm_factors; Robinson & Oshlack (2010)), the upper quartile (upper_quartile_factors; Bullard et al. (2010)) or the relative log expression (rle_factors; Anders & Huber (2010)). DESeq2 replaces library sizes by size factors, the median of each sample’s ratios to the geometric mean of every gene (median_ratio_size_factors; Anders & Huber (2010)). cpm and log_cpm give counts per million, and filter_by_expression keeps genes with enough counts in enough samples for the design Chen et al., 2016.

The negative binomial model¶

For the count $y_{gi}$ of gene $g$ in sample $i$,

with $\phi$ the dispersion and the offset $o_{gi}$ usually the log of the effective library size McCarthy et al., 2012. $\sqrt{\phi}$ is the biological coefficient of variation. The coefficients minimize the deviance

which fit_negative_binomial_model finds for every gene, as edgeR’s glmFit does (NegativeBinomialFit), and which negative_binomial_deviance and negative_binomial_unit_deviance compute exactly Dunn & Smyth, 2018.

edgeR’s quasi-likelihood F-tests¶

The quasi-likelihood model gives each gene a second, gene-wise dispersion $\sigma_g^2$ on top of a negative binomial dispersion $\phi$ shared by all genes Lund et al., 2012:

fit_quasi_likelihood follows edgeR 4’s defaults Chen et al., 2025:

1. $\phi$ is the Cox-Reid common dispersion of the most highly expressed genes (common_dispersion): the maximum over $\phi$ of the summed adjusted profile log-likelihoods $\ell_g(\phi) - \tfrac{1}{2}\log\det(X' W_g X)$ Cox & Reid, 1987McCarthy et al., 2012.

2. Each gene’s quasi-dispersion is its deviance over its residual degrees of freedom, both adjusted for small counts: each unit deviance is approximated by a scaled chi-squared variable with the exact mean and variance of the unit deviance (negative_binomial_unit_deviance_moments, which gives UnitDevianceMoments), where edgeR approximates them.

3. The quasi-dispersions are moderated by empirical Bayes towards a trend in the average log CPM (average_log_cpm), with limma’s estimator for unequal degrees of freedom (Empirical Bayes moderation).

The fit is a QuasiLikelihoodFit. quasi_likelihood_f_test fits the model without the tested coefficients or contrasts, and divides the likelihood ratio statistic by its degrees of freedom and by the posterior quasi-dispersion; its QuasiLikelihoodTest has edgeR’s columns: log fold change, log CPM, F, p-value and FDR. The edgeR User’s Guide (section 1.4) recommends these F-tests for bulk RNA-seq because they account for the uncertainty of the dispersion estimates, which gives stricter error rate control.

DESeq2¶

DESeq2 models the same distribution on the log2 scale with size factors $s_i$, $\log_2 \mu_{gi} = \log_2 s_i + x_i' \beta_g$, and shrinks the dispersions in four steps Love et al., 2014, which deseq2_dispersions reproduces (DESeq2Dispersions):

1. gene-wise estimates, the maxima of the Cox-Reid adjusted likelihood;

2. a trend in the mean normalized count $m$, $\alpha_{tr}(m) = a_0 + a_1 / m$, fitted by a gamma regression;

3. a normal prior on the log dispersions around the trend, with its variance from the spread of the gene-wise estimates;

4. maximum a posteriori estimates, keeping the gene-wise estimate of genes far above the trend (dispersion outliers).

fit_deseq2 fits every gene at those dispersions (DESeq2Fit). deseq2_wald_test tests a coefficient or contrast $c$ by $c'\hat\beta / \sqrt{c'\Sigma c}$ against the normal distribution; deseq2_likelihood_ratio_test tests a reduced design by the difference of the deviances against a chi-squared distribution, which tests several coefficients at once, such as every level of a factor or every time point. Both report DESeq2’s results() columns (DESeq2Results), filter genes with low mean counts before adjusting the p-values (independent filtering; Bourgon et al. (2010)) and remove the p-values of genes with an influential outlier count (Cook’s distance; Cook (1977)).

What differs from the R packages¶

CytePool’s fits converge to the maximum where edgeR’s and DESeq2’s stop at a tolerance, and its deviances and unit deviance moments are exact where edgeR approximates them; the differences are reported, with their sizes, on the Validation page. DESeq2’s replacement of outlier counts, shrunken fold changes (lfcShrink) and edgeR’s legacy = TRUE methods are not implemented.