A bulk RNA-seq experiment is a matrix of counts, samples x genes. CytePool tests it with the methods of limma, edgeR and DESeq2, each checked against the R package (Validation). This page runs the three on one simulated experiment whose true changes are known, and compares their calls with the truth. Choosing a method explains how they differ, and the Quick start walks through limma-voom step by step.
The experiment is the simulated bulk dataset (Datasets):
three treatments (control, drug A and drug B) with four replicates each, sequenced in
two batches; drug A changes 500 of the 10,000 genes and drug B 300, and one sample,
drug_b_3, is noisier than the others. python docs/datasets.py prepare simulated_bulk writes it, as a cube with one group, bulk.
import numpy as np
import cytepool
from cytepool.design import Design
from cytepool.differential_expression import limma_voom
from cytepool.io import read_cube
from cytepool.normalization import filter_by_expression, tmm_factors
from cytepool.stats.deseq2 import (
deseq2_likelihood_ratio_test,
deseq2_wald_test,
fit_deseq2,
)
from cytepool.stats.quasi_likelihood import (
fit_quasi_likelihood,
quasi_likelihood_f_test,
)
cube = read_cube("simulated_bulk.cube.zarr")
truth = np.load("simulated_bulk_truth.npz")
print(f"CytePool {cytepool.__version__}, numpy {np.__version__}")
print(cube)CytePool 0.1.dev42+g6991e63f8, numpy 2.5.4
Cube with 12 samples x 1 groups x 10,000 features
layers: counts (uint32, lazy)
cells per entry: not recorded
samples: sample_id, condition, batch
groups: group_id
features: feature_id, symbol
The design¶
The batch enters the model so that its shift of every gene does not count as noise. Both drugs are compared with the control, the reference level.
design = Design(
cube.samples,
"~ batch + condition",
levels={"condition": ["control", "drug_a", "drug_b"]},
contrasts={
"drug_a_vs_control": "conditiondrug_a",
"drug_b_vs_control": "conditiondrug_b",
},
)
print(design)
print(design.check())Design: 12 samples x 4 columns, formula '~ batch + condition'
columns: (Intercept), batchbatch2, conditiondrug_a, conditiondrug_b
contrasts: drug_a_vs_control, drug_b_vs_control
Design check: 12 samples, 4 columns, rank 4, 8 residual degrees of freedom; no problems found.
limma-voom, with sample quality weights¶
limma_voom runs on the cube’s one group: genes are filtered by expression, TMM
factors normalize the libraries, voom gives every observation a precision weight,
and the contrasts are moderated by empirical Bayes. quality_weights=True also
estimates a weight for every sample from how closely it follows the model
(Precision weights). The noisy sample gets a
weight well below the others, so it counts less instead of being removed.
voom_result = limma_voom(cube, design, quality_weights=True)
print(voom_result)
samples = voom_result.samples
for sample_id, weight in zip(samples["sample_id"], samples["quality_weight"]):
print(f" {sample_id:<10} {weight:5.2f} {'#' * round(20 * weight)}")DifferentialExpressionResult (limma_voom): 1 group, 1 fitted; 2 contrasts: drug_a_vs_control, drug_b_vs_control
bulk: 12 samples, 9593 features
control_1 1.12 ######################
control_2 1.13 #######################
control_3 1.10 ######################
control_4 1.11 ######################
drug_a_1 1.14 #######################
drug_a_2 1.12 ######################
drug_a_3 1.20 ########################
drug_a_4 1.08 ######################
drug_b_1 1.05 #####################
drug_b_2 1.06 #####################
drug_b_3 0.30 ######
drug_b_4 1.18 ########################
edgeR’s quasi-likelihood F-tests and DESeq2¶
fit_quasi_likelihood and quasi_likelihood_f_test take the counts of the
genes that pass edgeR’s filter (filter_by_expression) with TMM factors
(tmm_factors). fit_deseq2 takes every gene and filters on the mean
count when it adjusts the p-values; deseq2_wald_test tests one contrast, and
deseq2_likelihood_ratio_test tests both drugs at once against the design with
the batch only.
counts = cube.layer("counts", "bulk")
keep = filter_by_expression(counts, design=design)
quasi = fit_quasi_likelihood(
counts[:, keep], design, norm_factors=tmm_factors(counts[:, keep])
)
deseq = fit_deseq2(counts, design)
either_drug = deseq2_likelihood_ratio_test(deseq, Design(cube.samples, "~ batch"))
print(f"{keep.sum():,} of {keep.size:,} genes pass the expression filter")
print(
f"quasi-likelihood: common dispersion {quasi.dispersion:.3f}, prior df {np.median(quasi.prior_df):.1f}"
)
print(
f"DESeq2 likelihood ratio test of both drugs: {np.sum(either_drug.adjusted_p_value < 0.05):,} genes at FDR 5%"
)9,593 of 10,000 genes pass the expression filter
quasi-likelihood: common dispersion 0.072, prior df 9.9
DESeq2 likelihood ratio test of both drugs: 563 genes at FDR 5%
Calls compared with the truth¶
At a false discovery rate of 5%, each method’s calls are compared with the genes the simulation changed: how many it called, how many of those truly changed, and the share of its calls that are false.
genes = np.asarray(cube.features["feature_id"])
truly_changed = {
"drug_a_vs_control": truth["true_log2_fold_change_drug_a"] != 0,
"drug_b_vs_control": truth["true_log2_fold_change_drug_b"] != 0,
}
print(f"{'contrast':<19}{'method':<18}{'called':>7}{'true':>7}{'false share':>13}")
for name in design.contrasts.names:
voom_table = voom_result.table("bulk", name)
voom_called = np.isin(
genes,
np.asarray(voom_table["feature_id"])[
np.asarray(voom_table["adjusted_p_value"]) < 0.05
],
)
quasi_called = np.zeros(genes.size, dtype=bool)
quasi_test = quasi_likelihood_f_test(quasi, contrasts=design.contrasts.select(name))
quasi_called[keep] = quasi_test.adjusted_p_value < 0.05
deseq_called = (
deseq2_wald_test(deseq, design.contrasts.select(name)).adjusted_p_value < 0.05
)
for method, called in [
("limma-voom", voom_called),
("edgeR QL", quasi_called),
("DESeq2 Wald", deseq_called),
]:
true_calls = np.sum(called & truly_changed[name])
share = 1 - true_calls / max(called.sum(), 1)
print(f"{name:<19}{method:<18}{called.sum():>7,}{true_calls:>7,}{share:>13.1%}")
print(
f"{'':<19}{'called by all':<18}{np.sum(voom_called & quasi_called & deseq_called):>7,}"
)contrast method called true false share
drug_a_vs_control limma-voom 355 337 5.1%
drug_a_vs_control edgeR QL 320 312 2.5%
drug_a_vs_control DESeq2 Wald 360 338 6.1%
called by all 320
drug_b_vs_control limma-voom 204 192 5.9%
drug_b_vs_control edgeR QL 193 184 4.7%
drug_b_vs_control DESeq2 Wald 238 203 14.7%
called by all 188
The three methods find mostly the same genes: 320 genes are called by all three for drug A. They miss some of the changed genes, as the smallest simulated changes, a fold change of about 1.4, are hard to detect with four replicates. In this simulation DESeq2’s Wald test made more false calls for drug B than the 5% that the adjustment allows on average, also when the noisy sample is left out; one simulated data set is no benchmark, and the Literature lists comparisons on many. Which method to prefer depends mostly on the design and on what comes next: limma’s tables carry moderated standard errors that a meta-analysis combines across studies, and DESeq2’s likelihood ratio test asks about several coefficients at once.