Read the dose response per feature, and find the window it happens in

docs/api.md

@@ -139,6 +139,9 @@ Writing to a layer suffixes the column, so `key_added="sphered"` flags `var["deg
     tl.edistance
     tl.transport
     tl.dose_response
+    tl.dose_features
+    tl.dose_direction
+    tl.dose_trajectory
     tl.nn_moa_classify
     tl.moa_enrichment
     tl.feature_sets
@@ -170,6 +173,9 @@ Writing to a layer suffixes the column, so `key_added="sphered"` flags `var["deg
 | `tl.edistance` | `uns["mantispy"][key_added]`, or `..._pairwise` when `reference=None` |
 | `tl.transport` | `uns["mantispy"][key_added]` and `..._units`, `obs[key_added + "_agreement"]` |
 | `tl.dose_response` | `uns["mantispy"][key_added]` |
+| `tl.dose_features` | `uns["mantispy"][key_added]`, one row per compound and feature |
+| `tl.dose_direction` | `uns["mantispy"][key_added]`, one row per compound and concentration, `obs[key_added + "_phase"]` |
+| `tl.dose_trajectory` | returns a new object: compounds by features-and-positions at `"perturbation"` resolution |
 | `tl.nn_moa_classify` | `obs[key_added + "_predicted"]`, `uns["mantispy"][key_added]` and `..._confusion` |
 | `tl.moa_enrichment` | `uns["mantispy"][key_added]` |
 | `tl.feature_sets` | returns a decoupler network; stores nothing |
@@ -265,6 +271,7 @@ pycytominer and similar tools expect.
     pl.effect_sizes
     pl.feature_volcano
     pl.dose_response
+    pl.dose_direction
     pl.moa_confusion
     pl.moa_enrichment
     pl.distance_heatmap

spec/schema-1.0.json

@@ -21,7 +21,7 @@
         "Metadata_Perturbation",
         "Metadata_Compound",
         "Metadata_Concentration",
-        "Metadata_ConcentrationNominal",
+        "Metadata_ConcentrationRecorded",
         "Metadata_MOA",
         "Metadata_CellLine",
         "Metadata_Control",

src/mantispy/_core/_reduce.py

@@ -38,13 +38,20 @@
     "get_matrix",
     "group_codes",
     "group_offsets",
+    "group_rows",
     "iter_groups",
     "reduce_grouped",
     "representation",
     "transform_grouped",
 ]
 
 
+def group_rows(codes: np.ndarray, n_groups: int) -> list[np.ndarray]:
+    """Row indices of each group, taken from one stable ordering rather than by scanning the codes per group."""
+    order, offsets = group_offsets(codes, n_groups)
+    return [order[offsets[group] : offsets[group + 1]] for group in range(n_groups)]
+
+
 def get_matrix(adata: AnnData, layer: str | None = None, rows: np.ndarray | None = None) -> np.ndarray:
     """Return the requested matrix as a dense ``float32`` array.
 

src/mantispy/_core/schema.py

@@ -47,8 +47,9 @@
     "Metadata_Perturbation",
     "Metadata_Compound",
     "Metadata_Concentration",
-    # The dose a well was meant to get, where a plate map records the same one to several precisions.
-    "Metadata_ConcentrationNominal",
+    # What the plate map wrote, where it records one dose to several precisions and Metadata_Concentration
+    # holds the one the well was meant to get.
+    "Metadata_ConcentrationRecorded",
     "Metadata_MOA",
     "Metadata_CellLine",
     "Metadata_Control",

src/mantispy/ds/_datasets.py

@@ -421,26 +421,45 @@ def chroma(cache_dir: str | Path | None = None, **kwargs: Any) -> AnnData:
 }
 
 
-#: Doses closer than this are one nominal concentration a plate map wrote to two precisions.
-_DOSE_TOLERANCE = 0.01
+def _decimals(value: float) -> int:
+    """How many decimal places a level was written with, from its shortest exact repr."""
+    return len(np.format_float_positional(value, trim="-").partition(".")[2])
 
 
-def _aligned_doses(doses: pd.Series) -> pd.Series:
+def _aligned_doses(doses: pd.Series, compounds: pd.Series) -> pd.Series:
     """The dose each well was meant to get, where plate maps record one concentration to several precisions.
 
-    Two of the OASIS plate maps write the dose to three decimals and the rest to four, so one concentration
-    arrives as both ``3.704`` and ``3.7037``. Anything grouping by dose then splits a compound's replicates
-    across two levels, and :func:`~mantispy.tl.dose_response` reports more distinct doses than were plated.
-    Levels within ``_DOSE_TOLERANCE`` of each other are one level, named by the value the most wells carry.
-    Rounding cannot do this: no fixed precision separates the pairs that differ from those that do not.
+    The OASIS plate maps disagree on precision rather than on value: one batch writes the concentration the
+    dilution actually produced, ``0.0152416`` uM, and another writes it rounded, ``0.015``. A level written to
+    fewer decimals is therefore the same dose as the finer level that rounds to it, which is a statement about
+    how the number was recorded rather than a tolerance fitted to the data.
+
+    Matching runs within a compound, because a compound's levels are one dilution series and cannot collide.
+    Across the whole plate map they can: berberine's 25 uM and the main ladder's 33.3 uM are a third apart and
+    genuinely different doses, closer together than ``0.000762`` and ``0.001`` are, which are one dose.
+
+    Rounding every level to a fixed precision cannot do this, and neither can a relative tolerance. The one case
+    it would get wrong is a ladder with two rungs inside a rounding step of each other, which a series coarser
+    than two-fold never has.
     """
-    counts = doses.value_counts()
-    levels = np.sort(counts.index.to_numpy(dtype=float))
-    # A relative tolerance says nothing about an undosed well, so zero is its own level and passes through.
-    levels = levels[levels > 0]
-    runs = np.split(levels, np.flatnonzero(levels[1:] / levels[:-1] - 1 >= _DOSE_TOLERANCE) + 1)
-    lookup = {level: max(run, key=lambda value: (counts[value], value)) for run in runs for level in run}
-    return doses.replace(lookup)
+
+    def align(block: pd.Series) -> pd.Series:
+        levels = np.sort(block.dropna().unique())
+        levels = levels[levels > 0]
+        places = [_decimals(level) for level in levels]
+        lookup = {}
+        for level, digits in zip(levels, places, strict=True):
+            finer = [
+                other
+                for other, deeper in zip(levels, places, strict=True)
+                if deeper > digits and round(other, digits) == level
+            ]
+            # The nearest one: rounding 0.006 up to 0.01 must not claim a level that was written as 0.01.
+            lookup[level] = min(finer, key=lambda other: abs(other - level), default=level)
+        return block.replace(lookup)
+
+    # dropna=False, so a well whose compound is blank keeps the dose it was recorded with.
+    return doses.groupby(compounds, observed=True, dropna=False).transform(align)
 
 
 def _oasis_platemaps(cache_dir: str | Path | None) -> pd.DataFrame:
@@ -468,8 +487,9 @@ def _oasis_platemaps(cache_dir: str | Path | None) -> pd.DataFrame:
             kept["Metadata_Compound"] = kept["Metadata_Compound"].fillna(frame["BROAD_ID"])
         frames.append(kept)
     platemap = pd.concat(frames, ignore_index=True)
-    platemap["Metadata_Concentration"] = pd.to_numeric(platemap["Metadata_Concentration"], errors="coerce")
-    platemap["Metadata_ConcentrationNominal"] = _aligned_doses(platemap["Metadata_Concentration"])
+    recorded = pd.to_numeric(platemap["Metadata_Concentration"], errors="coerce")
+    platemap["Metadata_ConcentrationRecorded"] = recorded
+    platemap["Metadata_Concentration"] = _aligned_doses(recorded, platemap["Metadata_Compound"])
     # One batch writes the line as HepRG and the others as HepaRG; two spellings would split every per-line grouping.
     platemap["Metadata_CellLine"] = platemap["Metadata_CellLine"].replace({"HepRG": "HepaRG"})
     return platemap
@@ -501,10 +521,12 @@ def oasis_pilot(annotate: bool = True, cache_dir: str | Path | None = None, **kw
         The assay-development batch doses DMSO itself, so a control well there carries a concentration.
         ``Metadata_Control`` marks the compound, not the dose.
 
-        Two of the plate maps write the dose to three decimals and the rest to four, so one concentration is
-        recorded as both ``3.704`` and ``3.7037``. ``Metadata_Concentration`` keeps what was recorded;
-        ``Metadata_ConcentrationNominal`` puts the levels that agree to within 1% onto the value the most wells
-        carry, and names the replicate groups. Group by the raw column and a treatment's wells split in two.
+        The plate maps disagree on precision: one batch writes the concentration the dilution produced,
+        ``0.0152416`` uM, and another writes it rounded, ``0.015``. ``Metadata_Concentration`` is the dose the
+        well was meant to get, reading a coarser spelling as the finer level it rounds to within each compound,
+        and it names the replicate groups. ``Metadata_ConcentrationRecorded`` keeps what the plate map wrote:
+        read against that column, every dosed compound here carries eighteen levels where ten were plated, and a
+        treatment's wells split across two spellings.
     """
     adata = _profiles("oasis_pilot", cache_dir, select=lambda name: name.endswith(".csv.gz"), **kwargs)
     if not annotate:
@@ -517,7 +539,7 @@ def oasis_pilot(annotate: bool = True, cache_dir: str | Path | None = None, **kw
         get_logger().warning("oasis_pilot: %d of %d wells have no plate-map row", unmatched, len(merged))
     adata.obs["Metadata_Compound"] = merged["Metadata_Compound"].to_numpy()
     adata.obs["Metadata_Concentration"] = merged["Metadata_Concentration"].to_numpy(dtype=float)
-    adata.obs["Metadata_ConcentrationNominal"] = merged["Metadata_ConcentrationNominal"].to_numpy(dtype=float)
+    adata.obs["Metadata_ConcentrationRecorded"] = merged["Metadata_ConcentrationRecorded"].to_numpy(dtype=float)
     adata.obs["Metadata_CellLine"] = merged["Metadata_CellLine"].to_numpy()
     adata.obs["Metadata_Control"] = merged["Metadata_Compound"].astype(str).str.upper().eq("DMSO").to_numpy()
     # Replicates share a compound at a concentration, which is what the mode= shorthands of mt.tl.map compare.
@@ -526,15 +548,15 @@ def oasis_pilot(annotate: bool = True, cache_dir: str | Path | None = None, **kw
         np.where(
             is_control,
             "DMSO",
-            merged["Metadata_Compound"].astype(str) + "@" + merged["Metadata_ConcentrationNominal"].astype(str),
+            merged["Metadata_Compound"].astype(str) + "@" + merged["Metadata_Concentration"].astype(str),
         )
     )
     get_logger().info(
         "OASIS pilot: %d wells x %d features, %d compounds over %d concentrations, %d control wells",
         adata.n_obs,
         adata.n_vars,
         int(merged.loc[~is_control, "Metadata_Compound"].nunique()),
-        int(merged["Metadata_ConcentrationNominal"].nunique()),
+        int(merged["Metadata_Concentration"].nunique()),
         int(is_control.sum()),
     )
     return adata

src/mantispy/pl/__init__.py

@@ -9,7 +9,7 @@
 from mantispy.pl._evaluation import batch_variance, map, metrics, replicate_correlation, similarity
 from mantispy.pl._features import feature_correlation, feature_groups
 from mantispy.pl._heterogeneity import cell_cycle, cluster_composition, density, subpopulation_hits
-from mantispy.pl._hits import dose_response, effect_sizes, feature_volcano, hits
+from mantispy.pl._hits import dose_direction, dose_response, effect_sizes, feature_volcano, hits
 from mantispy.pl._moa import distance_heatmap, moa_confusion, moa_enrichment, pathway_coherence, sets_heatmap
 from mantispy.pl._plate import plate
 from mantispy.pl._qc import cell_counts, cytotoxicity, feature_distributions, nan_matrix, qc, replicate_saturation
@@ -26,6 +26,7 @@
     "cytotoxicity",
     "density",
     "distance_heatmap",
+    "dose_direction",
     "dose_response",
     "effect_sizes",
     "feature_correlation",

src/mantispy/pl/_hits.py

@@ -9,6 +9,7 @@
 from mantispy._core.frames import as_frame
 from mantispy.pl._common import axes as _axes
 from mantispy.pl._common import table as _table
+from mantispy.tl._dose import DOSE_PHASES
 
 if TYPE_CHECKING:
     import pandas as pd
@@ -75,11 +76,17 @@ def hits(adata: AnnData, key: str = "hits", label_top: int = 10, ax: Axes | None
     return ax
 
 
+def _rows_for(table: pd.DataFrame, column: str, value: str, key: str) -> pd.DataFrame:
+    """The table's rows for one group, or a KeyError naming the groups it does hold."""
+    selected = table[table[column].astype(str) == str(value)]
+    if selected.empty:
+        raise KeyError(f"no {column} {value!r} in uns['mantispy'][{key!r}]; it holds {sorted(set(table[column]))[:5]}")
+    return selected
+
+
 def _effects(adata: AnnData, group: str, key: str) -> tuple[pd.DataFrame, pd.Series | None]:
     table = _table(adata, key, "mt.tl.effect_size")
-    selected = table[table["group"].astype(str) == str(group)]
-    if selected.empty:
-        raise KeyError(f"no group {group!r} in uns['mantispy'][{key!r}]; it holds {sorted(set(table['group']))[:5]}")
+    selected = _rows_for(table, "group", group, key)
     var = as_frame(adata.var)
     families = var["feature_group"].astype(str) if "feature_group" in var else None
     return selected, families
@@ -210,9 +217,7 @@ def dose_response(
     from mantispy.tl._dose import four_parameter_logistic
 
     table = _table(adata, key, "mt.tl.dose_response")
-    row = table[table["compound"].astype(str) == str(compound)]
-    if row.empty:
-        raise KeyError(f"no compound {compound!r} in uns['mantispy'][{key!r}]")
+    row = _rows_for(table, "compound", compound, key)
 
     if response is None:
         response = _recorded_response(adata, "dose_response", "hits_row_distance")
@@ -248,3 +253,74 @@ def dose_response(
     ax.set_title(f"{compound}  (spearman {float(fitted['spearman']):.2f})", fontsize=9)
     ax.legend(fontsize=7)
     return ax
+
+
+#: Background colour of each phase, keyed by :data:`~mantispy.tl._dose.DOSE_PHASES` so the two cannot drift.
+#: Grey where nothing happens, warm where it does, green where it has arrived, red where the cells are gone.
+DOSE_PHASE_COLOURS = dict(zip(DOSE_PHASES, ("#f2f2f2", "#fde6c4", "#dbe8d4", "#f6d2d2"), strict=True))
+
+
+def dose_direction(
+    adata: AnnData,
+    compound: str,
+    key: str = "dose_direction",
+    ax: Axes | None = None,
+) -> Axes:
+    """One compound's ladder, with the background banded by what each concentration is doing.
+
+    The solid line is how far the profile sits from the controls, and the dashed line is how far it moved from
+    the concentration below it. The second is what says where the action is: a response that is still changing
+    has a large step, one that has arrived has a small one however high the solid line sits. The two dotted
+    horizontals are the floors those lines are read against, which control wells laid out the same way reach.
+
+    Args:
+        adata: Object holding the table :func:`~mantispy.tl.dose_direction` wrote.
+        compound: Which compound of that table to draw.
+        key: Name of that table in ``uns["mantispy"]``.
+        ax: Axes to draw on, or ``None`` for a new figure.
+
+    Returns:
+        The axes drawn on.
+
+    Raises:
+        KeyError: There is no such table, or it holds no such compound.
+    """
+    table = _table(adata, key, "mt.tl.dose_direction")
+    block = _rows_for(table, "compound", compound, key).sort_values("dose")
+
+    ax = _axes(ax, (5.2, 3.6))
+    doses = block["dose"].to_numpy(dtype=float)
+    # Each concentration owns the ladder up to halfway to its neighbours, measured in log dose, and half a step
+    # past the two ends.
+    log_dose = np.log10(doses)
+    gaps = np.diff(log_dose) if len(doses) > 1 else np.array([0.6])
+    padded = np.concatenate([[log_dose[0] - gaps[0]], log_dose, [log_dose[-1] + gaps[-1]]])
+    edges = 10.0 ** ((padded[:-1] + padded[1:]) / 2)
+    for left, right, phase in zip(edges[:-1], edges[1:], block["phase"], strict=True):
+        ax.axvspan(left, right, color=DOSE_PHASE_COLOURS.get(str(phase), "#ffffff"), lw=0, zorder=0)
+
+    floor = float(np.nanmedian(block["amplitude_null"].to_numpy(dtype=float)))
+    ax.axhline(floor, ls=":", lw=1, color="0.45", zorder=1)
+    ax.axhline(floor * np.sqrt(2.0), ls=":", lw=1, color="0.65", zorder=1)
+    ax.plot(
+        doses, block["amplitude"], marker="o", ms=4, lw=1.6, color="#2a4d69", label="distance from controls", zorder=3
+    )
+    ax.plot(
+        doses,
+        block["step_amplitude"],
+        marker="s",
+        ms=3,
+        lw=1.3,
+        ls="--",
+        color="#c1611f",
+        label="moved since the last",
+        zorder=3,
+    )
+
+    ax.set_xscale("log")
+    ax.set_xlim(edges[0], edges[-1])
+    ax.set_xlabel("concentration")
+    ax.set_ylabel("MADs per feature")
+    ax.set_title(str(compound), fontsize=10)
+    ax.legend(fontsize=7, frameon=False, loc="upper left")
+    return ax

src/mantispy/pp/_batch.py

@@ -32,7 +32,7 @@ def _median_polish_stack(grids: np.ndarray, max_iter: int, tol: float) -> tuple[
     """Median polish every feature of a ``(rows, columns, features)`` stack.
 
     Each feature's grid is independent, so the stack goes to one numba kernel that polishes a plane per thread.
-    A per-feature Python loop took 20 minutes on 132 JUMP plates, almost all of it interpreter and pandas overhead.
+    A per-feature Python loop spends almost all of its time in the interpreter and in pandas instead.
 
     Returns the fitted ``(row_effects, column_effects)``, both ``(positions, features)``.
     The grand level is not included, so subtracting the effects keeps each feature's level, as :func:`regress_out` does.

src/mantispy/tl/__init__.py

@@ -5,7 +5,7 @@
 from mantispy.tl._design import cytotoxicity, replicate_saturation
 from mantispy.tl._differential import differential_features
 from mantispy.tl._distance import edistance
-from mantispy.tl._dose import dose_response
+from mantispy.tl._dose import dose_direction, dose_features, dose_response, dose_trajectory
 from mantispy.tl._effect import effect_size, wasserstein_features
 from mantispy.tl._enrich import enrich, feature_sets, rank_features, rank_sets
 from mantispy.tl._heterogeneity import (
@@ -28,7 +28,10 @@
     "cluster_composition",
     "consensus",
     "cytotoxicity",
+    "dose_direction",
+    "dose_features",
     "dose_response",
+    "dose_trajectory",
     "differential_features",
     "edistance",
     "effect_size",

src/mantispy/tl/_aggregate.py

@@ -59,7 +59,7 @@ def aggregate(
         ValueError: ``func`` is not one of ``FUNCTIONS``.
 
     Notes:
-        This uses mantispy's own NaN-skipping kernel rather than :func:`scanpy.get.aggregate`, which propagates NaN and is measurably slower on both mean and median.
+        This uses mantispy's own NaN-skipping kernel rather than :func:`scanpy.get.aggregate`, which propagates NaN and is slower on both mean and median.
     """
     if func not in FUNCTIONS:
         raise ValueError(f"func must be one of {tuple(FUNCTIONS)}, got {func!r}")
@@ -120,7 +120,7 @@ def aggregate(
 
 def _site_counts(frame: pd.DataFrame, codes: np.ndarray, n_groups: int) -> np.ndarray:
     """Distinct fields of view among each group's cells, where site 1 of one well and of the next are different fields."""
-    # One integer per field, built from per-column codes rather than a MultiIndex of tuples, which is 7x slower.
+    # One integer per field, built from per-column codes rather than a MultiIndex, which would build a tuple per cell.
     field = np.zeros(len(frame), dtype=np.int64)
     for column in ("Metadata_Plate", "Metadata_Well", "Metadata_Site"):
         if column in frame: