from __future__ import annotations
import numpy as np
import pandas as pd
from anndata import AnnData
from mantispy._core._reduce import group_codes, group_offsets, representation
from mantispy._core._stats import 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.mutation import inplace_or_copy
WEIGHTS = ("activity", "equal")
def _effects(
values: np.ndarray, units: np.ndarray, codes: np.ndarray, keys: pd.Index, is_control: np.ndarray
) -> tuple[dict[str, np.ndarray], dict[str, np.ndarray], list[str]]:
"""Per setting, a groups-by-features effect matrix and each group's effect magnitude.
Groups missing from a setting get NaN rows.
"""
effect: dict[str, np.ndarray] = {}
activity: dict[str, np.ndarray] = {}
usable = []
order, offsets = group_offsets(codes, len(keys))
for unit in sorted(set(units.tolist())):
here = units == unit
if (here & is_control).sum() < 2:
continue
usable.append(unit)
centre = np.nanmedian(values[here & is_control], axis=0)
block = np.full((len(keys), values.shape[1]), np.nan)
for index in range(len(keys)):
group = order[offsets[index] : offsets[index + 1]]
rows = group[here[group] & ~is_control[group]]
if rows.size:
block[index] = np.nanmedian(values[rows], axis=0) - centre
effect[unit] = block
activity[unit] = np.linalg.norm(np.nan_to_num(block), axis=1)
activity[unit][~np.isfinite(block).any(axis=1)] = np.nan
return effect, activity, usable
def _quiet(block: np.ndarray) -> np.ndarray:
"""``block`` with all-NaN rows replaced by zeros, so ``nanmean`` does not warn."""
empty = ~np.isfinite(block).any(axis=1)
if not empty.any():
return block
out = block.copy()
out[empty] = 0.0
return out
def _standardize(block: np.ndarray) -> np.ndarray:
"""Rows centered and scaled to unit norm, so the dot product of two rows is a correlation.
Missing values are filled with zero after centering, as :func:`~mantispy.tl.hit_calling` does.
Without missing values the result is the Pearson correlation.
"""
with np.errstate(invalid="ignore"):
centre = np.where(
np.isfinite(block).any(axis=1, keepdims=True), np.nanmean(_quiet(block), axis=1, keepdims=True), 0.0
)
centred = np.nan_to_num(block - centre)
norm = np.sqrt(np.einsum("ij,ij->i", centred, centred))[:, None]
return np.divide(centred, norm, out=np.zeros_like(centred), where=norm > 0)
def _agreement_matrix(left: np.ndarray, right: np.ndarray) -> np.ndarray:
"""Correlation of every group's effect in one setting with every group's effect in another.
Entry ``[i, j]`` correlates group ``i``'s effect at the left setting with group ``j``'s at the right.
The diagonal holds each group's statistic and the off-diagonal entries form the null.
"""
return _standardize(left) @ _standardize(right).T
def _level_pairs(frame: pd.DataFrame, levels: list[str], unit_key: str) -> dict[str, list[tuple[str, str]]]:
"""Unit pairs assigned to each level, keyed in the order of ``levels``.
Each pair of units belongs only to the coarsest level at which the two differ.
"""
# drop=False: with a single level, the unit column is also the level column.
lookup = frame.drop_duplicates(unit_key).set_index(unit_key, drop=False)
units = list(lookup.index)
buckets: dict[str, list[tuple[str, str]]] = {level: [] for level in levels}
for i, one in enumerate(units):
for other in units[:i]:
for level in levels: # coarsest first, so the first difference wins
if lookup.loc[one, level] != lookup.loc[other, level]:
buckets[level].append((one, other))
break
return buckets
[docs]
@inplace_or_copy(expects=("well", "perturbation"))
def transport(
adata: AnnData,
by: str | list[str] = "Metadata_Plate",
groupby: str = "Metadata_Perturbation",
reference: str | None = "negcon",
use_rep: str | None = None,
weight: str = "activity",
min_shared: int = 3,
threshold: float = 0.05,
key_added: str = "transport",
copy: bool = False,
) -> AnnData | None:
"""Test whether each perturbation's effect reproduces across settings.
Args:
adata: Well-level profiles.
Each setting needs its own reference wells, since effects are measured against them; settings with fewer than two are left out.
by: ``obs`` column defining the setting, or a list of columns from coarsest to finest.
``["Metadata_Source", "Metadata_Plate"]`` reports agreement between plates of one source separately from agreement between sources, and the difference shows what a change of laboratory costs beyond a change of plate.
The finest level defines the units that are compared.
groupby: The perturbation column.
reference: Which rows are the negative controls, per setting.
use_rep: Score ``obsm[use_rep]`` instead of ``X``.
weight: ``"activity"`` weights each comparison by the smaller of the two effect magnitudes, since the correlation of an inactive perturbation is noise.
``"equal"`` weights all comparisons the same.
min_shared: Minimum number of shared perturbations for a pair of units to be compared.
threshold: q-value cutoff for ``transports``.
The null uses every mismatched pair of perturbations, so there is no null size or seed to set.
key_added: Name for the outputs.
copy: Return a modified copy instead of mutating in place.
Returns:
``None``, or the modified copy.
Writes ``uns["mantispy"][key_added]`` with one row per ``group`` and ``level`` and the columns ``n_pairs``, ``agreement``, ``pvalue``, ``qvalue`` and ``transports``.
Writes ``uns["mantispy"][key_added + "_units"]``, the units-by-units agreement matrix drawn by :func:`~mantispy.pl.setting_agreement`, and joins the finest level's agreement back onto the rows as ``obs[key_added + "_agreement"]``.
Raises:
ValueError: ``weight`` is not one of ``WEIGHTS``, the unit column has a single level, or fewer than two settings have at least two reference rows.
KeyError: ``obs`` has no column named by ``by``.
Notes:
The null pairs a perturbation at one unit with a different perturbation at another, so a screen in which all effects look alike does not count as reproducing.
This matters on screens with few mechanisms.
The mismatched null centers at +0.26 on BBBC021 and at +0.03 on JUMP.
The null is computed in full rather than sampled.
Standardizing the effect vectors turns all group-against-group correlations into one matrix product, so 301 groups over 48 unit pairs give 90,300 null values in about a second.
A sampled null would floor the p-values at ``1/(n + 1)``, too coarse for Benjamini-Hochberg over every ``(group, level)`` row, and the number of groups called would depend on the number of draws.
This is an observational measure.
It shows whether an effect reproduced at another site, not what the effect would have been there.
"""
if weight not in WEIGHTS:
raise ValueError(f"weight must be one of {WEIGHTS}, got {weight!r}")
levels = [by] if isinstance(by, str) else list(by)
obs = as_frame(adata.obs)
if missing := [level for level in levels if level not in obs.columns]:
raise KeyError(f"obs has no column(s) {missing}; by= names the setting(s) to compare across")
unit_key = levels[-1]
units = obs[unit_key].astype(str).to_numpy()
if len(set(units.tolist())) < 2:
raise ValueError(
f"transport compares settings, and obs[{unit_key!r}] has one level. Pass by= a column that "
"varies, such as the plate, batch or laboratory."
)
values = representation(adata, use_rep)
is_control = reference_mask(adata, reference)
codes, keys = group_codes(adata, groupby)
effect, activity, usable = _effects(values, units, codes, keys, is_control)
if len(usable) < 2:
raise ValueError(
f"only {len(usable)} level(s) of {unit_key!r} have at least two reference rows. Effects are "
"measured against each setting's own controls, so at least two settings need them."
)
frame = obs.loc[:, levels].astype(str).assign(**{unit_key: units})
buckets = _level_pairs(frame, levels, unit_key)
buckets = {
level: [pair for pair in pairs if pair[0] in usable and pair[1] in usable] for level, pairs in buckets.items()
}
records: list[dict[str, object]] = []
matrix_rows: list[dict[str, object]] = []
for level, pairs in buckets.items():
if not pairs:
get_logger().info("transport has no unit pair separating at %r; skipping that level", level)
continue
total = np.zeros((len(keys), len(keys)))
weights = np.zeros_like(total)
counted = np.zeros(len(keys), dtype=int)
for one, other in pairs:
left, right = effect[one], effect[other]
here = np.isfinite(left).any(axis=1)
there = np.isfinite(right).any(axis=1)
if int((here & there).sum()) < min_shared:
continue
block = _agreement_matrix(left, right)
if weight == "activity":
pair_weight = np.minimum(activity[one][:, None], activity[other][None, :])
else:
pair_weight = np.ones_like(block)
pair_weight = np.where(here[:, None] & there[None, :], np.nan_to_num(pair_weight), 0.0)
total += np.nan_to_num(block) * pair_weight
weights += pair_weight
counted += (here & there).astype(int)
shared = np.flatnonzero(here & there)
matrix_rows.extend(
{
"left": one,
"right": other,
"group": str(keys[int(index)]),
"agreement": float(block[index, index]),
"weight": float(pair_weight[index, index]),
}
for index in shared
)
with np.errstate(invalid="ignore", divide="ignore"):
combined = np.where(weights > 0, total / np.where(weights > 0, weights, 1.0), np.nan)
observed = np.diag(combined)
off_diagonal = combined[~np.eye(len(keys), dtype=bool)]
null = off_diagonal[np.isfinite(off_diagonal)]
for index in range(len(keys)):
if not counted[index] or not np.isfinite(observed[index]):
continue
pvalue = float((np.sum(null >= observed[index]) + 1) / (null.size + 1)) if null.size else np.nan
records.append(
{
"group": str(keys[int(index)]),
"level": level,
"n_pairs": int(counted[index]),
"agreement": float(observed[index]),
"pvalue": pvalue,
}
)
table = pd.DataFrame(records, columns=["group", "level", "n_pairs", "agreement", "pvalue"])
if len(table):
table["qvalue"] = benjamini_hochberg(table["pvalue"].to_numpy())
table["transports"] = table["qvalue"] < threshold
else:
table["qvalue"] = []
table["transports"] = []
adata.uns.setdefault("mantispy", {})[key_added] = table
pairs_frame = pd.DataFrame(matrix_rows, columns=["left", "right", "group", "agreement", "weight"])
adata.uns["mantispy"][f"{key_added}_units"] = _unit_matrix(pairs_frame, usable, weight)
finest = table[table["level"] == unit_key] if unit_key in set(table["level"]) else table
lookup = finest.drop_duplicates("group").set_index("group")["agreement"]
labels = obs[groupby].astype(str)
adata.obs[f"{key_added}_agreement"] = lookup.reindex(labels).to_numpy()
for level in buckets:
called = table[(table["level"] == level) & table["transports"]]
get_logger().info(
"transport(%s): %d of %d group(s) reproduce", level, len(called), int((table["level"] == level).sum())
)
return None
def _unit_matrix(pairs: pd.DataFrame, units: list[str], weight: str) -> pd.DataFrame:
"""Units by units: how well two settings agree over the perturbations they share."""
values = np.full((len(units), len(units)), np.nan)
np.fill_diagonal(values, 1.0)
matrix = pd.DataFrame(values, index=units, columns=units, dtype=float)
if not len(pairs):
return matrix
for (one, other), block in pairs.groupby(["left", "right"], observed=True):
weights = block["weight"].to_numpy() if weight == "activity" else np.ones(len(block))
if weights.sum() <= 0:
weights = np.ones(len(block))
value = float(np.average(block["agreement"].to_numpy(), weights=weights))
matrix.loc[one, other] = matrix.loc[other, one] = value
return matrix
