The Differential-Expression Test Suite — R Seurat vs Truecell (Python)¶

Wave 3's first side-by-side, and the last large untested surface in the library: find_markers offers eight statistical tests and none of them had ever been compared to R. Their unit tests assert self-consistency on synthetic fixtures — the same shape of coverage that let the CLR and SCTransform defects survive.

Dataset: pbmc3k — 2,700 PBMCs, 10x Genomics (2016). The comparison runs on clusters 0 and 1 (692 and 515 cells, 13,714 genes). R reference: Seurat 5.5.1 · MAST 1.38.0 · DESeq2 1.52.0 · Python: Truecell

Seurat	Truecell
`FindMarkers(obj, test.use = "wilcox")`	`find_markers(obj, test_use="wilcox")`
`test.use = "t"` · `"bimod"` · `"LR"`	`test_use="t"` · `"bimod"` · `"LR"`
`test.use = "negbinom"` · `"roc"`	`test_use="negbinom"` · `"roc"`
`test.use = "MAST"`	`test_use="mast"`
`test.use = "DESeq2"`	`test_use="deseq2"`

This tutorial found and fixed two defects. avg_log2FC put Seurat's pseudocount on the group mean rather than the group sum, which floored every fold change and — because logfc_threshold filters on that value — changed which genes were returned at all. And negbinom ran a moment-dispersion likelihood-ratio test where Seurat runs an ML-dispersion Wald test. Both are written up in what this tutorial found.

Both tools test the same cells: Python clusters pbmc3k and writes the assignment to figures_de/groups.csv, which the R side reads. Louvain numbering is not guaranteed to agree across implementations, and a clustering difference would look exactly like a DE difference.

Headline¶

Metric	Result
`avg_log2FC` vs Seurat, all 13,712 shared genes	max abs diff 6.44e-15
Tests reproducing Seurat's top 50 genes	7 of 7 per-cell tests (`roc` scores AUC, not p)
`wilcox` · `t` · `bimod` · `LR` — p-value Spearman	1.000000 · 0.999980 · 0.999994 · 0.999975
`mast` — Spearman (all genes / detected >5%)	0.9471 / 0.9979
`negbinom` — Spearman (all genes / detected >5%)	0.6943 / 0.9165
`roc` — max abs AUC difference	5.0e-04, which is Seurat's own 3-dp rounding
Before the fix — genes returned at `logfc_threshold=0.25`	truecell 2,298 vs Seurat 11,931 (Jaccard 0.193)

Setup¶

R (Seurat) Python (Truecell)

library(Seurat)

obj <- CreateSeuratObject(
  Read10X(".../hg19"), project = "pbmc3k_de",
  min.cells = 3, min.features = 200)
obj[["percent.mt"]] <- PercentageFeatureSet(
  obj, pattern = "^MT-")
obj <- subset(obj, subset =
  nFeature_RNA > 200 & nFeature_RNA < 2500 &
  percent.mt < 5)
obj <- NormalizeData(obj, verbose = FALSE)

# groups written by the Python side
g <- read.csv("figures_de/groups.csv", row.names = 1)
obj <- subset(obj, cells = rownames(g))
Idents(obj) <- factor(g[colnames(obj), "group"])

from truecell.datasets import pbmc3k
from truecell import create_truecell_object
from truecell.preprocessing import normalize_data
from truecell.markers import find_markers

obj = build()          # same QC, then cluster
groups = shared_groups(obj)
groups.to_csv("figures_de/groups.csv")

sub = obj.subset(cells=list(groups.index))
sub.idents = list(groups.values)

Both sides then run with logfc.threshold = 0 and min.pct = 0, so the comparison sees every gene rather than only those surviving a filter whose input is one of the numbers under test. Seurat's defaults would have hidden exactly the disagreement worth seeing.

Running the tests¶

R (Seurat) Python (Truecell)

res <- FindMarkers(obj, ident.1 = "0", ident.2 = "1",
                   test.use = "wilcox",
                   logfc.threshold = 0, min.pct = 0)
# R underflows 172 p-values to exactly 0; rank among the ones it scored.
res <- res[res$p_val > 0, ]
head(res[order(res$p_val), c("p_val","avg_log2FC")], 3)
#>                p_val avg_log2FC
#> TYROBP  4.593065e-219  -6.277211
#> CST3    1.282760e-214  -6.130653
#> S100A9  4.216270e-206  -7.446613

res = find_markers(sub, "0", "1", test_use="wilcox",
                   logfc_threshold=0, min_pct=0)
res.head(3)[["p_val", "avg_log2FC"]]
#>                p_val  avg_log2FC
#> TYROBP  4.593065e-219   -6.277211
#> CST3    1.282760e-214   -6.130653
#> S100A9  4.216270e-206   -7.446613

Both columns are identical to seven significant figures — the same p-values and the same fold changes, on the same cells.

All eight tests against Seurat

What this tutorial found¶

1. `avg_log2FC` put the pseudocount in the wrong place¶

Seurat 5's log1pdata.mean.fxn is, verbatim:

log(x = (rowSums(x = expm1(x = x)) + pseudocount.use) / NCOL(x), base = base)

One pseudocount added to the group's sum, then divided by n — so on the mean scale it is worth 1/n, not 1. truecell computed log2(mean(expm1(x)) + 1), adding a whole count to the mean. That is Seurat 4's formula; the repo targets Seurat 5.

The effect is to floor every fold change toward zero. A gene detected in 0 % of cluster 0 and 24 % of cluster 1 read −1.26 where Seurat reads −9.92.

avg_log2FC before and after

Why this is a defect and not a cosmetic difference. logfc_threshold filters on this value, so the error did not merely misreport fold changes — it changed which genes came back:

`logfc_threshold`	truecell (before)	Seurat	Jaccard
0.1 (Seurat's default)	4,903	13,009	0.377
0.25	2,298	11,931	0.193
0.5	981	10,228	0.096
1.0	335	6,928	0.048

At a common 0.25 threshold, fewer than one gene in five agreed.

What the threshold let through

The most telling part: where both groups express a gene, the two formulas nearly agree (Spearman 0.990 on the 1,362 genes with pct > 0.1 in both). The error was concentrated in sparse, marker-like genes — precisely what differential expression exists to find. After the fix, 6.44e-15 across all 13,712 genes.

There was already a test for this. test_avg_log2fc_matches_seurat_formula re-implemented the same wrong formula and checked that truecell agreed with itself. It was green throughout while carrying a name that claimed Seurat parity — the same shape as #48's test_fetch_data. It is corrected here.

2. `negbinom` was running a different test¶

Seurat's GLMDETest fits MASS::glm.nb — which estimates the dispersion by maximum likelihood — and reads the Wald p-value off the group coefficient (summary(...)$coef[2, 4]). truecell used a fixed method-of-moments dispersion and a likelihood-ratio test: a different estimator and a different statistic. HLA-DRA read 5.5e-128 against R's 1.1e-321.

After the fix the p-values agree exactly for every gene anyone would look at:

detection (max of the two groups)	genes	median \|log10 ratio\|	Spearman
> 25 %	980	0.000	0.988
10 – 25 %	1,519	0.000	0.933
5 – 10 %	1,919	0.000	0.773
1 – 5 %	5,120	0.063	0.285
< 1 %	1,928	0.144	0.085

What disagreement remains sits below 5 % detection, where the negative-binomial GLM is fitting almost-empty rows and neither tool is estimating anything meaningful. Seurat agrees: its min.cells.feature default drops those genes, and in this run R returned 11,466 genes against truecell's 13,714 — every one of the 2,248 it dropped was below 1 % detection in both groups (the highest reached 0.4 %). The headline Spearman of 0.69 is that tail; on genes Seurat would actually have tested, it is 0.92.

Differences left standing, and why¶

deseq2 is not Seurat's DESeq2. DESeq2DETest builds a DESeqDataSet with one column per cell and tests cells as replicates. truecell sums counts per sample and tests at the sample level. Treating cells as replicates is the practice Squair et al. (2021) showed inflates false positives, so pseudobulk is the better statistics — and because it requires sample_col, it cannot be silently mistaken for the per-cell test: it raises. Reported rather than changed in either direction.

mast is a hand-rolled hurdle model, not a call to the MAST package, which has no Python equivalent to depend on. Spearman 0.947 across all genes, 0.9979 on genes detected above 5 %, and the same top 50. Worth knowing: Seurat's MASTDETest fits ~ condition alone — it adds no cellular detection rate term unless you pass one. truecell's docstring previously advised passing CDR "to match Seurat's default CDR covariate", which had it backwards; that is corrected.

Seurat rounds myAUC to three decimals inside DifferentialAUC, so the ROC comparison cannot be tighter than 5e-4 however correct both sides are. That is R's rounding, not a divergence — stated because it looks like one.

R's wilcox returns NaN for the 365 genes with no expression in either group; truecell returns p = 1. A test that cannot be run has no evidence against the null, so 1 is the more useful answer, and R's NaN set is a subset of truecell's.

R also returns exactly 0 for 172 genes, and that is not double underflow: truecell scores 92 of them above 1e-50, the largest — SEPT1 — at 9.3e-17. It is not Seurat's wrapper either. Calling base R's wilcox.test on that gene's normalised row directly, outside Seurat, returns 0 as well (W = 212196.5 on 692 vs 515 cells). Both R's NaN rows and its zero rows are excluded from every correlation reported here, since neither carries a rank — but the zeros are worth knowing about, because they are not the harmless precision limit the previous version of this note called them.

Parity — verified against R Seurat¶

Test	Genes	max \|Δlog2FC\|	p Spearman (all)	p Spearman (detected >5 %)	Top 50
`wilcox`	13,712	6.2e-15	1.000000	1.0000	50/50
`t`	13,712	6.2e-15	0.999980	1.0000	50/50
`bimod`	13,712	6.2e-15	0.999994	1.0000	50/50
`LR`	13,712	6.2e-15	0.999975	1.0000	50/50
`negbinom`	11,466	6.2e-15	0.694340	0.9165	50/50
`roc`	13,712	6.2e-15	AUC 5.0e-04	—	—
`mast`	13,712	6.2e-15	0.947038	0.9979	50/50
`deseq2`	13,712	3.47	0.476942	0.1959	22/50

The last digits of these moved slightly when the CSV round-trip was fixed (see The two columns a person actually reads, below): they had been read back through a misparsing float reader. The change is at the ULP level and no band moved, but the table is the measured one, not the previous one.

deseq2's row is the pseudobulk-vs-per-cell divergence described above, not a defect; its fold change differs too because a pseudobulk fold change is computed on summed counts.

The two columns a person actually reads¶

The table above scores the statistic. It does not score the two numbers someone annotating clusters looks at: the fold change they sort by, and the adjusted p-value they threshold on. Neither was compared until an expert reviewer pointed out that the max-difference bound and a set overlap answer neither question.

Test	logFC Spearman	logFC Kendall	Top 50 by \|logFC\|	`p_val_adj` to 6 s.f.	Same call at 0.05	Genes differing
`wilcox`	1.000000	1.000000	50/50	0.9933	1.0000	0
`t`	1.000000	1.000000	50/50	0.9999	1.0000	0
`bimod`	1.000000	1.000000	50/50	0.9997	0.9999	2
`LR`	1.000000	1.000000	50/50	0.9915	1.0000	0
`negbinom`	1.000000	1.000000	50/50	0.8547	0.9958	48
`roc`	1.000000	1.000000	50/50	—	—	—
`mast`	1.000000	1.000000	50/50	0.7984	0.9895	144
`deseq2`	0.966510	0.847811	34/50	0.4314	0.8938	1,401

Rank correlation is reported alongside the max-difference bound rather than instead of it, because the two fail differently. A uniform scale error leaves every rank perfect and blows up the max; a handful of swapped mid-table genes leaves the max tiny and moves the ranks. Here both are clean: fold-change order is preserved exactly for all seven cell-level tests.

Do identical adjusted p-values occur? Mostly not, and the reason is worth stating rather than the rate. The correction is identical — both tools compute min(p × 13,714, 1), which holds bit-exactly on both sides — but the raw p-values feeding it differ by up to 0.54 % relative on wilcox, a real difference between SciPy's Wilcoxon and Seurat's, so the product rarely lands on the same double. What survives that is what matters: the ordering is exact, and every gene falls on the same side of 0.05 for wilcox, t and LR.

Two traps in measuring this, both of which had to be fixed before the numbers above meant anything:

Seurat clamps p_val_adj at 1, and 11,858 of 13,712 genes land there. A bare "fraction identical" therefore reads ≈0.87 before a single interesting gene is considered. --report scores the unclamped subset separately.
Neither side's CSV round-tripped a float64. R's write.csv renders 15 significant digits, and raising it does not help — R's own sprintf("%.17g") is not correctly rounded and emits digits denoting a different double. Pandas' read_csv misparses about a third of random doubles by an ULP at its default float_precision. So the R side now also writes r_<test>_exact.csv in C99 hex float (%a, a transcription of the IEEE-754 bits rather than a decimal approximation of them), and the Python side reads with float_precision="round_trip". Before that, an "is it identical" comparison was measuring the two languages' text formatters.

These numbers are asserted, not just printed¶

Every row above used to live only in this file, so nothing failed when one moved. deseq2's overlap was written here as 25/50 and had drifted to 22 — the clusters this tutorial tests are found by find_clusters, and the graph fixes in #67-#71 moved a few cells between them. That is a legitimate reason for the number to change, which is exactly why it needed a band rather than a sentence: a regression landing on 22 would have read the same way.

--report now checks each number against a declared range and exits non-zero if one falls outside:

band	range	why
top 50, the seven cell-level tests	= 50	Same statistic, same cells. One dropped gene is a regression.
top 50, `deseq2`	15 – 32	A divergence measurement. 20–26 over 20 resampled replicate splits; 25 on the previous clustering. Bounded well below 50 — reaching parity would mean `sample_col` had stopped being honoured.
p Spearman >5 %, `wilcox`/`t`/`bimod`/`LR`	≥ 0.9999	Measured at exactly 1.0.
p Spearman >5 %, `negbinom`	≥ 0.88	Same model, different optimiser: 0.9165.
p Spearman >5 %, `mast`	≥ 0.99	A hand-rolled hurdle model, not the MAST package: 0.9979.
p Spearman >5 %, `deseq2`	0.12 – 0.30	Pseudobulk against per-cell; a high value here would be the surprise.
max \|Δlog2FC\|, cell-level tests	≤ 1e-12	Arithmetic on the shared matrix. `deseq2` is excluded by name, not by threshold — its 3.47 is correct and would otherwise set everyone else's tolerance.
max \|ΔAUC\|, `roc`	≤ 5e-4	Half a unit in Seurat's third decimal. Measured 4.9986e-4, i.e. on the boundary.

The reference has to be the one the handoff asked for¶

--report also refuses to compare against an R run that predates the groups.csv it was supposed to answer. This is not hypothetical: on the working copy where these bands were written, the Python tables were from 25 July and the R tables from 19 July, taken on a different cluster assignment — and the report printed a full parity table showing wilcox at 48/50 and a Spearman of 0.907. Nothing in it said "stale file"; it read as a regression in the port.

The check is pct.1 and pct.2. They are counts of detected cells per group with no statistics in the way, so two runs over the same handoff agree to Seurat's three-decimal rounding and no worse. They differed for 12,491 of 13,712 genes.

Runtime, for scale: Seurat's slowest test here is negbinom at 91.8 s (MAST 60.5 s, DESeq2 60.3 s); truecell's are 36.0 s, 28.5 s and 3.4 s.

Running it¶

# 1. Python side — writes figures_de/groups.csv and py_<test>.csv
python tutorials/pbmc3k_de_tutorial.py

# 2. R side — writes figures_de/r_<test>.csv (needs MAST + DESeq2)
Rscript tutorials/pbmc3k_de_verify.R

# 3. Compare
python tutorials/pbmc3k_de_tutorial.py --report

# 4. Figures
python tutorials/generate_de_plots.py

The R side needs MAST and DESeq2:

BiocManager::install(c("MAST", "DESeq2"))

Not glmGamPoi. It is a Suggests of DESeq2 rather than an Imports, FindMarkers never calls it, and installing it flips sctransform's vst onto a different backend — which would move the SCTransform R reference that sctransform_vignette.md is pinned against. Installing with default dependencies leaves Suggests alone, which is what you want. This was verified rather than assumed: an SCTransform fingerprint taken before and after the install is byte-identical.