"""Feature-set enrichment over the parsed feature annotation.
``var`` records the object, measurement family and channel of every feature, which defines a
set-membership table. Scoring those sets reports which kinds of measurement changed, such as
the mitochondrial texture features, instead of a list of individual columns.
The sets are passed to decoupler, so the same call works for sets built from the feature
names and for sets from prior knowledge.
"""
from __future__ import annotations
import warnings
from collections.abc import Sequence
from typing import Any
import anndata as ad
import numpy as np
import pandas as pd
from anndata import AnnData
from mantispy._core._reduce import get_matrix
from mantispy._core._utils import as_frame, get_logger, inplace_or_copy
METHODS = ("ulm", "mlm", "ora")
#: Named shorthands for composite groupings.
COMPOSITES = {"group_by_channel": ["feature_group", "channel"]}
[docs]
def feature_sets(adata: AnnData, by: str | Sequence[str] = "feature_group") -> pd.DataFrame:
"""Build a decoupler network from the parsed feature annotation.
Args:
adata: Object whose ``var`` carries the parser's columns.
by: A ``var`` column, several of them (joined with ``|``), or ``"group_by_channel"``
for the ``feature_group``-and-``channel`` combination.
Returns:
A frame with ``source``, ``target`` and ``weight``, ready for :func:`enrich` or for
decoupler directly. Features with a missing annotation in any of the columns are left
out.
"""
requested = COMPOSITES.get(by, by) if isinstance(by, str) else list(by)
columns = [requested] if isinstance(requested, str) else list(requested)
var = as_frame(adata.var)
missing = [column for column in columns if column not in var]
if missing:
raise KeyError(f"var has no column(s) {missing}; available: {sorted(var.columns)}")
known = var[columns].notna().all(axis=1).to_numpy()
labels = var.loc[known, columns].astype(str).agg("|".join, axis=1)
usable = ~labels.str.contains("nan").to_numpy()
return pd.DataFrame(
{
"source": labels[usable].to_numpy(),
"target": var.index[known][usable].to_numpy(),
"weight": 1.0,
}
)
[docs]
@inplace_or_copy()
def enrich(
adata: AnnData,
net: pd.DataFrame | None = None,
by: str | Sequence[str] = "feature_group",
method: str = "ulm",
top_fraction: float = 0.05,
copy: bool = False,
**decoupler_kwargs: Any,
) -> AnnData | None:
"""Score every profile against every feature set.
Args:
adata: Profiles to score. Normalize first, since the methods use the values as given.
net: A decoupler network with ``source``, ``target`` and ``weight``, for example
prior-knowledge sets. Built from ``by`` when omitted.
by: Passed to :func:`feature_sets` when ``net`` is not given.
method: ``"ulm"`` fits a linear model per set and is the usual choice; ``"mlm"`` fits all
sets jointly, which handles overlapping sets; ``"ora"`` is an over-representation
test on the extremes.
top_fraction: Fraction of features, ranked by value, that ``method="ora"`` counts as
extreme. The default 0.05 tests the top twentieth against the rest. Ignored when
``n_up`` is passed, and by ``"ulm"`` and ``"mlm"``, which use every feature.
copy: Return a modified copy instead of mutating in place.
decoupler_kwargs: Passed through to decoupler, e.g. ``tmin`` for the smallest usable set.
Returns:
``None``, or the modified copy. decoupler writes ``obsm["score_<method>"]`` and
``obsm["padj_<method>"]``, both frames indexed by set name.
"""
import decoupler as dc
if method not in METHODS:
raise ValueError(f"method must be one of {METHODS}, got {method!r}")
network = feature_sets(adata, by) if net is None else net
if not len(network):
raise ValueError(f"no feature sets built from by={by!r}: every feature's annotation is missing")
if method == "ora":
if not 0.0 < top_fraction < 1.0:
raise ValueError(f"top_fraction must be in (0, 1), got {top_fraction!r}")
# decoupler keeps features ranked above n_up, so the top k of n needs n_up = n - k.
# Its own default selects the bottom 95%.
decoupler_kwargs.setdefault("n_up", int(round(adata.n_vars * (1.0 - top_fraction))))
getattr(dc.mt, method)(adata, network, **decoupler_kwargs)
get_logger().info("enrich(%s) scored %d set(s)", method, network["source"].nunique())
return None
[docs]
@inplace_or_copy()
def rank_features(
adata: AnnData,
groupby: str,
method: str = "wilcoxon",
key_added: str = "rank_features",
copy: bool = False,
) -> AnnData | None:
"""Rank features by how well they separate each group, with the annotation attached.
Wraps :func:`scanpy.tl.rank_genes_groups` and joins the parsed ``var`` annotation onto the
result, so each ranked feature carries its object, feature group and channel.
Args:
adata: Object to rank.
groupby: ``obs`` column defining the groups.
method: Passed to scanpy: ``"wilcoxon"``, ``"t-test"``, ``"logreg"``.
key_added: Name for the output table.
copy: Return a modified copy instead of mutating in place.
Returns:
``None``, or the modified copy. Writes ``uns["mantispy"][key_added]`` with ``group``,
``feature``, ``score``, ``pvalue``, ``qvalue`` and the ``object``, ``feature_group``
and ``channel`` the feature belongs to.
Notes:
scanpy's ``logfoldchanges`` column is dropped. Normalized morphology features are
signed, so the ratio is often negative and its log is NaN or meaningless. Rank by
``score``.
"""
import scanpy as sc
var = as_frame(adata.var)
# A shallow object: X is shared, so ranking a large matrix does not double memory.
scratch = ad.AnnData(X=get_matrix(adata), obs=as_frame(adata.obs)[[groupby]].astype("category"), var=var[[]])
with warnings.catch_warnings():
# scanpy always computes log fold changes, which warn on the negative ratios of signed
# features. The column is dropped below.
warnings.filterwarnings("ignore", "invalid value encountered in log2", RuntimeWarning)
sc.tl.rank_genes_groups(scratch, groupby=groupby, method=method)
table = sc.get.rank_genes_groups_df(scratch, group=None).rename(
columns={"names": "feature", "scores": "score", "pvals": "pvalue", "pvals_adj": "qvalue"}
)
table = table.drop(columns=[column for column in ("logfoldchanges",) if column in table])
annotation = [column for column in ("object", "feature_group", "channel") if column in var]
table = table.merge(var[annotation].astype(str), left_on="feature", right_index=True, how="left")
adata.uns.setdefault("mantispy", {})[key_added] = table
return None
[docs]
@inplace_or_copy()
def rank_sets(
adata: AnnData,
groupby: str,
score_key: str = "score_ulm",
key_added: str = "rank_sets",
copy: bool = False,
) -> AnnData | None:
"""Rank feature sets by how far each group's score sits from the rest.
Args:
adata: Object already scored by :func:`enrich`.
groupby: ``obs`` column defining the groups.
score_key: The ``obsm`` key written by :func:`enrich`.
key_added: Name for the output table.
copy: Return a modified copy instead of mutating in place.
Returns:
``None``, or the modified copy. Writes ``uns["mantispy"][key_added]`` with ``group``,
``set`` and ``score``, the group's mean enrichment score minus the mean over all other rows.
"""
if score_key not in adata.obsm:
raise KeyError(f"obsm has no {score_key!r}; run mt.tl.enrich first, which writes it")
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()
records = []
for group in pd.unique(groups):
inside = groups == group
if inside.all():
raise ValueError(f"{groupby!r} has a single group, so there is nothing to rank it against")
difference = np.nanmean(scores[inside], axis=0) - np.nanmean(scores[~inside], axis=0)
records.append(pd.DataFrame({"group": str(group), "set": names, "score": difference}))
adata.uns.setdefault("mantispy", {})[key_added] = pd.concat(records, ignore_index=True)
return None