"""Per-metacell quality metrics for CAMP partitions."""
from __future__ import annotations
import numpy as np
import pandas as pd
import scanpy as sc
from anndata import AnnData
from scipy import sparse
from scipy.stats import entropy as scipy_entropy
from sklearn.metrics import pairwise_distances
from campmc._logging import get_logger
logger = get_logger(__name__)
def _purity(cell_labels: pd.Series) -> float:
counts = cell_labels.value_counts(dropna=False)
if len(counts) == 0:
return np.nan
return float(counts.max()) / float(counts.sum())
def _label_entropy(cell_labels: pd.Series) -> float:
counts = cell_labels.value_counts(dropna=False).values.astype(float)
if counts.sum() == 0:
return np.nan
p = counts / counts.sum()
return float(scipy_entropy(p, base=2))
def _compactness(diff_coords: np.ndarray) -> float:
return float(-np.var(diff_coords, axis=0).mean())
def _inv(expr_block) -> float:
x = expr_block.toarray() if sparse.issparse(expr_block) else np.asarray(expr_block)
gene_var = np.var(x, axis=0)
gene_mean = np.mean(x, axis=0)
return float(np.percentile(gene_var / (gene_mean + 1e-8), 95))
def _ensure_diffmap(adata: AnnData, use_rep: str) -> None:
if "X_diffmap" in adata.obsm:
return
if use_rep not in adata.obsm:
raise KeyError(f"Missing obsm['{use_rep}'] for diffusion map computation.")
sc.pp.neighbors(adata, use_rep=use_rep)
sc.tl.diffmap(adata)
[docs]
def quality(
adata: AnnData,
*,
partition_key: str = "camp",
label_key: str | None = None,
use_rep: str = "X_pca",
) -> pd.DataFrame:
"""Compute per-metacell quality metrics.
Metrics include compactness and separation on the diffusion map,
an INV (intra-metacell variability) score, and optionally purity /
label entropy when ``label_key`` is provided.
Parameters
----------
adata
Annotated data with metacell labels in ``obs[partition_key]``.
partition_key
Column in ``adata.obs`` holding metacell assignments.
label_key
Optional cell-type (or other) labels for purity and entropy.
use_rep
Representation used to compute neighbors / diffusion map when
``X_diffmap`` is not already present.
Returns
-------
DataFrame
One row per metacell with columns such as ``size``, ``compactness``,
``separation``, ``sc_ratio``, ``INV``, and optionally ``purity`` /
``label_entropy``.
Examples
--------
.. exec-jupyter::
import campmc as cp
import scanpy as sc
adata = sc.datasets.pbmc68k_reduced()
cp.partition(adata, method="camp3", gamma=50, random_state=0)
qc = cp.mt.quality(adata, partition_key="camp", label_key="louvain")
qc[["metacell_id", "size", "purity", "sc_ratio"]].head()
"""
if partition_key not in adata.obs:
raise KeyError(f"obs missing '{partition_key}'")
work = adata.copy()
_ensure_diffmap(work, use_rep)
x_diff = work.obsm["X_diffmap"]
x_diff_dense = x_diff if not sparse.issparse(x_diff) else x_diff.toarray()
work.obs["metacell"] = work.obs[partition_key].astype(str)
grouped = work.obs.groupby("metacell", observed=True)
centroids: dict[str, np.ndarray] = {}
for metacell_id, group in grouped:
ix = work.obs.index.get_indexer(group.index)
centroids[metacell_id] = x_diff_dense[ix].mean(axis=0)
mc_ids = list(centroids.keys())
if not mc_ids:
return pd.DataFrame()
c_mat = np.vstack([centroids[mid] for mid in mc_ids])
dist = pairwise_distances(c_mat, c_mat)
np.fill_diagonal(dist, np.inf)
nearest_sep = np.min(dist, axis=1)
sep_map = {mid: float(nearest_sep[j]) for j, mid in enumerate(mc_ids)}
rows: list[dict] = []
for metacell_id, group in grouped:
ix = work.obs.index.get_indexer(group.index)
diff_coords = x_diff_dense[ix]
compact = _compactness(diff_coords)
separation = sep_map[metacell_id]
row = {
"metacell_id": metacell_id,
"size": len(group),
"compactness": compact,
"separation": separation,
"sc_ratio": separation / np.sqrt(max(1e-12, -compact)),
"INV": _inv(work.X[ix, :]),
}
if label_key is not None and label_key in work.obs:
labels = work.obs.loc[group.index, label_key]
row["purity"] = _purity(labels.astype(str))
row["label_entropy"] = _label_entropy(labels.astype(str))
rows.append(row)
logger.info("quality metrics for %d metacells", len(rows))
return pd.DataFrame(rows)