diff --git a/docs/api/preprocessing.md b/docs/api/preprocessing.md index 3b6c684deb..774e56062b 100644 --- a/docs/api/preprocessing.md +++ b/docs/api/preprocessing.md @@ -76,6 +76,16 @@ Also see {ref}`data integration tools ` and external {ref}`ext pp.scrublet_simulate_doublets ``` +## Sample demultiplexing + +```{eval-rst} +.. autosummary:: + :nosignatures: + :toctree: generated/ + + pp.hashsolo +``` + ## Neighbors ```{eval-rst} diff --git a/docs/conf.py b/docs/conf.py index 2b88b5b724..fb76d2c373 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -184,6 +184,7 @@ "pp.filter_cells": (["np", "sp", "da"], []), "pp.filter_genes": (["np", "sp", "da"], []), "pp.harmony_integrate": (["np"], []), + "pp.hashsolo": (["np", "sp"], []), "pp.highly_variable_genes": (["np", "sp", "da"], ["da[sp[csc]]"]), "pp.log1p": (["np", "sp", "da"], []), "pp.neighbors": (["np", "sp"], []), diff --git a/docs/external/preprocessing.md b/docs/external/preprocessing.md index 312f7bb121..9f830481b8 100644 --- a/docs/external/preprocessing.md +++ b/docs/external/preprocessing.md @@ -7,6 +7,7 @@ Previously found here, but now part of scanpy’s main API: - {func}`scanpy.pp.harmony_integrate` +- {func}`scanpy.pp.hashsolo` - {func}`scanpy.pp.scrublet` - {func}`scanpy.pp.scrublet_simulate_doublets` diff --git a/docs/release-notes/4303.feat.md b/docs/release-notes/4303.feat.md new file mode 100644 index 0000000000..4278c13e08 --- /dev/null +++ b/docs/release-notes/4303.feat.md @@ -0,0 +1,3 @@ +Move {func}`scanpy.external.pp.hashsolo` into scanpy’s main API as {func}`scanpy.pp.hashsolo` and deprecate the old location. +The new function takes {mod}`anndata.acc` references, so the hashing counts can live in {attr}`~anndata.AnnData.obs`, in the data matrix (e.g. `A.X[:, ["Hash1", "Hash2"]]`), or all together in {attr}`~anndata.AnnData.obsm` (e.g. `A.obsm["hto"]`), +and writes its results to `.obs[key_added]` and `.obsm[key_added]` instead of six unnamespaced `.obs` columns {smaller}`P Angerer` diff --git a/src/scanpy/external/pp/_hashsolo.py b/src/scanpy/external/pp/_hashsolo.py index 410a19417a..d98497a31f 100644 --- a/src/scanpy/external/pp/_hashsolo.py +++ b/src/scanpy/external/pp/_hashsolo.py @@ -1,296 +1,25 @@ -"""A probabilistic cell hashing demultiplexing method. +"""Deprecated location of the HashSolo cell hashing demultiplexing method. -HashSolo generates a noise distribution and signal distribution -for each hashing barcode from empirically observed counts. -These distributions are updates from the global signal and noise barcode distributions, -which helps in the setting where not many cells are observed. -For a hashing barcode: - -Signal distributions - are estimated from samples where that hashing barcode has the highest count. - -Noise distributions - are estimated from samples where that hashing barcode is one the k-2 lowest barcodes, - where k is the number of barcodes. - -We test each of the following hypotheses in a bayesian fashion, -and select the most probable hypothesis. - -A doublet - should have its two highest barcode counts most likely - coming from a signal distribution for those barcodes. - -A singlet - should have its highest barcode from a signal distribution, - and its second highest barcode from a noise distribution. - -A negative two highest barcodes - should come from noise distributions. +See :func:`scanpy.pp.hashsolo`. """ from __future__ import annotations -from itertools import product from typing import TYPE_CHECKING import numpy as np -import pandas as pd -from scipy.stats import norm +from scverse_misc import Deprecation, deprecated -from ..._utils import check_nonnegative_integers from ..._utils._doctests import doctest_skipif +from ...preprocessing._hashsolo import _hashsolo if TYPE_CHECKING: from collections.abc import Sequence from anndata import AnnData - from numpy.typing import ArrayLike, NDArray - - -def _calculate_log_likelihoods( # noqa: PLR0915 - data: np.ndarray, number_of_noise_barcodes: int -) -> tuple[NDArray[np.float64], NDArray[np.float64], dict[int, str]]: - """Calculate log likelihoods for each hypothesis, negative, singlet, doublet. - - Parameters - ---------- - data - cells by hashing counts matrix - number_of_noise_barcodes - number of barcodes to used to calculated noise distribution - - Returns - ------- - log_likelihoods_for_each_hypothesis - a 2d np.array log likelihood of each hypothesis - all_indices - counter_to_barcode_combo - - """ - - def gaussian_updates( - data: np.ndarray, mu_o: float, std_o: float - ) -> tuple[float, float]: - """Update parameters of your gaussian. - - See . - - Parameters - ---------- - data - 1-d array of counts - mu_o - global mean for hashing count distribution - std_o - global std for hashing count distribution - - Returns - ------- - mean - of gaussian - std - of gaussian - - """ - lam_o = 1 / (std_o**2) - n = len(data) - lam = 1 / np.var(data) if len(data) > 1 else lam_o - lam_n = lam_o + n * lam - mu_n = ( - (np.mean(data) * n * lam + mu_o * lam_o) / lam_n if len(data) > 0 else mu_o - ) - return mu_n, (1 / (lam_n / (n + 1))) ** (1 / 2) - - eps = 1e-15 - # probabilites for negative, singlet, doublets - log_likelihoods_for_each_hypothesis = np.zeros((data.shape[0], 3)) - - all_indices = np.empty(data.shape[0]) - num_of_barcodes = data.shape[1] - number_of_non_noise_barcodes = ( - num_of_barcodes - number_of_noise_barcodes - if number_of_noise_barcodes is not None - else 2 - ) - - num_of_noise_barcodes = num_of_barcodes - number_of_non_noise_barcodes - - # assume log normal - data = np.log(data + 1) - data_arg = np.argsort(data, axis=1) - data_sort = np.sort(data, axis=1) - - # global signal and noise counts useful for when we have few cells - # barcodes with the highest number of counts are assumed to be a true signal - # barcodes with rank < k are considered to be noise - global_signal_counts = np.ravel(data_sort[:, -1]) - global_noise_counts = np.ravel(data_sort[:, :-number_of_non_noise_barcodes]) - global_mu_signal_o, global_sigma_signal_o = ( - np.mean(global_signal_counts), - np.std(global_signal_counts), - ) - global_mu_noise_o, global_sigma_noise_o = ( - np.mean(global_noise_counts), - np.std(global_noise_counts), - ) - - noise_params_dict = {} - signal_params_dict = {} - - # for each barcode get empirical noise and signal distribution parameterization - for x in np.arange(num_of_barcodes): - sample_barcodes = data[:, x] - sample_barcodes_noise_idx = np.where(data_arg[:, :num_of_noise_barcodes] == x)[ - 0 - ] - sample_barcodes_signal_idx = np.where(data_arg[:, -1] == x) - - # get noise and signal counts - noise_counts = sample_barcodes[sample_barcodes_noise_idx] - signal_counts = sample_barcodes[sample_barcodes_signal_idx] - - # get parameters of distribution, assuming lognormal do update from global values - noise_param = gaussian_updates( - noise_counts, global_mu_noise_o, global_sigma_noise_o - ) - signal_param = gaussian_updates( - signal_counts, global_mu_signal_o, global_sigma_signal_o - ) - noise_params_dict[x] = noise_param - signal_params_dict[x] = signal_param - - counter_to_barcode_combo: dict[int, str] = {} - counter = 0 - - # for each combination of noise and signal barcode calculate probiltiy of in silico and real cell hypotheses - for noise_sample_idx, signal_sample_idx in product( - np.arange(num_of_barcodes), np.arange(num_of_barcodes) - ): - signal_subset = data_arg[:, -1] == signal_sample_idx - noise_subset = data_arg[:, -2] == noise_sample_idx - subset = signal_subset & noise_subset - if sum(subset) == 0: - continue - - indices = np.where(subset)[0] - barcode_combo = "_".join([str(noise_sample_idx), str(signal_sample_idx)]) - all_indices[np.where(subset)[0]] = counter - counter_to_barcode_combo[counter] = barcode_combo - counter += 1 - noise_params = noise_params_dict[noise_sample_idx] - signal_params = signal_params_dict[signal_sample_idx] - - # calculate probabilties for each hypothesis for each cell - data_subset = data[subset] - log_signal_signal_probs = np.log( - norm.pdf( - data_subset[:, signal_sample_idx], - *signal_params[:-2], - loc=signal_params[-2], - scale=signal_params[-1], - ) - + eps - ) - signal_noise_params = signal_params_dict[noise_sample_idx] - log_noise_signal_probs = np.log( - norm.pdf( - data_subset[:, noise_sample_idx], - loc=signal_noise_params[-2], - scale=signal_noise_params[-1], - ) - + eps - ) - - log_noise_noise_probs = np.log( - norm.pdf( - data_subset[:, noise_sample_idx], - loc=noise_params[-2], - scale=noise_params[-1], - ) - + eps - ) - log_signal_noise_probs = np.log( - norm.pdf( - data_subset[:, signal_sample_idx], - loc=noise_params[-2], - scale=noise_params[-1], - ) - + eps - ) - - probs_of_negative = np.sum( - [log_noise_noise_probs, log_signal_noise_probs], axis=0 - ) - probs_of_singlet = np.sum( - [log_noise_noise_probs, log_signal_signal_probs], axis=0 - ) - probs_of_doublet = np.sum( - [log_noise_signal_probs, log_signal_signal_probs], axis=0 - ) - log_probs_list = [probs_of_negative, probs_of_singlet, probs_of_doublet] - - # each cell and each hypothesis probability - for prob_idx, log_prob in enumerate(log_probs_list): - log_likelihoods_for_each_hypothesis[indices, prob_idx] = log_prob - return ( - log_likelihoods_for_each_hypothesis, - all_indices, - counter_to_barcode_combo, - ) - - -def _calculate_bayes_rule( - data: np.ndarray, priors: ArrayLike, number_of_noise_barcodes: int -) -> dict[str, np.ndarray]: - """Calculate bayes rule from log likelihoods. - - Parameters - ---------- - data - Anndata object filled only with hashing counts - priors - a list of your prior for each hypothesis - first element is your prior for the negative hypothesis - second element is your prior for the singlet hypothesis - third element is your prior for the doublet hypothesis - We use [0.01, 0.8, 0.19] by default because we assume the barcodes - in your cell hashing matrix are those cells which have passed QC - in the transcriptome space, e.g. UMI counts, pct mito reads, etc. - number_of_noise_barcodes - number of barcodes to used to calculated noise distribution - - Returns - ------- - A dict of bayes key results with the following entries: - - `"most_likely_hypothesis"` - A 1d np.array of the most likely hypothesis - `"probs_hypotheses"` - A 2d np.array probability of each hypothesis - `"log_likelihoods_for_each_hypothesis"` - A 2d np.array log likelihood of each hypothesis - - """ - priors = np.array(priors) - log_likelihoods_for_each_hypothesis, _, _ = _calculate_log_likelihoods( - data, number_of_noise_barcodes - ) - probs_hypotheses = ( - np.exp(log_likelihoods_for_each_hypothesis) - * priors - / np.sum( - np.multiply(np.exp(log_likelihoods_for_each_hypothesis), priors), - axis=1, - )[:, None] - ) - most_likely_hypothesis = np.argmax(probs_hypotheses, axis=1) - return { - "most_likely_hypothesis": most_likely_hypothesis, - "probs_hypotheses": probs_hypotheses, - "log_likelihoods_for_each_hypothesis": log_likelihoods_for_each_hypothesis, - } +@deprecated(Deprecation("1.13.0", "Use :func:`scanpy.pp.hashsolo` instead.")) @doctest_skipif(reason="Illustrative but not runnable doctest code") def hashsolo( adata: AnnData, @@ -303,8 +32,9 @@ def hashsolo( ) -> AnnData | None: """Probabilistic demultiplexing of cell hashing data using HashSolo :cite:p:`Bernstein2020`. - .. note:: - More information and bug reports `here `__. + :func:`scanpy.pp.hashsolo` accepts :mod:`anndata.acc` references, + so the hashing counts may live in :attr:`~anndata.AnnData.obs` *or* in the data matrix, + and it writes namespaced result fields. Parameters ---------- @@ -352,109 +82,41 @@ def hashsolo( Examples -------- - >>> import anndata - >>> import scanpy.external as sce - >>> adata = anndata.read_h5ad("data.h5ad") - >>> sce.pp.hashsolo(adata, ["Hash1", "Hash2", "Hash3"]) + >>> import scanpy as sc + >>> adata = sc.io.read_h5ad("data.h5ad") + >>> sc.external.pp.hashsolo(adata, ["Hash1", "Hash2", "Hash3"]) >>> adata.obs.head() """ print( "Please cite HashSolo paper:\nhttps://www.cell.com/cell-systems/fulltext/S2405-4712(20)30195-2" ) - adata = adata.copy() if not inplace else adata + adata = adata if inplace else adata.copy() + cell_hashing_columns = list(cell_hashing_columns) data = adata.obs[cell_hashing_columns].to_numpy() - if not check_nonnegative_integers(data): - msg = "Cell hashing counts must be non-negative" - raise ValueError(msg) - if (number_of_noise_barcodes is not None) and ( - number_of_noise_barcodes >= len(cell_hashing_columns) - ): - msg = "number_of_noise_barcodes must be at least one less \ - than the number of samples you have as determined by the number of \ - cell_hashing_columns you've given as input " - raise ValueError(msg) - num_of_cells = adata.shape[0] - results = pd.DataFrame( - np.zeros((num_of_cells, 6)), - columns=[ - "most_likely_hypothesis", - "probs_hypotheses", - "cluster_feature", - "negative_hypothesis_probability", - "singlet_hypothesis_probability", - "doublet_hypothesis_probability", - ], - index=adata.obs_names, + clusters = ( + None + if pre_existing_clusters is None + else adata.obs[pre_existing_clusters].to_numpy() ) - if pre_existing_clusters is not None: - cluster_features = pre_existing_clusters - unique_cluster_features = np.unique(adata.obs[cluster_features]) - for cluster_feature in unique_cluster_features: - cluster_feature_bool_vector = adata.obs[cluster_features] == cluster_feature - posterior_dict = _calculate_bayes_rule( - data[cluster_feature_bool_vector], - priors, - number_of_noise_barcodes, - ) - results.loc[cluster_feature_bool_vector, "most_likely_hypothesis"] = ( - posterior_dict["most_likely_hypothesis"] - ) - results.loc[cluster_feature_bool_vector, "cluster_feature"] = ( - cluster_feature - ) - results.loc[ - cluster_feature_bool_vector, "negative_hypothesis_probability" - ] = posterior_dict["probs_hypotheses"][:, 0] - results.loc[ - cluster_feature_bool_vector, "singlet_hypothesis_probability" - ] = posterior_dict["probs_hypotheses"][:, 1] - results.loc[ - cluster_feature_bool_vector, "doublet_hypothesis_probability" - ] = posterior_dict["probs_hypotheses"][:, 2] - else: - posterior_dict = _calculate_bayes_rule(data, priors, number_of_noise_barcodes) - results.loc[:, "most_likely_hypothesis"] = posterior_dict[ - "most_likely_hypothesis" - ] - results.loc[:, "cluster_feature"] = 0 - results.loc[:, "negative_hypothesis_probability"] = posterior_dict[ - "probs_hypotheses" - ][:, 0] - results.loc[:, "singlet_hypothesis_probability"] = posterior_dict[ - "probs_hypotheses" - ][:, 1] - results.loc[:, "doublet_hypothesis_probability"] = posterior_dict[ - "probs_hypotheses" - ][:, 2] + probs = _hashsolo( + data, + priors=priors, + clusters=clusters, + n_barcodes_noise=number_of_noise_barcodes, + ) + most_likely_hypothesis = np.argmax(probs, axis=1) - adata.obs["most_likely_hypothesis"] = results.loc[ - adata.obs_names, "most_likely_hypothesis" - ] - adata.obs["cluster_feature"] = results.loc[adata.obs_names, "cluster_feature"] - adata.obs["negative_hypothesis_probability"] = results.loc[ - adata.obs_names, "negative_hypothesis_probability" - ] - adata.obs["singlet_hypothesis_probability"] = results.loc[ - adata.obs_names, "singlet_hypothesis_probability" - ] - adata.obs["doublet_hypothesis_probability"] = results.loc[ - adata.obs_names, "doublet_hypothesis_probability" - ] + adata.obs["most_likely_hypothesis"] = most_likely_hypothesis.astype(float) + adata.obs["cluster_feature"] = 0.0 if clusters is None else clusters + for i, hypothesis in enumerate(["negative", "singlet", "doublet"]): + adata.obs[f"{hypothesis}_hypothesis_probability"] = probs[:, i] - adata.obs["Classification"] = None - adata.obs.loc[adata.obs["most_likely_hypothesis"] == 2, "Classification"] = ( - "Doublet" - ) - adata.obs.loc[adata.obs["most_likely_hypothesis"] == 0, "Classification"] = ( - "Negative" - ) - all_sings = adata.obs["most_likely_hypothesis"] == 1 - singlet_sample_index = np.argmax( - adata.obs.loc[all_sings, cell_hashing_columns].values, axis=1 + classification = np.asarray( + np.array(cell_hashing_columns)[np.argmax(data, axis=1)], dtype=object ) - adata.obs.loc[all_sings, "Classification"] = adata.obs[ - cell_hashing_columns - ].columns[singlet_sample_index] + classification[most_likely_hypothesis == 0] = "Negative" + classification[most_likely_hypothesis == 2] = "Doublet" + adata.obs["Classification"] = classification return adata if not inplace else None diff --git a/src/scanpy/get/__init__.py b/src/scanpy/get/__init__.py index f4a3f6467f..652b2da946 100644 --- a/src/scanpy/get/__init__.py +++ b/src/scanpy/get/__init__.py @@ -7,6 +7,7 @@ _check_mask, _get_arr, _get_vec, + _get_vec_compat, _Rep, _set_obs_rep, obs_df, @@ -20,6 +21,7 @@ "_check_mask", "_get_arr", "_get_vec", + "_get_vec_compat", "_set_obs_rep", "aggregate", "obs_df", diff --git a/src/scanpy/get/_aggregated.py b/src/scanpy/get/_aggregated.py index 20ec02272c..4a17068bc5 100644 --- a/src/scanpy/get/_aggregated.py +++ b/src/scanpy/get/_aggregated.py @@ -23,7 +23,7 @@ mean_var_csr, mean_var_dense, ) -from .get import _check_mask, _get_arr, _get_vec, _refs_dim, _resolve_ref +from .get import _check_mask, _get_arr, _get_vec_compat, _refs_dim, _resolve_ref if TYPE_CHECKING: from collections.abc import Collection, Iterable @@ -363,7 +363,7 @@ def aggregate( acc = A.X data = _get_arr(adata, acc, dim=dim, layer=layer, obsm=obsm, varm=varm) - values = _get_vec(adata, by, dim=dim) + values = _get_vec_compat(adata, by, dim=dim) dim_df = pd.DataFrame({ ( ref diff --git a/src/scanpy/get/get.py b/src/scanpy/get/get.py index c3b62ce70d..80e38249b1 100644 --- a/src/scanpy/get/get.py +++ b/src/scanpy/get/get.py @@ -635,7 +635,7 @@ def _check_mask[M: NDArray[np.bool] | NDArray[np.floating] | pd.Series | None]( msg = f"Cannot use refererence for {desc} without providing anndata object as argument" raise ValueError(msg) try: - mask_array = np.asarray(_get_vec(data, mask, dim=dim)) + mask_array = np.asarray(_get_vec_compat(data, mask, dim=dim)) except KeyError: if isinstance(mask, AdRef): msg = ( @@ -725,26 +725,68 @@ def _fetch_vec(adata: AnnData, ref: AdRef | str, *, dim: Literal["obs", "var"]) @overload -def _get_vec( +def _get_vec_compat( adata: AnnData, ref: Collection[AdRef] | Collection[str], *, dim: Literal["obs", "var"] | None = None, ) -> list[Any]: ... @overload -def _get_vec( +def _get_vec_compat( adata: AnnData, ref: AdRef | str, *, dim: Literal["obs", "var"] | None = None ) -> Any: ... -def _get_vec( +def _get_vec_compat( adata: AnnData, ref: AdRef | str | Collection[AdRef] | Collection[str], *, dim: Literal["obs", "var"] | None = None, ) -> Any: - """Get the 1D array a `ref`erence points to, resolving plain strings first.""" + """Get the 1D array(s) one or more `ref`erences point to. + + Treats strings as `obs` columns instead of `anndata.acc` specs when the preset is v1. + """ if _collection_of(ref, (AdRef, str)): dim = _refs_dim(ref, dim=dim) return [_fetch_vec(adata, r, dim=dim) for r in ref] + if isinstance(ref, Collection) and not isinstance(ref, str): + msg = f"Expected a single ref or collection of refs, got {ref!r}" + raise TypeError(msg) + dim = _ref_dim(ref, dim=dim) return _fetch_vec(adata, ref, dim=dim) + + +@overload +def _get_vec( + adata: AnnData, + ref: Collection[AdRef] | Collection[str], + *, + dim: Literal["obs", "var"] | None = None, +) -> list[Any]: ... +@overload +def _get_vec( + adata: AnnData, ref: AdRef | str, *, dim: Literal["obs", "var"] | None = None +) -> Any: ... +def _get_vec( + adata: AnnData, + ref: AdRef | str | Collection[AdRef] | Collection[str], + *, + dim: Literal["obs", "var"] | None = None, +) -> Any: + """Get the 1D array(s) one or more `ref`erences point to, using `anndata.acc`.""" + from anndata.acc import A, AdRef + + if _collection_of(ref, (AdRef, str)): + refs = [A.resolve(r, vec=True) if isinstance(r, str) else r for r in ref] + _refs_dim(refs, dim=dim) + return [adata[r] for r in refs] + + if isinstance(ref, Collection) and not isinstance(ref, str): + msg = f"Expected a single ref or collection of refs, got {ref!r}" + raise TypeError(msg) + + if isinstance(ref, str): + ref = A.resolve(ref, vec=True) + _ref_dim(ref, dim=dim) + return adata[ref] diff --git a/src/scanpy/preprocessing/__init__.py b/src/scanpy/preprocessing/__init__.py index cb1aedb90c..3e1480c3df 100644 --- a/src/scanpy/preprocessing/__init__.py +++ b/src/scanpy/preprocessing/__init__.py @@ -6,6 +6,7 @@ from ._combat import combat from ._deprecated.sampling import subsample from ._harmony import harmony_integrate +from ._hashsolo import hashsolo from ._highly_variable_genes import highly_variable_genes from ._normalization import normalize_total from ._pca import pca @@ -30,6 +31,7 @@ "filter_cells", "filter_genes", "harmony_integrate", + "hashsolo", "highly_variable_genes", "log1p", "neighbors", diff --git a/src/scanpy/preprocessing/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py new file mode 100644 index 0000000000..6e91143be2 --- /dev/null +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -0,0 +1,347 @@ +"""A probabilistic cell hashing demultiplexing method. + +HashSolo generates a noise distribution and signal distribution +for each hashing barcode from empirically observed counts. +These distributions are updates from the global signal and noise barcode distributions, +which helps in the setting where not many cells are observed. +For a hashing barcode: + +Signal distributions + are estimated from samples where that hashing barcode has the highest count. + +Noise distributions + are estimated from samples where that hashing barcode is one the k-2 lowest barcodes, + where k is the number of barcodes. + +We test each of the following hypotheses in a bayesian fashion, +and select the most probable hypothesis. + +A doublet + should have its two highest barcode counts most likely + coming from a signal distribution for those barcodes. + +A singlet + should have its highest barcode from a signal distribution, + and its second highest barcode from a noise distribution. + +A negative two highest barcodes + should come from noise distributions. +""" + +from __future__ import annotations + +from itertools import product +from typing import TYPE_CHECKING, NamedTuple + +import numpy as np +import pandas as pd +from fast_array_utils.conv import to_dense +from scipy.stats import norm + +from .._compat import CSBase +from .._utils import check_nonnegative_integers +from .._utils._doctests import doctest_needs +from ..get.get import MultiAcc, _get_arr, _get_vec + +if TYPE_CHECKING: + from collections.abc import Collection + + from anndata import AnnData + from anndata.acc import AdRef + from numpy.typing import ArrayLike, NDArray + + +# smaller than what log-counts can resolve +# large enough that its inverse doesn’t get close to overflowing +_VAR_EPS = 1e-8 +"""Variance floor, so degenerate spreads (e.g. a one-cell cluster) stay finite.""" + + +class Gaussian(NamedTuple): + """A gaussian, its fields named after :func:`scipy.stats.norm.pdf`’s parameters.""" + + loc: float | np.floating + """Mean of the gaussian.""" + scale: float | np.floating + """Standard deviation of the gaussian.""" + + def log_pdf(self, counts: NDArray[np.floating]) -> NDArray[np.float64]: + """Log of the probability density at each of `counts`.""" + eps = 1e-15 # avoid log(0) + return np.log(norm.pdf(counts, loc=self.loc, scale=self.scale) + eps) + + def posterior(self, data: NDArray[np.floating]) -> Gaussian: + """Update this gaussian, used as a prior, with the observed 1-d `data`. + + See . + """ + n = len(data) + lam_o = 1 / max(self.scale**2, _VAR_EPS) + lam = 1 / max(np.var(data), _VAR_EPS) if n > 1 else lam_o + lam_n = lam_o + n * lam + mu_n = (np.mean(data) * n * lam + self.loc * lam_o) / lam_n if n else self.loc + return Gaussian(mu_n, np.sqrt((n + 1) / lam_n)) + + +def _calculate_log_likelihoods( + data: NDArray[np.integer], n_barcodes_noise: int +) -> NDArray[np.float64]: + """Calculate log likelihoods for each hypothesis, negative, singlet, doublet. + + Parameters + ---------- + data + cells by hashing counts matrix + n_barcodes_noise + number of barcodes to used to calculated noise distribution + + Returns + ------- + A 2d array of shape `(n_cells, 3)` with the log likelihood of each hypothesis. + + """ + # probabilites for negative, singlet, doublets + log_likelihoods = np.zeros((data.shape[0], 3)) + + n_barcodes = data.shape[1] + + # assume log normal + data: NDArray[np.floating] = np.log(data + 1) + # per cell, the barcode indices ordered by ascending count + data_arg = np.argsort(data, axis=1) + data_sort = np.sort(data, axis=1) + + # global signal and noise counts useful for when we have few cells + # barcodes with the highest number of counts are assumed to be a true signal + # barcodes with rank < k are considered to be noise + global_signal = data_sort[:, -1] + global_noise = data_sort[:, :n_barcodes_noise] + prior_sig = Gaussian(np.mean(global_signal), np.std(global_signal)) + prior_noise = Gaussian(np.mean(global_noise), np.std(global_noise)) + + # for each barcode get empirical noise and signal distribution parameterization, + # assuming lognormal, as an update from the global values + p_noise: list[Gaussian] = [] + p_sig: list[Gaussian] = [] + for x in range(n_barcodes): + is_noise = (data_arg[:, :n_barcodes_noise] == x).any(axis=1) + is_sig = data_arg[:, -1] == x + p_noise.append(prior_noise.posterior(data[is_noise, x])) + p_sig.append(prior_sig.posterior(data[is_sig, x])) + + # for each combination of noise and signal barcode calculate probiltiy of in silico and real cell hypotheses + for i_noise, i_sig in product(range(n_barcodes), repeat=2): + mask = (data_arg[:, -1] == i_sig) & (data_arg[:, -2] == i_noise) + if not mask.any(): + continue + + # the distributions the 2nd-highest and highest barcode’s counts + # are assumed to come from under each hypothesis + hypotheses = [ + (p_noise[i_noise], p_noise[i_noise]), # negative: neither barcode is signal + (p_noise[i_noise], p_sig[i_sig]), # singlet: only the highest barcode is + (p_sig[i_noise], p_sig[i_sig]), # doublet: both are + ] + for i_prob, (p_2nd, p_top) in enumerate(hypotheses): + log_likelihoods[mask, i_prob] = p_2nd.log_pdf( + data[mask, i_noise] + ) + p_top.log_pdf(data[mask, i_sig]) + return log_likelihoods + + +def _calculate_bayes_rule( + data: NDArray[np.integer], priors: ArrayLike, n_barcodes_noise: int +) -> NDArray[np.float64]: + """Calculate the posterior probability of each hypothesis from log likelihoods. + + Parameters + ---------- + data + cells by hashing counts matrix + priors + prior for each hypothesis, in the order `[negative, singlet, doublet]` + n_barcodes_noise + number of barcodes to used to calculated noise distribution + + Returns + ------- + A 2d array of shape `(n_cells, 3)` with the probability of each hypothesis. + + """ + log_likelihoods = _calculate_log_likelihoods(data, n_barcodes_noise) + likelihoods = np.exp(log_likelihoods) * np.asarray(priors) + return likelihoods / likelihoods.sum(axis=1)[:, None] + + +def _hashsolo( + data: NDArray[np.integer], + *, + priors: ArrayLike, + clusters: NDArray | None, + n_barcodes_noise: int | None, +) -> NDArray[np.float64]: + """Validate counts and run the bayes rule, optionally per cluster.""" + if not check_nonnegative_integers(data): + msg = "Cell hashing counts must be non-negative integers" + raise ValueError(msg) + n_barcodes = data.shape[1] + if n_barcodes_noise is None: + n_barcodes_noise = n_barcodes - 2 + if not 1 <= n_barcodes_noise < n_barcodes: + msg = ( + f"The number of noise barcodes ({n_barcodes_noise}) must be at least 1 and smaller " + f"than the number of hashing barcodes ({n_barcodes}). " + f"Pass at least 3 `hashes` or set `n_noise_barcodes` explicitly." + ) + raise ValueError(msg) + + if clusters is None: + return _calculate_bayes_rule(data, priors, n_barcodes_noise) + + probs = np.zeros((data.shape[0], 3)) + # `factorize` then `np.unique` gives cells with no cluster the code -1 + codes = pd.factorize(clusters)[0] + for code in np.unique(codes): + mask = codes == code + probs[mask] = _calculate_bayes_rule(data[mask], priors, n_barcodes_noise) + return probs + + +def _ref_name(ref: AdRef | str) -> str: + """Get the name of the barcode a `ref` points to, e.g. `A.X[:, "Hash1"]` → `Hash1`.""" + if isinstance(ref, str): + return ref + idx = ref.idx[-1] if isinstance(ref.idx, tuple) else ref.idx + return str(idx) + + +def _get_hashes( + adata: AnnData, hashes: MultiAcc | Collection[AdRef] | Collection[str] +) -> tuple[NDArray, NDArray[np.str_]]: + """Get the hashing count matrix and each hashing barcode’s name.""" + if not isinstance(hashes, MultiAcc): + data = np.column_stack(_get_vec(adata, hashes, dim="obs")) + return data, np.array([_ref_name(ref) for ref in hashes]) + # a `MultiAcc` such as `A.obsm["hto"]` refers to all of its columns + data = _get_arr(adata, hashes, dim="obs") + if isinstance(data, pd.DataFrame): + return data.to_numpy(), data.columns.to_numpy(dtype=str) + if isinstance(data, CSBase): + data = to_dense(data) + return data, np.arange(data.shape[1]).astype(str) + + +@doctest_needs("anndata_acc") +def hashsolo( + adata: AnnData, + hashes: MultiAcc | Collection[AdRef] | Collection[str], + *, + priors: tuple[float, float, float] = (0.01, 0.8, 0.19), + pre_existing_clusters: AdRef | str | None = None, + n_noise_barcodes: int | None = None, + key_added: str = "hashsolo", + copy: bool = False, +) -> AnnData | None: + r"""Probabilistic demultiplexing of cell hashing data using HashSolo :cite:p:`Bernstein2020`. + + .. array-support:: pp.hashsolo + + Parameters + ---------- + adata + The (annotated) data matrix of shape `n_obs` × `n_vars`. + Rows correspond to cells and columns to genes. + hashes + References to the vectors holding the cell hashing counts, + e.g. `A.obs[["Hash1", "Hash2"]]` for columns in :attr:`~anndata.AnnData.obs` + or `A.X[:, ["Hash1", "Hash2"]]` for hashing barcodes among the + :attr:`~anndata.AnnData.var_names`. + A :class:`~anndata.acc.MultiAcc` such as `A.obsm["hto"]` refers to + *all* columns of what it points to, + named after the :class:`~pandas.DataFrame`’s columns or, + for a plain array, its column positions. + priors + Prior probabilities of each hypothesis, in + the order `[negative, singlet, doublet]`. The default is set to + `[0.01, 0.8, 0.19]` assuming barcode counts are from cells that + have passed QC in the transcriptome space, e.g. UMI counts, pct + mito reads, etc. + pre_existing_clusters + Reference to a vector of pre-existing cluster assignments\ [#ref]_ + (e.g. Leiden clusters or cell types, but not batch assignments). + If provided, demultiplexing is performed separately for each cluster. + n_noise_barcodes + The number of barcodes used to create the noise distribution. + Defaults to `len(hashes) - 2`. + key_added + Key under which to add the demultiplexing results. + copy + Whether to modify a copy of `adata` instead of `adata` itself. + + Returns + ------- + Returns `None` if `copy=False`, else the modified `adata`. + Sets the following fields: + + `adata.obs[key_added]` : :class:`~pandas.Categorical` (shape `(n_obs,)`) + Classification of each cell: the name of one of the `hashes`, + `"Negative"`, or `"Doublet"`. + `adata.obsm[key_added]` : :class:`~pandas.DataFrame` (shape `(n_obs, 3)`) + Probability of the `negative`, `singlet`, and `doublet` hypothesis. + + Examples + -------- + Simulate 300 cells, each carrying one of 3 hashtag oligos: + + >>> import numpy as np + >>> import scanpy as sc + >>> from anndata import AnnData + >>> from anndata.acc import A + >>> + >>> rng = np.random.default_rng(0) + >>> hto = rng.poisson(20, size=(300, 3)) # ambient background + >>> hto[np.arange(300), np.arange(300) % 3] = rng.poisson(1000, 300) # signal + >>> adata = AnnData(rng.poisson(1, (300, 5)).astype("f4"), obsm=dict(hto=hto)) + + A :class:`~anndata.acc.MultiAcc` demultiplexes using every column it points to. + A plain array has no column names, so the barcodes are named after their positions: + + >>> sc.pp.hashsolo(adata, A.obsm["hto"]) + >>> adata.obs["hashsolo"].cat.categories.astype("string") + Index(['0', '1', '2'], dtype='string') + >>> adata.obs["hashsolo"].value_counts().tolist() + [100, 100, 100] + + The same counts as :attr:`~anndata.AnnData.obs` columns + (as in a Cell Ranger run’s “Multiplexing Capture” features), + referenced one by one: + + >>> adata.obs[["Hash1", "Hash2", "Hash3"]] = hto + >>> sc.pp.hashsolo(adata, A.obs[["Hash1", "Hash2", "Hash3"]], key_added="hs2") + >>> adata.obs["hs2"].cat.categories.astype("string") + Index(['Hash1', 'Hash2', 'Hash3'], dtype='string') + + """ + adata = adata.copy() if copy else adata + data, names = _get_hashes(adata, hashes) + clusters = ( + None + if pre_existing_clusters is None + else np.asarray(_get_vec(adata, pre_existing_clusters, dim="obs")) + ) + probs = _hashsolo( + data, priors=priors, clusters=clusters, n_barcodes_noise=n_noise_barcodes + ) + + most_likely_hypothesis = np.argmax(probs, axis=1) + classification = pd.Series( + names[np.argmax(data, axis=1)], index=adata.obs_names, dtype="string" + ) + classification[most_likely_hypothesis == 0] = "Negative" + classification[most_likely_hypothesis == 2] = "Doublet" + + adata.obs[key_added] = classification.astype("category") + adata.obsm[key_added] = pd.DataFrame( + probs, columns=["negative", "singlet", "doublet"], index=adata.obs_names + ) + return adata if copy else None diff --git a/tests/external/test_hashsolo.py b/tests/external/test_hashsolo.py deleted file mode 100644 index 3779fddb9a..0000000000 --- a/tests/external/test_hashsolo.py +++ /dev/null @@ -1,42 +0,0 @@ -from __future__ import annotations - -import warnings - -import numpy as np -import pandas as pd -from anndata import AnnData, ImplicitModificationWarning - -import scanpy.external as sce - - -def test_cell_demultiplexing(): - from scipy import stats - - rng = np.random.default_rng() - - signal = stats.poisson.rvs(1000, 1, 990, random_state=rng) - doublet_signal = stats.poisson.rvs(1000, 1, 10, random_state=rng) - x = np.reshape(stats.poisson.rvs(500, 1, 10000, random_state=rng), (1000, 10)) - for idx, signal_count in enumerate(signal): - col_pos = idx % 10 - x[idx, col_pos] = signal_count - - for idx, signal_count in enumerate(doublet_signal): - col_pos = (idx % 10) - 1 - x[idx, col_pos] = signal_count - - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", category=ImplicitModificationWarning) - test_data = AnnData(rng.integers(0, 100, size=x.shape), obs=pd.DataFrame(x)) - sce.pp.hashsolo(test_data, test_data.obs.columns) - - doublets = ["Doublet"] * 10 - classes = np.repeat(np.arange(10), 98).reshape(98, 10, order="F").ravel().tolist() - negatives = ["Negative"] * 10 - expected = pd.array(doublets + classes + negatives, dtype="string") - classification = test_data.obs["Classification"].array.astype("string") - # This is a bit flaky, so allow some mismatches: - # (Series() because of https://github.com/pandas-dev/pandas/issues/63458) - if pd.Series(expected != classification).sum() > 3: - # Compare lists for better error message - assert classification.tolist() == expected.tolist() diff --git a/tests/test_hashsolo.py b/tests/test_hashsolo.py new file mode 100644 index 0000000000..9501f32573 --- /dev/null +++ b/tests/test_hashsolo.py @@ -0,0 +1,184 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING, NamedTuple + +import anndata +import numpy as np +import pandas as pd +import pytest +from anndata import AnnData +from scipy import sparse, stats + +import scanpy as sc +import scanpy.external as sce +from testing.scanpy._pytest.marks import needs + +if TYPE_CHECKING: + from collections.abc import Sequence + + from anndata.acc import AdRef, MultiAcc + +if TYPE_CHECKING or hasattr(anndata, "acc"): + from anndata.acc import A + + +HASHES = [f"Hash{i}" for i in range(10)] + + +class Hashed(NamedTuple): + """An `AnnData` with hashing counts, how to reference them, and their names.""" + + adata: AnnData + refs: MultiAcc | Sequence[AdRef | str] + names: Sequence[str] + + +@pytest.fixture +def counts() -> pd.DataFrame: + """Hashing counts: cell `i` is a singlet for barcode `i % 10`. + + Except for the first 10 cells (doublets) and the last 10 (negatives). + """ + rng = np.random.default_rng(0) + signal = stats.poisson.rvs(1000, 1, 990, random_state=rng) + doublet_signal = stats.poisson.rvs(1000, 1, 10, random_state=rng) + x = np.reshape(stats.poisson.rvs(500, 1, 10000, random_state=rng), (1000, 10)) + for idx, signal_count in enumerate(signal): + x[idx, idx % 10] = signal_count + for idx, signal_count in enumerate(doublet_signal): + x[idx, (idx % 10) - 1] = signal_count + return pd.DataFrame(x, columns=HASHES, index=[f"cell{i}" for i in range(len(x))]) + + +@pytest.fixture +def adata(counts: pd.DataFrame) -> AnnData: + """Hashing counts in `.obs`, with unrelated expression data in `X`.""" + rng = np.random.default_rng(0) + return AnnData(rng.integers(0, 100, size=counts.shape), obs=counts.copy()) + + +@pytest.fixture(params=["obs", "X", "X-sparse", "obsm-df", "obsm-array"]) +def hashed( + request: pytest.FixtureRequest, adata: AnnData, counts: pd.DataFrame +) -> Hashed: + """Build an `AnnData` holding the hashing counts, plus references to them.""" + match request.param: + case "obs": + return Hashed(adata, A.obs[HASHES], HASHES) + case "X" | "X-sparse": + x = counts.to_numpy() + if request.param == "X-sparse": + x = sparse.csr_matrix(x) # noqa: TID251 + adata = AnnData(x, var=pd.DataFrame(index=HASHES)) + return Hashed(adata, A.X[:, HASHES], HASHES) + case "obsm-df" | "obsm-array": # a `MultiAcc` means “all of its columns” + df = request.param == "obsm-df" + adata.obsm["hto"] = counts.copy() if df else counts.to_numpy() + # a plain array has no column names, so they fall back to positions + names = HASHES if df else [str(i) for i in range(len(HASHES))] + return Hashed(adata, A.obsm["hto"], names) + case _: + pytest.fail(f"Unknown param {request.param!r}") + + +@pytest.mark.parametrize( + ("n_hashes", "n_noise"), + [ + pytest.param(1, None, id="1-hash"), + pytest.param(2, None, id="2-hashes-default"), + pytest.param(len(HASHES), 0, id="explicit-0"), + ], +) +@needs.anndata_acc +def test_too_few_noise_barcodes( + adata: AnnData, n_hashes: int, n_noise: int | None +) -> None: + """Without a noise barcode, the noise distribution is undefined.""" + with pytest.raises(ValueError, match=r"noise barcodes?"): + sc.pp.hashsolo(adata, A.obs[HASHES[:n_hashes]], n_noise_barcodes=n_noise) + + +@needs.anndata_acc +def test_cell_demultiplexing(hashed: Hashed) -> None: + adata, refs, names = hashed + sc.pp.hashsolo(adata, refs) + + expected = pd.array( + ["Doublet"] * 10 + + np.repeat(names, 98).reshape(98, 10, order="F").ravel().tolist() + + ["Negative"] * 10, + dtype="string", + ) + classification = adata.obs["hashsolo"].array.astype("string") + # This is a bit flaky, so allow some mismatches: + # (Series() because of https://github.com/pandas-dev/pandas/issues/63458) + if pd.Series(expected != classification).sum() > 3: + # Compare lists for better error message + assert classification.tolist() == expected.tolist() + + probs = adata.obsm["hashsolo"] + assert isinstance(probs, pd.DataFrame) + assert list(probs.columns) == ["negative", "singlet", "doublet"] + np.testing.assert_allclose(probs.to_numpy().sum(axis=1), 1) + + +@needs.anndata_acc +def test_copy(adata: AnnData) -> None: + copied = sc.pp.hashsolo(adata, A.obs[HASHES], copy=True) + assert "hashsolo" not in adata.obs + assert "hashsolo" in copied.obs + + +@needs.anndata_acc +def test_pre_existing_clusters(adata: AnnData) -> None: + """Clustered demultiplexing equals demultiplexing each cluster on its own.""" + adata.obs["cl"] = np.where(np.arange(adata.n_obs) % 2, "a", "b") + sc.pp.hashsolo(adata, A.obs[HASHES], pre_existing_clusters=A.obs["cl"]) + + np.testing.assert_allclose(adata.obsm["hashsolo"].to_numpy().sum(axis=1), 1) + for cluster in ("a", "b"): + sub = adata[adata.obs["cl"] == cluster].copy() + sc.pp.hashsolo(sub, A.obs[HASHES]) + np.testing.assert_allclose( + sub.obsm["hashsolo"].to_numpy(), + adata.obsm["hashsolo"].loc[sub.obs_names].to_numpy(), + ) + + +@needs.anndata_acc +def test_pre_existing_clusters_missing_label(adata: AnnData) -> None: + """Unlabeled cells must not silently turn into confident negatives.""" + adata.obs["cl"] = pd.Categorical(["a"] * (adata.n_obs - 10) + [None] * 10) + sc.pp.hashsolo(adata, A.obs[HASHES], pre_existing_clusters=A.obs["cl"]) + np.testing.assert_allclose(adata.obsm["hashsolo"].to_numpy().sum(axis=1), 1) + + +@needs.anndata_acc +def test_pre_existing_clusters_singleton(adata: AnnData) -> None: + """A one-cell cluster has zero variance; it must not yield NaN probabilities.""" + adata.obs["cl"] = ["a"] * (adata.n_obs - 1) + ["b"] + sc.pp.hashsolo(adata, A.obs[HASHES], pre_existing_clusters=A.obs["cl"]) + assert np.isfinite(adata.obsm["hashsolo"].to_numpy()).all() + + +def test_legacy_api_needs_no_acc(subtests: pytest.Subtests, adata: AnnData) -> None: + with pytest.warns(FutureWarning, match=r"scanpy\.pp\.hashsolo"): + sce.pp.hashsolo(adata, HASHES) + + +@needs.anndata_acc +def test_legacy_api_matches(subtests: pytest.Subtests, adata: AnnData) -> None: + with pytest.warns(FutureWarning, match=r"scanpy\.pp\.hashsolo"): + sce.pp.hashsolo(adata, HASHES) + + expected_new = adata.copy() + sc.pp.hashsolo(expected_new, A.obs[HASHES]) + pd.testing.assert_series_equal( + adata.obs["Classification"].astype("string"), + expected_new.obs["hashsolo"].astype("string"), + check_names=False, + ) + np.testing.assert_array_equal( + adata.obs["most_likely_hypothesis"], + np.argmax(expected_new.obsm["hashsolo"].to_numpy(), axis=1), + )