"""Plots for mechanism retrieval and feature-set enrichment."""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
import pandas as pd
from mantispy._core.frames import as_frame
from mantispy.pl._common import axes as _axes
from mantispy.pl._common import maybe_interactive as _maybe_interactive
from mantispy.pl._common import returned as _returned
from mantispy.pl._common import table as _table
if TYPE_CHECKING:
from collections.abc import Sequence
from anndata import AnnData
from matplotlib.axes import Axes
def _heatmap(
ax: Axes,
values: np.ndarray,
rows: Sequence[str],
columns: Sequence[str],
cmap: str,
label: str,
fmt: str | None = None,
) -> Axes:
image = ax.imshow(values, aspect="auto", cmap=cmap)
ax.set_xticks(np.arange(len(columns)))
ax.set_xticklabels(columns, rotation=90, fontsize=6)
ax.set_yticks(np.arange(len(rows)))
ax.set_yticklabels(rows, fontsize=6)
if fmt is not None and values.size <= 400:
for i in range(values.shape[0]):
for j in range(values.shape[1]):
ax.text(j, i, format(values[i, j], fmt), ha="center", va="center", fontsize=5)
ax.figure.colorbar(image, ax=ax, label=label)
return ax
[docs]
def moa_confusion(adata: AnnData, key: str = "moa", normalize: bool = True, ax: Axes | None = None) -> Axes | None:
"""The confusion matrix of :func:`~mantispy.tl.nn_moa_classify`, as a heatmap.
With row normalization the diagonal is per-mechanism recall, and an off-diagonal block marks a pair of mechanisms the morphology does not separate.
Such pairs usually have similar phenotypes.
Args:
adata: Object holding the confusion table :func:`~mantispy.tl.nn_moa_classify` wrote.
key: Name that run's outputs were stored under, whose confusion table is ``key + "_confusion"``.
normalize: Divide each row by its total, which turns the counts into per-mechanism recall.
ax: Axes to draw on, or ``None`` for a new figure.
Returns:
The axes when the caller passed ``ax``, else ``None`` because the plot then owns the figure it created.
When returned they hold true against predicted mechanisms with the values printed when there are at most 400 cells, and the run's scheme and accuracy in the title.
Raises:
KeyError: ``uns["mantispy"]`` holds no ``key + "_confusion"`` table.
"""
table = _table(adata, f"{key}_confusion", "mt.tl.nn_moa_classify")
matrix = table.pivot_table(index="true", columns="predicted", values="count", aggfunc="sum", fill_value=0)
labels = sorted(set(matrix.index) | set(matrix.columns))
counts = matrix.reindex(index=labels, columns=labels, fill_value=0).to_numpy(dtype=float)
if normalize:
totals = counts.sum(axis=1, keepdims=True)
counts = np.divide(counts, totals, out=np.zeros_like(counts), where=totals > 0)
ax = _axes(ax, (0.45 * len(labels) + 3, 0.4 * len(labels) + 2.5))
_heatmap(ax, counts, labels, labels, "Blues", "fraction" if normalize else "count", ".2f" if normalize else ".0f")
ax.set_xlabel("predicted")
ax.set_ylabel("true")
summary = adata.uns.get("mantispy", {}).get(key, {})
title = f"{summary.get('scheme', '')} accuracy {float(summary.get('accuracy', float('nan'))):.1%}"
ax.set_title(title, fontsize=9)
_maybe_interactive(
"heatmap",
ax=ax,
matrix=counts,
rows=labels,
columns=labels,
value_label="fraction" if normalize else "count",
title=title,
)
return _returned(ax)
[docs]
def moa_enrichment(
adata: AnnData, group: str, key: str = "moa_enrichment", top: int = 10, ax: Axes | None = None
) -> Axes | None:
"""Which mechanisms one profile's neighborhood is enriched for.
Args:
adata: Object holding the table :func:`~mantispy.tl.moa_enrichment` wrote.
group: Which group of that table to draw.
key: Name of that table in ``uns["mantispy"]``.
top: How many mechanisms to draw, taken by p-value.
ax: Axes to draw on, or ``None`` for a new figure.
Returns:
The axes when the caller passed ``ax``, else ``None`` because the plot then owns the figure it created.
When returned they hold one bar of ``-log10`` q per mechanism, labeled by how many of the neighbors carried it, and a reference line at ``q = 0.05``.
Raises:
KeyError: There is no such table, or it holds no such group.
"""
table = _table(adata, key, "mt.tl.moa_enrichment")
selected = table[table["group"].astype(str) == str(group)]
if selected.empty:
raise KeyError(f"no group {group!r} in uns['mantispy'][{key!r}]")
best = selected.nsmallest(top, "pvalue")[::-1]
labels = [f"{moa} ({n})" for moa, n in zip(best["moa"], best["n_neighbours"], strict=True)]
neg_log_q = -np.log10(np.clip(best["qvalue"].to_numpy(dtype=float), 1e-12, None))
ax = _axes(ax, (5.5, 0.3 * len(best) + 1.5))
ax.barh(np.arange(len(best)), neg_log_q)
ax.set_yticks(np.arange(len(best)))
ax.set_yticklabels(labels, fontsize=6)
ax.axvline(-np.log10(0.05), color="grey", ls="--", lw=1)
ax.set_xlabel("-log10 q")
ax.set_title(str(group), fontsize=9)
tidy = pd.DataFrame({"mechanism": labels, "-log10 q": neg_log_q})
_maybe_interactive("barh", ax=ax, data=tidy, x="-log10 q", y="mechanism", title=str(group))
return _returned(ax)
[docs]
def distance_heatmap(
adata: AnnData, key: str = "edistance", groupby: str | None = "Metadata_MOA", ax: Axes | None = None
) -> Axes | None:
"""The group-by-group distance matrix, ordered so related groups sit together.
Args:
adata: Object holding the pairwise matrix :func:`~mantispy.tl.edistance` wrote with ``reference=None``.
key: Name that run's outputs were stored under, whose matrix is ``key + "_pairwise"``.
groupby: ``obs`` column to order the groups by, or ``None`` to keep the matrix's own order.
Ordering also needs an ``obs`` column naming the groups of the matrix, and without one the matrix is drawn unordered rather than refused.
ax: Axes to draw on, or ``None`` for a new figure.
Returns:
The axes when the caller passed ``ax``, else ``None`` because the plot then owns the figure it created.
When returned they hold the distance matrix with a white line at each ``groupby`` boundary.
Raises:
KeyError: ``uns["mantispy"]`` holds no ``key + "_pairwise"`` matrix.
"""
matrix = _table(adata, f"{key}_pairwise", "mt.tl.edistance(reference=None)")
labels = list(matrix.columns)
matrix.index = pd.Index(labels)
obs = as_frame(adata.obs)
annotation = None
if groupby and groupby in obs:
naming = next((column for column in obs.columns if set(obs[column].astype(str)) >= set(labels)), None)
if naming is not None:
lookup = obs.groupby(obs[naming].astype(str), observed=True)[groupby].first().astype(str)
annotation = [str(lookup.get(label, "")) for label in labels]
order = np.argsort(annotation, kind="stable")
labels = [labels[index] for index in order]
annotation = [annotation[index] for index in order]
matrix = matrix.loc[labels, labels]
ax = _axes(ax, (0.25 * len(labels) + 3, 0.22 * len(labels) + 2.5))
_heatmap(ax, matrix.to_numpy(dtype=float), labels, labels, "magma", "energy distance")
if annotation:
for position in np.flatnonzero(np.asarray(annotation[1:]) != np.asarray(annotation[:-1])) + 1:
ax.axhline(position - 0.5, color="white", lw=0.8)
ax.axvline(position - 0.5, color="white", lw=0.8)
_maybe_interactive(
"heatmap",
ax=ax,
matrix=matrix.to_numpy(dtype=float),
rows=labels,
columns=labels,
value_label="energy distance",
title="distance heatmap",
)
return _returned(ax)
[docs]
def sets_heatmap(
adata: AnnData, groupby: str, score_key: str = "score_ulm", top: int = 30, ax: Axes | None = None
) -> Axes | None:
"""Mean enrichment score per group per feature set.
Args:
adata: Object :func:`~mantispy.tl.enrich` has scored, holding the per-row scores in ``obsm``.
groupby: ``obs`` column whose groups become the rows.
score_key: ``obsm`` key holding those scores, named after the method that wrote them.
top: How many feature sets to draw, taken by their largest absolute mean score.
ax: Axes to draw on, or ``None`` for a new figure.
Returns:
The axes when the caller passed ``ax``, else ``None`` because the plot then owns the figure it created.
When returned they hold groups against feature sets, the sets it kept sorted by name, on a diverging scale.
Raises:
KeyError: ``obsm`` holds nothing under ``score_key``.
"""
if score_key not in adata.obsm:
raise KeyError(f"obsm has no {score_key!r}; run mt.tl.enrich first")
stored = adata.obsm[score_key]
names = list(stored.columns) if hasattr(stored, "columns") else [str(i) for i in range(np.shape(stored)[1])]
scores = np.asarray(stored, dtype=float)
groups = as_frame(adata.obs)[groupby].astype(str).to_numpy()
labels = list(dict.fromkeys(groups))
means = np.stack([np.nanmean(scores[groups == group], axis=0) for group in labels])
keep = np.argsort(-np.nanmax(np.abs(means), axis=0))[:top]
keep = keep[np.argsort([names[index] for index in keep])]
kept_names = [names[index] for index in keep]
ax = _axes(ax, (0.3 * len(keep) + 3, 0.28 * len(labels) + 2))
_heatmap(ax, means[:, keep], labels, kept_names, "coolwarm", "mean score")
_maybe_interactive(
"heatmap",
ax=ax,
matrix=means[:, keep],
rows=labels,
columns=kept_names,
value_label="mean score",
title="feature-set scores",
)
return _returned(ax)
[docs]
def pathway_coherence(
adata: AnnData, key: str = "pathway_coherence", top: int = 15, ax: Axes | None = None
) -> Axes | None:
"""Coherence per gene set, the significant ones marked.
Sets are ordered by coherence, as in the table.
Under a permutation null every coherent set ties at the p-value floor, so the q-value marks significance and coherence ranks the sets.
Args:
adata: Object holding the table :func:`~mantispy.tl.pathway_coherence` wrote.
key: Name of that table in ``uns["mantispy"]``.
top: How many sets to draw, taken by coherence.
ax: Axes to draw on, or ``None`` for a new figure.
Returns:
The axes when the caller passed ``ax``, else ``None`` because the plot then owns the figure it created.
When returned they hold one bar per set labeled by how many of its genes were in the screen, colored by whether its q-value is below 0.05.
Raises:
KeyError: ``uns["mantispy"]`` holds no table under ``key``.
ValueError: That table is empty, which is what happens when no set had enough of its genes in the screen.
"""
import matplotlib.pyplot as plt
table = _table(adata, key, "mt.tl.pathway_coherence", "no set had enough of its genes in the screen")
best = table.nlargest(min(top, len(table)), "coherence")[::-1]
significant = best["qvalue"].to_numpy(dtype=float) < 0.05
labels = [f"{name} ({n})" for name, n in zip(best["set"], best["n_genes"], strict=True)]
ax = _axes(ax, (6, 0.3 * len(best) + 1.5))
ax.barh(np.arange(len(best)), best["coherence"], color=np.where(significant, "crimson", "lightgrey"))
ax.set_yticks(np.arange(len(best)))
ax.set_yticklabels(labels, fontsize=6)
ax.set_xlabel("mean similarity among the set's genes")
ax.legend(
handles=[
plt.Line2D([], [], color="crimson", lw=6, label="q < 0.05"),
plt.Line2D([], [], color="lightgrey", lw=6, label="not significant"),
],
fontsize=6,
loc="lower right",
)
tidy = pd.DataFrame(
{
"set": labels,
"coherence": best["coherence"].to_numpy(dtype=float),
"significant": np.where(significant, "q < 0.05", "not significant"),
}
)
_maybe_interactive("barh", ax=ax, data=tidy, x="coherence", y="set", color="significant", title="pathway coherence")
return _returned(ax)
