From 34060bb61451d0bdf4da8a1a2528246d4bfef7c5 Mon Sep 17 00:00:00 2001 From: "Philipp A." Date: Tue, 18 Aug 2026 11:28:25 +0200 Subject: [PATCH 01/10] feat: vendor hashsolo --- docs/api/preprocessing.md | 10 + docs/conf.py | 1 + docs/external/preprocessing.md | 1 + docs/release-notes/4302.feat.md | 3 + src/scanpy/external/pp/_hashsolo.py | 393 +--------------------- src/scanpy/preprocessing/__init__.py | 2 + src/scanpy/preprocessing/_hashsolo.py | 449 ++++++++++++++++++++++++++ tests/external/test_hashsolo.py | 42 --- tests/test_hashsolo.py | 105 ++++++ 9 files changed, 588 insertions(+), 418 deletions(-) create mode 100644 docs/release-notes/4302.feat.md create mode 100644 src/scanpy/preprocessing/_hashsolo.py delete mode 100644 tests/external/test_hashsolo.py create mode 100644 tests/test_hashsolo.py 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/4302.feat.md b/docs/release-notes/4302.feat.md new file mode 100644 index 0000000000..b80405a421 --- /dev/null +++ b/docs/release-notes/4302.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` *or* in the data matrix (e.g. `A.X[:, ["Hash1", "Hash2"]]`), +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..b303d706a5 100644 --- a/src/scanpy/external/pp/_hashsolo.py +++ b/src/scanpy/external/pp/_hashsolo.py @@ -1,296 +1,24 @@ -"""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 _legacy_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 `scanpy.pp.hashsolo` instead.")) @doctest_skipif(reason="Illustrative but not runnable doctest code") def hashsolo( adata: AnnData, @@ -303,8 +31,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,9 +81,9 @@ def hashsolo( Examples -------- - >>> import anndata + >>> import scanpy as sc >>> import scanpy.external as sce - >>> adata = anndata.read_h5ad("data.h5ad") + >>> adata = sc.io.read_h5ad("data.h5ad") >>> sce.pp.hashsolo(adata, ["Hash1", "Hash2", "Hash3"]) >>> adata.obs.head() @@ -362,99 +91,11 @@ def hashsolo( 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 - 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, - ) - 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] - - 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["Classification"] = None - adata.obs.loc[adata.obs["most_likely_hypothesis"] == 2, "Classification"] = ( - "Doublet" + return _legacy_hashsolo( + adata, + cell_hashing_columns, + priors=priors, + pre_existing_clusters=pre_existing_clusters, + number_of_noise_barcodes=number_of_noise_barcodes, + inplace=inplace, ) - 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 - ) - adata.obs.loc[all_sings, "Classification"] = adata.obs[ - cell_hashing_columns - ].columns[singlet_sample_index] - - return adata if not inplace else None 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..e2ed974942 --- /dev/null +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -0,0 +1,449 @@ +"""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 + +import numpy as np +import pandas as pd +from scipy.stats import norm + +from .._utils import check_nonnegative_integers +from .._utils._doctests import doctest_skipif +from ..get.get import _get_vec + +if TYPE_CHECKING: + from collections.abc import Collection, Sequence + + from anndata import AnnData + from anndata.acc import AdRef + 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], + 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 | None +) -> 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]` + number_of_noise_barcodes + 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, number_of_noise_barcodes) + likelihoods = np.exp(log_likelihoods) * np.asarray(priors) + return likelihoods / likelihoods.sum(axis=1)[:, None] + + +def _hashsolo( + data: NDArray, + *, + priors: ArrayLike, + clusters: NDArray | None, + number_of_noise_barcodes: 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" + raise ValueError(msg) + if number_of_noise_barcodes is not None and number_of_noise_barcodes >= ( + n_barcodes := data.shape[1] + ): + msg = ( + f"`number_of_noise_barcodes` ({number_of_noise_barcodes}) must be smaller " + f"than the number of hashing barcodes ({n_barcodes})." + ) + raise ValueError(msg) + + if clusters is None: + return _calculate_bayes_rule(data, priors, number_of_noise_barcodes) + + probs = np.zeros((data.shape[0], 3)) + for cluster in pd.unique(clusters): + mask = clusters == cluster + probs[mask] = _calculate_bayes_rule( + data[mask], priors, number_of_noise_barcodes + ) + 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) + + +@doctest_skipif(reason="Illustrative but not runnable doctest code") +def hashsolo( + adata: AnnData, + cell_hashing_columns: Collection[AdRef | str], + *, + priors: tuple[float, float, float] = (0.01, 0.8, 0.19), + pre_existing_clusters: AdRef | str | None = None, + number_of_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. + cell_hashing_columns + References to the vectors holding the cell hashing counts\ [#ref]_, + 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`. + 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. + number_of_noise_barcodes + The number of barcodes used to create the noise distribution. + Defaults to `len(cell_hashing_columns) - 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 `cell_hashing_columns`, + `"Negative"`, or `"Doublet"`. + `adata.obsm[key_added]` : :class:`~pandas.DataFrame` (shape `(n_obs, 3)`) + Probability of the `negative`, `singlet`, and `doublet` hypothesis. + + Examples + -------- + Hashing counts stored as :attr:`~anndata.AnnData.var_names` + (e.g. the “Multiplexing Capture” features of a Cell Ranger run): + + >>> import scanpy as sc + >>> from anndata.acc import A + >>> adata = sc.io.read_h5ad("data.h5ad") + >>> sc.pp.hashsolo(adata, A.X[:, ["Hash1", "Hash2", "Hash3"]]) + >>> adata.obs["hashsolo"].value_counts() + + Hashing counts stored as :attr:`~anndata.AnnData.obs` columns: + + >>> sc.pp.hashsolo(adata, A.obs[["Hash1", "Hash2", "Hash3"]]) + + .. [#ref] If :attr:`scanpy.settings.preset` is :attr:`~scanpy.Preset.ScanpyV2Preview`, + :class:`str`\ s are :meth:`anndata.acc.AdAcc.resolve`\ d to :class:`~anndata.acc.AdRef`\ s, + otherwise interpreted as :attr:`anndata.AnnData.obs` columns. + + """ + adata = adata.copy() if copy else adata + refs = list(cell_hashing_columns) + data = np.column_stack(_get_vec(adata, refs, dim="obs")) + 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, + number_of_noise_barcodes=number_of_noise_barcodes, + ) + + names = np.array([_ref_name(ref) for ref in refs]) + 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 + + +def _legacy_hashsolo( + adata: AnnData, + cell_hashing_columns: Sequence[str], + *, + priors: tuple[float, float, float], + pre_existing_clusters: str | None, + number_of_noise_barcodes: int | None, + inplace: bool, +) -> AnnData | None: + """Implement the pre-1.13 `scanpy.external.pp.hashsolo` API.""" + adata = adata if inplace else adata.copy() + cell_hashing_columns = list(cell_hashing_columns) + data = adata.obs[cell_hashing_columns].to_numpy() + clusters = ( + None + if pre_existing_clusters is None + else adata.obs[pre_existing_clusters].to_numpy() + ) + probs = _hashsolo( + data, + priors=priors, + clusters=clusters, + number_of_noise_barcodes=number_of_noise_barcodes, + ) + most_likely_hypothesis = np.argmax(probs, axis=1) + + 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] + + classification = np.asarray( + np.array(cell_hashing_columns)[np.argmax(data, axis=1)], dtype=object + ) + 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/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..d3e758809d --- /dev/null +++ b/tests/test_hashsolo.py @@ -0,0 +1,105 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +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 + +if TYPE_CHECKING: + from collections.abc import Callable + +HASHES = [f"Hash{i}" for i in range(10)] + + +@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))]) + + +def _in_obs(counts: pd.DataFrame) -> tuple[AnnData, list]: + from anndata.acc import A + + rng = np.random.default_rng(0) + adata = AnnData(rng.integers(0, 100, size=counts.shape), obs=counts.copy()) + return adata, A.obs[HASHES] + + +def _in_x(counts: pd.DataFrame) -> tuple[AnnData, list]: + from anndata.acc import A + + adata = AnnData(counts.to_numpy(), var=pd.DataFrame(index=HASHES)) + return adata, A.X[:, HASHES] + + +def _in_x_sparse(counts: pd.DataFrame) -> tuple[AnnData, list]: + adata, refs = _in_x(counts) + adata.X = sparse.csr_matrix(adata.X) # noqa: TID251 + return adata, refs + + +@pytest.mark.parametrize( + "make", [_in_obs, _in_x, _in_x_sparse], ids=["obs", "X", "X-sparse"] +) +def test_cell_demultiplexing(counts: pd.DataFrame, make: Callable) -> None: + adata, refs = make(counts) + sc.pp.hashsolo(adata, refs) + + expected = pd.array( + ["Doublet"] * 10 + + np.repeat(HASHES, 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 list(probs.columns) == ["negative", "singlet", "doublet"] + np.testing.assert_allclose(probs.to_numpy().sum(axis=1), 1) + + +def test_copy(counts: pd.DataFrame) -> None: + adata, refs = _in_obs(counts) + copied = sc.pp.hashsolo(adata, refs, copy=True) + assert "hashsolo" not in adata.obs + assert "hashsolo" in copied.obs + + +def test_legacy_api(counts: pd.DataFrame) -> None: + adata, _ = _in_obs(counts) + with pytest.warns(FutureWarning, match=r"scanpy\.pp\.hashsolo"): + sce.pp.hashsolo(adata, HASHES) + + expected_new = adata.copy() + sc.pp.hashsolo(expected_new, 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), + ) From 6308164af820d79e3d558720c996a45537251573 Mon Sep 17 00:00:00 2001 From: "Philipp A." Date: Tue, 18 Aug 2026 11:30:39 +0200 Subject: [PATCH 02/10] fix: PR number --- docs/release-notes/{4302.feat.md => 4303.feat.md} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/release-notes/{4302.feat.md => 4303.feat.md} (100%) diff --git a/docs/release-notes/4302.feat.md b/docs/release-notes/4303.feat.md similarity index 100% rename from docs/release-notes/4302.feat.md rename to docs/release-notes/4303.feat.md From 40b277e438bf4c20494256088574035c839b10bb Mon Sep 17 00:00:00 2001 From: "Philipp A." Date: Tue, 18 Aug 2026 11:37:34 +0200 Subject: [PATCH 03/10] test both --- tests/test_hashsolo.py | 72 ++++++++++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 28 deletions(-) diff --git a/tests/test_hashsolo.py b/tests/test_hashsolo.py index d3e758809d..dbc0015a38 100644 --- a/tests/test_hashsolo.py +++ b/tests/test_hashsolo.py @@ -10,9 +10,12 @@ import scanpy as sc import scanpy.external as sce +from testing.scanpy._pytest.marks import needs if TYPE_CHECKING: - from collections.abc import Callable + from collections.abc import Sequence + + from anndata.acc import AdRef HASHES = [f"Hash{i}" for i in range(10)] @@ -34,32 +37,45 @@ def counts() -> pd.DataFrame: return pd.DataFrame(x, columns=HASHES, index=[f"cell{i}" for i in range(len(x))]) -def _in_obs(counts: pd.DataFrame) -> tuple[AnnData, list]: - from anndata.acc import A - +@pytest.fixture +def adata(counts: pd.DataFrame) -> AnnData: + """Hashing counts in `.obs`, with unrelated expression data in `X`.""" rng = np.random.default_rng(0) - adata = AnnData(rng.integers(0, 100, size=counts.shape), obs=counts.copy()) - return adata, A.obs[HASHES] - - -def _in_x(counts: pd.DataFrame) -> tuple[AnnData, list]: - from anndata.acc import A - - adata = AnnData(counts.to_numpy(), var=pd.DataFrame(index=HASHES)) - return adata, A.X[:, HASHES] - - -def _in_x_sparse(counts: pd.DataFrame) -> tuple[AnnData, list]: - adata, refs = _in_x(counts) - adata.X = sparse.csr_matrix(adata.X) # noqa: TID251 - return adata, refs + return AnnData(rng.integers(0, 100, size=counts.shape), obs=counts.copy()) -@pytest.mark.parametrize( - "make", [_in_obs, _in_x, _in_x_sparse], ids=["obs", "X", "X-sparse"] +@pytest.fixture( + params=[ + "obs-str", + pytest.param("obs", marks=needs.anndata_acc), + pytest.param("X", marks=needs.anndata_acc), + pytest.param("X-sparse", marks=needs.anndata_acc), + ] ) -def test_cell_demultiplexing(counts: pd.DataFrame, make: Callable) -> None: - adata, refs = make(counts) +def hashed( + request: pytest.FixtureRequest, adata: AnnData, counts: pd.DataFrame +) -> tuple[AnnData, Sequence[AdRef | str]]: + """Build an `AnnData` holding the hashing counts, plus references to them.""" + match request.param: + case "obs-str": # plain `.obs` column names work without `anndata.acc` + return adata, HASHES + case "obs": + from anndata.acc import A + + return adata, A.obs[HASHES] + case "X" | "X-sparse": + from anndata.acc import A + + x = counts.to_numpy() + if request.param == "X-sparse": + x = sparse.csr_matrix(x) # noqa: TID251 + return AnnData(x, var=pd.DataFrame(index=HASHES)), A.X[:, HASHES] + case _: + pytest.fail(f"Unknown param {request.param!r}") + + +def test_cell_demultiplexing(hashed: tuple[AnnData, Sequence[AdRef | str]]) -> None: + adata, refs = hashed sc.pp.hashsolo(adata, refs) expected = pd.array( @@ -76,19 +92,19 @@ def test_cell_demultiplexing(counts: pd.DataFrame, make: Callable) -> None: 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) -def test_copy(counts: pd.DataFrame) -> None: - adata, refs = _in_obs(counts) - copied = sc.pp.hashsolo(adata, refs, copy=True) +def test_copy(adata: AnnData) -> None: + copied = sc.pp.hashsolo(adata, HASHES, copy=True) assert "hashsolo" not in adata.obs assert "hashsolo" in copied.obs -def test_legacy_api(counts: pd.DataFrame) -> None: - adata, _ = _in_obs(counts) +def test_legacy_api(adata: AnnData) -> None: + """The deprecated `sce.pp.hashsolo` matches `sc.pp.hashsolo` and needs no `anndata.acc`.""" with pytest.warns(FutureWarning, match=r"scanpy\.pp\.hashsolo"): sce.pp.hashsolo(adata, HASHES) From 912c488a8acb6d9fc0fd5ca2922673b24b677544 Mon Sep 17 00:00:00 2001 From: "Philipp A." Date: Tue, 18 Aug 2026 12:12:51 +0200 Subject: [PATCH 04/10] some cleanup --- src/scanpy/external/pp/_hashsolo.py | 3 +- src/scanpy/preprocessing/_hashsolo.py | 108 ++++++++++---------------- 2 files changed, 41 insertions(+), 70 deletions(-) diff --git a/src/scanpy/external/pp/_hashsolo.py b/src/scanpy/external/pp/_hashsolo.py index b303d706a5..e88a4efd82 100644 --- a/src/scanpy/external/pp/_hashsolo.py +++ b/src/scanpy/external/pp/_hashsolo.py @@ -82,9 +82,8 @@ def hashsolo( Examples -------- >>> import scanpy as sc - >>> import scanpy.external as sce >>> adata = sc.io.read_h5ad("data.h5ad") - >>> sce.pp.hashsolo(adata, ["Hash1", "Hash2", "Hash3"]) + >>> sc.external.pp.hashsolo(adata, ["Hash1", "Hash2", "Hash3"]) >>> adata.obs.head() """ diff --git a/src/scanpy/preprocessing/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py index e2ed974942..fcf7307129 100644 --- a/src/scanpy/preprocessing/_hashsolo.py +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -31,7 +31,7 @@ from __future__ import annotations from itertools import product -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, NamedTuple import numpy as np import pandas as pd @@ -49,8 +49,17 @@ from numpy.typing import ArrayLike, NDArray -def _calculate_log_likelihoods( # noqa: PLR0915 - data: np.ndarray, number_of_noise_barcodes: int +class GaussianParams(NamedTuple): + """Parameters of a gaussian, named after :func:`scipy.stats.norm.pdf`’s parameters.""" + + loc: float + """Mean of the gaussian.""" + scale: float + """Standard deviation of the gaussian.""" + + +def _calculate_log_likelihoods( + data: np.ndarray, number_of_noise_barcodes: int | None ) -> tuple[NDArray[np.float64], NDArray[np.float64], dict[int, str]]: """Calculate log likelihoods for each hypothesis, negative, singlet, doublet. @@ -70,9 +79,7 @@ def _calculate_log_likelihoods( # noqa: PLR0915 """ - def gaussian_updates( - data: np.ndarray, mu_o: float, std_o: float - ) -> tuple[float, float]: + def gaussian_updates(data: np.ndarray, mu_o: float, std_o: float) -> GaussianParams: """Update parameters of your gaussian. See . @@ -88,10 +95,7 @@ def gaussian_updates( Returns ------- - mean - of gaussian - std - of gaussian + The updated parameters of the gaussian. """ lam_o = 1 / (std_o**2) @@ -101,7 +105,7 @@ def gaussian_updates( 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) + return GaussianParams(mu_n, (1 / (lam_n / (n + 1))) ** (1 / 2)) eps = 1e-15 # probabilites for negative, singlet, doublets @@ -136,8 +140,8 @@ def gaussian_updates( np.std(global_noise_counts), ) - noise_params_dict = {} - signal_params_dict = {} + noise_params_dict: dict[int, GaussianParams] = {} + signal_params_dict: dict[int, GaussianParams] = {} # for each barcode get empirical noise and signal distribution parameterization for x in np.arange(num_of_barcodes): @@ -152,14 +156,12 @@ def gaussian_updates( signal_counts = sample_barcodes[sample_barcodes_signal_idx] # get parameters of distribution, assuming lognormal do update from global values - noise_param = gaussian_updates( + noise_params_dict[x] = gaussian_updates( noise_counts, global_mu_noise_o, global_sigma_noise_o ) - signal_param = gaussian_updates( + signal_params_dict[x] = 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 @@ -184,54 +186,24 @@ def gaussian_updates( # 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], - 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] + data_noise = data_subset[:, noise_sample_idx] + data_signal = data_subset[:, signal_sample_idx] + # the distributions the 2nd-highest and highest barcode’s counts + # are assumed to come from under each hypothesis + hypotheses = [ + (noise_params, noise_params), # negative: neither barcode is signal + (noise_params, signal_params), # singlet: only the highest barcode is + (signal_params_dict[noise_sample_idx], signal_params), # doublet: both are + ] # each cell and each hypothesis probability - for prob_idx, log_prob in enumerate(log_probs_list): + for prob_idx, log_prob in enumerate( + ( + np.log(norm.pdf(data_noise, loc=d_noise.loc, scale=d_noise.scale) + eps) + + np.log(norm.pdf(data_signal, loc=d_sig.loc, scale=d_sig.scale) + eps) + for d_noise, d_sig in hypotheses + ) + ): log_likelihoods_for_each_hypothesis[indices, prob_idx] = log_prob return ( log_likelihoods_for_each_hypothesis, @@ -307,7 +279,7 @@ def _ref_name(ref: AdRef | str) -> str: @doctest_skipif(reason="Illustrative but not runnable doctest code") def hashsolo( adata: AnnData, - cell_hashing_columns: Collection[AdRef | str], + hashes: Collection[AdRef | str], *, priors: tuple[float, float, float] = (0.01, 0.8, 0.19), pre_existing_clusters: AdRef | str | None = None, @@ -324,7 +296,7 @@ def hashsolo( adata The (annotated) data matrix of shape `n_obs` × `n_vars`. Rows correspond to cells and columns to genes. - cell_hashing_columns + hashes References to the vectors holding the cell hashing counts\ [#ref]_, e.g. `A.obs[["Hash1", "Hash2"]]` for columns in :attr:`~anndata.AnnData.obs` or `A.X[:, ["Hash1", "Hash2"]]` for hashing barcodes among the @@ -341,7 +313,7 @@ def hashsolo( If provided, demultiplexing is performed separately for each cluster. number_of_noise_barcodes The number of barcodes used to create the noise distribution. - Defaults to `len(cell_hashing_columns) - 2`. + Defaults to `len(hashes) - 2`. key_added Key under which to add the demultiplexing results. copy @@ -353,7 +325,7 @@ def hashsolo( Sets the following fields: `adata.obs[key_added]` : :class:`~pandas.Categorical` (shape `(n_obs,)`) - Classification of each cell: the name of one of the `cell_hashing_columns`, + 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. @@ -379,7 +351,7 @@ def hashsolo( """ adata = adata.copy() if copy else adata - refs = list(cell_hashing_columns) + refs = list(hashes) data = np.column_stack(_get_vec(adata, refs, dim="obs")) clusters = ( None From 02e0ca89866d2fdf8b2299c70112eae67e8d8d51 Mon Sep 17 00:00:00 2001 From: "Philipp A." Date: Tue, 18 Aug 2026 13:15:11 +0200 Subject: [PATCH 05/10] multiacc support --- docs/release-notes/4303.feat.md | 2 +- src/scanpy/external/pp/_hashsolo.py | 34 +++++++++-- src/scanpy/preprocessing/_hashsolo.py | 82 +++++++++++---------------- tests/test_hashsolo.py | 37 +++++++++--- 4 files changed, 90 insertions(+), 65 deletions(-) diff --git a/docs/release-notes/4303.feat.md b/docs/release-notes/4303.feat.md index b80405a421..4278c13e08 100644 --- a/docs/release-notes/4303.feat.md +++ b/docs/release-notes/4303.feat.md @@ -1,3 +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` *or* in the data matrix (e.g. `A.X[:, ["Hash1", "Hash2"]]`), +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 e88a4efd82..59780456ef 100644 --- a/src/scanpy/external/pp/_hashsolo.py +++ b/src/scanpy/external/pp/_hashsolo.py @@ -7,10 +7,11 @@ from typing import TYPE_CHECKING +import numpy as np from scverse_misc import Deprecation, deprecated from ..._utils._doctests import doctest_skipif -from ...preprocessing._hashsolo import _legacy_hashsolo +from ...preprocessing._hashsolo import _hashsolo if TYPE_CHECKING: from collections.abc import Sequence @@ -90,11 +91,32 @@ def hashsolo( print( "Please cite HashSolo paper:\nhttps://www.cell.com/cell-systems/fulltext/S2405-4712(20)30195-2" ) - return _legacy_hashsolo( - adata, - cell_hashing_columns, + adata = adata if inplace else adata.copy() + cell_hashing_columns = list(cell_hashing_columns) + data = adata.obs[cell_hashing_columns].to_numpy() + clusters = ( + None + if pre_existing_clusters is None + else adata.obs[pre_existing_clusters].to_numpy() + ) + probs = _hashsolo( + data, priors=priors, - pre_existing_clusters=pre_existing_clusters, + clusters=clusters, number_of_noise_barcodes=number_of_noise_barcodes, - inplace=inplace, ) + most_likely_hypothesis = np.argmax(probs, axis=1) + + 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] + + classification = np.asarray( + np.array(cell_hashing_columns)[np.argmax(data, axis=1)], dtype=object + ) + 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/preprocessing/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py index fcf7307129..44c5cd0d6f 100644 --- a/src/scanpy/preprocessing/_hashsolo.py +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -35,14 +35,16 @@ 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_skipif -from ..get.get import _get_vec +from ..get.get import MultiAcc, _get_arr, _get_vec if TYPE_CHECKING: - from collections.abc import Collection, Sequence + from collections.abc import Collection from anndata import AnnData from anndata.acc import AdRef @@ -144,7 +146,7 @@ def gaussian_updates(data: np.ndarray, mu_o: float, std_o: float) -> GaussianPar signal_params_dict: dict[int, GaussianParams] = {} # for each barcode get empirical noise and signal distribution parameterization - for x in np.arange(num_of_barcodes): + for x in range(num_of_barcodes): sample_barcodes = data[:, x] sample_barcodes_noise_idx = np.where(data_arg[:, :num_of_noise_barcodes] == x)[ 0 @@ -168,7 +170,7 @@ def gaussian_updates(data: np.ndarray, mu_o: float, std_o: float) -> GaussianPar # 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) + range(num_of_barcodes), repeat=2 ): signal_subset = data_arg[:, -1] == signal_sample_idx noise_subset = data_arg[:, -2] == noise_sample_idx @@ -276,10 +278,27 @@ def _ref_name(ref: AdRef | str) -> str: 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): + refs = list(hashes) + data = np.column_stack(_get_vec(adata, refs, dim="obs")) + return data, np.array([_ref_name(ref) for ref in refs]) + # 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_skipif(reason="Illustrative but not runnable doctest code") def hashsolo( adata: AnnData, - hashes: Collection[AdRef | str], + hashes: MultiAcc | Collection[AdRef] | Collection[str], *, priors: tuple[float, float, float] = (0.01, 0.8, 0.19), pre_existing_clusters: AdRef | str | None = None, @@ -301,6 +320,10 @@ def hashsolo( 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 @@ -345,14 +368,17 @@ def hashsolo( >>> sc.pp.hashsolo(adata, A.obs[["Hash1", "Hash2", "Hash3"]]) + All hashing counts stored together in :attr:`~anndata.AnnData.obsm`: + + >>> sc.pp.hashsolo(adata, A.obsm["hto"]) + .. [#ref] If :attr:`scanpy.settings.preset` is :attr:`~scanpy.Preset.ScanpyV2Preview`, :class:`str`\ s are :meth:`anndata.acc.AdAcc.resolve`\ d to :class:`~anndata.acc.AdRef`\ s, otherwise interpreted as :attr:`anndata.AnnData.obs` columns. """ adata = adata.copy() if copy else adata - refs = list(hashes) - data = np.column_stack(_get_vec(adata, refs, dim="obs")) + data, names = _get_hashes(adata, hashes) clusters = ( None if pre_existing_clusters is None @@ -365,7 +391,6 @@ def hashsolo( number_of_noise_barcodes=number_of_noise_barcodes, ) - names = np.array([_ref_name(ref) for ref in refs]) most_likely_hypothesis = np.argmax(probs, axis=1) classification = pd.Series( names[np.argmax(data, axis=1)], index=adata.obs_names, dtype="string" @@ -378,44 +403,3 @@ def hashsolo( probs, columns=["negative", "singlet", "doublet"], index=adata.obs_names ) return adata if copy else None - - -def _legacy_hashsolo( - adata: AnnData, - cell_hashing_columns: Sequence[str], - *, - priors: tuple[float, float, float], - pre_existing_clusters: str | None, - number_of_noise_barcodes: int | None, - inplace: bool, -) -> AnnData | None: - """Implement the pre-1.13 `scanpy.external.pp.hashsolo` API.""" - adata = adata if inplace else adata.copy() - cell_hashing_columns = list(cell_hashing_columns) - data = adata.obs[cell_hashing_columns].to_numpy() - clusters = ( - None - if pre_existing_clusters is None - else adata.obs[pre_existing_clusters].to_numpy() - ) - probs = _hashsolo( - data, - priors=priors, - clusters=clusters, - number_of_noise_barcodes=number_of_noise_barcodes, - ) - most_likely_hypothesis = np.argmax(probs, axis=1) - - 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] - - classification = np.asarray( - np.array(cell_hashing_columns)[np.argmax(data, axis=1)], dtype=object - ) - 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/tests/test_hashsolo.py b/tests/test_hashsolo.py index dbc0015a38..98fc1f8782 100644 --- a/tests/test_hashsolo.py +++ b/tests/test_hashsolo.py @@ -1,6 +1,6 @@ from __future__ import annotations -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, NamedTuple import numpy as np import pandas as pd @@ -15,11 +15,19 @@ if TYPE_CHECKING: from collections.abc import Sequence - from anndata.acc import AdRef + from anndata.acc import AdRef, MultiAcc 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`. @@ -50,37 +58,48 @@ def adata(counts: pd.DataFrame) -> AnnData: pytest.param("obs", marks=needs.anndata_acc), pytest.param("X", marks=needs.anndata_acc), pytest.param("X-sparse", marks=needs.anndata_acc), + pytest.param("obsm-df", marks=needs.anndata_acc), + pytest.param("obsm-array", marks=needs.anndata_acc), ] ) def hashed( request: pytest.FixtureRequest, adata: AnnData, counts: pd.DataFrame -) -> tuple[AnnData, Sequence[AdRef | str]]: +) -> Hashed: """Build an `AnnData` holding the hashing counts, plus references to them.""" match request.param: case "obs-str": # plain `.obs` column names work without `anndata.acc` - return adata, HASHES + return Hashed(adata, HASHES, HASHES) case "obs": from anndata.acc import A - return adata, A.obs[HASHES] + return Hashed(adata, A.obs[HASHES], HASHES) case "X" | "X-sparse": from anndata.acc import A x = counts.to_numpy() if request.param == "X-sparse": x = sparse.csr_matrix(x) # noqa: TID251 - return AnnData(x, var=pd.DataFrame(index=HASHES)), A.X[:, HASHES] + 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” + from anndata.acc import A + + 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}") -def test_cell_demultiplexing(hashed: tuple[AnnData, Sequence[AdRef | str]]) -> None: - adata, refs = hashed +def test_cell_demultiplexing(hashed: Hashed) -> None: + adata, refs, names = hashed sc.pp.hashsolo(adata, refs) expected = pd.array( ["Doublet"] * 10 - + np.repeat(HASHES, 98).reshape(98, 10, order="F").ravel().tolist() + + np.repeat(names, 98).reshape(98, 10, order="F").ravel().tolist() + ["Negative"] * 10, dtype="string", ) From 79ac83297dadb321493ca73b5cd0c6a617db6725 Mon Sep 17 00:00:00 2001 From: "Philipp A." Date: Tue, 18 Aug 2026 13:59:00 +0200 Subject: [PATCH 06/10] run doctest --- src/scanpy/preprocessing/_hashsolo.py | 37 ++++++++++++++++++--------- 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/src/scanpy/preprocessing/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py index 44c5cd0d6f..6b87243799 100644 --- a/src/scanpy/preprocessing/_hashsolo.py +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -40,7 +40,7 @@ from .._compat import CSBase from .._utils import check_nonnegative_integers -from .._utils._doctests import doctest_skipif +from .._utils._doctests import doctest_needs from ..get.get import MultiAcc, _get_arr, _get_vec if TYPE_CHECKING: @@ -295,7 +295,7 @@ def _get_hashes( return data, np.arange(data.shape[1]).astype(str) -@doctest_skipif(reason="Illustrative but not runnable doctest code") +@doctest_needs("anndata_acc") def hashsolo( adata: AnnData, hashes: MultiAcc | Collection[AdRef] | Collection[str], @@ -355,22 +355,35 @@ def hashsolo( Examples -------- - Hashing counts stored as :attr:`~anndata.AnnData.var_names` - (e.g. the “Multiplexing Capture” features of a Cell Ranger run): + 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 - >>> adata = sc.io.read_h5ad("data.h5ad") - >>> sc.pp.hashsolo(adata, A.X[:, ["Hash1", "Hash2", "Hash3"]]) - >>> adata.obs["hashsolo"].value_counts() + >>> + >>> 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)) - Hashing counts stored as :attr:`~anndata.AnnData.obs` columns: - - >>> sc.pp.hashsolo(adata, A.obs[["Hash1", "Hash2", "Hash3"]]) - - All hashing counts stored together in :attr:`~anndata.AnnData.obsm`: + 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') .. [#ref] If :attr:`scanpy.settings.preset` is :attr:`~scanpy.Preset.ScanpyV2Preview`, :class:`str`\ s are :meth:`anndata.acc.AdAcc.resolve`\ d to :class:`~anndata.acc.AdRef`\ s, From 874a498911d3979ccb87e1b17743019c558a424a Mon Sep 17 00:00:00 2001 From: Phil Schaf Date: Thu, 20 Aug 2026 13:59:48 +0200 Subject: [PATCH 07/10] simplify --- src/scanpy/external/pp/_hashsolo.py | 2 +- src/scanpy/preprocessing/_hashsolo.py | 180 ++++++++------------------ 2 files changed, 58 insertions(+), 124 deletions(-) diff --git a/src/scanpy/external/pp/_hashsolo.py b/src/scanpy/external/pp/_hashsolo.py index 59780456ef..cf07d820a6 100644 --- a/src/scanpy/external/pp/_hashsolo.py +++ b/src/scanpy/external/pp/_hashsolo.py @@ -19,7 +19,7 @@ from anndata import AnnData -@deprecated(Deprecation("1.13.0", "Use `scanpy.pp.hashsolo` instead.")) +@deprecated(Deprecation("1.13.0", "Use :func:`scanpy.pp.hashsolo` instead.")) @doctest_skipif(reason="Illustrative but not runnable doctest code") def hashsolo( adata: AnnData, diff --git a/src/scanpy/preprocessing/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py index 6b87243799..34e40b98e8 100644 --- a/src/scanpy/preprocessing/_hashsolo.py +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -51,18 +51,35 @@ from numpy.typing import ArrayLike, NDArray -class GaussianParams(NamedTuple): - """Parameters of a gaussian, named after :func:`scipy.stats.norm.pdf`’s parameters.""" +class Gaussian(NamedTuple): + """A gaussian, its fields named after :func:`scipy.stats.norm.pdf`’s parameters.""" - loc: float + loc: float | np.floating """Mean of the gaussian.""" - scale: float + 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 / (self.scale**2) + lam = 1 / np.var(data) 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: np.ndarray, number_of_noise_barcodes: int | None -) -> tuple[NDArray[np.float64], NDArray[np.float64], dict[int, str]]: + data: NDArray[np.floating], number_of_noise_barcodes: int | None +) -> NDArray[np.float64]: """Calculate log likelihoods for each hypothesis, negative, singlet, doublet. Parameters @@ -74,144 +91,61 @@ def _calculate_log_likelihoods( Returns ------- - log_likelihoods_for_each_hypothesis - a 2d np.array log likelihood of each hypothesis - all_indices - counter_to_barcode_combo + A 2d array of shape `(n_cells, 3)` with the log likelihood of each hypothesis. """ - - def gaussian_updates(data: np.ndarray, mu_o: float, std_o: float) -> GaussianParams: - """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 - ------- - The updated parameters of the 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 GaussianParams(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)) + log_likelihoods = 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 + n_barcodes = data.shape[1] + n_barcodes_noise = ( + number_of_noise_barcodes if number_of_noise_barcodes is not None - else 2 + else n_barcodes - 2 ) - num_of_noise_barcodes = num_of_barcodes - number_of_non_noise_barcodes - # assume log normal - data = np.log(data + 1) + 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_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: dict[int, GaussianParams] = {} - signal_params_dict: dict[int, GaussianParams] = {} - - # for each barcode get empirical noise and signal distribution parameterization - for x in range(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_params_dict[x] = gaussian_updates( - noise_counts, global_mu_noise_o, global_sigma_noise_o - ) - signal_params_dict[x] = gaussian_updates( - signal_counts, global_mu_signal_o, global_sigma_signal_o - ) - - counter_to_barcode_combo: dict[int, str] = {} - counter = 0 + 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 noise_sample_idx, signal_sample_idx in product( - range(num_of_barcodes), repeat=2 - ): - 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: + 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 - 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] - data_noise = data_subset[:, noise_sample_idx] - data_signal = data_subset[:, signal_sample_idx] # the distributions the 2nd-highest and highest barcode’s counts # are assumed to come from under each hypothesis hypotheses = [ - (noise_params, noise_params), # negative: neither barcode is signal - (noise_params, signal_params), # singlet: only the highest barcode is - (signal_params_dict[noise_sample_idx], signal_params), # doublet: both are + (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 ] - - # each cell and each hypothesis probability - for prob_idx, log_prob in enumerate( - ( - np.log(norm.pdf(data_noise, loc=d_noise.loc, scale=d_noise.scale) + eps) - + np.log(norm.pdf(data_signal, loc=d_sig.loc, scale=d_sig.scale) + eps) - for d_noise, d_sig in hypotheses - ) - ): - log_likelihoods_for_each_hypothesis[indices, prob_idx] = log_prob - return ( - log_likelihoods_for_each_hypothesis, - all_indices, - counter_to_barcode_combo, - ) + 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( @@ -233,7 +167,7 @@ def _calculate_bayes_rule( A 2d array of shape `(n_cells, 3)` with the probability of each hypothesis. """ - log_likelihoods, _, _ = _calculate_log_likelihoods(data, number_of_noise_barcodes) + log_likelihoods = _calculate_log_likelihoods(data, number_of_noise_barcodes) likelihoods = np.exp(log_likelihoods) * np.asarray(priors) return likelihoods / likelihoods.sum(axis=1)[:, None] From 243d0418a89d83a6d5e23e666b2b6fe390dc7a55 Mon Sep 17 00:00:00 2001 From: Phil Schaf Date: Thu, 20 Aug 2026 15:02:28 +0200 Subject: [PATCH 08/10] no legacy mode --- src/scanpy/get/__init__.py | 2 ++ src/scanpy/get/_aggregated.py | 4 +-- src/scanpy/get/get.py | 52 ++++++++++++++++++++++++--- src/scanpy/preprocessing/_hashsolo.py | 17 ++++----- tests/test_hashsolo.py | 32 ++++++++--------- 5 files changed, 73 insertions(+), 34 deletions(-) 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/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py index 34e40b98e8..a602af7dc6 100644 --- a/src/scanpy/preprocessing/_hashsolo.py +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -78,7 +78,7 @@ def posterior(self, data: NDArray[np.floating]) -> Gaussian: def _calculate_log_likelihoods( - data: NDArray[np.floating], number_of_noise_barcodes: int | None + data: NDArray[np.integer], number_of_noise_barcodes: int | None ) -> NDArray[np.float64]: """Calculate log likelihoods for each hypothesis, negative, singlet, doublet. @@ -149,7 +149,7 @@ def _calculate_log_likelihoods( def _calculate_bayes_rule( - data: np.ndarray, priors: ArrayLike, number_of_noise_barcodes: int | None + data: NDArray[np.integer], priors: ArrayLike, number_of_noise_barcodes: int | None ) -> NDArray[np.float64]: """Calculate the posterior probability of each hypothesis from log likelihoods. @@ -173,7 +173,7 @@ def _calculate_bayes_rule( def _hashsolo( - data: NDArray, + data: NDArray[np.integer], *, priors: ArrayLike, clusters: NDArray | None, @@ -217,9 +217,8 @@ def _get_hashes( ) -> tuple[NDArray, NDArray[np.str_]]: """Get the hashing count matrix and each hashing barcode’s name.""" if not isinstance(hashes, MultiAcc): - refs = list(hashes) - data = np.column_stack(_get_vec(adata, refs, dim="obs")) - return data, np.array([_ref_name(ref) for ref in refs]) + 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): @@ -250,7 +249,7 @@ def hashsolo( 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\ [#ref]_, + 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`. @@ -319,10 +318,6 @@ def hashsolo( >>> adata.obs["hs2"].cat.categories.astype("string") Index(['Hash1', 'Hash2', 'Hash3'], dtype='string') - .. [#ref] If :attr:`scanpy.settings.preset` is :attr:`~scanpy.Preset.ScanpyV2Preview`, - :class:`str`\ s are :meth:`anndata.acc.AdAcc.resolve`\ d to :class:`~anndata.acc.AdRef`\ s, - otherwise interpreted as :attr:`anndata.AnnData.obs` columns. - """ adata = adata.copy() if copy else adata data, names = _get_hashes(adata, hashes) diff --git a/tests/test_hashsolo.py b/tests/test_hashsolo.py index 98fc1f8782..1c12bec3ab 100644 --- a/tests/test_hashsolo.py +++ b/tests/test_hashsolo.py @@ -52,23 +52,12 @@ def adata(counts: pd.DataFrame) -> AnnData: return AnnData(rng.integers(0, 100, size=counts.shape), obs=counts.copy()) -@pytest.fixture( - params=[ - "obs-str", - pytest.param("obs", marks=needs.anndata_acc), - pytest.param("X", marks=needs.anndata_acc), - pytest.param("X-sparse", marks=needs.anndata_acc), - pytest.param("obsm-df", marks=needs.anndata_acc), - pytest.param("obsm-array", marks=needs.anndata_acc), - ] -) +@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-str": # plain `.obs` column names work without `anndata.acc` - return Hashed(adata, HASHES, HASHES) case "obs": from anndata.acc import A @@ -93,6 +82,7 @@ def hashed( pytest.fail(f"Unknown param {request.param!r}") +@needs.anndata_acc def test_cell_demultiplexing(hashed: Hashed) -> None: adata, refs, names = hashed sc.pp.hashsolo(adata, refs) @@ -116,19 +106,29 @@ def test_cell_demultiplexing(hashed: Hashed) -> None: 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, HASHES, copy=True) + from anndata.acc import A + + copied = sc.pp.hashsolo(adata, A.obs[HASHES], copy=True) assert "hashsolo" not in adata.obs assert "hashsolo" in copied.obs -def test_legacy_api(adata: AnnData) -> None: - """The deprecated `sce.pp.hashsolo` matches `sc.pp.hashsolo` and needs no `anndata.acc`.""" +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: + from anndata.acc import A + with pytest.warns(FutureWarning, match=r"scanpy\.pp\.hashsolo"): sce.pp.hashsolo(adata, HASHES) expected_new = adata.copy() - sc.pp.hashsolo(expected_new, HASHES) + sc.pp.hashsolo(expected_new, A.obs[HASHES]) pd.testing.assert_series_equal( adata.obs["Classification"].astype("string"), expected_new.obs["hashsolo"].astype("string"), From 0d6c2efe1a799bfec207047a3f7f47ea45482d18 Mon Sep 17 00:00:00 2001 From: Phil Schaf Date: Thu, 20 Aug 2026 15:17:06 +0200 Subject: [PATCH 09/10] better error when too few noise barcodes --- src/scanpy/external/pp/_hashsolo.py | 2 +- src/scanpy/preprocessing/_hashsolo.py | 44 +++++++++++---------------- tests/test_hashsolo.py | 19 ++++++++++++ 3 files changed, 38 insertions(+), 27 deletions(-) diff --git a/src/scanpy/external/pp/_hashsolo.py b/src/scanpy/external/pp/_hashsolo.py index cf07d820a6..d98497a31f 100644 --- a/src/scanpy/external/pp/_hashsolo.py +++ b/src/scanpy/external/pp/_hashsolo.py @@ -103,7 +103,7 @@ def hashsolo( data, priors=priors, clusters=clusters, - number_of_noise_barcodes=number_of_noise_barcodes, + n_barcodes_noise=number_of_noise_barcodes, ) most_likely_hypothesis = np.argmax(probs, axis=1) diff --git a/src/scanpy/preprocessing/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py index a602af7dc6..3d8f99126b 100644 --- a/src/scanpy/preprocessing/_hashsolo.py +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -78,7 +78,7 @@ def posterior(self, data: NDArray[np.floating]) -> Gaussian: def _calculate_log_likelihoods( - data: NDArray[np.integer], number_of_noise_barcodes: int | None + data: NDArray[np.integer], n_barcodes_noise: int ) -> NDArray[np.float64]: """Calculate log likelihoods for each hypothesis, negative, singlet, doublet. @@ -86,7 +86,7 @@ def _calculate_log_likelihoods( ---------- data cells by hashing counts matrix - number_of_noise_barcodes + n_barcodes_noise number of barcodes to used to calculated noise distribution Returns @@ -98,11 +98,6 @@ def _calculate_log_likelihoods( log_likelihoods = np.zeros((data.shape[0], 3)) n_barcodes = data.shape[1] - n_barcodes_noise = ( - number_of_noise_barcodes - if number_of_noise_barcodes is not None - else n_barcodes - 2 - ) # assume log normal data: NDArray[np.floating] = np.log(data + 1) @@ -149,7 +144,7 @@ def _calculate_log_likelihoods( def _calculate_bayes_rule( - data: NDArray[np.integer], priors: ArrayLike, number_of_noise_barcodes: int | None + data: NDArray[np.integer], priors: ArrayLike, n_barcodes_noise: int ) -> NDArray[np.float64]: """Calculate the posterior probability of each hypothesis from log likelihoods. @@ -159,7 +154,7 @@ def _calculate_bayes_rule( cells by hashing counts matrix priors prior for each hypothesis, in the order `[negative, singlet, doublet]` - number_of_noise_barcodes + n_barcodes_noise number of barcodes to used to calculated noise distribution Returns @@ -167,7 +162,7 @@ def _calculate_bayes_rule( A 2d array of shape `(n_cells, 3)` with the probability of each hypothesis. """ - log_likelihoods = _calculate_log_likelihoods(data, number_of_noise_barcodes) + log_likelihoods = _calculate_log_likelihoods(data, n_barcodes_noise) likelihoods = np.exp(log_likelihoods) * np.asarray(priors) return likelihoods / likelihoods.sum(axis=1)[:, None] @@ -177,30 +172,30 @@ def _hashsolo( *, priors: ArrayLike, clusters: NDArray | None, - number_of_noise_barcodes: int | 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" raise ValueError(msg) - if number_of_noise_barcodes is not None and number_of_noise_barcodes >= ( - n_barcodes := data.shape[1] - ): + 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"`number_of_noise_barcodes` ({number_of_noise_barcodes}) must be smaller " - f"than the number of hashing barcodes ({n_barcodes})." + 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, number_of_noise_barcodes) + return _calculate_bayes_rule(data, priors, n_barcodes_noise) probs = np.zeros((data.shape[0], 3)) for cluster in pd.unique(clusters): mask = clusters == cluster - probs[mask] = _calculate_bayes_rule( - data[mask], priors, number_of_noise_barcodes - ) + probs[mask] = _calculate_bayes_rule(data[mask], priors, n_barcodes_noise) return probs @@ -235,7 +230,7 @@ def hashsolo( *, priors: tuple[float, float, float] = (0.01, 0.8, 0.19), pre_existing_clusters: AdRef | str | None = None, - number_of_noise_barcodes: int | None = None, + n_noise_barcodes: int | None = None, key_added: str = "hashsolo", copy: bool = False, ) -> AnnData | None: @@ -267,7 +262,7 @@ def hashsolo( 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. - number_of_noise_barcodes + n_noise_barcodes The number of barcodes used to create the noise distribution. Defaults to `len(hashes) - 2`. key_added @@ -327,10 +322,7 @@ def hashsolo( else np.asarray(_get_vec(adata, pre_existing_clusters, dim="obs")) ) probs = _hashsolo( - data, - priors=priors, - clusters=clusters, - number_of_noise_barcodes=number_of_noise_barcodes, + data, priors=priors, clusters=clusters, n_barcodes_noise=n_noise_barcodes ) most_likely_hypothesis = np.argmax(probs, axis=1) diff --git a/tests/test_hashsolo.py b/tests/test_hashsolo.py index 1c12bec3ab..0706f209e0 100644 --- a/tests/test_hashsolo.py +++ b/tests/test_hashsolo.py @@ -82,6 +82,25 @@ def hashed( 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.""" + from anndata.acc import A + + 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 From ece9b6df278ccf3fafb14355078f8b6a53e3f035 Mon Sep 17 00:00:00 2001 From: Phil Schaf Date: Thu, 20 Aug 2026 15:44:11 +0200 Subject: [PATCH 10/10] address comments --- src/scanpy/preprocessing/_hashsolo.py | 18 +++++++--- tests/test_hashsolo.py | 49 ++++++++++++++++++++------- 2 files changed, 50 insertions(+), 17 deletions(-) diff --git a/src/scanpy/preprocessing/_hashsolo.py b/src/scanpy/preprocessing/_hashsolo.py index 3d8f99126b..6e91143be2 100644 --- a/src/scanpy/preprocessing/_hashsolo.py +++ b/src/scanpy/preprocessing/_hashsolo.py @@ -51,6 +51,12 @@ 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.""" @@ -70,8 +76,8 @@ def posterior(self, data: NDArray[np.floating]) -> Gaussian: See . """ n = len(data) - lam_o = 1 / (self.scale**2) - lam = 1 / np.var(data) if n > 1 else lam_o + 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)) @@ -176,7 +182,7 @@ def _hashsolo( ) -> 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" + msg = "Cell hashing counts must be non-negative integers" raise ValueError(msg) n_barcodes = data.shape[1] if n_barcodes_noise is None: @@ -193,8 +199,10 @@ def _hashsolo( return _calculate_bayes_rule(data, priors, n_barcodes_noise) probs = np.zeros((data.shape[0], 3)) - for cluster in pd.unique(clusters): - mask = clusters == cluster + # `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 diff --git a/tests/test_hashsolo.py b/tests/test_hashsolo.py index 0706f209e0..9501f32573 100644 --- a/tests/test_hashsolo.py +++ b/tests/test_hashsolo.py @@ -2,6 +2,7 @@ from typing import TYPE_CHECKING, NamedTuple +import anndata import numpy as np import pandas as pd import pytest @@ -17,6 +18,10 @@ 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)] @@ -59,20 +64,14 @@ def hashed( """Build an `AnnData` holding the hashing counts, plus references to them.""" match request.param: case "obs": - from anndata.acc import A - return Hashed(adata, A.obs[HASHES], HASHES) case "X" | "X-sparse": - from anndata.acc import A - 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” - from anndata.acc import A - 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 @@ -95,8 +94,6 @@ 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.""" - from anndata.acc import A - with pytest.raises(ValueError, match=r"noise barcodes?"): sc.pp.hashsolo(adata, A.obs[HASHES[:n_hashes]], n_noise_barcodes=n_noise) @@ -127,13 +124,43 @@ def test_cell_demultiplexing(hashed: Hashed) -> None: @needs.anndata_acc def test_copy(adata: AnnData) -> None: - from anndata.acc import A - 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) @@ -141,8 +168,6 @@ def test_legacy_api_needs_no_acc(subtests: pytest.Subtests, adata: AnnData) -> N @needs.anndata_acc def test_legacy_api_matches(subtests: pytest.Subtests, adata: AnnData) -> None: - from anndata.acc import A - with pytest.warns(FutureWarning, match=r"scanpy\.pp\.hashsolo"): sce.pp.hashsolo(adata, HASHES)