"""Feature selection, with pycytominer's semantics.
Every operation returns a boolean keep mask over ``var``.
"""
from __future__ import annotations
import warnings
from collections.abc import Callable, Mapping, Sequence
from typing import Any
import numpy as np
import pandas as pd
from anndata import AnnData
from mantispy._core._corr import _blockwise, correlated_pairs, rank_revealing_subset
from mantispy._core._numba import column_max_robust_z
from mantispy._core._reduce import get_matrix, group_codes, group_offsets
from mantispy._core._stats import MAD_TO_SIGMA, nanvar
from mantispy._core.features import feature_kind_of, load_blocklist, var_blocklist_hits
from mantispy._core.frames import as_frame
from mantispy._core.logging import get_logger
from mantispy._core.mutation import inplace_or_copy
#: Operations, in the order pycytominer applies them by default, after mantispy's own ``drop_degenerate``.
OPERATIONS = (
"drop_degenerate",
"variance_threshold",
"frequency_threshold",
"correlation_threshold",
"drop_na_columns",
"blocklist",
"drop_outliers",
"noise_removal",
)
#: The default operations, cheapest first: removing features by name, by missingness and by variance before the
#: pairwise correlation means the correlation, the slowest step, compares only the features that are left.
#: Pass pycytominer's order, ``("variance_threshold", "correlation_threshold", "drop_na_columns", "blocklist")``,
#: to reproduce its selection.
DEFAULT_OPERATIONS = (
"drop_degenerate",
"blocklist",
"drop_na_columns",
"variance_threshold",
"correlation_threshold",
)
def feature_blocklist(name: str = "default") -> list[str]:
"""Feature names a blocklist drops.
:func:`feature_select` applies the blocklist for you (pass ``blocklist=...``); this is for inspecting the list itself.
Args:
name: Which blocklist to load; ``"default"`` is the bundled pycytominer list.
Returns:
The blocked feature names.
Raises:
ValueError: `name` names a blocklist that does not exist.
"""
return load_blocklist(name)
def _op_drop_degenerate(adata: AnnData) -> np.ndarray:
"""Drop the features :func:`~mantispy.pp.normalize` flagged in ``var["degenerate_scale"]`` because it could not scale them.
An object :func:`~mantispy.pp.normalize` never flagged has no such column and keeps every feature.
"""
if "degenerate_scale" not in adata.var:
return np.ones(adata.n_vars, dtype=bool)
return ~as_frame(adata.var)["degenerate_scale"].to_numpy(dtype=bool)
def _op_variance_threshold(X: np.ndarray, min_variance: float = 1e-6) -> np.ndarray:
"""Drop features whose variance is at or below ``min_variance``.
sklearn's ``VarianceThreshold`` keeps ``variance > threshold`` and uses the population variance, so ``ddof=0``.
"""
return np.nan_to_num(nanvar(X), nan=0.0, posinf=0.0) > min_variance
def _op_frequency_threshold(X: np.ndarray, freq_cut: float = 0.05, unique_cut: float = 0.01) -> np.ndarray:
"""Drop near-constant features, by value frequency or by uniqueness.
The two rules are independent and combined with OR, both with a strict ``<``.
"""
n_obs, n_vars = X.shape
keep = np.ones(n_vars, dtype=bool)
for j in range(n_vars):
column = X[:, j]
counts = np.sort(np.unique(column[~np.isnan(column)], return_counts=True)[1])[::-1]
frequency_ratio = counts[1] / counts[0] if counts.size >= 2 else 0.0
unique_ratio = counts.size / n_obs
keep[j] = not (frequency_ratio < freq_cut or unique_ratio < unique_cut)
return keep
def _greedy_keep(pairs: np.ndarray, total: np.ndarray, n_vars: int) -> np.ndarray:
"""Keep mask that drops the more-connected member of every over-threshold pair."""
ranking = np.argsort(total, kind="stable")
rank = np.empty(n_vars, dtype=np.int64)
rank[ranking] = np.arange(n_vars)
keep = np.ones(n_vars, dtype=bool)
if pairs.size:
first, second = pairs[:, 0], pairs[:, 1]
keep[np.where(rank[first] > rank[second], first, second)] = False
return keep
def _op_correlation_threshold(
X: np.ndarray,
threshold: float = 0.9,
method: str = "pearson",
*,
absolute: bool = False,
iterative: bool = False,
order: np.ndarray | None = None,
window: int | None = None,
stride: int | None = None,
) -> np.ndarray:
"""Drop one feature from every pair correlated above ``threshold``.
With the defaults this matches pycytominer: each over-threshold pair is judged on its own against a ranking of total absolute correlation, the member ranked as more correlated overall is dropped, the comparison is on the signed correlation, and there is no re-sweep.
"""
n_vars = X.shape[1]
reduce = _iterative_correlation_drop if iterative else _greedy_keep
if window is None:
pairs, total = correlated_pairs(X, threshold, method=method, absolute=absolute)
return reduce(pairs, total, n_vars)
pairs, total = correlated_pairs(
X, threshold, method=method, absolute=absolute, order=order, window=window, stride=stride
)
keep = reduce(pairs, total, n_vars)
survivors = np.flatnonzero(keep)
if survivors.size > 1:
sub_pairs, sub_total = correlated_pairs(X[:, survivors], threshold, method=method, absolute=absolute)
keep[survivors] = reduce(sub_pairs, sub_total, survivors.size)
return keep
def _iterative_correlation_drop(pairs: np.ndarray, total: np.ndarray, n_vars: int) -> np.ndarray:
"""Remove the most-connected feature until no over-threshold pair is left.
At each step the still-connected feature belonging to the most surviving pairs is dropped, ties broken by the largest total ``|r|`` and then the lowest index.
"""
keep = np.ones(n_vars, dtype=bool)
active = np.ones(len(pairs), dtype=bool)
while active.any():
connected, degree = np.unique(pairs[active].ravel(), return_counts=True)
# lexsort reads keys last-first: degree, then total, then index.
worst = connected[np.lexsort((connected, -total[connected], -degree))[0]]
keep[worst] = False
active &= (pairs[:, 0] != worst) & (pairs[:, 1] != worst)
return keep
def _op_drop_na_columns(X: np.ndarray, cutoff: float = 0.05) -> np.ndarray:
"""Drop features missing in more than ``cutoff`` of rows."""
return np.isnan(X).mean(axis=0) <= cutoff
def _op_drop_outliers(X: np.ndarray, outlier_cutoff: float = 500.0) -> np.ndarray:
"""Drop features whose largest absolute value exceeds ``outlier_cutoff``."""
# One column block at a time (in _blockwise), so np.abs never copies more than one block.
largest = _blockwise(X, lambda block: np.nanmax(np.abs(block), axis=0), np.float64)
return ~(np.nan_to_num(largest, nan=0.0) > outlier_cutoff)
def _op_blocklist(adata: AnnData, blocklist: str | Sequence[str] = "default") -> np.ndarray:
"""Drop features named in the blocklist.
Matched against the current names and, like :func:`~mantispy.pp.filter_features`, against ``var["original_name"]``, and for names without an object prefix, as cp_measure writes them, on the parsed feature identity.
"""
return ~var_blocklist_hits(adata, blocklist)
def _op_noise_removal(
X: np.ndarray, codes: np.ndarray, stdev_cutoff: float = 0.8, *, mode: str = "absolute"
) -> np.ndarray:
"""Drop features that vary too much within a perturbation group.
The within-group spread of a feature is the mean, over groups, of each group's population standard deviation (``ddof=0``).
``mode`` sets what ``stdev_cutoff`` is compared against:
- ``"absolute"``: the spread itself, pycytominer's behaviour, meaningful only on the scale the normalization left.
- ``"ratio"``: the spread divided by the feature's total standard deviation, so the cutoff means the same whichever normalization preceded it.
- ``"quantile"``: ``stdev_cutoff`` is read as a quantile of the observed within-group spreads, keeping the features below it.
"""
if mode not in ("absolute", "ratio", "quantile"):
raise ValueError(f"noise_removal_cutoff_mode must be 'absolute', 'ratio' or 'quantile', got {mode!r}")
n_groups = int(codes.max()) + 1
order, offsets = group_offsets(codes, n_groups)
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning) # "Degrees of freedom <= 0 for slice", once per all-NaN group
deviations = np.stack(
[np.nanstd(X[order[offsets[group] : offsets[group + 1]]], axis=0, ddof=0) for group in range(n_groups)]
)
within = np.nan_to_num(deviations.mean(axis=0), nan=0.0)
cutoff = float(np.nanquantile(within, stdev_cutoff)) if mode == "quantile" else stdev_cutoff
if mode == "ratio":
total = np.nan_to_num(np.nanstd(X, axis=0, ddof=0), nan=0.0)
# A feature with no total spread has no noise to remove, so its ratio is zero and it is kept.
within = np.divide(within, total, out=np.zeros_like(within), where=total > 0)
return ~(within > cutoff)
[docs]
@inplace_or_copy()
def feature_select(
adata: AnnData,
operations: Sequence[str] = DEFAULT_OPERATIONS,
min_variance: float = 1e-6,
freq_cut: float = 0.05,
unique_cut: float = 0.01,
corr_threshold: float = 0.9,
corr_method: str = "pearson",
corr_window: int | None = None,
corr_stride: int | None = None,
na_cutoff: float = 0.05,
outlier_cutoff: float = 500.0,
blocklist: str | Sequence[str] = "default",
noise_removal_perturb_groups: str = "Metadata_Perturbation",
noise_removal_stdev_cutoff: float = 0.8,
corr_absolute: bool = False,
corr_iterative: bool = False,
decorrelate: bool = False,
decorr_threshold: float = 0.99,
decorr_method: str = "pearson",
key_added: str = "selected",
copy: bool = False,
*,
noise_removal_cutoff_mode: str = "absolute",
subset: bool = False,
) -> AnnData | None:
"""Flag the features worth keeping.
Args:
adata: Object to select features on.
Usually well-level profiles.
operations: Which operations to run, from ``OPERATIONS``, in the order they run.
The default, ``DEFAULT_OPERATIONS``, holds pycytominer's usual operations plus ``drop_degenerate``, which removes nothing from an object :func:`~mantispy.pp.normalize` did not flag, ordered cheapest first so ``correlation_threshold`` runs last on what is left.
Passing pycytominer's order, ``["variance_threshold", "correlation_threshold", "drop_na_columns", "blocklist"]``, gives pycytominer's selection.
Without ``drop_outliers``, a warning names the kept features whose values sit more than ``outlier_cutoff`` robust standard deviations from their own median.
min_variance: ``variance_threshold``: keep features with variance above this.
freq_cut: ``frequency_threshold``: drop a feature when the count of its second most common value divided by the count of its most common is below this.
Either this rule or ``unique_cut`` drops a feature.
unique_cut: ``frequency_threshold``: drop a feature when its share of distinct values is below this.
corr_threshold: ``correlation_threshold``: drop one member of every pair correlated above this.
corr_method: ``correlation_threshold``: ``"pearson"`` or ``"spearman"``.
corr_window: ``correlation_threshold``: ``None`` runs the exact pass (the default, matching pycytominer).
An int switches to the two-pass fast path (prune redundancy within name-sorted windows of that size, then run the exact pass on the survivors); 500 is a good default, several times faster on large screens.
It keeps a different set of features (a different member of each correlated group) but preserves the information and the downstream signal.
corr_stride: ``correlation_threshold``: step between windows, default half the window.
corr_absolute: ``correlation_threshold``: threshold ``|r|`` rather than the signed correlation, so a strongly anti-correlated pair is also reduced.
Off by default, matching pycytominer.
corr_iterative: ``correlation_threshold``: drop the most-connected feature one at a time until no pair is left, keeping at least as many features as the single pass.
Off by default.
With ``corr_absolute`` this approximates cytominer's R path (``caret::findCorrelation``).
decorrelate: Run an extra, experimental redundancy step after the operations, on the features they keep.
Unlike ``correlation_threshold`` it removes features that are a linear combination of several others, not just pairwise duplicates, by a rank-revealing QR.
Off by default, and not part of pycytominer.
decorr_threshold: ``decorrelate``: drop a feature once its multiple correlation with the kept set reaches this.
The default 0.99 removes only near-collinear features, so the kept set spans almost the same space; lower it towards ``corr_threshold`` for a smaller, more aggressive set.
decorr_method: ``decorrelate``: ``"pearson"`` or ``"spearman"``.
na_cutoff: ``drop_na_columns``: drop features missing in more than this fraction of rows.
outlier_cutoff: ``drop_outliers``: drop features whose absolute value exceeds this.
blocklist: ``blocklist``: ``"default"`` for the bundled list, or explicit names.
Matched against the current names and against ``var["original_name"]``, so it works either side of :func:`~mantispy.pp.standardize_feature_names`.
A name without an object prefix, as cp_measure writes them, is matched on its parsed group, feature, channels and scale, so the default list also applies to cp_measure features.
noise_removal_perturb_groups: ``noise_removal``: ``obs`` column grouping replicates.
noise_removal_stdev_cutoff: ``noise_removal``: drop features whose within-group standard deviation, averaged over groups, is above this.
An absolute threshold on the scale :func:`~mantispy.pp.normalize` left the values on, so it is only meaningful next to the normalization that produced them; pycytominer's default assumes whole-plate standardization.
key_added: Name of the boolean ``var`` column to write.
copy: Return a modified copy instead of mutating in place.
noise_removal_cutoff_mode: ``noise_removal``: what ``noise_removal_stdev_cutoff`` is compared against.
``"absolute"`` (the default) keeps pycytominer's behaviour, comparing the raw within-group spread and so tied to the normalization scale.
``"ratio"`` compares the within-group spread to the feature's total spread, and ``"quantile"`` reads the cutoff as a quantile of the observed within-group spreads; both keep a useful subset under the control-based or MAD normalization the package recommends, where the absolute cutoff keeps everything or nothing.
subset: Inplace subset to the selected features, as well as flagging them in ``var[key_added]``.
With ``copy=True`` the returned copy is the subset.
Returns:
``None``, or the modified copy.
Writes ``var[key_added]`` and a per-operation count of removals to ``uns["mantispy"]["feature_select"]``, each count being what that operation removed from the features the earlier ones kept.
Nothing is dropped unless ``subset=True``; :func:`subset_features` does it afterwards.
Raises:
ValueError: If ``operations`` names an operation that is not in ``OPERATIONS``.
KeyError: If ``noise_removal`` is requested but ``noise_removal_perturb_groups`` is not an ``obs`` column.
Notes:
``drop_degenerate`` runs first: a feature :func:`~mantispy.pp.normalize` could not scale can hold values large enough to decide the correlation ranking of every feature it is compared with.
The other operations then run in the order ``operations`` lists them, each judging only the features the earlier ones kept, as pycytominer does.
The counts in ``uns["mantispy"]["feature_select"]`` therefore do not overlap, and they sum to the number of features removed.
Order matters only for ``correlation_threshold``, which ranks a feature against the other features still present; every other operation judges each feature on its own values.
So the default order, which runs the correlation last, can keep a different member of a correlated pair than pycytominer's order, which runs it before ``drop_na_columns`` and ``blocklist``; it is also much faster on data with many mostly-missing features.
``decorrelate`` runs last, on the features the operations kept.
``correlation_threshold`` is the most expensive operation.
pycytominer uses ``pandas.DataFrame.corr``, one Cython pass per column pair.
Here the pairs come from chunked matrix products over blocks of columns, so no ``n_vars ** 2`` array is held in memory.
"""
unknown = set(operations) - set(OPERATIONS)
if unknown:
raise ValueError(f"unknown operation(s) {sorted(unknown)}; choose from {OPERATIONS}")
# A non-empty not_measured implies the feature_kind column exists (feature_kind_of returns {} otherwise).
# Warn only when the non-measurement features are the majority: a few embedding columns among mostly real
# measurements still selects meaningfully on the measurements.
not_measured = feature_kind_of(adata) & {"embedding", "derived"}
n_not_measured = int(adata.var["feature_kind"].isin(not_measured).sum()) if not_measured else 0
if n_not_measured * 2 > adata.n_vars:
get_logger().warning(
"feature_select: %d of %d features are %s, not measurements, so the selection is unlikely to mean much",
n_not_measured,
adata.n_vars,
"/".join(sorted(not_measured)),
)
full = get_matrix(adata)
keep = np.ones(adata.n_vars, dtype=bool)
removed: dict[str, int] = {}
if "drop_degenerate" in operations:
keep = _op_drop_degenerate(adata)
removed["drop_degenerate"] = int((~keep).sum())
judged = np.flatnonzero(keep)
X = full[:, judged] if judged.size < adata.n_vars else full
# Each operation judges only the columns of X that the earlier ones kept, as pycytominer does.
alive = np.ones(judged.size, dtype=bool)
n_judged_by: dict[str, int] = {}
for operation in operations:
if operation == "drop_degenerate":
continue
n_judged_by[operation] = int(alive.sum())
if operation == "variance_threshold":
mask = _op_variance_threshold(X, min_variance)
elif operation == "frequency_threshold":
mask = np.ones(judged.size, dtype=bool)
mask[alive] = _op_frequency_threshold(X[:, alive], freq_cut, unique_cut)
elif operation == "correlation_threshold":
# The only operation whose verdict on a column depends on the other columns present.
columns = np.flatnonzero(alive)
order = None
if corr_window is not None:
# Names of the columns passed in, so order indexes them 0..columns.size-1.
names = np.asarray(adata.var_names[judged[columns]], dtype=str)
order = np.argsort(names, kind="stable")
mask = np.ones(judged.size, dtype=bool)
mask[columns] = _op_correlation_threshold(
X if columns.size == judged.size else X[:, columns],
corr_threshold,
corr_method,
absolute=corr_absolute,
iterative=corr_iterative,
order=order,
window=corr_window,
stride=corr_stride,
)
elif operation == "drop_na_columns":
mask = _op_drop_na_columns(X, na_cutoff)
elif operation == "blocklist":
mask = _op_blocklist(adata, blocklist)[judged]
elif operation == "drop_outliers":
mask = _op_drop_outliers(X, outlier_cutoff)
else:
# Unknown names were rejected above, so what is left is `noise_removal`.
if noise_removal_perturb_groups not in adata.obs:
raise KeyError(f"obs has no column {noise_removal_perturb_groups!r} to group replicates by")
codes, _ = group_codes(adata, noise_removal_perturb_groups)
mask = _op_noise_removal(X, codes, noise_removal_stdev_cutoff, mode=noise_removal_cutoff_mode)
removed[operation] = int((alive & ~mask).sum())
alive &= mask
keep[judged] = alive
if decorrelate:
survivors = np.flatnonzero(keep)
mask = rank_revealing_subset(full[:, survivors], decorr_threshold, decorr_method)
removed["decorrelate"] = int((~mask).sum())
keep[survivors[~mask]] = False
noise_stripped = removed.get("noise_removal", 0)
# Only the absolute cutoff is tied to the normalization scale, so only it earns the scale-relative advice.
n_noise_judged = n_judged_by.get("noise_removal", 0)
near_total_noise = noise_removal_cutoff_mode == "absolute" and noise_stripped > 0.95 * n_noise_judged
if not keep.any() and adata.n_vars:
get_logger().warning(
"feature_select: kept none of the %d features",
adata.n_vars,
deep="each cutoff is absolute on the scale pp.normalize left, and noise_removal's stdev_cutoff is the usual cause",
)
elif near_total_noise:
get_logger().warning(
"feature_select: noise_removal dropped %d of %d features",
noise_stripped,
n_noise_judged,
deep="stdev_cutoff is absolute, so set noise_removal_cutoff_mode='ratio' or 'quantile' for a relative cutoff",
)
# A feature kept with NaN in fewer than na_cutoff of its rows still carries those missing values, which
# pp.sphere and sc.pp.pca cannot take. Warn with the lever to remove them here, rather than let the matrix
# fail several steps later with an error that does not point back to feature selection.
kept = np.flatnonzero(keep)
# any(axis=0) over the full matrix, then index the kept columns: avoids copying the kept
# sub-matrix (fancy column indexing on a wide screen is a full-size float copy).
with_nan = kept[np.isnan(full).any(axis=0)[kept]] if kept.size else kept
if with_nan.size:
lever = (
"lower na_cutoff to 0.0" if "drop_na_columns" in operations else "add drop_na_columns with na_cutoff=0.0"
)
get_logger().warning(
"feature_select: %d selected feature(s) still hold NaN, which pp.sphere and PCA cannot take",
with_nan.size,
deep=f"{lever} and subset, or fill them with mt.pp.impute",
)
if "drop_outliers" not in operations:
_warn_extreme(adata, full, kept, outlier_cutoff)
adata.var[key_added] = keep
adata.uns.setdefault("mantispy", {})["feature_select"] = removed
get_logger().info(
"feature_select kept %d of %d features (%s)",
int(keep.sum()),
adata.n_vars,
", ".join(f"{name} -{count}" for name, count in removed.items()),
)
if subset:
adata._inplace_subset_var(keep)
return None
def _warn_extreme(adata: AnnData, full: np.ndarray, kept: np.ndarray, cutoff: float) -> None:
"""Warn when a kept feature holds values far outside its own spread, which the default operations let through.
Each feature is scaled by its own median and MAD, so the check means the same on raw and on normalized values, and ``outlier_cutoff`` is the bound ``drop_outliers`` would apply to a robustly scaled feature.
"""
if not kept.size:
return
# One parallel pass over the kept columns, each read straight from the matrix: no float64 copy of the
# kept block and no per-column Python dispatch, which is otherwise the whole cost of feature_select.
largest = column_max_robust_z(full, kept, MAD_TO_SIGMA)
extreme = kept[largest > cutoff]
if extreme.size:
names = [str(adata.var_names[index]) for index in extreme]
shown = ", ".join(names[:5]) + (f" and {len(names) - 5} more" if len(names) > 5 else "")
get_logger().warning(
"feature_select: %d selected feature(s) lie more than %g robust standard deviations from their median (%s)",
len(names),
cutoff,
shown,
deep="add 'drop_outliers' to operations, or check the single cells they come from",
)
[docs]
def subset_features(adata: AnnData, key: str = "selected") -> AnnData:
"""Return a new object holding only the features flagged by ``var[key]``.
Args:
adata: Object to subset.
Never modified.
key: Boolean ``var`` column naming the features to keep, as :func:`feature_select` writes.
Returns:
A new object holding only the flagged features.
Raises:
KeyError: If ``var`` has no column ``key``.
"""
if key not in adata.var:
raise KeyError(f"var has no column {key!r}; run mt.pp.feature_select first")
return adata[:, as_frame(adata.var)[key].to_numpy(dtype=bool)].copy()
[docs]
def decorr_threshold_sweep(
adata: AnnData,
thresholds: Sequence[float],
scorers: Mapping[str, Callable[[AnnData], float]] | Sequence[Callable[[AnnData], float]],
*,
key: str = "selected",
**feature_select_kwargs: Any,
) -> pd.DataFrame:
"""Score ``decorrelate`` at a range of thresholds so the smallest set that holds a metric can be chosen.
``decorr_threshold`` trades size for signal, and the right setting depends on the screen and on the metric
that matters, so there is no universal best. This runs :func:`feature_select` with ``decorrelate=True`` at each
threshold on a copy, scores the selected object with each callable, and returns one row per threshold.
Args:
adata: Object to select features on. Never modified.
thresholds: ``decorr_threshold`` values to try.
scorers: Scoring callables, each taking the selected :class:`~anndata.AnnData` and returning a scalar (for
example replicate or activity mAP). A mapping names the columns; a sequence names them by
``callable.__name__``.
key: Boolean ``var`` column :func:`feature_select` writes and :func:`subset_features` reads.
feature_select_kwargs: Passed through to :func:`feature_select` (``operations``, ``corr_threshold``, and so
on); ``decorrelate``, ``decorr_threshold``, ``key_added`` and ``copy`` are set here.
Returns:
A frame with a ``threshold`` column, an ``n_kept`` column, and one column per scorer.
"""
named = list(scorers.items() if isinstance(scorers, Mapping) else ((s.__name__, s) for s in scorers))
records = []
for threshold in thresholds:
selected = feature_select(
adata,
decorrelate=True,
decorr_threshold=float(threshold),
key_added=key,
copy=True,
**feature_select_kwargs,
)
assert selected is not None # feature_select(copy=True) returns the modified copy
kept = subset_features(selected, key)
row: dict[str, float] = {"threshold": float(threshold), "n_kept": int(kept.n_vars)}
for name, scorer in named:
row[name] = float(scorer(kept))
records.append(row)
return pd.DataFrame(records)
