Differential expression

limma_voom() and limma_trend() test the contrasts of a Design in every group (cell type) of a cube, one group at a time, as the limma and edgeR User’s Guides do for RNA-seq counts, and dream_voom() does so with the random intercepts of the design, such as a donor’s, as variancePartition’s dream does. The design is matched to the cube by sample_id; samples with too few cells are left out of a group’s fit, and a group whose design becomes rank deficient is fitted on its column space, with the contrasts its samples can estimate. Options for pseudobulk data, each off by default, add sample quality weights, prior weights from the number of cells, a minimum total count per sample, and p-values adjusted across groups. The result holds a table per group and contrast, a summary of every group, the samples used and left out with the weights of each fit, and a record that traces every number to the profiles of the cube that produced it.

Differential expression between samples, one group (cell type) of a cube at a time.

limma_voom() and limma_trend() run the limma workflows for RNA-seq counts (Ritchie et al. 2015) on the pseudobulk counts of each group of a cube: limma-voom, which gives every observation a precision weight from the mean-variance trend of the log counts per million (Law et al. 2014), and limma-trend, which lets the prior variance of the empirical Bayes moderation follow that trend instead. dream_voom() runs variancePartition’s dream (Hoffman and Roussos 2021), limma-voom with a linear mixed model in place of the linear model, for designs with random effects such as several samples of each donor. Each tests the named contrasts of a Design and returns, for every group and contrast, a table of the tested features with topTable’s statistics, together with a summary, the samples behind every fit and a provenance record.

Samples:

The design’s sample table must have a sample_id column. Its samples are matched to the cube’s by sample_id, never by position, and they may be a subset of the cube’s samples; the order of the design’s rows does not matter, as every group is fitted with its samples in the cube’s order. In each group, a sample is left out of the fit when it has fewer than min_cells cells (when the cube records cell numbers), no counts, or fewer than min_counts counts in all (its total over every feature of the cube), and the result records why. min_counts is 0 by default; the scverse single-cell best practices (Heumos et al. 2023) keep samples with at least 10 cells and 1,000 counts (min_cells=10, min_counts=1000), as decoupler’s filter_samples does, which keeps a sample with exactly as many (found by probing).

The design of a group:

The design matrix of a group is the rows of the full design matrix for the samples left in it, so every column keeps its meaning: factor levels, reference levels and covariate values are those of the full design. With fewer samples that matrix can be rank deficient, for example when every sample of one donor is left out. Such a group is fitted on the column space of its design matrix, keeping the columns that R’s qr() keeps, as limma’s lmFit does; each other column (aliased) equals a combination of the kept columns before it. A contrast c is estimable from the group’s samples when it is orthogonal to the null space of the group’s design matrix (Searle 1971): when its weight on every aliased column is the sum of its weights on the kept columns, weighted by that column’s combination. This is the rule of Design.check. c applied to the coefficients then equals its weights on the kept columns applied to the fit on the kept columns, whatever the coefficients that the samples cannot separate, so its estimate and its exact standard error are those of any full-rank parametrisation of the same model. A contrast that is not estimable is reported, with the relations it violates, and not computed.

Groups that are not fitted:

A group is not fitted, and its status in the summary says why, when no sample is left, the design leaves no residual degrees of freedom, no contrast is estimable, for dream the random effects cannot be estimated from the group’s samples, no feature passes the expression filter (or, for voom, fewer than two), a sample has no counts in the features that pass, or voom, the sample quality weights or the moderation fail on its data. A run never stops for one group.

Steps:

For each group, with the R functions they match in the limma and edgeR User’s Guides:

1. Leave out samples, as above; the design matrix is the rows of the rest.

2. Filter features by expression with the group’s design, by edgeR’s filterByExpr(y, design) (Chen, Lun and Smyth 2016) with its defaults (expression_filter), and recompute the library sizes over the features kept (y[keep, , keep.lib.sizes = FALSE]).

3. Normalization factors by TMM (Robinson and Oshlack 2010), edgeR’s calcNormFactors (normalization).

4. limma-voom: voom on the library sizes times the normalization factors, with its LOWESS span chosen from the number of features tested, as limma’s voom chooses it by default since limma 3.66 (span), then lmFit with its weights. limma-trend: log2 counts per million with a prior count of 3, cpm(y, log = TRUE, prior.count = 3) as the limma User’s Guide uses for limma-trend, then lmFit. Sample weights, when asked, enter here (below).

5. contrasts.fit, eBayes (Smyth 2004), with trend = TRUE for limma-trend and robust (Phipson et al. 2016) when asked, and topTable(sort.by = "none"): the Benjamini-Hochberg adjustment (Benjamini and Hochberg 1995) runs over the tested features of one group and contrast, and with adjust_across_groups also over every group (below).

Sample weights:

Both kinds are off by default, and they are alternatives: limma’s voomWithQualityWeights takes no prior weights, and edgeR’s voomLmFit refuses prior weights together with estimated sample weights (found by probing), so CytePool refuses the two together too. The sample table and the record of each fitted group hold the sample weights of its fit. With limma-voom, the quality weights recorded are the final ones; voomWithQualityWeights estimates them twice, and its first estimate only weights the fit of voom’s second trend. voom_with_quality_weights(), given the group’s counts, design rows and normalization factors, reproduces both.

• quality_weights: sample quality weights, estimated in each group from how closely each sample’s values follow the linear model (Ritchie et al. 2006; Liu et al. 2015), which a limma author recommends for pseudobulk data (Smyth 2025, on the Bioconductor support site). limma-voom runs limma’s voomWithQualityWeights (voom_with_quality_weights()), gene by gene by default, as limma does. limma-trend estimates them from the log2 counts per million as limma’s arrayWeights does (sample_quality_weights()), by REML by default, the method that arrayWeights chooses for such data, and passes them to lmFit as sample weights. Where a donor has two samples, as in ~ donor + condition, their residuals are equal and opposite (up to voom’s weights), so both get the low weight of an outlying sample. Where the weights cannot be estimated, as with fewer than two residual degrees of freedom, they are one, as limma’s are, with a warning that names the group and that its record keeps; where the estimate fails, as when a sample carries no information about its variance and prior_n is 0, the group is not fitted.

• prior_weights="n_cells": each sample’s number of cells in the group, as prior weights. A pseudobulk profile sums the counts of its cells, so its log2 counts per million estimate the mean expression of those cells, whose variance falls as one over their number: the precision of a profile is proportional to its number of cells, as the dreamlet paper argues (Hoffman et al. 2026, Methods), and dreamlet’s processAssays weights its fits by the number of cells of each sample, as its documentation describes. limma-voom passes them to voom as weights: they weight the fits behind voom’s trend, and the precision weights are the trend’s times them, as limma’s help page says voom’s weights modify prior weights, and as edgeR’s voomLmFit(prior.weights = ...) forms them (found by probing). limma-trend passes them to lmFit as sample weights. Only their ratios within a group matter: weights that are all ten times larger give the same results.

dream:

dream_voom() takes the steps of limma-voom up to the normalization, then those of variancePartition’s voomWithDreamWeights, dream, eBayes and topTable: precision weights from the mean-variance trend of the design’s mixed model, fitted to every feature by maximum likelihood (voom_with_dream_weights()); the mixed model of every feature fitted by REML with those weights (fit_mixed_model()); t-tests of the estimable contrasts (fixed_effect_tests()) with the degrees of freedom dream chooses by default (ddf = "adaptive"), Kenward and Roger’s (1997) in a group of fewer than 20 samples and Satterthwaite’s (1946) otherwise; F-tests of sets of contrasts; and their empirical Bayes moderation (moderated_fixed_effect_tests()). The random effects of a group are those of its samples: a grouping factor’s levels without samples in the group drop out. cytepool.stats.dream documents the conventions found by probing variancePartition and the deliberate differences from it. Prior weights from the number of cells are divided by their mean in each group, as variancePartition divides prior weights, which changes no test.

P-values adjusted across groups:

adjusted_p_value adjusts the p-values of one group and contrast, as limma’s topTable does. With adjust_across_groups, the column adjusted_p_value_across_groups adjusts each contrast’s p-values by Benjamini-Hochberg over every test of that contrast in the analysis: every feature tested in every group that was fitted and can estimate the contrast, stacked into one family. Each contrast is a family of its own, and groups left out of groups are not in it. muscat reports both adjustments (Crowell et al. 2020), whose comparison found the adjustment over all groups very conservative. It is not larger for every test: a family with more small p-values can lower an adjusted value.

Deliberate differences from limma:

validation/differential_expression.py reports each of them.

• In a rank-deficient group, a contrast that is estimable but weighs an aliased column is computed. limma’s contrasts.fit stops with an error for such a contrast, and with it for every contrast passed together with it (found by probing), so in limma it needs a full-rank parametrisation.

• A contrast that is not estimable, but weighs no aliased column, is reported as not estimable. limma computes it from the kept columns, whose meaning has then changed: with every sample of the reference donor D1 left out of ~ donor + condition, limma reports for the coefficient donorD2 (donor D2 against D1) the difference between donor D2 and the last donor.

• With contrast_standard_errors="exact", the default of limma_voom(), the standard errors of contrasts that combine several coefficients are exact; limma’s contrasts.fit approximates them when weights differ between features, and "limma" reproduces it. limma’s approximation, unlike the exact value, depends on the parametrisation, so only exact standard errors are the same in every parametrisation of a group’s design.

Memory:

One group is in memory at a time: its counts layer, read from the cube when the group is fitted, and the fit of its features. A cube opened with read_cube() therefore stays on disk except for the group being fitted. Results hold one table per group and contrast, and a record that identifies the cube by its fingerprint, which costs no read for a cube file (“Provenance” in limma_voom()).

References

Law CW, Chen Y, Shi W, Smyth GK (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15, R29. doi:10.1186/gb-2014-15-2-r29

Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43, e47. doi:10.1093/nar/gkv007

Chen Y, Lun ATL, Smyth GK (2016). From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research 5, 1438. doi:10.12688/f1000research.8987.2

Phipson B, Lee S, Majewski IJ, Alexander WS, Smyth GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Annals of Applied Statistics 10, 946-963. doi:10.1214/16-AOAS920

Searle SR (1971). Linear Models. Wiley, New York.

Ritchie ME, Diyagama D, Neilson J, van Laar R, Dobrovic A, Holloway A, Smyth GK (2006). Empirical array quality weights in the analysis of microarray data. BMC Bioinformatics 7, 261. doi:10.1186/1471-2105-7-261

Liu R, Holik AZ, Su S, Jansz N, Chen K, Leong HS, Blewitt ME, Asselin-Labat ML, Smyth GK, Ritchie ME (2015). Why weight? Modelling sample and observational level variability improves power in RNA-seq analyses. Nucleic Acids Research 43, e97. doi:10.1093/nar/gkv412

Hoffman GE, Lee D, Bendl J, Prashant NM, Hong A, Casey C, Alvia M, Shao Z, Argyriou S, Therrien K, Venkatesh S, Voloudakis G, Haroutunian V, Fullard JF, Roussos P (2026). Efficient differential expression analysis of large-scale single-cell transcriptomics data using Dreamlet. Nature Communications 17, 9597. doi:10.1038/s41467-026-75680-8

Crowell HL, Soneson C, Germain PL, Calini D, Collin L, Raposo C, Malhotra D, Robinson MD (2020). muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nature Communications 11, 6077. doi:10.1038/s41467-020-19894-4

Hoffman GE, Roussos P (2021). Dream: powerful differential expression analysis for repeated measures designs. Bioinformatics 37, 192-201. doi:10.1093/bioinformatics/btaa687

Smyth GK (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, article 3. doi:10.2202/1544-6115.1027

Robinson MD, Oshlack A (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology 11, R25. doi:10.1186/gb-2010-11-3-r25

Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57, 289-300. doi:10.1111/j.2517-6161.1995.tb02031.x

Satterthwaite FE (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin 2, 110-114. doi:10.2307/3002019

Kenward MG, Roger JH (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics 53, 983-997. doi:10.2307/2533558

Heumos L, Schaar AC, Lance C, Litinetskaya A, Drost F, Zappia L, Luecken MD, Strobl DC, Henao J, Curion F, Single-cell Best Practices Consortium, Schiller HB, Theis FJ (2023). Best practices for single-cell analysis across modalities. Nature Reviews Genetics 24, 550-572. doi:10.1038/s41576-023-00586-w

Smyth GK (2025). Answer to “Guidance on limma for pseudobulk scRNA-seq with multi-factor design”, 11 September 2025. Bioconductor support site, https://support.bioconductor.org/p/9162569/

class cytepool.differential_expression.DifferentialExpressionResult(method, contrasts, summary, samples, tables, provenance)

Differential expression of every group of a cube, from limma or dream.

Tables hold, for every group and estimable contrast, one row per tested feature, in the cube’s feature order (topTable(sort.by = "none")):

• feature_id and, when the cube’s feature table has one, symbol;

• log2_fold_change: the contrast’s estimate (limma’s logFC);

• standard_error: its moderated standard error, the contrast’s unscaled standard deviation times the feature’s posterior standard deviation (limma’s stdev.unscaled * sqrt(s2.post)), so that t is log2_fold_change / standard_error; it is what a meta-analysis of studies reads (differential_expression_meta_analysis());

• average_log2_cpm: the mean log2 count per million over the group’s samples (AveExpr);

• t and p_value: the moderated t-statistic and its two-sided p-value;

• adjusted_p_value: Benjamini-Hochberg over the group’s tested features (adj.P.Val);

• adjusted_p_value_across_groups, with adjust_across_groups only: Benjamini-Hochberg over the tested features of every group that tests the contrast;

• b_statistic: the log-odds that the feature is differentially expressed (B).

Tables from dream_voom() also give df, the degrees of freedom of each test (Satterthwaite’s or Kenward and Roger’s plus the prior’s, df.total), after t, and z_statistic, the z-statistic with the test’s p-value (topTable’s z.std), after b_statistic; a feature whose mixed model has no minimum has NaN statistics. Its F-tests have tables of their own, whose columns dream_voom() lists.

The summary has one row per group: group_id; status, "fitted" or why the group was not fitted ("no_samples", "no_residual_df", "no_estimable_contrast", "random_effects_not_estimable", "no_features", "too_few_features", "empty_library" or "failed"); reason, the explanation (empty when fitted); n_samples_used and n_samples_left_out; n_features_tested; df_residual, for dream the median over the features of their approximate residual degrees of freedom (in a group not fitted, the residual degrees of freedom of the fixed effects, as for limma); prior_df, the prior degrees of freedom (their median over features with robust moderation, where each feature has its own); for dream, df_method, the degrees of freedom of the group’s tests ("satterthwaite" or "kenward_roger"; empty where the group was not fitted); and not_estimable, the contrasts (and F-tests) the group’s samples cannot estimate, separated by ", ". Numbers are NaN, or 0, where a group was not fitted that far.

The sample table has one row per group and design sample, in the cube’s sample order: group_id, sample_id, n_cells (when the cube records cell numbers), used (in the group’s design), reason (why a sample was left out; empty when used), and the library_size (over the tested features) and norm_factor that entered the fit; with prior_weights, the prior_weight of each sample, and with quality_weights, its quality_weight. Numbers that did not enter a fit are NaN: those of samples left out, and of groups not fitted.

Built by limma_voom(), limma_trend() and dream_voom(), from the contrasts, tables and record they computed. Immutable.

Raises:

• TypeError – If method is not a string, contrasts not a Contrasts, summary or samples not a Table, tables not a mapping of (group id, contrast name) to Table, or provenance not a mapping that JSON can write.

• ValueError – If method is not "limma_voom", "limma_trend" or "dream_voom", or provenance holds a NaN or an infinity.

Explained in: Choosing a method

property method

"limma_voom", "limma_trend" or "dream_voom".

property groups

The ids of the groups analysed, in order.

property contrasts

The contrasts tested, over the columns of the full design.

property summary

The status, samples, features and degrees of freedom of every group.

property samples

Every design sample of every group, used in its fit or left out, and why.

property tables

The table of every group and contrast, or F-test, that was tested, by (group, contrast) or (group, F-test).

property f_tests

The F-tests of dream_voom(), each with the contrasts it tests at once, as a new dict; empty for the other methods or without F-tests.

property settings

The settings of the analysis, as a new dict (also in the provenance).

property provenance

A JSON-serialisable record of what produced every number, as a new dict.

It holds the method; the CytePool, numpy, numba and Python versions; the cube ("cube"): its fingerprint, the digests it is built from and its own provenance; the layer read; the design (formula, columns, levels, the values of its variables for every sample, and the weights of every contrast; for dream also its random effects, "random", and the values of the variables they group by); the settings; and, for every group, its status, the samples that entered its design, in the order they were fitted, those left out with the reason, the design columns its samples cannot separate, the contrasts it cannot estimate with the reason, the number of features tested, for limma-voom and dream the span of its trend (voom_span, null where voom did not run), for dream the number of features whose mixed model has no minimum (features_not_fitted) and of those whose REML fit stopped without converging (features_not_converged; their tests use the estimates where it stopped), the degrees of freedom of its tests (df_method), the number of features whose Kenward-Roger tests fell back to Satterthwaite’s (features_satterthwaite_fallback) and the numerator degrees of freedom of each F-test (f_test_df1), the sample weights its fit used, in the order of its samples (prior_weights with prior_weights, quality_weights with quality_weights; null where the group was not fitted), and the warnings that the estimation of its sample quality weights gave (warnings, such as weights that could not be estimated and are one). Each number of a table comes from the profiles of its group and of those samples in the cube with that fingerprint: the counts of every feature of these profiles enter its filter, normalization and trend. The fingerprint, and the digest of the group’s counts, say whether a cube at hand holds those profiles (Cube.verify_fingerprint).

table(group, contrast)

Return the table of one group and contrast, or F-test.

Parameters:

• group – A group id.

• contrast – A contrast name, or the name of an F-test of dream_voom().

Raises:

KeyError – If the group or contrast was not analysed, or the group was not fitted or cannot estimate the contrast; the message says why.

Explained in: Choosing a method

cytepool.differential_expression.dream_voom(cube, design, *, groups=None, contrasts=None, f_tests=None, min_cells=10, min_counts=0, expression_filter=True, normalization='tmm', span=None, prior_weights=None, df_method='adaptive', kenward_roger_weights='average', robust=False, trend=False, adjust_across_groups=False, n_threads=None)

Test contrasts in every group of a cube by dream, with random effects.

variancePartition’s dream workflow for repeated measures (Hoffman and Roussos 2021), such as several samples of each donor: as limma_voom(), but the design’s random intercepts, such as (1 | donor), enter every fit. For each group, the samples are left out, the features filtered and the counts normalized as for limma_voom; voom’s precision weights come from the mean-variance trend of mixed models fitted by maximum likelihood (voomWithDreamWeights, voom_with_dream_weights()); every feature’s mixed model is fitted by REML with those weights (dream, fit_mixed_model()); its estimable contrasts are tested (fixed_effect_tests()) with the degrees of freedom of df_method; and the tests are moderated by empirical Bayes as variancePartition’s eBayes moderates them (moderated_fixed_effect_tests()). The module documentation gives the steps and the differences from variancePartition.

Every group’s mixed model uses the random effects of the samples left in it: the levels of a grouping factor without samples there drop out, and a group whose samples cannot estimate a random effect, such as one donor left, or one sample per donor, is not fitted (status "random_effects_not_estimable"). A feature whose mixed model has no minimum, which the fixed and random effects reproduce, is tested with NaN statistics; its group’s record counts them.

Degrees of freedom:

As dream’s ddf, by default ("adaptive") Kenward and Roger’s method in a group with fewer than 20 samples and Satterthwaite’s otherwise. variancePartition documents “20 samples or fewer”; with 20 its dream uses Satterthwaite’s (found by probing variancePartition 1.36.2). By default (kenward_roger_weights="average"), Kenward and Roger’s adjustment is computed as dream computes it with pbkrtest, with every observation at its feature’s average weight; with voom’s weights that adjustment can make standard errors smaller than the unadjusted ones, which "precision", the method for the weighted model, never does (see fixed_effect_tests()). As dream does, a feature whose adjusted covariance matrix is not positive definite or has a condition number above 1e6 is tested by Satterthwaite’s method instead; the group’s record counts such features. That number depends on the units of continuous covariates: centring and scaling them avoids it.

F-tests:

f_tests names sets of contrasts to test at once, such as every difference among three time points: each is tested by dream’s moderated F-statistic, which topTable gives for several coefficients (see cytepool.stats.dream), in every group that can estimate all its contrasts. Its table has, for every tested feature, feature_id (and symbol), log2_fold_change_<contrast> for each contrast, average_log2_cpm, f_statistic, df (the denominator degrees of freedom), p_value, adjusted_p_value (and adjusted_p_value_across_groups), and standardized_f_statistic, the F-statistic on infinite denominator degrees of freedom with the same p-value (F.std). The numerator degrees of freedom are the number of linearly independent contrasts, which the group’s record gives.

Parameters:

• cube – A cube with a counts layer.

• design – A Design whose random effects are random intercepts, such as random="(1 | donor)", matched to the cube by sample_id as for limma_voom().

• groups – A group id or a list of them; all groups by default.

• contrasts – The contrasts to test, matched to the design’s columns by name; by default the design’s named contrasts.

• f_tests – A mapping of an F-test’s name to the names of the contrasts it tests at once, among those tested, such as {"time": ["t1_vs_t0", "t2_vs_t0"]}; None for none. A name must differ from every contrast’s.

• min_cells – The fewest cells a sample needs in a group to enter its fit; ignored for cubes without cell numbers.

• min_counts – The fewest counts a sample needs in a group, over every feature of the cube, to enter its fit; 0 by default.

• expression_filter – True for edgeR’s filterByExpr defaults with the group’s fixed effects, a mapping of changed settings, or False.

• normalization – "tmm" for TMM normalization factors, or None.

• span – The LOWESS span of the trend, in (0, 1], for every group; None chooses it in each group from the number of features tested, as limma_voom does. variancePartition’s voomWithDreamWeights uses 0.5.

• prior_weights – "n_cells" to give each sample its number of cells in the group as a prior weight, as dreamlet weights samples, or None for none. Needs a cube that records cell numbers.

• df_method – "adaptive" (dream’s default, above), "satterthwaite" or "kenward_roger" (with dream’s fallback, above) for every group.

• kenward_roger_weights – "average", as dream computes Kenward and Roger’s adjustment, or "precision", with each observation’s own weight; see above.

• robust – Whether to estimate the prior of the moderation robustly, as eBayes(robust = TRUE) does.

• trend – Whether the prior variance of the moderation follows the mean log2 CPM of the features, as eBayes(trend = TRUE) does. dream moderates the residual variance of an observation of average weight, which keeps much of the mean-variance trend that voom’s weights remove (see cytepool.stats.dream); variancePartition’s default is False.

• adjust_across_groups – Whether to add the column adjusted_p_value_across_groups, as for limma_voom().

• n_threads – The most threads for the mixed model fits of the features (see fit_mixed_model()), at least 1; None keeps numba’s setting. The results do not depend on it.

Returns:

The tables of every group and estimable contrast, and of every F-test, a summary of every group, the samples used and left out, the settings and a provenance record, as for limma_voom(); the tables of contrasts also give each test’s degrees of freedom (df) and z-statistic (z_statistic), and the summary the degrees of freedom of each group’s tests (df_method).

Raises:

• TypeError – If cube is not a Cube, design not a Design, an option has the wrong type, or contrasts, f_tests, groups, min_cells, expression_filter or n_threads is or holds a numpy masked array.

• ValueError – If the design has no random effects or one that is not a random intercept, has no sample_id column, names samples the cube does not have or names a sample twice, if there are no contrasts to test, an option is invalid, an F-test names a contrast that is not tested or shares a contrast’s name, or prior_weights is "n_cells" and the cube does not record cell numbers.

• KeyError – If a group is not in the cube.

Explained in: Choosing a method, Mixed models, REML and dream, Multiple testing

cytepool.differential_expression.limma_trend(cube, design, *, groups=None, contrasts=None, min_cells=10, min_counts=0, expression_filter=True, normalization='tmm', prior_count=3.0, prior_weights=None, quality_weights=False, robust=False, adjust_across_groups=False)

Test contrasts in every group of a cube by limma-trend.

As limma_voom(), but the counts become log2 counts per million with a prior count (edgeR’s cpm(y, log = TRUE, prior.count = 3), the limma User’s Guide’s choice for limma-trend), every feature is fitted without observation weights, and the prior variance of the empirical Bayes moderation follows the mean log2 CPM of the features (eBayes(trend = TRUE)). The standard errors of contrasts are exact, as limma’s are, also with sample weights. The result records the cube’s fingerprint and digests, at the cost “Provenance” in limma_voom() gives.

The options for pseudobulk data are those of limma_voom(); sample weights (prior or quality weights) enter lmFit as sample weights.

Parameters:

• cube – A cube with a counts layer.

• design – A Design matched to the cube by sample_id, as for limma_voom().

• groups – A group id or a list of them; all groups by default.

• contrasts – The contrasts to test; by default the design’s named contrasts.

• min_cells – The fewest cells a sample needs in a group to enter its fit; ignored for cubes without cell numbers.

• min_counts – The fewest counts a sample needs in a group, over every feature of the cube, to enter its fit; 0 by default.

• expression_filter – True for edgeR’s filterByExpr defaults with the group’s design, a mapping of changed settings, or False.

• normalization – "tmm" for TMM normalization factors, or None.

• prior_count – The average count added to every observation before taking logs, scaled by library size as edgeR does; positive. 3 by default, as in the limma User’s Guide.

• prior_weights – "n_cells" to weight each sample by its number of cells in the group, or None for none. Needs a cube that records cell numbers.

• quality_weights – True to estimate sample quality weights from the log2 counts per million of each group, as limma’s arrayWeights does, with its settings for such data (method="reml", prior_n=10, max_iterations=50, tolerance=1e-5; see sample_quality_weights()); a mapping with some of these settings to change them, such as {"method": "gene_by_gene"}; or False (or None) for none. Give at most one of prior_weights and quality_weights.

• robust – Whether to estimate the prior of the moderation robustly.

• adjust_across_groups – Whether to add the column adjusted_p_value_across_groups, as for limma_voom().

Returns:

The tables of every group and estimable contrast, a summary of every group, the samples used and left out, the settings and a provenance record.

Warns:

RuntimeWarning – As for limma_voom().

Raises:

• TypeError – If cube is not a Cube, design not a Design, an option has the wrong type, or contrasts, groups, min_cells, expression_filter or quality_weights is or holds a numpy masked array.

• ValueError – As for limma_voom(), or if prior_count is not positive.

• KeyError – If a group is not in the cube.

Explained in: Choosing a method, Meta-analysis, Multiple testing

cytepool.differential_expression.limma_voom(cube, design, *, groups=None, contrasts=None, min_cells=10, min_counts=0, expression_filter=True, normalization='tmm', span=None, prior_weights=None, quality_weights=False, robust=False, contrast_standard_errors='exact', adjust_across_groups=False, n_threads=None)

Test contrasts in every group of a cube by limma-voom.

For each group, the samples with fewer than min_cells cells or no counts are left out, features are filtered by expression, the counts are normalized, voom gives every observation a precision weight from the mean-variance trend of the log2 counts per million, a weighted linear model is fitted to every feature, and its contrasts are moderated by empirical Bayes: the steps of the limma User’s Guide for RNA-seq (voom, lmFit, contrasts.fit, eBayes, topTable) after edgeR’s filterByExpr and calcNormFactors. The module documentation describes the steps, rank-deficient groups and the deliberate differences from limma.

voom smooths its mean-variance trend with a LOWESS span that, by default, it chooses in each group from the number of features tested there, by limma’s chooseLowessSpan rule (choose_lowess_span()), as limma’s voom does by default since limma 3.66 (adaptive.span = TRUE). span fixes it instead; limma’s default before 3.66 was 0.5. The record of each group gives the span its trend used, as voom_span.

Options for pseudobulk data, each off by default so that the results are limma’s standard workflow unless asked: a minimum total count per sample (min_counts), sample quality weights (quality_weights, limma’s voomWithQualityWeights), prior weights from the number of cells behind each profile (prior_weights="n_cells", as dreamlet weights samples), and p-values adjusted across groups as well as within each (adjust_across_groups). The module documentation gives the rules and their sources.

Provenance:

The result records the cube’s fingerprint, the digests it is built from and the cube’s provenance (Cube.fingerprint, Cube.digests, Cube.provenance), computed after the fits. For a cube opened with read_cube() this reads nothing: the file stores the digests, and each group the fits read whole is checked against its digest. A cube in memory hashes its layers once, at about 1 s per GB of values where the processor has SHA instructions. A lazy cube derived from a cube file with select, regroup or collapse_samples reads, once, the groups whose values its operations changed and that the fits did not read whole: those of its other layers, and those of groups not analysed. A lazy layer of another array-like, which can change, is read whole once more.

Parameters:

• cube – A cube with a counts layer, such as one from from_cells() or read_cube().

• design – A Design whose sample table has a sample_id column with cube sample ids; its rows may be some of the cube’s samples, in any order. A design with random effects is refused: these fits have fixed effects only, and dream_voom() fits random intercepts.

• groups – A group id or a list of them; all groups by default.

• contrasts – The contrasts to test, as accepted by make_contrasts() and matched to the design’s columns by name, such as {"stim_vs_ctrl": "conditionstim"}. By default, the design’s named contrasts.

• min_cells – The fewest cells a sample needs in a group to enter its fit. Ignored when the cube does not record cell numbers (bulk data), where only samples without counts are left out.

• min_counts – The fewest counts a sample needs in a group, summed over every feature of the cube (its library size before the expression filter), to enter its fit; at least 0. The default, 0, leaves out only samples without counts.

• expression_filter – True for edgeR’s filterByExpr with its defaults (min_count=10, min_total_count=15, large_n=10, min_prop=0.7), using the group’s design, as the edgeR and limma User’s Guides do; a mapping with some of these settings to change them (see filter_by_expression()); or False to test every feature.

• normalization – "tmm" for edgeR’s TMM normalization factors, the default of calcNormFactors, or None for library sizes alone.

• span – The LOWESS span of voom’s trend, in (0, 1], for every group; None chooses it in each group from the number of features tested, as voom() does.

• prior_weights – "n_cells" to give each sample its number of cells in the group as a prior weight, which voom’s precision weights multiply, or None for none. Needs a cube that records cell numbers.

• quality_weights – True to estimate sample quality weights in each group and multiply them into voom’s precision weights, as limma’s voomWithQualityWeights does, with its settings (method="gene_by_gene", prior_n=10, max_iterations=50, tolerance=1e-5; see voom_with_quality_weights()); a mapping with some of these settings to change them, such as {"method": "reml"}; or False (or None) for none. Give at most one of prior_weights and quality_weights.

• robust – Whether to estimate the prior of the moderation robustly, as eBayes(robust = TRUE) does, which protects against hypervariable features.

• contrast_standard_errors – "exact" computes the standard error of every contrast from each feature’s own weighted fit; "limma" reproduces the approximation of limma’s contrasts.fit. They differ only for contrasts that combine correlated coefficients (see contrasts_fit()).

• adjust_across_groups – Whether to add the column adjusted_p_value_across_groups to every table: each contrast’s p-values adjusted by Benjamini-Hochberg over every group that tests it, besides adjusted_p_value, which adjusts them within the group.

• n_threads – The most threads for the weighted least-squares fits of the features, which run in parallel (numba’s threads; see numba.set_num_threads and fit_linear_model()), at least 1; None keeps numba’s setting. The results do not depend on it. With 1, no parallel region runs, so that forked processes and several Python threads can fit at once.

Returns:

The tables of every group and estimable contrast, a summary of every group, the samples used and left out, the settings and a provenance record.

Warns:

RuntimeWarning – If sample quality weights cannot be estimated in a group, such as with fewer than two residual degrees of freedom (they are then one), or their REML stops without converging; the message names the group, and its record keeps it.

Raises:

• TypeError – If cube is not a Cube, design not a Design, an option has the wrong type, or contrasts, groups, min_cells, expression_filter, quality_weights or n_threads is or holds a numpy masked array.

• ValueError – If the design has random effects or no sample_id column, names samples the cube does not have or names a sample twice, if there are no contrasts to test, an option is invalid, prior_weights and quality_weights are both given, or prior_weights is "n_cells" and the cube does not record cell numbers.

• KeyError – If a group is not in the cube.

Used in: Differential expression per cell type, Bulk RNA-seq

Explained in: Choosing a method, Designs, linear models and contrasts, Precision weights, Meta-analysis, Multiple testing, Cube files and provenance records