Which genes does a treatment change, and in which cell types? This page answers it for
the interferon-beta stimulation of Kang et al. Kang et al., 2018: blood cells of eight lupus
patients, each donor’s cells split into a culture stimulated with interferon-beta for
six hours and an untreated one. The cells of each sample and cell type are pooled into
one profile, so that every test compares donors, not cells (Choosing a
method says why). The cube is the kang2018
dataset (Datasets), a download of 77 MB that
python docs/datasets.py prepare kang2018 makes.
import matplotlib.pyplot as plt
import cytepool
from cytepool.design import Design
from cytepool.differential_expression import limma_trend, limma_voom
from cytepool.io import read_cube
from cytepool.plotting import (
plot_cell_numbers,
plot_md,
plot_significant_genes,
plot_volcano,
)
cube = read_cube("kang2018.cube.zarr")
print(f"CytePool {cytepool.__version__}")
print(cube)
plot_cell_numbers(cube)
plt.show()CytePool 0.1.dev54+g3f0cbb52e
Cube with 16 samples x 8 groups x 35,635 features
layers: counts (uint32, lazy), detected (uint32, lazy)
filled: 128 of 128 sample x group entries have cells (100.0%)
samples: sample_id, donor, condition
groups: group_id
features: feature_id, symbol
A profile pooled from a handful of cells is noisy, so limma_voom leaves the
profiles of fewer than 10 cells out of their cell type’s fit (min_cells=10); they
are pale in the heatmap. Two samples of the dendritic cells and four of the
megakaryocytes are left out, and every other profile enters.
A paired design¶
Every donor gives a control and a stimulated sample. The design ~ donor + condition
compares each donor’s two samples, so that the differences between donors do not count
as noise, and its coefficient conditionstim is the log2 fold change of the
stimulation. Name the levels of the condition: the first is the reference.
Design.check looks for designs that cannot answer the question before anything
is fitted (Designs, linear models and contrasts).
design = Design(
cube.samples,
"~ donor + condition",
levels={"condition": ["ctrl", "stim"]},
contrasts={"stim_vs_ctrl": "conditionstim"},
)
print(design)
print(design.check())Design: 16 samples x 9 columns, formula '~ donor + condition'
columns: (Intercept), donor1015, donor1016, donor1039, donor107, donor1244, donor1256, donor1488, conditionstim
contrasts: stim_vs_ctrl
Design check: 16 samples, 9 columns, rank 9, 7 residual degrees of freedom; no problems found.
Test every cell type¶
limma_voom runs limma-voom in each cell type in turn: it leaves out the profiles
with too few cells, keeps the genes with enough counts for the design, normalizes the
libraries by TMM, gives every observation a voom precision weight, fits the design to
every gene and moderates the contrast by empirical Bayes Law et al., 2014Ritchie et al., 2015. The
figure counts the genes up and down at an adjusted p-value below 0.05 in each cell
type, the first figure to look at in a study of many cell types.
result = limma_voom(cube, design)
print(result)
plot_significant_genes(result)
plt.show()DifferentialExpressionResult (limma_voom): 8 groups, 8 fitted; 1 contrast: stim_vs_ctrl
B cells: 16 samples, 2768 features
CD14+ Monocytes: 16 samples, 5125 features
CD4 T cells: 16 samples, 6295 features
CD8 T cells: 16 samples, 1453 features
Dendritic cells: 14 samples, 2281 features
FCGR3A+ Monocytes: 16 samples, 2947 features
Megakaryocytes: 12 samples, 704 features
NK cells: 16 samples, 2292 features
A cell type left with too few samples to estimate a variance is not fitted, and both
the printed summary and the figure say so (From cells to
results shows one); result.summary gives the
reason and result.samples lists every profile with the reason it was left out. Each fit keeps
the rows of the design for its own samples, so the megakaryocytes, which have no
profile of donor 107, are fitted without that donor’s column.
One cell type in detail¶
Each cell type and contrast has a table with the columns of limma’s topTable
(DifferentialExpressionResult.table). The mean-difference plot shows the log2
fold changes against the average expression, and the volcano plot the p-values against
the fold changes, with the genes of the smallest p-values named.
table = result.table("CD14+ Monocytes", "stim_vs_ctrl")
figure, (left, right) = plt.subplots(1, 2, figsize=(11, 4.2), layout="constrained")
plot_md(table, ax=left)
plot_volcano(table, ax=right)
plt.show()
The table converts to a pandas DataFrame (Table.to_pandas) for sorting and
export; here are the ten genes with the smallest p-values.
top = table.to_pandas().sort_values("p_value").head(10).set_index("symbol")
top[["log2_fold_change", "average_log2_cpm", "t", "adjusted_p_value"]]voom or trend¶
limma_trend fits the log-CPM values without precision weights and lets the prior
variance follow the expression level. It suits libraries of similar sizes, within
about threefold of each other by the limma User’s Guide (Choosing a
method). A pseudobulk library grows with the
number of cells pooled, and the libraries of each cell type here differ far more:
used = result.samples.filter(result.samples["used"])
for group in result.groups:
sizes = used.filter(used["group_id"] == group)["library_size"]
print(
f"{group:<18} largest library {sizes.max() / sizes.min():4.1f} times the smallest"
)
trend = limma_trend(cube, design)
figure, (left, right) = plt.subplots(
1, 2, figsize=(11, 4.5), layout="constrained", sharey=True
)
plot_significant_genes(result, ax=left)
left.set_title("limma-voom")
plot_significant_genes(trend, ax=right)
right.set_title("limma-trend")
plt.show()B cells largest library 14.5 times the smallest
CD14+ Monocytes largest library 5.5 times the smallest
CD4 T cells largest library 9.8 times the smallest
CD8 T cells largest library 31.2 times the smallest
Dendritic cells largest library 3.9 times the smallest
FCGR3A+ Monocytes largest library 7.1 times the smallest
Megakaryocytes largest library 8.5 times the smallest
NK cells largest library 10.1 times the smallest
limma-trend calls fewer genes than limma-voom in every cell type, and the gap is widest where the libraries differ most, as in the CD8 T cells; where they differ least, as in the dendritic cells and the CD14+ monocytes, the two agree more closely. voom’s weights give the small libraries less say.
Adjusting within or across cell types¶
The p-values of each cell type and contrast are adjusted over its own genes, as limma’s
topTable adjusts them. adjust_across_groups=True adds a column adjusted over the
genes of every cell type for the same contrast, one family of tests for the study
(Multiple testing).
across = limma_voom(cube, design, adjust_across_groups=True)
figure, (left, right) = plt.subplots(
1, 2, figsize=(11, 4.5), layout="constrained", sharey=True
)
plot_significant_genes(across, ax=left)
left.set_title("Adjusted within each cell type")
plot_significant_genes(across, across_groups=True, ax=right)
right.set_title("Adjusted across cell types")
plt.show()
Adjusting across cell types is not stricter for every cell type. The many small p-values of the monocytes lower the adjusted p-values of the other cell types’ genes, which gain calls, while the monocytes, now adjusted with genes that respond less, lose some.
Choices¶
voom or trend.
limma_voomis the default, and the safer choice when library sizes vary with the number of cells pooled;limma_trendis faster and suits libraries of similar sizes.Profiles with few cells or counts.
min_cells(10) andmin_counts(0) decide which profiles enter each fit; raise them when small pools look noisy on an MDS plot, or keep every profile and weight it by its number of cells withprior_weights="n_cells"Hoffman et al., 2026.quality_weights=Truedown-weights the samples that follow the model less closely instead. Normalization, filtering and sample quality shows both.Within or across cell types. Adjust within each cell type to rank its genes; adjust across cell types when the study’s list of genes, of every cell type, is the result.
A variable of the donors. The paired design cannot test a variable that is constant within each donor, such as sex or disease: Repeated measures and donors tests it with a random effect for the donor.
Other models. edgeR’s quasi-likelihood F-tests and DESeq2 run on one cell type’s counts at a time (
Cube.layer), as on bulk RNA-seq.
The responding genes feed the next questions: which pathways and regulators they belong to (Gene sets and pathways) and which regulators are active (Pathway and regulator activity).
- Kang, H. M., Subramaniam, M., Targ, S., Nguyen, M., Maliskova, L., McCarthy, E., Wan, E., Wong, S., Byrnes, L., Lanata, C. M., Gate, R. E., Mostafavi, S., Marson, A., Zaitlen, N., Criswell, L. A., & Ye, C. J. (2018). Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nature Biotechnology, 36(1), 89–94. 10.1038/nbt.4042
- Law, C. W., Chen, Y., Shi, W., & Smyth, G. K. (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology, 15(2), R29. 10.1186/gb-2014-15-2-r29
- Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., & Smyth, G. K. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47. 10.1093/nar/gkv007
- Hoffman, G. E., Lee, D., Bendl, J., Prashant, N. M., Hong, A., Casey, C., Alvia, M., Shao, Z., Argyriou, S., Therrien, K., Venkatesh, S., Voloudakis, G., Haroutunian, V., Fullard, J. F., & Roussos, P. (2026). Efficient differential expression analysis of large-scale single-cell transcriptomics data using Dreamlet. Nature Communications, 17(1), 9597. 10.1038/s41467-026-75680-8