Plotting

Every function here returns a matplotlib Figure, so you save it, compose it, or let a notebook display it:

fig = truecell.dim_plot(pbmc, reduction="umap", label=True)
fig.savefig("umap.png", dpi=150, bbox_inches="tight")

matplotlib is optional — it lives in the analysis extra, and import truecell works without it. That is why the Figure return annotation is an if TYPE_CHECKING: import: it documents the return type without making the import mandatory.

Theme

Every plot on this page scales its text from one base size and takes its group colours from one palette. Set them once instead of passing the same overrides to each call:

import truecell as tc

tc.set_theme(base_size=13, style="seurat")   # bigger text, cowplot's look
fig = tc.dim_plot(pbmc)

with tc.theme_context(base_size=8):          # scoped to the block
    panel = tc.feature_plot(pbmc, ["MS4A1", "LYZ"])

The default (base_size=10, no style) leaves matplotlib's rcParams untouched and reproduces the sizes the module used before the theme existed.

set_theme

set_theme(base_size: Optional[float] = None, palette: Optional[list] = None, font: Optional[str] = None, dpi: Optional[float] = None, style: Optional[str] = None) -> dict

Set the look of every plot this module draws.

The Python counterpart of applying a ggplot theme once and having Seurat's plots pick it up, rather than passing the same overrides to every call.

Only the arguments you pass are changed; the rest keep their current values. Use reset_theme to go back to the defaults, or theme_context to scope a change to a with block.

Parameters:

• base_size (Optional[float], default: None ) –

  float, optional Base font size in points (default 10.0). Every text element scales from it, so this is the one knob for "make the labels bigger" — a figure destined for a slide wants 13-14, a dense multi-panel figure wants 8.

• palette (Optional[list], default: None ) –

  list, optional Categorical colours to use instead of ggplot's hue_pal. Ignored for any plot with more groups than the list has entries, which falls back to hue_pal rather than repeating a colour.

• font (Optional[str], default: None ) –

  str, optional matplotlib font family, e.g. "Helvetica", "DejaVu Sans".

• dpi (Optional[float], default: None ) –

  float, optional Figure and savefig dpi.

• style (Optional[str], default: None ) –

  str, optional "seurat" (cowplot's look — no grid, black spines) or "minimal" (faint grid, grey spines). None leaves matplotlib's rcParams untouched.

Returns:

• dict –

  The theme as applied, so a caller can stash it and restore it later.

Examples:

>>> import truecell as tc
>>> tc.set_theme(base_size=13, style="seurat")
>>> tc.dim_plot(pbmc)

Source code in truecell/plotting.py

def set_theme(
    base_size: Optional[float] = None,
    palette: Optional[list] = None,
    font: Optional[str] = None,
    dpi: Optional[float] = None,
    style: Optional[str] = None,
) -> dict:
    """Set the look of every plot this module draws.

    The Python counterpart of applying a ggplot theme once and having Seurat's
    plots pick it up, rather than passing the same overrides to every call.

    Only the arguments you pass are changed; the rest keep their current values.
    Use :func:`reset_theme` to go back to the defaults, or :func:`theme_context`
    to scope a change to a ``with`` block.

    Parameters
    ----------
    base_size : float, optional
        Base font size in points (default 10.0). Every text element scales from
        it, so this is the one knob for "make the labels bigger" — a figure
        destined for a slide wants 13-14, a dense multi-panel figure wants 8.
    palette : list, optional
        Categorical colours to use instead of ggplot's ``hue_pal``. Ignored for
        any plot with more groups than the list has entries, which falls back to
        ``hue_pal`` rather than repeating a colour.
    font : str, optional
        matplotlib font family, e.g. ``"Helvetica"``, ``"DejaVu Sans"``.
    dpi : float, optional
        Figure and savefig dpi.
    style : str, optional
        ``"seurat"`` (cowplot's look — no grid, black spines) or ``"minimal"``
        (faint grid, grey spines). ``None`` leaves matplotlib's rcParams
        untouched.

    Returns
    -------
    dict
        The theme as applied, so a caller can stash it and restore it later.

    Examples
    --------
    >>> import truecell as tc
    >>> tc.set_theme(base_size=13, style="seurat")   # doctest: +SKIP
    >>> tc.dim_plot(pbmc)                            # doctest: +SKIP
    """
    if base_size is not None:
        if base_size <= 0:
            raise ValueError(f"base_size must be positive, got {base_size!r}")
        _THEME["base_size"] = float(base_size)
    if palette is not None:
        _THEME["palette"] = list(palette)
    if font is not None:
        _THEME["font"] = font
    if dpi is not None:
        _THEME["dpi"] = float(dpi)
    if style is not None:
        if style not in _STYLES:
            raise ValueError(
                f"Unknown style {style!r}. Available: {sorted(_STYLES)}"
            )
        _THEME["style"] = style
    _apply_rc()
    return dict(_THEME)

theme_context

theme_context(**kwargs)

Apply a theme for the duration of a with block, then restore it.

with theme_context(base_size=14, style="seurat"): # doctest: +SKIP ... fig = dim_plot(pbmc)

Source code in truecell/plotting.py

@contextmanager
def theme_context(**kwargs):
    """Apply a theme for the duration of a ``with`` block, then restore it.

    >>> with theme_context(base_size=14, style="seurat"):   # doctest: +SKIP
    ...     fig = dim_plot(pbmc)
    """
    previous = dict(_THEME)
    try:
        plt = _mpl()
        saved_rc = dict(plt.rcParams)
    except ImportError:
        saved_rc = None
    try:
        set_theme(**kwargs)
        yield get_theme()
    finally:
        _THEME.update(previous)
        if saved_rc is not None:
            # `rcParams.update` rather than assignment: a few keys are
            # read-only derived entries that raise if written back.
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                _mpl().rcParams.update(saved_rc)

get_theme

get_theme() -> dict

Return a copy of the active theme.

Source code in truecell/plotting.py

def get_theme() -> dict:
    """Return a copy of the active theme."""
    return dict(_THEME)

reset_theme

reset_theme() -> dict

Restore the default theme and matplotlib's own rcParams.

Source code in truecell/plotting.py

def reset_theme() -> dict:
    """Restore the default theme and matplotlib's own rcParams."""
    _THEME.update(_THEME_DEFAULTS)
    _THEME["palette"] = None
    try:
        plt = _mpl()
    except ImportError:
        return dict(_THEME)
    plt.rcdefaults()
    return dict(_THEME)

Colour

hue_pal

hue_pal(n: int, h: tuple = (15, 375), c: float = 100, lum: float = 65) -> list[str]

n evenly-spaced HCL hues — ggplot2's default discrete colour scale.

The Python counterpart of scales::hue_pal()(n), which is what Seurat's DimPlot/VlnPlot/DoHeatmap colour groups with. Verified exact against R for n = 1-6, 8 and 9 in tests/test_plotting_theme.py.

Computing the ramp beats carrying a fixed list, for two reasons. The palette depends on n — ggplot spreads the same hue circle across however many groups there are, so the colours for 9 clusters are not the colours for 8 plus one more — and a fixed list therefore only matches Seurat at the one length it was copied from. It also cannot run out or repeat.

Source code in truecell/plotting.py

def hue_pal(n: int, h: tuple = (15, 375), c: float = 100, lum: float = 65) -> list[str]:
    """``n`` evenly-spaced HCL hues — ggplot2's default discrete colour scale.

    The Python counterpart of ``scales::hue_pal()(n)``, which is what Seurat's
    ``DimPlot``/``VlnPlot``/``DoHeatmap`` colour groups with. Verified exact
    against R for n = 1-6, 8 and 9 in ``tests/test_plotting_theme.py``.

    Computing the ramp beats carrying a fixed list, for two reasons. The palette
    depends on ``n`` — ggplot spreads the same hue circle across however many
    groups there are, so the colours for 9 clusters are *not* the colours for 8
    plus one more — and a fixed list therefore only matches Seurat at the one
    length it was copied from. It also cannot run out or repeat.
    """
    if n < 1:
        return []
    return [_hcl_to_hex(x, c, lum) for x in np.linspace(h[0], h[1], n + 1)[:n]]

Plots

dim_plot

dim_plot(obj, reduction: str = 'umap', group_by: Optional[str] = None, label: bool = True, label_size: Optional[float] = None, pt_size: float = 4.0, alpha: float = 0.7, figsize: tuple = (7, 6), palette: Optional[list] = None, title: Optional[str] = None, raster: Optional[bool] = None, split_by: Optional[str] = None, ncol: Optional[int] = None, repel: bool = False) -> 'Figure'

Plot cells in a reduced-dimension embedding coloured by identity.

Mirrors R's DimPlot(pbmc, reduction = "umap", label = TRUE).

Parameters:

• reduction (str, default: 'umap' ) –

  which reduction to use ("umap", "pca", …)

• group_by (Optional[str], default: None ) –

  metadata column for colouring (default: active idents). A categorical column is drawn in its category order, as R draws a factor's levels.

• label (bool, default: True ) –

  label each group at the median of its cells' coordinates, in each panel, as Seurat's LabelClusters does

• label_size (Optional[float], default: None ) –

  font size for the group labels (default: scales with the theme)

• repel (bool, default: False ) –

  move the labels so that none overlaps another, runs off the panel, or covers another group's median or more than half of any group's cells, keeping each as close to its median as that allows. A line joins a label to its group when it has moved away, going round the other labels where there is room. Seurat's repel does this with ggrepel, which starts from random jitter; this placement is deterministic. It is worked out again whenever the figure is drawn, so it still holds after the figure is resized or saved at another size.

• pt_size (float, default: 4.0 ) –

  scatter point size

• alpha (float, default: 0.7 ) –

  point transparency

• raster (Optional[bool], default: None ) –

  draw the cells as a raster layer instead of vector paths. None (default) rasterises above 100,000 cells, the same rule and threshold Seurat uses.

• split_by (Optional[str], default: None ) –

  metadata column to facet on — one panel per level, holding only that level's cells. Colours and axis limits are shared across panels, matching facet_wrap's fixed scales: the panels are only comparable if a position and a colour mean the same thing in each.

• ncol (Optional[int], default: None ) –

  panels per row when split_by is set.

Source code in truecell/plotting.py

def dim_plot(
    obj,
    reduction: str = "umap",
    group_by: Optional[str] = None,
    label: bool = True,
    label_size: Optional[float] = None,
    pt_size: float = 4.0,
    alpha: float = 0.7,
    figsize: tuple = (7, 6),
    palette: Optional[list] = None,
    title: Optional[str] = None,
    raster: Optional[bool] = None,
    split_by: Optional[str] = None,
    ncol: Optional[int] = None,
    repel: bool = False,
) -> "Figure":
    """Plot cells in a reduced-dimension embedding coloured by identity.

    Mirrors R's ``DimPlot(pbmc, reduction = "umap", label = TRUE)``.

    Parameters
    ----------
    reduction : which reduction to use ("umap", "pca", …)
    group_by  : metadata column for colouring (default: active idents). A categorical
                column is drawn in its category order, as R draws a factor's levels.
    label     : label each group at the median of its cells' coordinates, in
                each panel, as Seurat's ``LabelClusters`` does
    label_size: font size for the group labels (default: scales with the theme)
    repel     : move the labels so that none overlaps another, runs off the
                panel, or covers another group's median or more than half of
                any group's cells, keeping each as close to its median as that
                allows. A line joins a label to its group when it has moved
                away, going round the other labels where there is room.
                Seurat's ``repel`` does this with ggrepel, which starts from
                random jitter; this placement is deterministic. It is worked
                out again whenever the figure is drawn, so it still holds after
                the figure is resized or saved at another size.
    pt_size   : scatter point size
    alpha     : point transparency
    raster    : draw the cells as a raster layer instead of vector paths.
                ``None`` (default) rasterises above 100,000 cells, the same rule
                and threshold Seurat uses.
    split_by  : metadata column to facet on — one panel per level, holding only
                that level's cells. Colours and axis limits are shared across
                panels, matching ``facet_wrap``'s fixed scales: the panels are
                only comparable if a position and a colour mean the same thing
                in each.
    ncol      : panels per row when ``split_by`` is set.
    """
    plt = _mpl()
    emb = _get_embedding(obj, reduction)
    groups = _get_groups(obj, group_by)
    unique = _group_levels(obj, group_by, groups)
    colors = palette or _palette(len(unique))
    # On the whole embedding, not the per-group slice: the decision is about the
    # figure's total path count, and a 200k-cell object split across 20 clusters
    # would otherwise never cross the threshold in any single call below.
    rast = _should_raster(raster, len(emb))

    ax1_label, ax2_label = _axis_labels(obj, reduction)
    colour_of = dict(zip(unique, colors))

    if split_by is None:
        panels: list = [(None, np.ones(len(emb), dtype=bool))]
        fig, axes = plt.subplots(figsize=figsize, squeeze=False)
        axes_flat = axes.flatten()
        xlim = ylim = None
    else:
        split_labels, levels = _split_levels(obj, split_by)
        panels = [(lv, split_labels == lv) for lv in levels]
        nrow, nc = _subplot_grid(len(levels), ncol)
        if figsize == (7, 6):  # untouched default — size to the grid instead
            figsize = (nc * 4.6, nrow * 4.2)
        fig, axes = plt.subplots(nrow, nc, figsize=figsize, squeeze=False)
        axes_flat = axes.flatten()
        xlim, ylim = _shared_limits(emb)

    for i, (level, in_panel) in enumerate(panels):
        ax = axes_flat[i]
        for g in unique:
            mask = (groups == g) & in_panel
            # Draw even when empty, so the legend carries every group in every
            # panel and the colours stay attached to the same labels.
            ax.scatter(emb[mask, 0], emb[mask, 1], s=pt_size, alpha=alpha,
                       color=colour_of[g], label=g if i == 0 else None,
                       linewidths=0, rasterized=rast)

        if label:
            texts, anchors, members = [], [], []
            for g in unique:
                mask = (groups == g) & in_panel
                if not mask.any():
                    continue
                # The median, as LabelClusters takes it. A mean is pulled towards
                # a group's outlying cells, off the cells the label names.
                cx, cy = np.nanmedian(emb[mask], axis=0)
                texts.append(ax.text(
                    cx, cy, g,
                    fontsize=_fs("label") if label_size is None else label_size,
                    fontweight="bold", ha="center", va="center",
                    bbox=dict(boxstyle="round,pad=0.25", fc="white",
                              alpha=0.75, ec="none")))
                anchors.append((cx, cy))
                members.append(emb[mask])
            if repel and texts:
                # Imported here: it subclasses a matplotlib artist, and
                # matplotlib is optional for the package as a whole.
                from ._repel import attach
                attach(ax, texts, anchors, members)

        if xlim is not None:
            ax.set_xlim(*xlim)
            ax.set_ylim(*ylim)
        ax.set_xlabel(ax1_label)
        ax.set_ylabel(ax2_label)
        if level is None:
            ax.set_title(title or f"{reduction.upper()} — coloured by "
                                  f"{'ident' if group_by is None else group_by}",
                         fontsize=_fs("large"))
        else:
            ax.set_title(str(level), fontsize=_fs("title"), fontweight="bold")
        _strip_axes(ax)

    for i in range(len(panels), len(axes_flat)):
        axes_flat[i].set_visible(False)

    # One legend for the figure. Anchoring it to the first panel would put it
    # inside the grid; on the figure it sits clear of every panel.
    handles, labels_ = axes_flat[0].get_legend_handles_labels()
    if split_by is None:
        axes_flat[0].legend(handles, labels_, markerscale=2.5,
                            fontsize=_fs("small"), bbox_to_anchor=(1.01, 1),
                            loc="upper left", frameon=False)
        fig.tight_layout()
    else:
        if title:
            fig.suptitle(title, fontsize=_fs("suptitle"), fontweight="bold")
        fig.tight_layout()
        fig.legend(handles, labels_, markerscale=2.5, fontsize=_fs("small"),
                   loc="center left", bbox_to_anchor=(1.0, 0.5), frameon=False)
    return fig

feature_plot

feature_plot(obj, features: Union[str, list[str]], reduction: str = 'umap', assay: Optional[str] = None, layer: Optional[str] = None, ncol: Optional[int] = None, order: bool = True, min_cutoff: float | str | None = None, max_cutoff: float | str | None = None, colormap: str = 'YlOrRd', pt_size: float = 3.0, figsize: Optional[tuple] = None, raster: Optional[bool] = None, split_by: Optional[str] = None) -> 'Figure'

Visualise feature expression on a dimensionality reduction embedding.

Mirrors R's FeaturePlot(pbmc, features = c("MS4A1", "LYZ")).

Parameters:

• features (Union[str, list[str]]) –

  gene name(s) or metadata column(s) to plot

• reduction (str, default: 'umap' ) –

  which reduction to use ("umap", "pca", …)

• order (bool, default: True ) –

  plot cells with highest expression on top

• min_cutoff (float | str | None, default: None ) –

  expression at or below this is drawn as not expressed (grey). A number, or a quantile such as "q05"; default 0.

• max_cutoff (float | str | None, default: None ) –

  expression above this takes the top of the colour scale. A number, or a quantile such as "q95"; default the feature's maximum. As in Seurat, a quantile is taken over the cells that express the feature (values above zero), so "q05"/"q95" suit a gene detected in a few percent of cells. If the two cutoffs leave the scale no range, a warning says so and the feature is drawn over its own range instead.

• colormap (str, default: 'YlOrRd' ) –

  matplotlib colormap name for expression

• pt_size (float, default: 3.0 ) –

  scatter point size

• raster (Optional[bool], default: None ) –

  draw the cells as a raster layer instead of vector paths. None (default) rasterises above 100,000 cells, the same rule and threshold Seurat uses.

• split_by (Optional[str], default: None ) –

  metadata column to facet on. Lays out one row per feature and one column per level, as Seurat does. The colour scale and the axis limits are computed once per feature over all cells and shared down the row, so a colour means the same expression in every panel — per-panel scales would make the levels look alike however different they are.

Source code in truecell/plotting.py

def feature_plot(
    obj,
    features: Union[str, list[str]],
    reduction: str = "umap",
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    ncol: Optional[int] = None,
    order: bool = True,
    min_cutoff: float | str | None = None,
    max_cutoff: float | str | None = None,
    colormap: str = "YlOrRd",
    pt_size: float = 3.0,
    figsize: Optional[tuple] = None,
    raster: Optional[bool] = None,
    split_by: Optional[str] = None,
) -> "Figure":
    """Visualise feature expression on a dimensionality reduction embedding.

    Mirrors R's ``FeaturePlot(pbmc, features = c("MS4A1", "LYZ"))``.

    Parameters
    ----------
    features    : gene name(s) or metadata column(s) to plot
    reduction   : which reduction to use ("umap", "pca", …)
    order       : plot cells with highest expression on top
    min_cutoff  : expression at or below this is drawn as not expressed (grey). A
                  number, or a quantile such as ``"q05"``; default 0.
    max_cutoff  : expression above this takes the top of the colour scale. A number,
                  or a quantile such as ``"q95"``; default the feature's maximum.
                  As in Seurat, a quantile is taken over the cells that express the
                  feature (values above zero), so ``"q05"``/``"q95"`` suit a gene
                  detected in a few percent of cells. If the two cutoffs leave the
                  scale no range, a warning says so and the feature is drawn over
                  its own range instead.
    colormap    : matplotlib colormap name for expression
    pt_size     : scatter point size
    raster      : draw the cells as a raster layer instead of vector paths.
                  ``None`` (default) rasterises above 100,000 cells, the same
                  rule and threshold Seurat uses.
    split_by    : metadata column to facet on. Lays out one row per feature and
                  one column per level, as Seurat does. The colour scale and the
                  axis limits are computed once per feature over *all* cells and
                  shared down the row, so a colour means the same expression in
                  every panel — per-panel scales would make the levels look
                  alike however different they are.
    """
    plt = _mpl()
    if isinstance(features, str):
        features = [features]

    emb = _get_embedding(obj, reduction)
    rast = _should_raster(raster, len(emb))
    ax1_label, ax2_label = _axis_labels(obj, reduction)

    if split_by is None:
        nrow, nc = _subplot_grid(len(features), ncol)
        if figsize is None:
            figsize = (nc * 4.5, nrow * 4)
        fig, axes = plt.subplots(nrow, nc, figsize=figsize, squeeze=False)
        axes_flat = axes.flatten()
        cells = [(None, np.ones(len(emb), dtype=bool))]
        xlim = ylim = None
    else:
        split_labels, levels = _split_levels(obj, split_by)
        cells = [(lv, split_labels == lv) for lv in levels]
        nrow, nc = len(features), len(levels)
        if figsize is None:
            figsize = (nc * 4.2, nrow * 3.9)
        fig, axes = plt.subplots(nrow, nc, figsize=figsize, squeeze=False)
        axes_flat = axes.flatten()
        xlim, ylim = _shared_limits(emb)

    for fi, feat in enumerate(features):
        expr = _get_expression(obj, feat, assay, layer)

        # Cutoffs from every cell, not the panel's subset — that is what makes
        # the row comparable.
        vmin = 0.0 if min_cutoff is None else _set_quantile(min_cutoff, expr)
        vmax = float(expr.max()) if max_cutoff is None else _set_quantile(max_cutoff, expr)
        if not vmax > vmin and expr.max() > 0:
            # The feature varies, so it is the cutoffs that left the scale no range:
            # one colour for every cell, and nothing to tell the plot is wrong.
            warnings.warn(
                f"feature_plot: for {feat!r}, min_cutoff={min_cutoff!r} and "
                f"max_cutoff={max_cutoff!r} give {vmin:.3g} and {vmax:.3g}, which leaves "
                f"the colour scale no range. Drawing it over the feature's own range, "
                f"0 to {expr.max():.3g}, instead.",
                stacklevel=2,
            )
            vmin, vmax = 0.0, float(expr.max())
        vmax = max(vmax, vmin + 1e-9)

        for ci, (level, in_panel) in enumerate(cells):
            ax = axes_flat[fi * len(cells) + ci]

            # R Seurat default: non-expressing (zero) cells rendered in light gray,
            # expressing cells drawn on top sorted by expression level.
            zero_mask = (expr <= vmin) & in_panel
            nonzero_mask = (expr > vmin) & in_panel
            ax.scatter(emb[zero_mask, 0], emb[zero_mask, 1],
                       c="#D3D3D3", s=pt_size, linewidths=0, rasterized=rast)

            expr_nz = expr[nonzero_mask]
            emb_nz = emb[nonzero_mask]
            if order:
                sort_idx = np.argsort(expr_nz)
                expr_nz = expr_nz[sort_idx]
                emb_nz = emb_nz[sort_idx]

            sc = ax.scatter(emb_nz[:, 0], emb_nz[:, 1], c=expr_nz, s=pt_size,
                            cmap=colormap, vmin=vmin, vmax=vmax, linewidths=0,
                            rasterized=rast)
            # One colourbar per row when split: it applies to the whole row, and
            # repeating it per panel would say otherwise.
            if level is None or ci == len(cells) - 1:
                plt.colorbar(sc, ax=ax, shrink=0.6, pad=0.01, aspect=25)
            if xlim is not None:
                ax.set_xlim(*xlim)
                ax.set_ylim(*ylim)
            ax.set_xlabel(ax1_label, fontsize=_fs("small"))
            ax.set_ylabel(ax2_label, fontsize=_fs("small"))
            ax.set_title(feat if level is None else f"{feat} — {level}",
                         fontsize=_fs("title"), fontweight="bold")
            ax.set_xticks([])
            ax.set_yticks([])
            for spine in ax.spines.values():
                spine.set_visible(False)

    for i in range(len(features) * len(cells), len(axes_flat)):
        axes_flat[i].set_visible(False)

    fig.tight_layout()
    return fig

vln_plot

vln_plot(obj, features: Union[str, list[str]], group_by: Optional[str] = None, assay: Optional[str] = None, layer: Optional[str] = None, 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, split_by: Optional[str] = None) -> 'Figure'

Violin plot of feature expression per cluster/identity.

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 (Union[str, list[str]]) –

  gene name(s) or metadata column(s)

• group_by (Optional[str], default: None ) –

  metadata column used for grouping (default: active idents). A categorical column is drawn in its category order, as R draws a factor's levels.

• pt_size (Optional[float], default: None ) –

  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 (Optional[int], default: None ) –

  number of columns in subplot grid

• figsize (Optional[tuple], default: None ) –

  figure size in inches; auto-computed if None

• palette (Optional[list], default: None ) –

  list of colours per group

• violin_width (float, default: 0.8 ) –

  width of a full violin in x-axis units

• jitter_seed (Optional[int], default: 0 ) –

  seed for the point jitter, so a figure redraws identically. None draws fresh jitter each call.

• raster (Optional[bool], default: None ) –

  rasterise the jittered points. None follows the same 100,000-point rule as the other plots.

• split_by (Optional[str], default: None ) –

  metadata column to split each group by. Matches Seurat's split.plot = FALSE default — the levels are dodged side-by-side within each group's x position and coloured by level, not drawn as one split violin and not faceted into separate panels.

Source code in truecell/plotting.py

def vln_plot(
    obj,
    features: Union[str, list[str]],
    group_by: Optional[str] = None,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    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,
    split_by: Optional[str] = None,
) -> "Figure":
    """Violin plot of feature expression per cluster/identity.

    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). A categorical
               column is drawn in its category order, as R draws a factor's levels.
    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.
    split_by : metadata column to split each group by. Matches Seurat's
               ``split.plot = FALSE`` default — the levels are *dodged*
               side-by-side within each group's x position and coloured by
               level, not drawn as one split violin and not faceted into
               separate panels.
    """
    from scipy.stats import gaussian_kde

    plt = _mpl()
    if isinstance(features, str):
        features = [features]

    groups = _get_groups(obj, group_by)
    unique = _group_levels(obj, group_by, groups)
    if pt_size is None:
        pt_size = _auto_point_size(len(groups))

    # Seurat's `split.by` on a violin is a *dodge*, not a facet and not a split
    # violin: `VlnPlot` has `split.plot = FALSE` by default, and its own message
    # says "Separate violin plots are now plotted side-by-side". So each ident
    # keeps one x position and the levels sit shoulder to shoulder within it,
    # coloured by level rather than by ident.
    if split_by is None:
        sub_levels: list = [None]
        split_labels = None
        colors = palette or _palette(len(unique))
    else:
        split_labels, sub_levels = _split_levels(obj, split_by)
        colors = palette or _palette(len(sub_levels))

    n_sub = len(sub_levels)
    slot = violin_width / n_sub  # ggplot dodges within one group's width

    nrow, nc = _subplot_grid(len(features), ncol)
    if figsize is None:
        figsize = (max(5, nc * 4 * max(1, n_sub / 2)), nrow * 3.5)

    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)

        for j, g in enumerate(unique):
          for k, lv in enumerate(sub_levels):
            in_cell = (groups == g)
            if lv is not None:
                in_cell = in_cell & (split_labels == lv)
            vals = expr[in_cell]
            colour = colors[j if lv is None else k]
            # Centre of this violin's slot within the group's width.
            centre = j + (0 if lv is None else (k - (n_sub - 1) / 2) * slot)
            body = violin_width if lv is None else slot * 0.9

            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.
                # 0.4375 of the body, which is exactly +/-0.35 at the default
                # width of 0.8 — the constant this used before it had to scale
                # with a dodge slot. Keeping it identical matters: a cluster with
                # no variance in a `counts` layer hits this branch, and the two
                # tutorial figures that draw one would otherwise shift by 0.002.
                ax.plot([centre - body * 0.4375, centre + body * 0.4375],
                        [lo, lo], color=colour, 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 * (body / 2)
            ax.fill_betweenx(grid, centre - half, centre + half,
                             facecolor=colour, alpha=0.8, linewidth=0,
                             zorder=2,
                             label=(str(lv) if (lv is not None and j == 0
                                                and i == 0) else None))
            med = float(np.median(vals))
            w_at_med = float(np.interp(med, grid, half))
            ax.plot([centre - w_at_med, centre + w_at_med], [med, med],
                    color="black", linewidth=1.5, zorder=4)

        if pt_size > 0:
            jitter_span = 0.2 * (1 if n_sub == 1 else slot / violin_width)
            for j, g in enumerate(unique):
                for k, lv in enumerate(sub_levels):
                    in_cell = (groups == g)
                    if lv is not None:
                        in_cell = in_cell & (split_labels == lv)
                    vals = expr[in_cell]
                    centre = j + (0 if lv is None
                                  else (k - (n_sub - 1) / 2) * slot)
                    jitter = rng.uniform(-jitter_span, jitter_span, size=vals.size)
                    ax.scatter(centre + jitter, vals, s=pt_size, alpha=0.4,
                               color=colors[j if lv is None else k], 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,
                           ha="right" if len(unique) > 6 else "center", fontsize=_fs("label"))
        ax.set_ylabel("Expression")
        ax.set_title(feat, fontsize=_fs("title"), fontweight="bold")
        _strip_axes(ax)

    for i in range(len(features), len(axes_flat)):
        axes_flat[i].set_visible(False)

    fig.tight_layout()
    if split_by is not None:
        handles, labels_ = axes_flat[0].get_legend_handles_labels()
        if handles:
            fig.legend(handles, labels_, title=split_by, fontsize=_fs("small"),
                       title_fontsize=_fs("small"), loc="center left",
                       bbox_to_anchor=(1.0, 0.5), frameon=False)
    return fig

ridge_plot

ridge_plot(obj, features: Union[str, list[str]], group_by: Optional[str] = None, assay: Optional[str] = None, layer: Optional[str] = None, ncol: Optional[int] = None, figsize: Optional[tuple] = None, palette: Optional[list] = None) -> 'Figure'

Ridgeline (joy) plots of feature expression per group.

Mirrors R's RidgePlot(pbmc, features = c("LYZ", "CD3D")).

Requires scipy for KDE smoothing.

Group labels sit outside the axes, so a panel needs more width than its axes alone. The default figsize allows for that. Passing one that is too small makes matplotlib draw the panels over each other instead of erroring, so an explicit figsize is checked and warns if the panels end up colliding.

Source code in truecell/plotting.py

def ridge_plot(
    obj,
    features: Union[str, list[str]],
    group_by: Optional[str] = None,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    ncol: Optional[int] = None,
    figsize: Optional[tuple] = None,
    palette: Optional[list] = None,
) -> "Figure":
    """Ridgeline (joy) plots of feature expression per group.

    Mirrors R's ``RidgePlot(pbmc, features = c("LYZ", "CD3D"))``.

    Requires ``scipy`` for KDE smoothing.

    Group labels sit outside the axes, so a panel needs more width than its axes
    alone. The default ``figsize`` allows for that. Passing one that is too small
    makes matplotlib draw the panels over each other instead of erroring, so an
    explicit ``figsize`` is checked and warns if the panels end up colliding.
    """
    from scipy.stats import gaussian_kde
    plt = _mpl()

    if isinstance(features, str):
        features = [features]

    groups = _get_groups(obj, group_by)
    unique = _group_levels(obj, group_by, groups)[::-1]
    colors = palette or _palette(len(unique))

    nrow, nc = _subplot_grid(len(features), ncol)
    user_figsize = figsize
    default_figsize = (nc * 5, nrow * max(3, len(unique) * 0.5))
    if figsize is None:
        figsize = default_figsize

    fig, axes = plt.subplots(nrow, nc, figsize=figsize, squeeze=False)
    axes_flat = axes.flatten()

    for i, feat in enumerate(features):
        ax = axes_flat[i]
        expr = _get_expression(obj, feat, assay, layer)
        xmin, xmax = expr.min(), expr.max()
        xs = np.linspace(xmin - 0.1, xmax + 0.1, 300)

        for j, (g, color) in enumerate(zip(unique, colors[::-1])):
            vals = expr[groups == g]
            if vals.std() < 1e-9 or len(vals) < 3:
                continue
            kde = gaussian_kde(vals, bw_method="scott")
            ys = kde(xs)
            ys = ys / ys.max() * 0.9  # normalise height
            ax.fill_between(xs, j + ys, j, alpha=0.8, color=color)
            ax.plot(xs, j + ys, color="white", linewidth=0.5)

        ax.set_yticks(range(len(unique)))
        # Row j holds the density for unique[j] (see loop above), so the tick
        # labels must be `unique`, not its reverse.
        ax.set_yticklabels(unique, fontsize=_fs("label"))
        ax.set_xlabel("Expression")
        ax.set_title(feat, fontsize=_fs("title"), fontweight="bold")
        _strip_axes(ax)

    for i in range(len(features), len(axes_flat)):
        axes_flat[i].set_visible(False)

    fig.tight_layout()
    if user_figsize is not None:
        _warn_if_panels_collide(fig, "ridge_plot", default_figsize)
    return fig

dot_plot

dot_plot(obj, features: Union[str, list[str]], group_by: Optional[str] = None, assay: Optional[str] = None, layer: str = 'data', cols: tuple = ('lightgrey', 'blue'), col_min: float = -2.5, col_max: float = 2.5, dot_min: float = 0.0, dot_scale: float = 6.0, scale: bool = True, figsize: Optional[tuple] = None) -> 'Figure'

Dot plot of feature expression across groups.

Mirrors R's DotPlot(pbmc, features = c("LYZ", "CD3D")). For each (feature, group) the dot size encodes the fraction of cells in the group expressing the feature (counts > 0) and the colour encodes the average expression (z-scored across groups when scale=True, matching Seurat).

Parameters:

• features (Union[str, list[str]]) –

  gene name(s) to plot (x-axis).

• group_by (Optional[str], default: None ) –

  metadata column for grouping (default: active idents, y-axis). A categorical column is drawn in its category order, as R draws a factor's levels.

• scale (bool, default: True ) –

  z-score each feature's average expression across groups.

• col_min (float, default: -2.5 ) –
• dot_min (float, default: 0.0 ) –
• dot_scale (float, default: 6.0 ) –

Source code in truecell/plotting.py

def dot_plot(
    obj,
    features: Union[str, list[str]],
    group_by: Optional[str] = None,
    assay: Optional[str] = None,
    layer: str = "data",
    cols: tuple = ("lightgrey", "blue"),
    col_min: float = -2.5,
    col_max: float = 2.5,
    dot_min: float = 0.0,
    dot_scale: float = 6.0,
    scale: bool = True,
    figsize: Optional[tuple] = None,
) -> "Figure":
    """Dot plot of feature expression across groups.

    Mirrors R's ``DotPlot(pbmc, features = c("LYZ", "CD3D"))``. For each
    (feature, group) the **dot size** encodes the fraction of cells in the group
    expressing the feature (counts > 0) and the **colour** encodes the average
    expression (z-scored across groups when ``scale=True``, matching Seurat).

    Parameters
    ----------
    features : gene name(s) to plot (x-axis).
    group_by : metadata column for grouping (default: active idents, y-axis). A
               categorical column is drawn in its category order, as R draws a factor's levels.
    scale    : z-score each feature's average expression across groups.
    col_min/col_max : colour scale limits for the scaled average expression.
    dot_min  : minimum fraction-expressing to draw a dot.
    dot_scale: scaling factor for dot area.
    """
    plt = _mpl()
    import matplotlib as mpl

    if isinstance(features, str):
        features = [features]

    groups = _get_groups(obj, group_by)
    unique = _group_levels(obj, group_by, groups)

    # Per (feature, group): average expression and fraction expressing.
    avg = np.zeros((len(features), len(unique)))
    pct = np.zeros((len(features), len(unique)))
    for fi, feat in enumerate(features):
        expr = _get_expression(obj, feat, assay, layer)
        for gi, g in enumerate(unique):
            cell_vals = expr[groups == g]
            if cell_vals.size == 0:
                continue
            avg[fi, gi] = cell_vals.mean()
            pct[fi, gi] = float((cell_vals > 0).mean())

    # Scale average expression across groups (z-score per feature), like Seurat.
    color_vals = avg.copy()
    if scale:
        mu = avg.mean(axis=1, keepdims=True)
        sd = avg.std(axis=1, ddof=1, keepdims=True)
        sd[sd == 0] = 1.0
        color_vals = np.clip((avg - mu) / sd, col_min, col_max)
        vmin, vmax = col_min, col_max
    else:
        # log1p of average expression, as Seurat does for the unscaled case.
        color_vals = np.log1p(avg)
        vmin, vmax = float(color_vals.min()), float(color_vals.max())

    cmap = mpl.colors.LinearSegmentedColormap.from_list("dotplot", list(cols))
    norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)

    if figsize is None:
        figsize = (max(5, 0.7 * len(features) + 2), max(3, 0.5 * len(unique) + 1.5))
    fig, ax = plt.subplots(figsize=figsize)

    for fi in range(len(features)):
        for gi in range(len(unique)):
            frac = pct[fi, gi]
            if frac < dot_min or frac == 0:
                continue
            ax.scatter(
                fi, gi, s=(frac * dot_scale) ** 2 * np.pi,
                c=[cmap(norm(color_vals[fi, gi]))],
                edgecolors="black", linewidths=0.3, zorder=3,
            )

    ax.set_xticks(range(len(features)))
    ax.set_xticklabels(features, rotation=90, fontsize=_fs("label"))
    ax.set_yticks(range(len(unique)))
    ax.set_yticklabels(unique, fontsize=_fs("label"))
    ax.set_xlim(-0.5, len(features) - 0.5)
    ax.set_ylim(-0.5, len(unique) - 0.5)
    ax.set_axisbelow(True)
    ax.grid(True, color="0.9", linewidth=0.6)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)

    # Colour bar (average expression) and a size legend (% expressed).
    sm = mpl.cm.ScalarMappable(norm=norm, cmap=cmap)
    sm.set_array([])
    cbar = fig.colorbar(sm, ax=ax, shrink=0.4, pad=0.02, aspect=12)
    cbar.set_label("Average expression" + (" (scaled)" if scale else ""), fontsize=_fs("small"))

    for frac in (0.25, 0.5, 0.75, 1.0):
        ax.scatter([], [], s=(frac * dot_scale) ** 2 * np.pi, c="grey",
                   edgecolors="black", linewidths=0.3, label=f"{int(frac * 100)}%")
    ax.legend(title="% expressed", bbox_to_anchor=(1.22, 0.0), loc="lower left",
              labelspacing=1.2, frameon=False, fontsize=_fs("small"), title_fontsize=_fs("small"))

    fig.tight_layout()
    return fig

feature_scatter

feature_scatter(obj, feature1: str, feature2: str, group_by: Optional[str] = None, assay: Optional[str] = None, layer: Optional[str] = None, pt_size: float = 4.0, alpha: float = 0.6, figsize: tuple = (6, 5), palette: Optional[list] = None, raster: Optional[bool] = None) -> 'Figure'

Scatter plot of two features coloured by identity.

Mirrors R's FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt").

Parameters:

• feature1 (str) –
• group_by (column for colouring; default: active idents. A categorical column is, default: None ) –

     coloured in its category order, as R colours a factor's levels.
  

• raster (Optional[bool], default: None ) –

     ``None`` (default) rasterises above 100,000 cells.
  

Source code in truecell/plotting.py

def feature_scatter(
    obj,
    feature1: str,
    feature2: str,
    group_by: Optional[str] = None,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    pt_size: float = 4.0,
    alpha: float = 0.6,
    figsize: tuple = (6, 5),
    palette: Optional[list] = None,
    raster: Optional[bool] = None,
) -> "Figure":
    """Scatter plot of two features coloured by identity.

    Mirrors R's ``FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")``.

    Parameters
    ----------
    feature1 / feature2 : gene names or metadata columns
    group_by : column for colouring; default: active idents. A categorical column is
               coloured in its category order, as R colours a factor's levels.
    raster   : draw the cells as a raster layer instead of vector paths.
               ``None`` (default) rasterises above 100,000 cells.
    """
    plt = _mpl()
    x = _get_expression(obj, feature1, assay, layer)
    y = _get_expression(obj, feature2, assay, layer)
    groups = _get_groups(obj, group_by)
    unique = _group_levels(obj, group_by, groups)
    colors = palette or _palette(len(unique))
    rast = _should_raster(raster, len(x))

    fig, ax = plt.subplots(figsize=figsize)
    for g, color in zip(unique, colors):
        mask = groups == g
        ax.scatter(x[mask], y[mask], s=pt_size, alpha=alpha,
                   color=color, label=g, linewidths=0, rasterized=rast)

    ax.set_xlabel(feature1)
    ax.set_ylabel(feature2)
    ax.set_title(f"{feature1} vs {feature2}", fontsize=_fs("large"))
    ax.legend(markerscale=2.5, fontsize=_fs("small"),
              bbox_to_anchor=(1.01, 1), loc="upper left", frameon=False)
    _strip_axes(ax)
    fig.tight_layout()
    return fig

elbow_plot

elbow_plot(obj, reduction: str = 'pca', ndims: int = 20, figsize: tuple = (7, 4)) -> 'Figure'

Rank principal components by standard deviation.

Mirrors R's ElbowPlot(pbmc).

Parameters:

• reduction (str, default: 'pca' ) –

  name of the PCA-like reduction

• ndims (int, default: 20 ) –

  number of PCs to show

Source code in truecell/plotting.py

def elbow_plot(
    obj,
    reduction: str = "pca",
    ndims: int = 20,
    figsize: tuple = (7, 4),
) -> "Figure":
    """Rank principal components by standard deviation.

    Mirrors R's ``ElbowPlot(pbmc)``.

    Parameters
    ----------
    reduction : name of the PCA-like reduction
    ndims     : number of PCs to show
    """
    plt = _mpl()
    if reduction not in obj.reductions:
        raise KeyError(f"Reduction '{reduction}' not found.")

    stdev = obj.reductions[reduction].stdev
    n = min(ndims, len(stdev))
    pcs = np.arange(1, n + 1)

    fig, ax = plt.subplots(figsize=figsize)
    ax.plot(pcs, stdev[:n], "o-", color="#F8766D", markersize=5, linewidth=1.8)
    ax.set_xlabel("Principal Component")
    ax.set_ylabel("Standard Deviation")
    ax.set_title("Elbow Plot", fontsize=_fs("large"))
    ax.set_xticks(pcs[::max(1, n // 10)])
    _strip_axes(ax)
    fig.tight_layout()
    return fig

variable_feature_plot

variable_feature_plot(obj, assay: Optional[str] = None, log: bool = True, label: bool = True, n_label: int = 10, figsize: tuple = (9, 5), raster: Optional[bool] = None, repel: bool = False) -> 'Figure'

Plot mean expression vs dispersion and highlight variable features.

Mirrors R's VariableFeaturePlot(pbmc).

Parameters:

• log (bool, default: True ) –

  use log10 axes

• label (bool, default: True ) –

  annotate the top n_label HVGs by name

• n_label (int, default: 10 ) –

  number of top HVGs to label

• raster (Optional[bool], default: None ) –

  draw the points as a raster layer instead of vector paths. None (default) rasterises above 100,000 points. The points here are genes, not cells, so the default rarely engages.

• repel (bool, default: False ) –

  move the names apart and off every named point, as Seurat's LabelPoints(repel = TRUE) does; a line joins a name that has moved away to its point. Without it each name sits just above and to the right of its point, and names in a dense region overlap.

Source code in truecell/plotting.py

def variable_feature_plot(
    obj,
    assay: Optional[str] = None,
    log: bool = True,
    label: bool = True,
    n_label: int = 10,
    figsize: tuple = (9, 5),
    raster: Optional[bool] = None,
    repel: bool = False,
) -> "Figure":
    """Plot mean expression vs dispersion and highlight variable features.

    Mirrors R's ``VariableFeaturePlot(pbmc)``.

    Parameters
    ----------
    log     : use log10 axes
    label   : annotate the top *n_label* HVGs by name
    n_label : number of top HVGs to label
    raster  : draw the points as a raster layer instead of vector paths.
              ``None`` (default) rasterises above 100,000 points. The points
              here are genes, not cells, so the default rarely engages.
    repel   : move the names apart and off every named point, as Seurat's
              ``LabelPoints(repel = TRUE)`` does; a line joins a name that has
              moved away to its point. Without it each name sits just above and
              to the right of its point, and names in a dense region overlap.
    """
    import scipy.sparse as sp
    plt = _mpl()

    assay_obj = _get_assay_obj(obj, assay)
    feat_names = assay_obj._all_feature_names
    hvg_set = set(assay_obj.variable_features)
    top_labeled = assay_obj.variable_features[:n_label]
    is_hvg = np.array([f in hvg_set for f in feat_names])

    # Prefer the statistics find_variable_features stored, whichever method
    # ran. `VariableFeaturePlot` picks its y axis off the column names the same
    # way: `variance.standardized` for vst, `mvp.dispersion` for the dispersion
    # methods -- and the axis label is the only thing on the figure that says
    # which one you are looking at.
    md = getattr(assay_obj, "meta_data", None)
    stored = set(md.columns) if md is not None else set()
    use_std_var = {"mean", "variance.standardized"} <= stored
    use_mvp = not use_std_var and {"mvp.mean", "mvp.dispersion"} <= stored

    if md is not None and use_std_var:
        means = md["mean"].values
        y_vals = md["variance.standardized"].values
        y_label = "Standardized Variance"
        # clip extreme standardized values for readability
        y_vals = np.clip(y_vals, 0, np.percentile(y_vals[y_vals > 0], 99.5))
        log = False  # standardized variance on linear scale matches R
    elif md is not None and use_mvp:
        # Seurat plots columns 1 and 2 of the mvp HVFInfo -- the mean and the
        # *raw* dispersion, not the scaled one that the cutoff tests -- and
        # leaves `log` off, both axes already being logarithms.
        means = md["mvp.mean"].values
        y_vals = md["mvp.dispersion"].values
        y_label = "Dispersion"
        log = False
    else:
        mat = _get_data_matrix(assay_obj)
        if sp.issparse(mat):
            means = np.array(mat.mean(axis=1)).flatten()
            sq_means = np.array(mat.power(2).mean(axis=1)).flatten()
        else:
            d = np.asarray(mat, dtype=float)
            means = d.mean(axis=1)
            sq_means = (d ** 2).mean(axis=1)
        y_vals = sq_means - means ** 2
        y_label = "Dispersion (log)" if log else "Dispersion"

    rast = _should_raster(raster, len(feat_names))
    fig, ax = plt.subplots(figsize=figsize)
    ax.scatter(means[~is_hvg], y_vals[~is_hvg], s=3, alpha=0.3,
               color="#AAAAAA", label="Other genes", linewidths=0,
               rasterized=rast)
    ax.scatter(means[is_hvg], y_vals[is_hvg], s=4, alpha=0.7,
               color="#F8766D", label="Variable features", linewidths=0,
               rasterized=rast)

    if label:
        named = []
        for gene in top_labeled:
            if gene in feat_names:
                i = feat_names.index(gene)
                named.append((gene, (means[i], y_vals[i])))
        if repel:
            # Each name starts on its point. Passing every named point as a
            # group of one keeps a name off its own point and off the others'.
            texts = [ax.text(x, y, gene, fontsize=_fs("tiny"), color="#333333",
                             ha="center", va="center")
                     for gene, (x, y) in named]
            from ._repel import attach
            attach(ax, texts, [xy for _, xy in named], [[xy] for _, xy in named])
        else:
            for gene, xy in named:
                ax.annotate(gene, xy, fontsize=_fs("tiny"), xytext=(3, 3),
                            textcoords="offset points", color="#333333")

    if log and not use_std_var:
        ax.set_xscale("log")
        ax.set_yscale("log")
        ax.set_xlabel("Average Expression (log)")
    else:
        ax.set_xlabel("Average Expression")

    ax.set_ylabel(y_label)
    ax.set_title(f"Highly Variable Features  ({len(hvg_set):,} selected)",
                 fontsize=_fs("large"))
    ax.legend(markerscale=3, fontsize=_fs("label"))
    _strip_axes(ax)
    fig.tight_layout()
    return fig

viz_dim_loadings

viz_dim_loadings(obj, reduction: str = 'pca', dims: Union[int, list[int]] = [1, 2], n_features: int = 15, ncol: Optional[int] = None, figsize: Optional[tuple] = None) -> 'Figure'

Horizontal bar charts of top positive and negative loading genes per PC.

Mirrors R's VizDimLoadings(pbmc, dims = 1:2, reduction = "pca").

Parameters:

• dims (Union[int, list[int]], default: [1, 2] ) –

  PC indices (1-based) to visualise

• n_features (int, default: 15 ) –

  number of top genes to show per direction (positive + negative)

Source code in truecell/plotting.py

def viz_dim_loadings(
    obj,
    reduction: str = "pca",
    dims: Union[int, list[int]] = [1, 2],
    n_features: int = 15,
    ncol: Optional[int] = None,
    figsize: Optional[tuple] = None,
) -> "Figure":
    """Horizontal bar charts of top positive and negative loading genes per PC.

    Mirrors R's ``VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")``.

    Parameters
    ----------
    dims       : PC indices (1-based) to visualise
    n_features : number of top genes to show per direction (positive + negative)
    """
    plt = _mpl()
    if isinstance(dims, int):
        dims = [dims]
    dims_0 = [d - 1 for d in dims]

    dr = obj.reductions[reduction]
    loadings = dr.feature_loadings
    feat_names = list(dr._feature_names) if hasattr(dr, "_feature_names") else []

    nrow, nc = _subplot_grid(len(dims), ncol)
    if figsize is None:
        figsize = (nc * 4.5, nrow * max(4, n_features * 0.35))

    fig, axes = plt.subplots(nrow, nc, figsize=figsize, squeeze=False,
                             layout="constrained")
    axes_flat = axes.flatten()

    for plot_i, (dim, d0) in enumerate(zip(dims, dims_0)):
        ax = axes_flat[plot_i]
        col = loadings[:, d0]

        top_pos_idx = np.argsort(col)[::-1][:n_features]
        top_neg_idx = np.argsort(col)[:n_features]
        # Combine: negatives at top (sorted most negative first), positives below
        idx = np.concatenate([top_pos_idx[::-1], top_neg_idx[::-1]])
        genes = [feat_names[i] if i < len(feat_names) else f"gene_{i}" for i in idx]
        vals  = col[idx]
        colors = ["#F8766D" if v > 0 else "#00BFC4" for v in vals]

        y_pos = np.arange(len(genes))
        ax.barh(y_pos, vals, color=colors, edgecolor="none", height=0.75)
        ax.set_yticks(y_pos)
        # Shrinks as rows are added so the labels stay legible, but the shrink
        # and its floor both scale with the theme rather than pinning to 9 and 6.
        ax.set_yticklabels(
            genes,
            fontsize=max(_fs("label") * 0.67, _fs("label") - n_features // 10),
        )
        ax.axvline(0, color="#333333", linewidth=0.8)
        ax.set_xlabel("Loading score")
        ax.set_title(f"{reduction.upper()} {dim}", fontsize=_fs("large"), fontweight="bold")
        _strip_axes(ax)

    for i in range(len(dims), len(axes_flat)):
        axes_flat[i].set_visible(False)

    # One legend under the panels. Inside each panel, at lower right, it sat on
    # the positive loadings' bars, and no corner stays clear at every size.
    from matplotlib.patches import Patch
    fig.legend(handles=[Patch(color="#F8766D", label="Positive"),
                        Patch(color="#00BFC4", label="Negative")],
               fontsize=_fs("small"), frameon=False, loc="outside lower center", ncol=2)
    return fig

dim_heatmap

dim_heatmap(obj, reduction: str = 'pca', dims: Union[int, list[int]] = 1, cells: int | None = None, balanced: bool = True, ncol: Optional[int] = None, figsize: Optional[tuple] = None) -> 'Figure'

Heatmap of gene loadings for selected principal components.

Shows cells in order of their score, highest on the left, against the top loading genes for each PC — mirrors R's DimHeatmap(pbmc, dims = 1:6).

Parameters:

• dims (Union[int, list[int]], default: 1 ) –

  PC index (1-based int) or list of indices

• cells (int | None, default: None ) –

  how many cells to show per PC, picked the way DimHeatmap picks them. None (the default) shows every cell, as Seurat does; the PBMC 3k vignette passes 500.

• balanced (bool, default: True ) –

  if True, take equal numbers from both extremes of the PC score

Source code in truecell/plotting.py

def dim_heatmap(
    obj,
    reduction: str = "pca",
    dims: Union[int, list[int]] = 1,
    cells: int | None = None,
    balanced: bool = True,
    ncol: Optional[int] = None,
    figsize: Optional[tuple] = None,
) -> "Figure":
    """Heatmap of gene loadings for selected principal components.

    Shows cells in order of their score, highest on the left, against the top
    loading genes for each PC — mirrors R's ``DimHeatmap(pbmc, dims = 1:6)``.

    Parameters
    ----------
    dims     : PC index (1-based int) or list of indices
    cells    : how many cells to show per PC, picked the way DimHeatmap picks
               them. ``None`` (the default) shows every cell, as Seurat does; the
               PBMC 3k vignette passes 500.
    balanced : if True, take equal numbers from both extremes of the PC score
    """
    import scipy.sparse as sp
    plt = _mpl()

    if isinstance(dims, int):
        dims = [dims]
    dims_0 = [d - 1 for d in dims]  # convert to 0-based

    dr = obj.reductions[reduction]
    emb = dr.cell_embeddings
    loadings = dr.feature_loadings
    feat_names = dr._feature_names if hasattr(dr, "_feature_names") else []

    n_top = 10
    nrow, nc = _subplot_grid(len(dims), ncol)
    if figsize is None:
        figsize = (nc * 5, nrow * 5)

    fig, axes = plt.subplots(nrow, nc, figsize=figsize, squeeze=False)
    axes_flat = axes.flatten()

    assay_obj = _get_assay_obj(obj, None)
    mat, scaled_feats = _resolve_layer(assay_obj, "scale.data")
    scaled_row = {g: i for i, g in enumerate(scaled_feats)}

    for plot_i, (dim, d0) in enumerate(zip(dims, dims_0)):
        ax = axes_flat[plot_i]
        scores = emb[:, d0]
        col_loads = loadings[:, d0]

        # DimHeatmap's cells, highest score on the left. `n_half` is where the
        # high end meets the low end, for the separator below.
        n_cells = len(scores) if cells is None else cells
        sel_cells = _dim_heatmap_cells(scores, n_cells, balanced)
        n_half = round(min(n_cells, len(scores)) / 2)

        # Select top-loading genes
        top_genes_idx = np.concatenate([
            np.argsort(col_loads)[::-1][:n_top],
            np.argsort(col_loads)[:n_top]
        ])
        top_genes_idx = np.unique(top_genes_idx)
        top_genes_idx = top_genes_idx[np.argsort(col_loads[top_genes_idx])[::-1]]

        # Build expression matrix: top genes × selected cells. A gene's row is its
        # place in scale.data's own feature list, which is narrower than the assay's.
        rows = []
        gene_labels = []
        for gi in top_genes_idx:
            if gi < len(feat_names):
                gname = feat_names[gi]
            else:
                continue
            if gname in scaled_row:
                row = mat[scaled_row[gname], :]
                row = np.asarray(row.todense()).flatten() if sp.issparse(row) else np.asarray(row).flatten()
                rows.append(row[sel_cells])
                gene_labels.append(gname)

        if not rows:
            ax.set_visible(False)
            continue

        mat_sub = np.array(rows)
        mat_sub = np.clip(mat_sub, -2.5, 2.5)

        im = ax.imshow(mat_sub, aspect="auto", cmap="RdBu_r",
                       vmin=-2.5, vmax=2.5, interpolation="none")
        ax.set_yticks(range(len(gene_labels)))
        ax.set_yticklabels(gene_labels, fontsize=_fs("tiny"))
        ax.set_xticks([])
        if balanced:
            ax.axvline(n_half - 0.5, color="white", linewidth=1.5)
        ax.set_title(f"PC {dim}", fontsize=_fs("title"))
        plt.colorbar(im, ax=ax, shrink=0.6, label="Scaled expr.")

    for i in range(len(dims), len(axes_flat)):
        axes_flat[i].set_visible(False)

    fig.suptitle(f"DimHeatmap — {reduction.upper()}", fontsize=_fs("suptitle"),
                 fontweight="bold", y=1.01)
    fig.tight_layout()
    return fig

do_heatmap

do_heatmap(obj, features: list[str], group_by: Optional[str] = None, assay: Optional[str] = None, layer: str = 'scale.data', label: bool = True, cells: Optional[list[str]] = None, figsize: Optional[tuple] = None, palette: Optional[list] = None) -> 'Figure'

Expression heatmap of selected genes across cells, sorted by cluster.

Mirrors R's DoHeatmap(pbmc, features = top10$gene).

Parameters:

• features (list[str]) –

  list of gene names to show as rows

• group_by (Optional[str], default: None ) –

  column used to sort and colour cells (default: active idents). A categorical column is ordered by its categories, as R orders a factor's levels.

• layer (str, default: 'scale.data' ) –

  which data layer to use (default: "scale.data")

• label (bool, default: True ) –

  annotate cluster boundaries with group names

• cells (Optional[list[str]], default: None ) –

  restrict to these cells and show them in this exact order, rather than the default group-sorted order. Mirrors R's cells argument; mixscape_heatmap uses it to order cells by knockout probability.

Source code in truecell/plotting.py

def do_heatmap(
    obj,
    features: list[str],
    group_by: Optional[str] = None,
    assay: Optional[str] = None,
    layer: str = "scale.data",
    label: bool = True,
    cells: Optional[list[str]] = None,
    figsize: Optional[tuple] = None,
    palette: Optional[list] = None,
) -> "Figure":
    """Expression heatmap of selected genes across cells, sorted by cluster.

    Mirrors R's ``DoHeatmap(pbmc, features = top10$gene)``.

    Parameters
    ----------
    features : list of gene names to show as rows
    group_by : column used to sort and colour cells (default: active idents). A
               categorical column is ordered by its categories, as R orders a factor's levels.
    layer    : which data layer to use (default: "scale.data")
    label    : annotate cluster boundaries with group names
    cells    : restrict to these cells *and* show them in this exact order,
               rather than the default group-sorted order. Mirrors R's ``cells``
               argument; :func:`mixscape_heatmap` uses it to order cells by
               knockout probability.
    """
    import scipy.sparse as sp
    plt = _mpl()

    assay_obj = _get_assay_obj(obj, assay)
    mat, all_feats = _resolve_layer(assay_obj, layer)

    groups = _get_groups(obj, group_by)
    unique = _group_levels(obj, group_by, groups)
    colors = palette or _palette(len(unique))

    if cells is None:
        cell_order = np.argsort([unique.index(g) for g in groups])
    else:
        pos = {c: i for i, c in enumerate(obj.cell_names())}
        missing = [c for c in cells if c not in pos]
        if missing:
            raise KeyError(f"Cells not in object: {missing[:5]}")
        cell_order = np.array([pos[c] for c in cells], dtype=int)
    sorted_groups = groups[cell_order]

    def _extract(gene):
        if gene not in all_feats:
            return np.zeros(len(groups))
        idx = all_feats.index(gene)
        row = mat[idx, :]
        arr = np.asarray(row.todense()).flatten() if sp.issparse(row) else np.asarray(row).flatten()
        return arr

    rows = [_extract(g)[cell_order] for g in features]
    matrix = np.clip(np.array(rows), -2.5, 2.5)

    if figsize is None:
        figsize = (max(10, len(groups) / 200), min(14, max(4, len(features) * 0.28)))

    fig, axes = plt.subplots(
        2, 1, figsize=figsize,
        gridspec_kw={"height_ratios": [0.05, 1], "hspace": 0.02}
    )

    # Top colour bar showing cluster membership
    colour_row = np.array([unique.index(g) for g in sorted_groups])[np.newaxis, :]
    axes[0].imshow(colour_row, aspect="auto",
                   cmap=plt.matplotlib.colors.ListedColormap(colors),
                   vmin=0, vmax=len(unique) - 1, interpolation="none")
    axes[0].set_xticks([])
    axes[0].set_yticks([])
    for spine in axes[0].spines.values():
        spine.set_visible(False)

    # Main heatmap
    im = axes[1].imshow(matrix, aspect="auto", cmap="RdBu_r",
                        vmin=-2.5, vmax=2.5, interpolation="none")
    axes[1].set_yticks(range(len(features)))
    axes[1].set_yticklabels(
        features,
        fontsize=max(_fs("label") * 0.56, _fs("label") - len(features) // 20),
    )
    axes[1].set_xticks([])

    # Cluster boundary lines + labels inside the colour bar row.
    # Walk contiguous runs rather than whole groups: with the default sort each
    # group is one run, but an explicit `cells` order (e.g. by knockout
    # probability) interleaves them, and a run-based pass labels both correctly.
    if label:
        boundaries = [0] + [
            i for i in range(1, len(sorted_groups))
            if sorted_groups[i] != sorted_groups[i - 1]
        ] + [len(sorted_groups)]
        # Interleaved orders produce many short runs; label only the ones wide
        # enough to read, so the colour bar does not fill with overlapping text.
        min_run = 0 if cells is None else max(1, len(sorted_groups) // 50)
        for start, end in zip(boundaries[:-1], boundaries[1:]):
            # Separator lines in both the colour bar and the heatmap
            axes[0].axvline(start - 0.5, color="white", linewidth=1.2)
            axes[1].axvline(start - 0.5, color="white", linewidth=0.8)
            if end - start >= min_run:
                # Labels sit inside axes[0] (the colour bar) — y=0.5 centres vertically
                axes[0].text((start + end) / 2, 0.5, sorted_groups[start],
                             ha="center", va="center",
                             fontsize=_fs("tiny"), color="white", fontweight="bold")

    plt.colorbar(im, ax=axes[1], shrink=0.4, pad=0.01, label="Scaled expression")
    axes[1].set_xlabel("Cells (sorted by cluster)")
    fig.suptitle("Expression Heatmap", fontsize=_fs("suptitle"), fontweight="bold", y=1.0)
    return fig

Spatial

image_dim_plot

image_dim_plot(obj, group_by: Optional[str] = None, image: Optional[Union[str, list]] = None, size: float = 1.0, cols: Optional[dict] = None, ncol: Optional[int] = None, flip_y: bool = True, figsize: Optional[tuple] = None, raster: Optional[bool] = None) -> 'Figure'

Plot cell centroids in physical space, coloured by a grouping variable.

Matplotlib equivalent of Seurat's ImageDimPlot — drawn directly from centroids (obj.get_tissue_coordinates), so it is immune to the ggplot2-4.x blank-render issue. One panel per image.

Parameters:

• group_by (Optional[str], default: None ) –

  metadata column (default: active idents) for point colour.

• image (Optional[Union[str, list]], default: None ) –

  image name(s) to draw (default: all).

• cols (Optional[dict], default: None ) –

  optional {group: colour} mapping.

• flip_y (bool, default: True ) –

  invert the y-axis so images match Seurat's orientation.

• raster (Optional[bool], default: None ) –

  draw the points as a raster layer instead of vector paths. None (default) rasterises above 100,000 points.

Source code in truecell/plotting.py

def image_dim_plot(
    obj,
    group_by: Optional[str] = None,
    image: Optional[Union[str, list]] = None,
    size: float = 1.0,
    cols: Optional[dict] = None,
    ncol: Optional[int] = None,
    flip_y: bool = True,
    figsize: Optional[tuple] = None,
    raster: Optional[bool] = None,
) -> "Figure":
    """Plot cell centroids in physical space, coloured by a grouping variable.

    Matplotlib equivalent of Seurat's ``ImageDimPlot`` — drawn directly from
    centroids (``obj.get_tissue_coordinates``), so it is immune to the
    ggplot2-4.x blank-render issue. One panel per image.

    Parameters
    ----------
    group_by : metadata column (default: active idents) for point colour.
    image    : image name(s) to draw (default: all).
    cols     : optional ``{group: colour}`` mapping.
    flip_y   : invert the y-axis so images match Seurat's orientation.
    raster   : draw the points as a raster layer instead of vector paths.
               ``None`` (default) rasterises above 100,000 points.
    """
    plt = _mpl()
    coords = obj.get_tissue_coordinates(image)
    if coords.empty:
        raise ValueError("Object has no spatial coordinates to plot.")

    labels = _get_groups(obj, group_by)
    lab_by_cell = dict(zip(obj.cell_names(), labels))
    coords = coords.assign(group=coords["cell"].map(lab_by_cell).astype(str))

    images = list(dict.fromkeys(coords["image"]))
    uniq = sorted(coords["group"].unique())
    if cols is None:
        cols = dict(zip(uniq, _palette(len(uniq))))

    nrow, ncol = _subplot_grid(len(images), ncol)
    if figsize is None:
        figsize = (5 * ncol, 4.5 * nrow)
    # Constrained, so the legend beside the panels gets the width it needs.
    fig, axes = plt.subplots(nrow, ncol, figsize=figsize, squeeze=False,
                             layout="constrained")
    axes_flat = axes.ravel()

    for ax, img in zip(axes_flat, images):
        d = coords[coords["image"] == img]
        for g, dd in d.groupby("group"):
            ax.scatter(dd["x"], dd["y"], s=size, c=[cols.get(g, "grey")],
                       rasterized=_should_raster(raster, len(coords)),
                       label=g, linewidths=0)
        ax.set_title(str(img), fontsize=_fs("label"), fontweight="bold")
        ax.set_aspect("equal")
        ax.axis("off")
        if flip_y:
            ax.invert_yaxis()
    for ax in axes_flat[len(images):]:
        ax.axis("off")

    handles = [plt.Line2D([], [], marker="o", linestyle="", markersize=6,
                          markerfacecolor=cols.get(g, "grey"), markeredgewidth=0)
               for g in uniq]
    # Outside the panels, where ggplot2 puts Seurat's legend. The fixed 12% strip
    # this had was narrower than a legend of cell-type names, which then lay
    # over the tissue: 4,151 cells in the Xenium tutorial's figure.
    fig.legend(handles, uniq, title=group_by or "ident", loc="outside right center",
               fontsize=_fs("small"), title_fontsize=_fs("small"), frameon=False)
    return fig

image_feature_plot

image_feature_plot(obj, feature: str, image: Optional[Union[str, list]] = None, size: float = 1.0, cmap: str = 'viridis', assay: Optional[str] = None, layer: Optional[str] = None, ncol: Optional[int] = None, flip_y: bool = True, figsize: Optional[tuple] = None, raster: Optional[bool] = None) -> 'Figure'

Plot cell centroids in physical space, coloured by a feature's expression.

Matplotlib equivalent of Seurat's ImageFeaturePlot. One panel per image. raster : draw the points as a raster layer instead of vector paths. None (default) rasterises above 100,000 points.

Source code in truecell/plotting.py

def image_feature_plot(
    obj,
    feature: str,
    image: Optional[Union[str, list]] = None,
    size: float = 1.0,
    cmap: str = "viridis",
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    ncol: Optional[int] = None,
    flip_y: bool = True,
    figsize: Optional[tuple] = None,
    raster: Optional[bool] = None,
) -> "Figure":
    """Plot cell centroids in physical space, coloured by a feature's expression.

    Matplotlib equivalent of Seurat's ``ImageFeaturePlot``. One panel per image.
    raster   : draw the points as a raster layer instead of vector paths.
               ``None`` (default) rasterises above 100,000 points.
    """
    plt = _mpl()
    coords = obj.get_tissue_coordinates(image)
    if coords.empty:
        raise ValueError("Object has no spatial coordinates to plot.")

    expr = _get_expression(obj, feature, assay, layer)
    expr_by_cell = dict(zip(obj.cell_names(), expr))
    coords = coords.assign(val=coords["cell"].map(expr_by_cell).astype(float))

    images = list(dict.fromkeys(coords["image"]))
    nrow, ncol = _subplot_grid(len(images), ncol)
    if figsize is None:
        figsize = (5 * ncol, 4.5 * nrow)
    fig, axes = plt.subplots(nrow, ncol, figsize=figsize, squeeze=False)
    axes_flat = axes.ravel()
    vmax = float(np.nanmax(coords["val"])) if len(coords) else 1.0

    sc = None
    for ax, img in zip(axes_flat, images):
        d = coords[coords["image"] == img]
        sc = ax.scatter(d["x"], d["y"], s=size, c=d["val"], cmap=cmap,
                        vmin=0, vmax=vmax, linewidths=0,
                        rasterized=_should_raster(raster, len(coords)))
        ax.set_title(str(img), fontsize=_fs("label"), fontweight="bold")
        ax.set_aspect("equal")
        ax.axis("off")
        if flip_y:
            ax.invert_yaxis()
    for ax in axes_flat[len(images):]:
        ax.axis("off")
    if sc is not None:
        fig.colorbar(sc, ax=axes_flat.tolist(), shrink=0.5, label=feature)
    return fig

spatial_dim_plot

spatial_dim_plot(obj, group_by: Optional[str] = None, image: Optional[Union[str, list]] = None, cols: Optional[dict] = None, pt_size_factor: float = 1.6, size: float = 12.0, alpha: float = 1.0, image_alpha: float = 1.0, resolution: Optional[str] = None, crop: bool = True, ncol: Optional[int] = None, figsize: Optional[tuple] = None, raster: Optional[bool] = None) -> 'Figure'

Plot spots on the tissue image, coloured by a grouping variable.

Matplotlib equivalent of Seurat's SpatialDimPlot. Each panel draws the H&E photo held by a VisiumV2 image and overlays the spots at their true diameter. An image slot with no photo (a plain FOV, or a Visium bundle loaded with image=False) degrades to a bare scatter of the same spots — the plot still works, it just has no tissue underneath.

Parameters:

• group_by (Optional[str], default: None ) –

  metadata column (default: active idents) for spot colour.

• image (Optional[Union[str, list]], default: None ) –

  image name(s) to draw (default: all).

• cols (Optional[dict], default: None ) –

  optional {group: colour} mapping.

• pt_size_factor (float, default: 1.6 ) –

  scales the spots relative to their real diameter, as in Seurat. Ignored when the image has no scalefactors_json.json, in which case size (a plain scatter point size) applies instead.

• image_alpha (float, default: 1.0 ) –

  opacity of the tissue photo — drop it to make spots pop.

• resolution (Optional[str], default: None ) –

  "hires" / "lowres"; defaults to whatever was loaded.

• crop (bool, default: True ) –

  zoom to the spots rather than showing the whole slide.

• raster (Optional[bool], default: None ) –

  draw the points as a raster layer instead of vector paths. None (default) rasterises above 100,000 points.

Source code in truecell/plotting.py

def spatial_dim_plot(
    obj,
    group_by: Optional[str] = None,
    image: Optional[Union[str, list]] = None,
    cols: Optional[dict] = None,
    pt_size_factor: float = 1.6,
    size: float = 12.0,
    alpha: float = 1.0,
    image_alpha: float = 1.0,
    resolution: Optional[str] = None,
    crop: bool = True,
    ncol: Optional[int] = None,
    figsize: Optional[tuple] = None,
    raster: Optional[bool] = None,
) -> "Figure":
    """Plot spots on the tissue image, coloured by a grouping variable.

    Matplotlib equivalent of Seurat's ``SpatialDimPlot``. Each panel draws the
    H&E photo held by a :class:`~truecell.spatial.visium.VisiumV2` image and
    overlays the spots at their true diameter. An image slot with no photo (a
    plain ``FOV``, or a Visium bundle loaded with ``image=False``) degrades to a
    bare scatter of the same spots — the plot still works, it just has no
    tissue underneath.

    Parameters
    ----------
    group_by : metadata column (default: active idents) for spot colour.
    image    : image name(s) to draw (default: all).
    cols     : optional ``{group: colour}`` mapping.
    pt_size_factor : scales the spots relative to their real diameter, as in
        Seurat. Ignored when the image has no ``scalefactors_json.json``, in
        which case ``size`` (a plain scatter point size) applies instead.
    image_alpha : opacity of the tissue photo — drop it to make spots pop.
    resolution : ``"hires"`` / ``"lowres"``; defaults to whatever was loaded.
    crop     : zoom to the spots rather than showing the whole slide.
    raster   : draw the points as a raster layer instead of vector paths.
               ``None`` (default) rasterises above 100,000 points.
    """
    plt = _mpl()
    fovs = _resolve_fovs(obj, image)

    labels = _get_groups(obj, group_by)
    lab_by_cell = dict(zip(obj.cell_names(), [str(v) for v in labels]))

    panels = {}
    for name, fov in fovs.items():
        coords, radius, img = _spatial_panel(fov, resolution)
        coords = coords.assign(group=[lab_by_cell.get(c) for c in coords.index])
        coords = coords[coords["group"].notna()]
        if not coords.empty:
            panels[name] = (coords, radius, img)
    if not panels:
        raise ValueError("No spots left to plot: the images share no cells with the object.")

    uniq = sorted({g for coords, _, _ in panels.values() for g in coords["group"]})
    if cols is None:
        cols = dict(zip(uniq, _palette(len(uniq))))

    nrow, ncol = _subplot_grid(len(panels), ncol)
    if figsize is None:
        figsize = (5 * ncol, 4.5 * nrow)
    # Constrained, so the legend beside the panels gets the width it needs.
    fig, axes = plt.subplots(nrow, ncol, figsize=figsize, squeeze=False,
                             layout="constrained")
    axes_flat = axes.ravel()

    for ax, (name, (coords, radius, img)) in zip(axes_flat, panels.items()):
        if img is not None:
            ax.imshow(img, alpha=image_alpha)
        face = [cols.get(g, "grey") for g in coords["group"]]
        coll = _spot_collection(ax, coords, radius, pt_size_factor)
        if coll is None:
            ax.scatter(coords["x"], coords["y"], s=size, c=face,
                       rasterized=_should_raster(raster, len(coords)),
                       alpha=alpha, linewidths=0)
        else:
            coll.set_facecolor(face)
            coll.set_alpha(alpha)
            ax.add_collection(coll)
        _spatial_limits(ax, coords, radius, img, crop)
        ax.set_title(str(name), fontsize=_fs("label"), fontweight="bold")
        ax.set_aspect("equal")
        ax.axis("off")
    for ax in axes_flat[len(panels):]:
        ax.axis("off")

    handles = [plt.Line2D([], [], marker="o", linestyle="", markersize=6,
                          markerfacecolor=cols.get(g, "grey"), markeredgewidth=0)
               for g in uniq]
    # Outside the panels, where ggplot2 puts Seurat's legend. The fixed 12% strip
    # this had was narrower than a legend of long group names, which then lay
    # over the spots.
    fig.legend(handles, uniq, title=group_by or "ident", loc="outside right center",
               fontsize=_fs("small"), title_fontsize=_fs("small"), frameon=False)
    return fig

spatial_feature_plot

spatial_feature_plot(obj, feature: str, image: Optional[Union[str, list]] = None, cmap: str = 'viridis', pt_size_factor: float = 1.6, size: float = 12.0, alpha: float = 1.0, image_alpha: float = 1.0, assay: Optional[str] = None, layer: Optional[str] = None, resolution: Optional[str] = None, crop: bool = True, ncol: Optional[int] = None, figsize: Optional[tuple] = None, raster: Optional[bool] = None) -> 'Figure'

Plot spots on the tissue image, coloured by a feature's expression.

Matplotlib equivalent of Seurat's SpatialFeaturePlot. Behaves exactly like spatial_dim_plot — same tissue background, same true-to-scale spots, same fallback when no image is stored — but colours the spots by a continuous value on a shared scale across panels. raster : draw the points as a raster layer instead of vector paths. None (default) rasterises above 100,000 points.

Source code in truecell/plotting.py

def spatial_feature_plot(
    obj,
    feature: str,
    image: Optional[Union[str, list]] = None,
    cmap: str = "viridis",
    pt_size_factor: float = 1.6,
    size: float = 12.0,
    alpha: float = 1.0,
    image_alpha: float = 1.0,
    assay: Optional[str] = None,
    layer: Optional[str] = None,
    resolution: Optional[str] = None,
    crop: bool = True,
    ncol: Optional[int] = None,
    figsize: Optional[tuple] = None,
    raster: Optional[bool] = None,
) -> "Figure":
    """Plot spots on the tissue image, coloured by a feature's expression.

    Matplotlib equivalent of Seurat's ``SpatialFeaturePlot``. Behaves exactly
    like :func:`spatial_dim_plot` — same tissue background, same true-to-scale
    spots, same fallback when no image is stored — but colours the spots by a
    continuous value on a shared scale across panels.
    raster   : draw the points as a raster layer instead of vector paths.
               ``None`` (default) rasterises above 100,000 points.
    """
    plt = _mpl()
    fovs = _resolve_fovs(obj, image)

    expr = _get_expression(obj, feature, assay, layer)
    expr_by_cell = dict(zip(obj.cell_names(), expr))

    panels = {}
    for name, fov in fovs.items():
        coords, radius, img = _spatial_panel(fov, resolution)
        coords = coords.assign(val=[expr_by_cell.get(c, np.nan) for c in coords.index])
        coords = coords[coords["val"].notna()]
        if not coords.empty:
            panels[name] = (coords, radius, img)
    if not panels:
        raise ValueError("No spots left to plot: the images share no cells with the object.")

    vals = np.concatenate([c["val"].to_numpy(dtype=float) for c, _, _ in panels.values()])
    vmax = float(np.nanmax(vals)) if vals.size else 1.0
    if not np.isfinite(vmax) or vmax <= 0:
        vmax = 1.0

    nrow, ncol = _subplot_grid(len(panels), ncol)
    if figsize is None:
        figsize = (5 * ncol, 4.5 * nrow)
    fig, axes = plt.subplots(nrow, ncol, figsize=figsize, squeeze=False)
    axes_flat = axes.ravel()

    mappable = None
    for ax, (name, (coords, radius, img)) in zip(axes_flat, panels.items()):
        if img is not None:
            ax.imshow(img, alpha=image_alpha)
        v = coords["val"].to_numpy(dtype=float)
        coll = _spot_collection(ax, coords, radius, pt_size_factor)
        if coll is None:
            mappable = ax.scatter(coords["x"], coords["y"], s=size, c=v, cmap=cmap,
                                  rasterized=_should_raster(raster, len(coords)),
                                  vmin=0, vmax=vmax, alpha=alpha, linewidths=0)
        else:
            coll.set_array(v)
            coll.set_cmap(cmap)
            coll.set_clim(0, vmax)
            coll.set_alpha(alpha)
            ax.add_collection(coll)
            mappable = coll
        _spatial_limits(ax, coords, radius, img, crop)
        ax.set_title(str(name), fontsize=_fs("label"), fontweight="bold")
        ax.set_aspect("equal")
        ax.axis("off")
    for ax in axes_flat[len(panels):]:
        ax.axis("off")

    if mappable is not None:
        fig.colorbar(mappable, ax=axes_flat.tolist(), shrink=0.5, label=feature)
    return fig

Mixscape diagnostics

plot_perturb_score

plot_perturb_score(obj, target_gene_ident: str, target_gene_class: str = 'gene', mixscape_class: str = 'mixscape_class', col: str = 'orange2', split_by: Optional[str] = None, before_mixscape: bool = False, prtb_type: str = 'KO', assay: str = 'PRTB', seed: int = 0, figsize: Optional[tuple] = None) -> 'Figure'

Density of one target gene's perturbation scores (Seurat's PlotPerturbScore).

Mirrors PlotPerturbScore(object, target.gene.ident = "IFNGR2", target.gene.class = "gene", mixscape.class = "mixscape_class"). This is the diagnostic that shows why mixscape split a guide the way it did: the score is each cell's projection onto the gene's perturbation axis, and the plot overlays the NT control density against the guide's own. A guide with a real effect is bimodal — one lobe sitting on the NT curve (the escapers) and one shifted away from it (the knockouts) — which is exactly the structure the mixture model is asked to find.

Two views, per R's before.mixscape:

• before_mixscape=False (default) — colour by mixscape_class, i.e. after the call: NT, "<gene> NP", and "<gene> KO" get their own curves, so you see where mixscape actually drew the line.
• before_mixscape=True — colour by the raw guide label only, the view you would have had without mixscape: NT against the whole guide population.

Cells are also drawn as a jittered strip — controls above the axis, the target gene below — so single-cell density is visible where the curves overlap.

Parameters:

• target_gene_ident (str) –

  the target gene to plot (must be one mixscape tested).

• target_gene_class (str, default: 'gene' ) –

  metadata column of per-cell guide class (default "gene").

• mixscape_class (str, default: 'mixscape_class' ) –

  metadata column of mixscape classifications.

• col (str, default: 'orange2' ) –

  colour for the target gene / knockout class. Controls and non-perturbed cells are fixed greys, as in R.

• split_by (Optional[str], default: None ) –

  metadata column to facet on, for screens spanning more than one cell type.

• before_mixscape (bool, default: False ) –

  colour by raw guide label instead of the mixscape call.

• prtb_type (str, default: 'KO' ) –

  perturbation label used by run_mixscape ("KO").

• assay (str, default: 'PRTB' ) –

  perturbation-signature assay the scores were stored under.

• seed (int, default: 0 ) –

  random state for the jitter strip (determinism).

Returns:

• matplotlib Figure –

Source code in truecell/plotting.py

def plot_perturb_score(
    obj,
    target_gene_ident: str,
    target_gene_class: str = "gene",
    mixscape_class: str = "mixscape_class",
    col: str = "orange2",
    split_by: Optional[str] = None,
    before_mixscape: bool = False,
    prtb_type: str = "KO",
    assay: str = "PRTB",
    seed: int = 0,
    figsize: Optional[tuple] = None,
) -> "Figure":
    """Density of one target gene's perturbation scores (Seurat's ``PlotPerturbScore``).

    Mirrors ``PlotPerturbScore(object, target.gene.ident = "IFNGR2",
    target.gene.class = "gene", mixscape.class = "mixscape_class")``. This is the
    diagnostic that shows *why* mixscape split a guide the way it did: the score
    is each cell's projection onto the gene's perturbation axis, and the plot
    overlays the NT control density against the guide's own. A guide with a real
    effect is bimodal — one lobe sitting on the NT curve (the escapers) and one
    shifted away from it (the knockouts) — which is exactly the structure the
    mixture model is asked to find.

    Two views, per R's ``before.mixscape``:

    * ``before_mixscape=False`` (default) — colour by ``mixscape_class``, i.e.
      *after* the call: NT, ``"<gene> NP"``, and ``"<gene> KO"`` get their own
      curves, so you see where mixscape actually drew the line.
    * ``before_mixscape=True`` — colour by the raw guide label only, the view you
      would have had without mixscape: NT against the whole guide population.

    Cells are also drawn as a jittered strip — controls above the axis, the target
    gene below — so single-cell density is visible where the curves overlap.

    Parameters
    ----------
    target_gene_ident : the target gene to plot (must be one mixscape tested).
    target_gene_class : metadata column of per-cell guide class (default ``"gene"``).
    mixscape_class    : metadata column of mixscape classifications.
    col               : colour for the target gene / knockout class. Controls and
                        non-perturbed cells are fixed greys, as in R.
    split_by          : metadata column to facet on, for screens spanning more
                        than one cell type.
    before_mixscape   : colour by raw guide label instead of the mixscape call.
    prtb_type         : perturbation label used by :func:`run_mixscape` (``"KO"``).
    assay             : perturbation-signature assay the scores were stored under.
    seed              : random state for the jitter strip (determinism).

    Returns
    -------
    matplotlib Figure
    """
    from scipy.stats import gaussian_kde
    plt = _mpl()

    scores = _mixscape_scores(obj, target_gene_ident, assay)
    pvec = scores["pvec"].to_numpy(dtype=float)
    cells = list(scores.index)
    # R indexes this column *by name* — `prtb_score[, target.gene.class]`, whose
    # default "gene" is `RunMixscape`'s default `labels`. Fall back to position
    # when the name is absent, because `run_mixscape(labels=...)` names the
    # column after whatever was passed and R would simply fail there; the frame
    # is always `["pvec", <labels>]`, so column 1 is the guide either way.
    if target_gene_class in scores.columns:
        guide_col = target_gene_class
    elif len(scores.columns) > 1:
        guide_col = scores.columns[1]
    else:
        raise KeyError(
            f"No guide-label column in the stored perturbation scores: wanted "
            f"{target_gene_class!r}, found {list(scores.columns)}. Pass "
            f"target_gene_class= the `labels` column run_mixscape was given."
        )
    guide = scores[guide_col].astype(str).to_numpy()

    nt_labels = sorted(set(guide) - {target_gene_ident})
    nt_name = nt_labels[0] if nt_labels else "NT"

    col = _r_colour(col)
    if before_mixscape:
        classes = [nt_name, target_gene_ident]
        colours = {nt_name: _r_colour("grey49"), target_gene_ident: col}
        group = guide
    else:
        if mixscape_class not in obj.meta_data.columns:
            raise KeyError(
                f"Column {mixscape_class!r} not in meta_data — run run_mixscape "
                "first, or pass before_mixscape=True."
            )
        group = obj.meta_data[mixscape_class].reindex(cells).astype(str).to_numpy()
        ko_name = f"{target_gene_ident} {prtb_type}"
        np_name = f"{target_gene_ident} NP"
        classes = [nt_name, np_name, ko_name]
        colours = {
            nt_name: _r_colour("grey49"),
            np_name: _r_colour("grey79"),
            ko_name: col,
        }

    if split_by is None:
        facets = [(None, np.arange(len(cells)))]
    else:
        if split_by not in obj.meta_data.columns:
            raise KeyError(f"split_by column {split_by!r} not found in meta_data.")
        split_vals = obj.meta_data[split_by].reindex(cells).astype(str).to_numpy()
        facets = [
            (lvl, np.where(split_vals == lvl)[0])
            for lvl in sorted(set(split_vals))
        ]

    nrow, ncol = _subplot_grid(len(facets))
    if figsize is None:
        figsize = (5.5 * ncol, 4.0 * nrow)
    fig, axes = plt.subplots(nrow, ncol, figsize=figsize, squeeze=False)
    axes_flat = axes.ravel()
    rng = np.random.default_rng(seed)

    for ax, (level, idx) in zip(axes_flat, facets):
        sub_p, sub_g = pvec[idx], group[idx]

        # Densities first — their peak sets the band the jitter strip occupies,
        # mirroring R's read of the built panel's y-range.
        top = 0.0
        curves = []
        for cls in classes:
            vals = sub_p[sub_g == cls]
            if vals.size < 2 or np.ptp(vals) <= 0:
                continue
            grid = np.linspace(sub_p.min(), sub_p.max(), 256)
            dens = gaussian_kde(vals)(grid)
            curves.append((cls, grid, dens))
            top = max(top, float(dens.max()))
        for cls, grid, dens in curves:
            ax.plot(grid, dens, color=colours[cls], linewidth=1.5, label=cls)

        if top <= 0:
            top = 1.0
        band = top / 10.0
        for cls in classes:
            vals = sub_p[sub_g == cls]
            if vals.size == 0:
                continue
            # Controls sit above the axis, the target gene's cells below.
            if cls == nt_name:
                y = rng.uniform(0.001, band, size=vals.size)
            else:
                y = rng.uniform(-band, 0.0, size=vals.size)
            ax.scatter(vals, y, color=colours[cls], s=1.0, linewidths=0)

        ax.axhline(0.0, color="black", linewidth=0.5)
        ax.set_xlabel("perturbation score")
        ax.set_ylabel("Cell density")
        if level is not None:
            ax.set_title(str(level), fontsize=_fs("title"))
        _strip_axes(ax)
        ax.legend(frameon=False, fontsize=_fs("label"))

    for ax in axes_flat[len(facets):]:
        ax.axis("off")

    fig.suptitle(f"Perturbation score — {target_gene_ident}",
                 fontsize=_fs("suptitle"), fontweight="bold")
    fig.tight_layout()
    return fig

mixscape_heatmap

mixscape_heatmap(obj, ident_1: str, ident_2: Optional[str] = None, balanced: bool = True, logfc_threshold: float = 0.25, assay: str = 'RNA', max_genes: int = 100, test_use: str = 'wilcox', max_cells_group: Optional[int] = None, order_by_prob: bool = True, mixscape_class: str = 'mixscape_class', prtb_type: str = 'KO', fc_name: str = 'avg_log2FC', pval_cutoff: float = 0.05, seed: int = 0, **kwargs) -> 'Figure'

DE heatmap with cells ordered by knockout probability (MixscapeHeatmap).

Mirrors MixscapeHeatmap(object, ident.1 = "NT", ident.2 = "IFNGR2 KO", balanced = TRUE, max.genes = 20). Where plot_perturb_score shows the one-dimensional score mixscape split on, this shows the genes underneath it: the DE genes between the two classes, with every cell ordered left-to-right by its knockout posterior. Read together with the class colour bar, a clean screen shows the expression block turning on in step with the probability — the escapers at the low-probability end still looking like control.

ident_1 / ident_2 are mixscape_class levels (e.g. "NT", "IFNGR2 KO", "IFNGR2 NP"), which run_mixscape also leaves as the active identity.

Parameters:

• ident_1 (str) –

  the two classes to contrast (ident_2=None → all others).

• ident_2 (str) –

  the two classes to contrast (ident_2=None → all others).

• balanced (bool, default: True ) –

  take up to max_genes genes from each direction of the fold change; otherwise only up-regulated ones.

• max_genes (int, default: 100 ) –

  cap on DE genes per direction.

• max_cells_group (Optional[int], default: None ) –

  downsample each class to this many cells.

• order_by_prob (bool, default: True ) –

  order cells by <mixscape_class>_p_<type>, highest first. If False, cells are shuffled (as R does).

• mixscape_class (str, default: 'mixscape_class' ) –

  metadata column of mixscape classifications; also names the posterior column read for the ordering.

• prtb_type (str, default: 'KO' ) –

  perturbation label used by run_mixscape ("KO").

• fc_name (str, default: 'avg_log2FC' ) –

  fold-change column in the DE table ("avg_log2FC").

• seed (int, default: 0 ) –

  random state for downsampling / shuffling (determinism).

• **kwargs –

  forwarded to do_heatmap.

Returns:

• matplotlib Figure –

Source code in truecell/plotting.py

def mixscape_heatmap(
    obj,
    ident_1: str,
    ident_2: Optional[str] = None,
    balanced: bool = True,
    logfc_threshold: float = 0.25,
    assay: str = "RNA",
    max_genes: int = 100,
    test_use: str = "wilcox",
    max_cells_group: Optional[int] = None,
    order_by_prob: bool = True,
    mixscape_class: str = "mixscape_class",
    prtb_type: str = "KO",
    fc_name: str = "avg_log2FC",
    pval_cutoff: float = 0.05,
    seed: int = 0,
    **kwargs,
) -> "Figure":
    """DE heatmap with cells ordered by knockout probability (``MixscapeHeatmap``).

    Mirrors ``MixscapeHeatmap(object, ident.1 = "NT", ident.2 = "IFNGR2 KO",
    balanced = TRUE, max.genes = 20)``. Where :func:`plot_perturb_score` shows the
    one-dimensional score mixscape split on, this shows the genes underneath it:
    the DE genes between the two classes, with every cell ordered left-to-right by
    its knockout posterior. Read together with the class colour bar, a clean
    screen shows the expression block turning on in step with the probability —
    the escapers at the low-probability end still looking like control.

    ``ident_1`` / ``ident_2`` are ``mixscape_class`` levels (e.g. ``"NT"``,
    ``"IFNGR2 KO"``, ``"IFNGR2 NP"``), which :func:`run_mixscape` also leaves as
    the active identity.

    Parameters
    ----------
    ident_1, ident_2 : the two classes to contrast (``ident_2=None`` → all others).
    balanced         : take up to ``max_genes`` genes from *each* direction of the
                       fold change; otherwise only up-regulated ones.
    max_genes        : cap on DE genes per direction.
    max_cells_group  : downsample each class to this many cells.
    order_by_prob    : order cells by ``<mixscape_class>_p_<type>``, highest first.
                       If False, cells are shuffled (as R does).
    mixscape_class   : metadata column of mixscape classifications; also names the
                       posterior column read for the ordering.
    prtb_type        : perturbation label used by :func:`run_mixscape` (``"KO"``).
    fc_name          : fold-change column in the DE table (``"avg_log2FC"``).
    seed             : random state for downsampling / shuffling (determinism).
    **kwargs         : forwarded to :func:`do_heatmap`.

    Returns
    -------
    matplotlib Figure
    """
    from .markers import find_markers
    from .preprocessing import scale_data

    if mixscape_class not in obj.meta_data.columns:
        raise KeyError(
            f"Column {mixscape_class!r} not in meta_data — run run_mixscape first."
        )

    # find_markers resolves ident_1/ident_2 against the active identity, so drive
    # the DE off the mixscape classes — restored afterwards.
    saved_ident = pd.Categorical(list(obj.idents))
    obj.idents = list(obj.meta_data[mixscape_class].astype(str))
    try:
        markers = find_markers(
            obj,
            ident_1=ident_1,
            ident_2=ident_2,
            assay=assay,
            test_use=test_use,
            only_pos=False,
            logfc_threshold=logfc_threshold,
        )
    finally:
        obj.idents = saved_ident

    if markers.empty or fc_name not in markers.columns:
        raise ValueError(
            f"No DE genes between {ident_1!r} and {ident_2 or 'the rest'!r}; "
            "nothing to plot."
        )

    sig = markers[markers["p_val"] < pval_cutoff]
    pos = list(sig.index[sig[fc_name] > logfc_threshold])[:max_genes]
    neg = list(sig.index[sig[fc_name] < -logfc_threshold])[:max_genes] if balanced else []
    marker_list = pos + neg
    if not marker_list:
        raise ValueError(
            f"No DE genes passed p_val < {pval_cutoff} and "
            f"|{fc_name}| > {logfc_threshold}; nothing to plot."
        )

    idents = [ident_1] + ([ident_2] if ident_2 is not None else [])
    klass = obj.meta_data[mixscape_class].astype(str)
    if ident_2 is None:
        keep = list(obj.meta_data.index[klass == ident_1])
        keep += list(obj.meta_data.index[klass != ident_1])
    else:
        keep = list(obj.meta_data.index[klass.isin(idents)])
    if not keep:
        raise ValueError(f"No cells in classes {idents}.")

    sub = obj.subset(cells=keep)
    rng = np.random.default_rng(seed)

    if max_cells_group is not None:
        sub_klass = sub.meta_data[mixscape_class].astype(str)
        picked: list[str] = []
        for lvl in sorted(set(sub_klass)):
            lvl_cells = list(sub.meta_data.index[sub_klass == lvl])
            if len(lvl_cells) > max_cells_group:
                lvl_cells = list(rng.choice(lvl_cells, max_cells_group, replace=False))
            picked += lvl_cells
        sub = sub.subset(cells=picked)

    scale_data(sub, features=marker_list, assay=assay)

    if order_by_prob:
        p_col = f"{mixscape_class}_p_{prtb_type.lower()}"
        if p_col not in sub.meta_data.columns:
            raise KeyError(
                f"Posterior column {p_col!r} not in meta_data — run run_mixscape "
                f"with prtb_type={prtb_type!r}, or pass order_by_prob=False."
            )
        # NT cells carry no posterior; R initialises the column to 0, so they sort
        # to the low-probability end rather than dropping out.
        probs = sub.meta_data[p_col].to_numpy(dtype=float)
        probs = np.nan_to_num(probs, nan=0.0)
        order = np.argsort(-probs, kind="stable")
    else:
        order = rng.permutation(len(sub.cell_names()))
    ordered_cells = [sub.cell_names()[i] for i in order]

    return do_heatmap(
        sub,
        features=marker_list,
        group_by=mixscape_class,
        assay=assay,
        cells=ordered_cells,
        **kwargs,
    )