"""Plots for hit calling, effect sizes and dose response."""
from __future__ import annotations
from typing import TYPE_CHECKING
import matplotlib.pyplot as plt
import numpy as np
from mantispy._core._utils import as_frame
from mantispy.pl._common import axes as _axes
from mantispy.pl._common import table as _table
if TYPE_CHECKING:
from anndata import AnnData
#: Values below this are clipped so they stay on the plot.
_FLOOR = 1e-12
def _significance(values) -> np.ndarray:
return -np.log10(np.clip(np.asarray(values, dtype=float), _FLOOR, None))
def _threshold(adata: AnnData, function: str, default: float = 0.05) -> float:
"""The threshold the run used, so the drawn line matches the table."""
recorded = adata.uns.get("mantispy", {}).get("params", {}).get(function, {})
return float(recorded.get("threshold", default))
[docs]
def hits(adata: AnnData, key: str = "hits", label_top: int = 10, ax: plt.Axes | None = None):
"""Distance from the controls against significance, with the most distant groups labeled.
A point in the upper right moved far from the controls and is significant under the
permutation null. The dashed line is the q-value threshold the run used.
"""
table = _table(adata, key, "mt.tl.hit_calling")
ax = _axes(ax, (5.5, 4.5))
significance = _significance(table["qvalue"])
called = table["is_hit"].to_numpy(dtype=bool)
ax.scatter(table["distance"][~called], significance[~called], s=14, color="lightgrey", label="not called")
ax.scatter(table["distance"][called], significance[called], s=14, color="crimson", label="hit")
threshold = _threshold(adata, key)
ax.axhline(-np.log10(threshold), color="grey", ls="--", lw=1, label=f"q = {threshold}")
for _, row in table.nlargest(label_top, "distance").iterrows():
ax.annotate(str(row["group"]), (row["distance"], -np.log10(max(float(row["qvalue"]), _FLOOR))), fontsize=6)
ax.set_xlabel("distance from the controls")
ax.set_ylabel("-log10 q")
ax.legend(fontsize=7)
return ax
def _effects(adata: AnnData, group: str, key: str):
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]}")
var = as_frame(adata.var)
families = var["feature_group"].astype(str) if "feature_group" in var else None
return selected, families
def _family_colours(families, names) -> tuple[list, dict]:
"""One color per feature family, assigned in sorted order of the families present."""
if families is None:
return ["tab:blue"] * len(names), {}
labels = [str(families.get(name, "unknown")) for name in names]
palette = {label: plt.get_cmap("tab20")(index % 20) for index, label in enumerate(sorted(set(labels)))}
return [palette[label] for label in labels], palette
[docs]
def effect_sizes(adata: AnnData, group: str, key: str = "effect", top: int = 30, ax: plt.Axes | None = None):
"""The largest effects for one group, colored by feature family."""
selected, families = _effects(adata, group, key)
strongest = selected.reindex(selected["effect"].abs().sort_values(ascending=False).index).head(top)[::-1]
ax = _axes(ax, (6, 0.22 * len(strongest) + 1.5))
colours, palette = _family_colours(families, strongest["feature"])
ax.barh(np.arange(len(strongest)), strongest["effect"], color=colours)
ax.set_yticks(np.arange(len(strongest)))
ax.set_yticklabels(strongest["feature"], fontsize=6)
ax.axvline(0, color="black", lw=0.8)
ax.set_xlabel("effect size")
ax.set_title(str(group), fontsize=9)
if palette:
ax.legend(
handles=[plt.Line2D([], [], color=colour, lw=6, label=label) for label, colour in palette.items()],
fontsize=6,
loc="lower right",
)
return ax
[docs]
def feature_volcano(adata: AnnData, group: str, key: str = "effect", label_top: int = 8, ax: plt.Axes | None = None):
"""Effect against significance, per feature, for one group."""
selected, families = _effects(adata, group, key)
ax = _axes(ax, (5.5, 4.5))
colours, palette = _family_colours(families, selected["feature"])
significance = _significance(selected["qvalue"])
ax.scatter(selected["effect"], significance, s=12, color=colours)
for _, row in (
selected.reindex(selected["effect"].abs().sort_values(ascending=False).index).head(label_top).iterrows()
):
ax.annotate(str(row["feature"]), (row["effect"], -np.log10(max(float(row["qvalue"]), _FLOOR))), fontsize=5)
ax.axhline(-np.log10(0.05), color="grey", ls="--", lw=1)
ax.axvline(0, color="black", lw=0.6)
ax.set_xlabel("effect size")
ax.set_ylabel("-log10 q")
ax.set_title(str(group), fontsize=9)
if palette:
ax.legend(
handles=[plt.Line2D([], [], color=colour, lw=6, label=label) for label, colour in palette.items()],
fontsize=5,
loc="upper left",
)
return ax
[docs]
def dose_response(
adata: AnnData,
compound: str,
key: str = "dose_response",
compound_key: str = "Metadata_Compound",
dose_key: str = "Metadata_Concentration",
response: str = "hits_distance",
ax: plt.Axes | None = None,
):
"""One compound's response against dose, with the fitted curve when there is one."""
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}]")
obs = as_frame(adata.obs)
selected = obs[obs[compound_key].astype(str) == str(compound)]
doses = selected[dose_key].to_numpy(dtype=float)
values = selected[response].to_numpy(dtype=float)
usable = np.isfinite(doses) & np.isfinite(values) & (doses > 0)
ax = _axes(ax, (5, 4))
ax.scatter(doses[usable], values[usable], s=18, label="wells")
ax.set_xscale("log")
fitted = row.iloc[0]
if bool(fitted["fit_ok"]):
grid = np.log10(np.geomspace(doses[usable].min(), doses[usable].max(), 100))
curve = four_parameter_logistic(
grid,
float(fitted["bottom"]),
float(fitted["top"]),
np.log10(float(fitted["ec50"])),
float(fitted["hill_slope"]),
)
ax.plot(10.0**grid, curve, color="crimson", lw=1.5, label=f"EC50 = {float(fitted['ec50']):.3g}")
ax.set_xlabel(dose_key.replace("Metadata_", ""))
ax.set_ylabel(response)
ax.set_title(f"{compound} (spearman {float(fitted['spearman']):.2f})", fontsize=9)
ax.legend(fontsize=7)
return ax