diff --git a/CHANGELOG.md b/CHANGELOG.md index 5aada97..ec4dc3e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,48 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Fixed + +- **`vln_plot` drew the wrong violin, three ways.** All three were checked + against R 4.6.1 / Seurat 5.5.1 running locally, not against recollection of + what Seurat does. + + **The bandwidth was 2.3x too wide.** The density used scipy's `"scott"` rule, + which scales the sample standard deviation. R's `stats::density` — and so + `geom_violin`, and so Seurat — uses `bw.nrd0`, which takes + `min(sd, IQR/1.34)`. Expression is zero-inflated, so the IQR term is much the + smaller, and measured against R on zero-inflated draws scipy's bandwidth came + out 2.0-2.5x wider. The effect is over-smoothing that flattens the spike at + zero, which is the shape of the distribution. `_bw_nrd0` now reproduces R to + 5e-13 across seven cases including both zero-spread fallback branches. + + (`bw.nrd0` divides the IQR by **1.34**. The neighbouring rule `bw.nrd`, which + R also ships, uses 1.349 — taking that one is a silent 0.67% error wherever + the IQR term wins. The first version of this fix had exactly that bug, caught + by diffing against R.) + + **The violin was not trimmed.** `geom_violin(trim = TRUE)` limits the density + to the observed range; an untrimmed gaussian KDE tails off past it, so every + violin extended below zero where expression cannot go. + + **Points were off by default.** `pt_size` defaulted to `0`. Seurat's `VlnPlot` + passes `pt.size = NULL`, which `ExIPlot` resolves through `AutoPointSize` — + `min(1583/n, 1)` — so points are shown and shrink as the cell count grows. + `pt_size=None` is now the default and follows that rule; `pt_size=0` still + omits them. + + The violin outline now matches `geom_violin(scale = "width", trim = TRUE)`'s + own computed polygon to 0.2% of full width, with the support equal to the data + range exactly. Nine tutorial figures are regenerated; they are the nine drawn + by the six generators that call `vln_plot`, and no others moved. + + Two additions fall out of the rewrite: `violin_width` sets the width of a full + violin, and `jitter_seed` (default `0`) makes the point jitter reproducible so + a committed figure redraws identically. + + Still deliberately different from Seurat: the median bar. `geom_violin` draws + none, and this keeps drawing one. + ### Added - **`ridge_plot` warns when an explicit `figsize` is too small for the group diff --git a/tests/test_plotting_violin.py b/tests/test_plotting_violin.py new file mode 100644 index 0000000..2c22167 --- /dev/null +++ b/tests/test_plotting_violin.py @@ -0,0 +1,232 @@ +"""Violin geometry: bandwidth, trim, width scaling, and the jittered points. + +A Seurat violin is `geom_violin(scale = "width", adjust = adjust, trim = TRUE)` +smoothed with R's `nrd0` bandwidth, plus jittered points whose size comes from +`AutoPointSize`. Each of those four is asserted here against values taken from +R 4.6.1 / Seurat 5.5.1 rather than from a reimplementation of the same idea. +""" +import numpy as np +import pytest +import scipy.sparse as sp + +import truecell as tc +from truecell.plotting import _auto_point_size, _bw_nrd0 + +# `stats::bw.nrd0(x)` in R 4.6.1, printed to 12 places. +R_BW_NRD0 = { + "arange10": (np.arange(10.0), 1.719286404692), + "zeros_then_tail": (np.array([0.0] * 10 + [1, 2, 3, 8.0]), 0.297148380291), + "mostly_zero": (np.array([0, 0, 0, 0, 1.0]), 0.291718187405), + "zero_inflated": ( + np.array([0.1, 0.4, 0.9, 1.6, 2.5, 3.6, 4.9, 6.4, 8.1, 10.0, 0, 0, 0, 0, 0]), + 1.660765786387, + ), + "two_points": (np.array([1.0, 2.0]), 0.292349069764), + # Zero spread takes bw.nrd0's fallback chain: sd and IQR are both 0, so it + # falls to abs(x[1]) and then to 1. Both branches are exercised. + "all_same_nonzero": (np.array([5.0, 5, 5, 5]), 3.410362274648), + "all_zero": (np.array([0.0, 0, 0, 0]), 0.682072454930), +} + + +@pytest.mark.parametrize("name", sorted(R_BW_NRD0)) +def test_bw_nrd0_matches_r(name): + x, expected = R_BW_NRD0[name] + assert _bw_nrd0(x) == pytest.approx(expected, abs=1e-11) + + +def test_bw_nrd0_uses_the_right_divisor(): + """`bw.nrd0` divides the IQR by 1.34; the neighbouring rule `bw.nrd` uses + 1.349. Taking the wrong one is a silent 0.67% error in every bandwidth where + the IQR term wins — which, on zero-inflated expression, is most of them.""" + # The IQR term has to win *and* be non-zero. Too much zero inflation puts + # both quartiles at zero, which sends bw.nrd0 down its fallback chain + # instead — so a tight bulk with a heavy tail, not 80% zeros. + x = np.array([0.0] * 20 + [1.0] * 20 + [50.0] * 5) + hi = float(np.std(x, ddof=1)) + q75, q25 = np.percentile(x, [75, 25]) + iqr_term = (q75 - q25) / 1.34 + assert 0 < iqr_term < hi, "fixture does not exercise the IQR branch" + assert _bw_nrd0(x) == pytest.approx(0.9 * iqr_term * x.size ** -0.2) + # And the wrong divisor would be visibly different, not a rounding matter. + assert _bw_nrd0(x) != pytest.approx( + 0.9 * (q75 - q25) / 1.349 * x.size ** -0.2, rel=1e-4 + ) + + +def test_bw_nrd0_is_narrower_than_scipy_scott_on_expression_like_data(): + """The reason this exists. Scott scales the sd; nrd0 takes min(sd, IQR/1.34), + and zero inflation makes the IQR term much the smaller. scipy's default + over-smooths, flattening the spike at zero that is the shape of the data.""" + rng = np.random.default_rng(0) + x = np.concatenate([np.zeros(1400), rng.gamma(2, 1, 600)]) + scott = float(np.std(x, ddof=1)) * x.size ** -0.2 + assert scott > 2 * _bw_nrd0(x) + + +def test_bw_nrd0_rejects_a_single_point(): + with pytest.raises(ValueError, match="at least 2"): + _bw_nrd0(np.array([1.0])) + + +def test_auto_point_size_follows_seurats_rule(): + """Seurat: min(1583 / n, 1), capped so small objects get full-size points.""" + assert _auto_point_size(100) == pytest.approx(_auto_point_size(1583)) + assert _auto_point_size(1583) > _auto_point_size(2638) > _auto_point_size(10_000) + # The cap, and the 1/n^2 falloff in matplotlib's area units. + ratio = _auto_point_size(20_000) / _auto_point_size(10_000) + assert ratio == pytest.approx(0.25, rel=1e-6) + + +@pytest.fixture(scope="module") +def obj(): + pytest.importorskip("matplotlib") + rng = np.random.default_rng(0) + n_genes, n_cells = 12, 400 + counts = rng.poisson(0.4, size=(n_genes, n_cells)).astype(float) + counts[0, :120] = rng.poisson(30, size=120) # one clearly expressed gene + genes = [f"G{i:02d}" for i in range(n_genes)] + o = tc.create_truecell_object( + counts=sp.csc_matrix(counts), assay="RNA", feature_names=genes, + cell_names=[f"C{i:03d}" for i in range(n_cells)], project="vln", + ) + tc.normalize_data(o) + o.meta_data["grp"] = [["a", "b", "c", "d"][i % 4] for i in range(n_cells)] + return o + + +def _violin_paths(ax): + from matplotlib.collections import PolyCollection + return [c for c in ax.collections if isinstance(c, PolyCollection)] + + +def test_violins_are_trimmed_to_the_observed_range(obj): + """trim = TRUE. Expression cannot be negative, so a violin must not tail off + below the data — an untrimmed gaussian KDE does exactly that.""" + plt = pytest.importorskip("matplotlib.pyplot") + fig = tc.vln_plot(obj, ["G00"], group_by="grp", pt_size=0) + ax = fig.axes[0] + data_min = tc.plotting._get_expression(obj, "G00").min() + for coll in _violin_paths(ax): + ys = np.concatenate([p[:, 1] for p in coll.get_paths()[0].to_polygons()]) + assert ys.min() >= data_min - 1e-9, "violin extends below the data" + assert data_min >= 0 + plt.close(fig) + + +def test_every_violin_reaches_the_same_maximum_width(obj): + """scale = "width" — groups are compared on shape, not on cell count.""" + plt = pytest.importorskip("matplotlib.pyplot") + width = 0.8 + fig = tc.vln_plot(obj, ["G00"], group_by="grp", pt_size=0, violin_width=width) + widths = [] + for coll in _violin_paths(fig.axes[0]): + poly = coll.get_paths()[0].to_polygons()[0] + centre = np.round(poly[:, 0].mean()) + widths.append(2 * np.abs(poly[:, 0] - centre).max()) + assert len(widths) == 4 + assert all(w == pytest.approx(width, rel=1e-6) for w in widths) + plt.close(fig) + + +def test_the_rendered_violin_uses_nrd0_not_scipys_default(obj): + """The one that matters, and the one component-level bandwidth tests miss. + + Asserting `_bw_nrd0` matches R says nothing about whether `vln_plot` calls + it: swapping the KDE to `bw_method="scott"` left every other test in this + file green. So compare the *drawn* outline against both candidate densities + and require it to be the nrd0 one. + """ + from scipy.stats import gaussian_kde + plt = pytest.importorskip("matplotlib.pyplot") + + # Zero-inflated, where the two rules diverge most. + rng = np.random.default_rng(3) + n = 600 + vals = np.concatenate([np.zeros(420), rng.gamma(2.0, 1.0, n - 420)]) + o = tc.create_truecell_object( + counts=sp.csc_matrix(np.vstack([vals, vals])), + assay="RNA", feature_names=["F0", "F1"], + cell_names=[f"C{i:03d}" for i in range(n)], project="bw", + ) + o.assays["RNA"].layers["data"] = o.assays["RNA"].layers["counts"].copy() + o.meta_data["one"] = ["g"] * n + + width = 0.8 + fig = tc.vln_plot(o, ["F0"], group_by="one", pt_size=0, violin_width=width) + poly = _violin_paths(fig.axes[0])[0].get_paths()[0].to_polygons()[0] + plt.close(fig) + + # Right edge of the drawn outline, as half-width against y. + right = poly[poly[:, 0] >= poly[:, 0].mean()] + order = np.argsort(right[:, 1]) + ys, half = right[order, 1], right[order, 0] - poly[:, 0].mean() + + grid = np.linspace(vals.min(), vals.max(), 512) + sd = float(np.std(vals, ddof=1)) + nrd0 = gaussian_kde(vals, bw_method=_bw_nrd0(vals) / sd)(grid) + scott = gaussian_kde(vals, bw_method="scott")(grid) + pred_nrd0 = np.interp(ys, grid, nrd0 / nrd0.max()) * width / 2 + pred_scott = np.interp(ys, grid, scott / scott.max()) * width / 2 + + err_nrd0 = np.abs(half - pred_nrd0).max() + err_scott = np.abs(half - pred_scott).max() + + # Measured separation is ~590x (nrd0 4e-4 against scott 2.3e-1). The residual + # on the nrd0 side is polygon-extraction noise — the closing edges of the + # filled path and interpolation onto its vertices — not a bandwidth error, + # so the bound is loose in absolute terms and tight relative to the gap. + assert err_nrd0 < 5e-3, f"outline does not match nrd0 (err {err_nrd0:.2e})" + assert err_scott > 50 * err_nrd0, ( + f"nrd0 (err {err_nrd0:.2e}) and scott (err {err_scott:.2e}) are not " + "distinguishable here — the fixture is too weak to prove anything" + ) + + +def test_points_are_drawn_by_default(obj): + """Seurat's VlnPlot passes pt.size = NULL, which ExIPlot resolves through + AutoPointSize — so points are on unless asked otherwise.""" + plt = pytest.importorskip("matplotlib.pyplot") + fig = tc.vln_plot(obj, ["G00"], group_by="grp") + from matplotlib.collections import PathCollection + scatters = [c for c in fig.axes[0].collections if isinstance(c, PathCollection)] + assert scatters, "no jittered points drawn" + assert sum(len(c.get_offsets()) for c in scatters) == obj.meta_data.shape[0] + plt.close(fig) + + +def test_pt_size_zero_still_suppresses_points(obj): + plt = pytest.importorskip("matplotlib.pyplot") + fig = tc.vln_plot(obj, ["G00"], group_by="grp", pt_size=0) + from matplotlib.collections import PathCollection + assert not [c for c in fig.axes[0].collections + if isinstance(c, PathCollection)] + plt.close(fig) + + +def test_jitter_is_reproducible_by_default(obj): + """A figure that redraws differently every call cannot be diffed, and the + tutorial figures are committed.""" + plt = pytest.importorskip("matplotlib.pyplot") + from matplotlib.collections import PathCollection + + def offsets(**kw): + fig = tc.vln_plot(obj, ["G00"], group_by="grp", **kw) + got = np.concatenate([c.get_offsets() for c in fig.axes[0].collections + if isinstance(c, PathCollection)]) + plt.close(fig) + return got + + assert np.array_equal(offsets(), offsets()) + assert not np.array_equal(offsets(jitter_seed=1), offsets(jitter_seed=2)) + + +def test_a_constant_group_still_appears(obj): + """A group whose values are all identical has no density. It must not vanish + from the panel — that would read as "no cells" rather than "no spread".""" + plt = pytest.importorskip("matplotlib.pyplot") + obj.meta_data["const"] = ["only"] * obj.meta_data.shape[0] + fig = tc.vln_plot(obj, ["G11"], group_by="const", pt_size=0) + ax = fig.axes[0] + assert _violin_paths(ax) or ax.lines, "constant group drew nothing" + plt.close(fig) diff --git a/truecell/plotting.py b/truecell/plotting.py index a12fccb..d8512c0 100644 --- a/truecell/plotting.py +++ b/truecell/plotting.py @@ -484,6 +484,47 @@ def _subplot_grid(n: int, ncol: Optional[int] = None): return nrow, ncol +def _bw_nrd0(x: np.ndarray) -> float: + """R's ``bw.nrd0`` — the bandwidth ``stats::density`` uses, and so ggplot. + + Not the same as scipy's ``"scott"``, and the difference is large on exactly + the data this package plots. Scott's rule scales the sample standard + deviation; ``nrd0`` takes ``min(sd, IQR/1.349)``, and expression is + zero-inflated, so the IQR term is much the smaller of the two. Measured + against R on zero-inflated draws, scipy's bandwidth comes out 2.0-2.5x + wider, which visibly over-smooths a violin — the low-expressing groups lose + the spike at zero that is the whole shape of the distribution. + """ + x = np.asarray(x, dtype=float) + if x.size < 2: + raise ValueError("need at least 2 data points") + hi = float(np.std(x, ddof=1)) + q75, q25 = np.percentile(x, [75, 25]) + # 1.34, which is `bw.nrd0`'s divisor. `bw.nrd` — a different rule R also + # ships — uses 1.349, and taking that one instead is a 0.67% error in every + # bandwidth where the IQR term wins, which on zero-inflated data is most. + lo = min(hi, float(q75 - q25) / 1.34) + # R's fallback chain when the spread is exactly zero, in its own order. + if lo == 0: + lo = hi or abs(float(x[0])) or 1.0 + return 0.9 * lo * x.size ** -0.2 + + +# ggplot sizes a point by diameter in millimetres; matplotlib's `s` is area in +# points squared. 1 mm is 72.27/25.4 points, so a ggplot size of s becomes +# pi * (s * 2.845 / 2)**2 in matplotlib's units. +_GGPLOT_SIZE_TO_MPL_AREA = np.pi * (72.27 / 25.4 / 2) ** 2 + + +def _auto_point_size(n_cells: int) -> float: + """Seurat's ``AutoPointSize``, converted to matplotlib's area units. + + R: ``min(1583 / n, 1)``. The cap means small objects draw full-size points + and large ones shrink, so a 100k-cell violin does not become a solid block. + """ + return float(min(1583.0 / max(n_cells, 1), 1.0) ** 2 * _GGPLOT_SIZE_TO_MPL_AREA) + + # --------------------------------------------------------------------------- # 1. vln_plot — VlnPlot # --------------------------------------------------------------------------- @@ -494,25 +535,41 @@ def vln_plot( group_by: Optional[str] = None, assay: Optional[str] = None, layer: Optional[str] = None, - pt_size: float = 0.0, + pt_size: Optional[float] = None, ncol: Optional[int] = None, figsize: Optional[tuple] = None, palette: Optional[list] = None, + violin_width: float = 0.8, + jitter_seed: Optional[int] = 0, + raster: Optional[bool] = None, ) -> "Figure": """Violin plot of feature expression per cluster/identity. - Mirrors R's ``VlnPlot(pbmc, features = c("LYZ", "CD3D"))``. + Mirrors R's ``VlnPlot(pbmc, features = c("LYZ", "CD3D"))``, including the + three things that make a Seurat violin the shape it is: the density is + trimmed to the observed range, smoothed with R's ``nrd0`` bandwidth, and + scaled so every group reaches the same maximum width. Parameters ---------- obj : Truecell object features : gene name(s) or metadata column(s) group_by : metadata column used for grouping (default: active idents) - pt_size : size of individual data points overlaid on violins (0 = none) + pt_size : marker area for the jittered points, in matplotlib's units. + ``None`` (default) follows Seurat's ``AutoPointSize``, which + shows points and shrinks them as the cell count grows; ``0`` + omits them. ncol : number of columns in subplot grid figsize : figure size in inches; auto-computed if None palette : list of colours per group + violin_width : width of a full violin in x-axis units + jitter_seed : seed for the point jitter, so a figure redraws identically. + ``None`` draws fresh jitter each call. + raster : rasterise the jittered points. ``None`` follows the same + 100,000-point rule as the other plots. """ + from scipy.stats import gaussian_kde + plt = _mpl() if isinstance(features, str): features = [features] @@ -520,6 +577,8 @@ def vln_plot( groups = _get_groups(obj, group_by) unique = sorted(set(groups), key=lambda x: (int(x) if x.isdigit() else x)) colors = palette or _palette(len(unique)) + if pt_size is None: + pt_size = _auto_point_size(len(groups)) nrow, nc = _subplot_grid(len(features), ncol) if figsize is None: @@ -528,24 +587,54 @@ def vln_plot( fig, axes = plt.subplots(nrow, nc, figsize=figsize, squeeze=False) axes_flat = axes.flatten() + rng = np.random.default_rng(jitter_seed) + for i, feat in enumerate(features): ax = axes_flat[i] expr = _get_expression(obj, feat, assay, layer) - grp_data = [expr[groups == g] for g in unique] - parts = ax.violinplot(grp_data, positions=range(len(unique)), - showmedians=True, showextrema=False) - for j, pc in enumerate(parts["bodies"]): - pc.set_facecolor(colors[j]) - pc.set_alpha(0.8) - parts["cmedians"].set_color("black") - parts["cmedians"].set_linewidth(1.5) + for j, g in enumerate(unique): + vals = expr[groups == g] + if vals.size < 2: + continue + lo, hi = float(vals.min()), float(vals.max()) + if hi - lo < 1e-12: + # Constant within the group: a density is undefined, but the + # group still exists. Draw a flat marker at the value so it is + # not silently missing from the panel. + ax.plot([j - 0.35, j + 0.35], [lo, lo], + color=colors[j], linewidth=1.5, zorder=2) + continue + + # trim = TRUE: the support is the observed range, so the violin does + # not tail off below zero where expression cannot go. + grid = np.linspace(lo, hi, 512) + bw = _bw_nrd0(vals) + sd = float(np.std(vals, ddof=1)) + # scipy takes a factor multiplying the sample sd, not a bandwidth. + dens = gaussian_kde(vals, bw_method=bw / sd if sd > 0 else None)(grid) + + # scale = "width": every violin reaches the same maximum width, so + # groups are comparable in shape rather than in cell count. + peak = dens.max() + if peak <= 0: + continue + half = dens / peak * (violin_width / 2) + ax.fill_betweenx(grid, j - half, j + half, + facecolor=colors[j], alpha=0.8, linewidth=0, + zorder=2) + med = float(np.median(vals)) + w_at_med = float(np.interp(med, grid, half)) + ax.plot([j - w_at_med, j + w_at_med], [med, med], + color="black", linewidth=1.5, zorder=4) if pt_size > 0: for j, g in enumerate(unique): - jitter = np.random.uniform(-0.2, 0.2, size=(groups == g).sum()) - ax.scatter(j + jitter, expr[groups == g], - s=pt_size, alpha=0.4, color=colors[j], zorder=3) + vals = expr[groups == g] + jitter = rng.uniform(-0.2, 0.2, size=vals.size) + ax.scatter(j + jitter, vals, s=pt_size, alpha=0.4, + color=colors[j], zorder=3, + rasterized=_should_raster(raster, len(expr))) ax.set_xticks(range(len(unique))) ax.set_xticklabels(unique, rotation=45 if len(unique) > 6 else 0, diff --git a/tutorials/figures/01_qc_violin.png b/tutorials/figures/01_qc_violin.png index e0655bd..55e436b 100644 Binary files a/tutorials/figures/01_qc_violin.png and b/tutorials/figures/01_qc_violin.png differ diff --git a/tutorials/figures/09_marker_violins.png b/tutorials/figures/09_marker_violins.png index 21f0182..278639c 100644 Binary files a/tutorials/figures/09_marker_violins.png and b/tutorials/figures/09_marker_violins.png differ diff --git a/tutorials/figures/09b_marker_violins_counts.png b/tutorials/figures/09b_marker_violins_counts.png index 2d84757..c04f1c6 100644 Binary files a/tutorials/figures/09b_marker_violins_counts.png and b/tutorials/figures/09b_marker_violins_counts.png differ diff --git a/tutorials/figures_advanced/01_qc_violin.png b/tutorials/figures_advanced/01_qc_violin.png index 4143bc3..0e02865 100644 Binary files a/tutorials/figures_advanced/01_qc_violin.png and b/tutorials/figures_advanced/01_qc_violin.png differ diff --git a/tutorials/figures_advanced/10_tnk_subset_violins.png b/tutorials/figures_advanced/10_tnk_subset_violins.png index ca6a699..bf93e22 100644 Binary files a/tutorials/figures_advanced/10_tnk_subset_violins.png and b/tutorials/figures_advanced/10_tnk_subset_violins.png differ diff --git a/tutorials/figures_hashing/py_03_ncount_violin.png b/tutorials/figures_hashing/py_03_ncount_violin.png index aff2952..f85b0c1 100644 Binary files a/tutorials/figures_hashing/py_03_ncount_violin.png and b/tutorials/figures_hashing/py_03_ncount_violin.png differ diff --git a/tutorials/figures_multimodal/10_adt_weight_by_celltype.png b/tutorials/figures_multimodal/10_adt_weight_by_celltype.png index db26640..975da20 100644 Binary files a/tutorials/figures_multimodal/10_adt_weight_by_celltype.png and b/tutorials/figures_multimodal/10_adt_weight_by_celltype.png differ diff --git a/tutorials/figures_sctransform/05_sct_violins.png b/tutorials/figures_sctransform/05_sct_violins.png index 04b9db3..bd3a872 100644 Binary files a/tutorials/figures_sctransform/05_sct_violins.png and b/tutorials/figures_sctransform/05_sct_violins.png differ diff --git a/tutorials/figures_spatial/01_qc_violin.png b/tutorials/figures_spatial/01_qc_violin.png index 401d318..d6323ab 100644 Binary files a/tutorials/figures_spatial/01_qc_violin.png and b/tutorials/figures_spatial/01_qc_violin.png differ