"""Feature selection, with pycytominer's semantics.
Every operation returns a boolean keep mask over ``var``. :func:`feature_select` combines
them into one boolean column without dropping anything, and :func:`subset_features` does
the subsetting.
Operation names and behavior follow pycytominer. ``variance_threshold`` is an
sklearn-style variance cut, and the frequency and uniqueness rules are in the separate
``frequency_threshold``. ``correlation_threshold`` judges each pair against a ranking
computed once from the full matrix instead of sweeping greedily, and thresholds the signed
correlation, so two features correlated at -1.0 are both kept.
"""
from __future__ import annotations
from collections.abc import Sequence
import numpy as np
from anndata import AnnData
from mantispy._core._corr import correlated_pairs
from mantispy._core._reduce import get_matrix, group_codes
from mantispy._core._utils import as_frame, get_logger, inplace_or_copy
from mantispy._core.features import blocklist_hits
#: Operations, in the order pycytominer applies them by default.
OPERATIONS = (
"variance_threshold",
"frequency_threshold",
"correlation_threshold",
"drop_na_columns",
"blocklist",
"drop_outliers",
"noise_removal",
)
DEFAULT_OPERATIONS = (
"variance_threshold",
"correlation_threshold",
"drop_na_columns",
"blocklist",
)
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``.
"""
with np.errstate(invalid="ignore"):
variance = np.nanvar(X, axis=0, ddof=0)
return np.nan_to_num(variance, 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 _op_correlation_threshold(X: np.ndarray, threshold: float = 0.9, method: str = "pearson") -> np.ndarray:
"""Drop one feature from every pair correlated above ``threshold``.
Each over-threshold pair is judged on its own against a ranking of total absolute
correlation computed once from the full matrix; the member ranked as more correlated
overall is dropped. There is no iterative sweep, so a feature already dropped by one
pair does not spare its partner in another.
"""
n_vars = X.shape[1]
# Signed correlation, as in pycytominer, which keeps a pair correlated at -1.0.
# correlated_pairs never builds the full matrix, so this scales to tens of thousands of features.
pairs, total = correlated_pairs(X, threshold, method=method)
# Rank features by how correlated they are with everything else, ascending.
order = np.argsort(total, kind="stable")
rank = np.empty(n_vars, dtype=np.int64)
rank[order] = 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_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``.
Ratios with a near-zero denominator blow up like this.
"""
with np.errstate(invalid="ignore"):
largest = np.nanmax(np.abs(X), axis=0)
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"]``, so the blocklist still applies after
:func:`~mantispy.pp.standardize_feature_names` has renamed the features.
"""
names = [adata.var_names.to_numpy()]
if "original_name" in adata.var:
names.append(adata.var["original_name"].astype(str).to_numpy())
return ~blocklist_hits(names, blocklist)
def _op_noise_removal(X: np.ndarray, codes: np.ndarray, stdev_cutoff: float = 0.8) -> np.ndarray:
"""Drop features that vary too much within a perturbation group.
The statistic is the mean, over groups, of each group's population standard
deviation (``ddof=0``).
"""
n_groups = int(codes.max()) + 1
with np.errstate(invalid="ignore"):
deviations = np.stack([np.nanstd(X[codes == group], axis=0, ddof=0) for group in range(n_groups)])
return ~(np.nan_to_num(deviations.mean(axis=0), nan=0.0) > stdev_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",
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,
key_added: str = "selected",
copy: 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``. The default omits
``frequency_threshold``, ``drop_outliers`` and ``noise_removal``, matching
pycytominer's own default.
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"``.
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`.
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.
key_added: Name of the boolean ``var`` column to write.
copy: Return a modified copy instead of mutating in place.
Returns:
``None``, or the modified copy. Writes ``var[key_added]`` and a per-operation count
of removals to ``uns["mantispy"]["feature_select"]``. Nothing is dropped; use
:func:`subset_features` for that.
Notes:
``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}")
X = get_matrix(adata)
keep = np.ones(adata.n_vars, dtype=bool)
removed: dict[str, int] = {}
for operation in operations:
before = keep.copy()
if operation == "variance_threshold":
keep &= _op_variance_threshold(X, min_variance)
elif operation == "frequency_threshold":
keep &= _op_frequency_threshold(X, freq_cut, unique_cut)
elif operation == "correlation_threshold":
keep &= _op_correlation_threshold(X, corr_threshold, corr_method)
elif operation == "drop_na_columns":
keep &= _op_drop_na_columns(X, na_cutoff)
elif operation == "blocklist":
keep &= _op_blocklist(adata, blocklist)
elif operation == "drop_outliers":
keep &= _op_drop_outliers(X, outlier_cutoff)
elif operation == "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)
keep &= _op_noise_removal(X, codes, noise_removal_stdev_cutoff)
removed[operation] = int((before & ~keep).sum())
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()),
)
return None
[docs]
def subset_features(adata: AnnData, key: str = "selected") -> AnnData:
"""Return a new object holding only the features flagged by ``var[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()