Read the control scale off the controls that did not fit the transform

docs/api.md

@@ -166,7 +166,7 @@ Writing to a layer suffixes the column, so `key_added="sphered"` flags `var["deg
 | `tl.consensus` | returns a new object at `"perturbation"` resolution; `obs["Metadata_ReplicateCount"]`, `uns["mantispy"]["consensus_weights"]` |
 | `tl.effect_size` | `varm[key_added]`, `uns["mantispy"][key_added]` and `..._groups` |
 | `tl.wasserstein_features` | `varm[key_added]`, `uns["mantispy"][key_added]` and `..._groups` |
-| `tl.hit_calling` | `uns["mantispy"][key_added]`, `obs[key_added + "_distance"]`, `obs[key_added + "_row_distance"]`, `obs[key_added + "_qvalue"]` |
+| `tl.hit_calling` | `uns["mantispy"][key_added]`, `obs[key_added + "_distance"]`, `obs[key_added + "_row_distance"]`, `obs[key_added + "_qvalue"]`, `obs[key_added + "_reference_held_out"]` |
 | `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]` |

spec/schema-1.0.json

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

src/mantispy/_core/masks.py

@@ -89,3 +89,17 @@ def reference_mask(adata: AnnData, reference: str | None) -> np.ndarray:
         "a control row",
         " Fill them, or check that the platemap covers every well.",
     )
+
+
+def held_out_reference(adata: AnnData, is_control: np.ndarray, distance_key: str) -> np.ndarray:
+    """Narrow a control mask to the rows that did not fit the transform ``distance_key`` is measured in.
+
+    A row that fitted the centroid and the covariance sits closer to the centroid than one that did not, so a scale read off every control comes out low.
+    :func:`~mantispy.tl.hit_calling` records the half it held out under the ``key_added`` the distance column already carries.
+    A missing column is not an error, on the same terms :func:`feature_mask` applies: the distance came from somewhere else, and every control is kept.
+    """
+    column = f"{distance_key.removesuffix('_row_distance')}_reference_held_out"
+    if column not in adata.obs:
+        get_logger().debug("obs has no column %r; reading the control scale off every control row", column)
+        return is_control
+    return is_control & reference_mask(adata, column)

src/mantispy/_core/schema.py

@@ -47,6 +47,8 @@
     "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",
     "Metadata_MOA",
     "Metadata_CellLine",
     "Metadata_Control",

src/mantispy/ds/_datasets.py

@@ -421,6 +421,28 @@ 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 _aligned_doses(doses: 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.
+    """
+    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 _oasis_platemaps(cache_dir: str | Path | None) -> pd.DataFrame:
     """The plate maps of every OASIS batch, read down to plate, well, compound and concentration."""
     frames = []
@@ -447,6 +469,7 @@ def _oasis_platemaps(cache_dir: str | Path | None) -> pd.DataFrame:
         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"])
     # 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
@@ -477,6 +500,11 @@ 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.
     """
     adata = _profiles("oasis_pilot", cache_dir, select=lambda name: name.endswith(".csv.gz"), **kwargs)
     if not annotate:
@@ -489,6 +517,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_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.
@@ -497,15 +526,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_Concentration"].astype(str),
+            merged["Metadata_Compound"].astype(str) + "@" + merged["Metadata_ConcentrationNominal"].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_Concentration"].nunique()),
+        int(merged["Metadata_ConcentrationNominal"].nunique()),
         int(is_control.sum()),
     )
     return adata

src/mantispy/pl/_hits.py

@@ -201,7 +201,8 @@ def dose_response(
         ax: Axes to draw on, or ``None`` for a new figure.
 
     Returns:
-        The axes drawn on.
+        The axes drawn on, with a log ``y`` scale where every drawn response is positive, since a distance from the controls has a long right tail and one stray well would otherwise flatten the rest onto the baseline.
+        The scale is decided per axes, so panels drawn side by side can differ; set it on the returned axes to compare them.
 
     Raises:
         KeyError: There is no such table, or it holds no such compound.
@@ -225,6 +226,10 @@ def dose_response(
     ax = _axes(ax, (5, 4))
     ax.scatter(doses[usable], values[usable], s=18, label="wells")
     ax.set_xscale("log")
+    # A distance from the controls has a long right tail, and one stray well is enough to flatten every
+    # other point onto the baseline. A response that reaches zero or below is drawn on a linear scale.
+    if usable.any() and values[usable].min() > 0:
+        ax.set_yscale("log")
 
     fitted = row.iloc[0]
     if bool(fitted["fit_ok"]):

src/mantispy/tl/_design.py

@@ -12,7 +12,7 @@
 from mantispy._core._reduce import group_codes, group_offsets, representation
 from mantispy._core.frames import as_frame
 from mantispy._core.logging import get_logger
-from mantispy._core.masks import reference_mask
+from mantispy._core.masks import held_out_reference, reference_mask
 from mantispy._core.mutation import inplace_or_copy
 
 
@@ -177,7 +177,7 @@ def cytotoxicity(
     reference: str | None = "negcon",
     count_key: str = "Metadata_CellCount",
     site_key: str | None = "Metadata_SiteCount",
-    distance_key: str = "hits_distance",
+    distance_key: str = "hits_row_distance",
     min_viability: float = 0.7,
     key_added: str = "cytotoxicity",
     copy: bool = False,
@@ -191,7 +191,7 @@ def cytotoxicity(
         count_key: ``obs`` column holding the cell count.
         site_key: ``obs`` column holding the number of fields of view that count covers.
             Where present, viability compares cells per field, so a well missing a field does not read as cell loss. ``None`` compares the counts as they are.
-        distance_key: ``obs`` column holding the distance from the controls, as written by :func:`~mantispy.tl.hit_calling`.
+        distance_key: ``obs`` column holding the per-row distance from the controls, as written by :func:`~mantispy.tl.hit_calling`. Its group-level sibling ``hits_distance`` is one number repeated over each group's rows, so the median below would return the value it was handed.
         min_viability: Fraction of the control cell count below which a group counts as having lost cells.
         key_added: Name for the outputs.
         copy: Return a modified copy instead of mutating in place.
@@ -232,7 +232,7 @@ def cytotoxicity(
     if distance_key not in obs:
         raise KeyError(
             f"obs has no column {distance_key!r}; run mt.tl.hit_calling first, which writes "
-            "obs['hits_distance'], or name another column"
+            "obs['hits_row_distance'], or name another column"
         )
 
     is_control = reference_mask(adata, reference)
@@ -241,7 +241,7 @@ def cytotoxicity(
         counts = counts / np.maximum(obs[site_key].to_numpy(dtype=float), 1)
     distances = obs[distance_key].to_numpy(dtype=float)
     control_count = float(np.nanmedian(counts[is_control]))
-    control_distance = float(np.nanmedian(distances[is_control]))
+    control_distance = float(np.nanmedian(distances[held_out_reference(adata, is_control, distance_key)]))
     if not np.isfinite(control_count) or control_count <= 0:
         raise ValueError(f"the reference rows have no usable {count_key!r} to normalize viability against")
 

src/mantispy/tl/_dose.py

@@ -13,7 +13,7 @@
 from mantispy._core._stats import MAD_TO_SIGMA, benjamini_hochberg
 from mantispy._core.frames import as_frame
 from mantispy._core.logging import get_logger
-from mantispy._core.masks import reference_mask
+from mantispy._core.masks import held_out_reference, reference_mask
 from mantispy._core.mutation import inplace_or_copy
 
 #: Column order of the output table, so an empty result still carries its columns.
@@ -278,7 +278,7 @@ def _baseline_and_cutoff(
     # and the trend need no controls, so only the hit call is left out.
     if reference is not None and not (reference == "negcon" and "Metadata_Control" not in adata.obs):
         values = as_frame(adata.obs)[response].to_numpy(dtype=float)
-        control = values[reference_mask(adata, reference)]
+        control = values[held_out_reference(adata, reference_mask(adata, reference), response)]
         control = control[np.isfinite(control)]
 
     if control.size < 2:
@@ -341,7 +341,7 @@ def dose_response(
         min_doses: Distinct doses below which the curve is skipped and only the trend is reported.
         min_r_squared: Coefficient of determination a fit needs before it is marked ok.
         reference: Rows that set the baseline the response is read against and the spread the cutoff comes from. ``None`` leaves the hit call out unless ``cutoff`` is given.
-        cutoff: Response a curve has to clear to count as active. The default takes three times the controls' MAD, the ToxCast pipeline's ``3 * bmad``.
+        cutoff: Response a curve has to clear to count as active. The default takes three times the controls' MAD, the ToxCast pipeline's ``3 * bmad``, over the controls :func:`~mantispy.tl.hit_calling` held out of its own fit rather than over all of them.
         key_added: Name for the output table.
         copy: Return a modified copy instead of mutating in place.
 

src/mantispy/tl/_hits.py

@@ -109,6 +109,7 @@ def hit_calling(
         Also writes ``obs[key_added + "_row_distance"]``, each row's own distance from the control centroid rather than its group's.
         That is the column a dose-response fit wants: the group statistic is one number repeated over the group's rows, so the controls show no spread and nothing downstream can read a scale off them.
         Rows that fitted the centroid sit a little closer to it than the held-out controls do, by the same split the Notes describe.
+        Which rows those are is written to ``obs[key_added + "_reference_held_out"]``, true for the controls that did not fit, so that a scale taken from the controls can be taken from the honest half.
 
     Raises:
         ValueError: ``method`` is not one of ``METHODS``, or ``reference`` selects fewer than four rows.
@@ -247,5 +248,6 @@ def hit_calling(
     # that reads a response per row, a dose-response fit above all, has no spread to work with. The row's own
     # distance is already computed here.
     adata.obs[f"{key_added}_row_distance"] = to_control
+    adata.obs[f"{key_added}_reference_held_out"] = held_out
     get_logger().info("hit_calling(%s) called %d of %d groups", method, int(table["is_hit"].sum()), len(table))
     return None

tests/test_datasets_oasis.py

@@ -4,6 +4,7 @@
 import pytest
 
 import mantispy as mt
+from mantispy.ds._datasets import _DOSE_TOLERANCE
 
 
 @pytest.fixture(scope="module")
@@ -38,6 +39,34 @@ def test_the_dose_range_survives_the_join(oasis):
     assert np.isfinite(treated["Metadata_Concentration"].to_numpy(dtype=float)).all()
 
 
+@pytest.mark.network
+def test_one_concentration_written_twice_is_one_dose(oasis):
+    """Two plate maps write the dose to three decimals and the rest to four, so 3.704 and 3.7037 are one level."""
+    obs = oasis.obs
+    raw = obs["Metadata_Concentration"].to_numpy(dtype=float)
+    aligned = obs["Metadata_ConcentrationNominal"].to_numpy(dtype=float)
+    dosed = raw > 0
+
+    levels = np.unique(aligned[dosed])
+    gaps = np.diff(levels) / levels[:-1]
+    assert (gaps >= _DOSE_TOLERANCE).all(), "two levels within the tolerance are one dose written twice"
+    assert len(levels) < len(np.unique(raw[dosed])), "nothing collapsed"
+    # The dose a well was meant to get is one it was recorded at, never an average of two.
+    assert set(levels) <= set(np.unique(raw[dosed]))
+    assert np.max(np.abs(aligned[dosed] / raw[dosed] - 1)) < _DOSE_TOLERANCE, "no well moves past the tolerance"
+
+
+@pytest.mark.network
+def test_replicates_of_one_treatment_land_in_one_group(oasis):
+    """Metadata_Perturbation is named by the aligned dose, so the two spellings do not halve a group."""
+    treated = oasis.obs[~oasis.obs["Metadata_Control"]]
+    sizes = treated.groupby(["Metadata_CellLine", "Metadata_Perturbation"], observed=True).size()
+    # HepaRG is eight plates, so its treatments carry eight wells each; U2OS is two and genuinely thin.
+    heparg = sizes.loc["HepaRG"]
+    assert heparg.median() >= 6, "the eight-plate line should keep its replicates together"
+    assert int((heparg < 3).sum()) < len(heparg) // 5
+
+
 @pytest.mark.network
 def test_annotate_false_leaves_the_profiles_alone(oasis):
     raw = mt.ds.oasis_pilot(annotate=False)

tests/test_pl_hits.py

@@ -150,3 +150,12 @@ def test_an_inhibitory_curve_is_drawn_the_way_the_data_runs(inhibitor_adata):
     line = next(line for line in ax.get_lines() if line.get_label().startswith("EC50"))
     drawn = line.get_ydata()
     assert drawn[0] > drawn[-1], "the curve runs uphill while the data runs downhill"
+
+
+def test_one_stray_well_does_not_flatten_the_rest_onto_the_baseline(inhibitor_adata):
+    """A distance from the controls has a long right tail, so a linear axis hides the response."""
+    assert mt.pl.dose_response(inhibitor_adata, compound="cpd").get_yscale() == "log"
+
+    # A response that reaches zero has no log scale to be drawn on.
+    inhibitor_adata.obs.loc[inhibitor_adata.obs.index[0], "hits_row_distance"] = 0.0
+    assert mt.pl.dose_response(inhibitor_adata, compound="cpd").get_yscale() == "linear"

tests/test_tl_design.py

@@ -47,7 +47,7 @@ def test_cytotoxicity_flags_cell_loss_not_morphology(profiles):
     toxic = (profiles.obs["Metadata_Perturbation"] == "pert00").to_numpy()
     counts[toxic] = counts[toxic] // 5
     profiles.obs["Metadata_CellCount"] = counts
-    profiles.obs["hits_distance"] = np.where(toxic, 10.0, 1.0)
+    profiles.obs["hits_row_distance"] = np.where(toxic, 10.0, 1.0)
 
     mt.tl.cytotoxicity(profiles)
     table = profiles.uns["mantispy"]["cytotoxicity"].set_index("group")
@@ -63,22 +63,45 @@ def test_cell_loss_alone_is_not_suspect(profiles):
     quiet = (profiles.obs["Metadata_Perturbation"] == "pert01").to_numpy()
     counts[quiet] = counts[quiet] // 5
     profiles.obs["Metadata_CellCount"] = counts
-    profiles.obs["hits_distance"] = 1.0
+    profiles.obs["hits_row_distance"] = 1.0
 
     mt.tl.cytotoxicity(profiles)
     table = profiles.uns["mantispy"]["cytotoxicity"].set_index("group")
     assert table.loc["pert01", "viability"] < 0.5
     assert not bool(table.loc["pert01", "suspect"])
 
 
+def test_cytotoxicity_medians_the_rows_rather_than_a_group_statistic(profiles):
+    """Regression for #84.
+
+    ``hits_distance`` is one number repeated over a group's rows, so the median this function
+    documents returned the value it was handed, and no genuinely per-row response could be given.
+    """
+    toxic = (profiles.obs["Metadata_Perturbation"] == "pert00").to_numpy()
+    counts = profiles.obs["Metadata_CellCount"].to_numpy().copy()
+    counts[toxic] = counts[toxic] // 5
+    profiles.obs["Metadata_CellCount"] = counts
+    # A quiet group with one well far out, and the group statistic that would hide the difference.
+    rows = np.ones(profiles.n_obs)
+    rows[np.flatnonzero(toxic)[0]] = 100.0
+    profiles.obs["hits_row_distance"] = rows
+    profiles.obs["hits_distance"] = np.where(toxic, 100.0, 1.0)
+
+    mt.tl.cytotoxicity(profiles)
+    table = profiles.uns["mantispy"]["cytotoxicity"].set_index("group")
+    assert table.loc["pert00", "viability"] < 0.5
+    assert table.loc["pert00", "distance"] == 1.0, "one well far out does not move the group's median"
+    assert not bool(table.loc["pert00", "suspect"])
+
+
 def test_cytotoxicity_says_what_to_run_first(profiles):
     with pytest.raises(KeyError, match="mt.tl.hit_calling"):
         mt.tl.cytotoxicity(profiles, distance_key="not_computed")
 
 
 def test_a_missing_count_names_where_to_get_one(profiles):
     del profiles.obs["Metadata_CellCount"]
-    profiles.obs["hits_distance"] = 1.0
+    profiles.obs["hits_row_distance"] = 1.0
     with pytest.raises(KeyError, match="mt.tl.aggregate.*mt.ds.*count_key="):
         mt.tl.cytotoxicity(profiles)
 
@@ -88,7 +111,7 @@ def test_viability_compares_cells_per_field(profiles):
     halved = (profiles.obs["Metadata_Perturbation"] == "pert01").to_numpy()
     profiles.obs["Metadata_SiteCount"] = np.where(halved, 2.0, 4.0)
     profiles.obs["Metadata_CellCount"] = np.where(halved, 10.0, 20.0)
-    profiles.obs["hits_distance"] = 1.0
+    profiles.obs["hits_row_distance"] = 1.0
 
     mt.tl.cytotoxicity(profiles)
     assert profiles.uns["mantispy"]["cytotoxicity"].set_index("group").loc["pert01", "viability"] == 1.0