diff --git a/HISTORY.rst b/HISTORY.rst index 5b279cfa..201e66c6 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -2,6 +2,31 @@ History ======= +1.5.1 (2026-11-07) +------------------ +* ``st.tl.cci.run`` now uses a spots-by-spots, sparse graph: + ``adata.obsp["spot_neighbours"]``. It replaces the comma-joined strings + ``adata.obsm["spot_neighbours"]`` and ``adata.obsm["spot_neigh_bcs"]``. +* With ``use_label``, ``run`` counts cell types over the same graph instead of + searching the neighbours a second time (``het.count`` takes it as ``neighbours=``). +* ``st.tl.cache.write_subset_h5ad`` takes ``obsp_keys``, and + ``merge_h5ad_into_adata`` copies ``obsp``, so a cached CCI run keeps its neighbours. + +API and Bug Fixes: + +* ``stlearn.tl.cci.neighbours.calc_neighbours`` returns that sparse matrix instead + of a numba ``List``. ``neighbourhood_matrix`` and + ``neighbourhoods_from_barcodes`` are removed. +* ``LRScorer`` takes a sparse matrix. ``get_scores``, ``adjacency`` and + ``flatten_neighbours`` are removed from ``stlearn.tl.cci.scoring``. +* ``st.tl.cci.run_cci`` and ``get_edges`` read the stored graph directly; + ``stlearn.tl.cci.neighbours.get_neighbourhoods`` and the edge helpers in + ``stlearn.tl.cci.edges`` are removed. ``get_interaction_matrix`` and + ``get_interaction_pvals`` take the graph's CSR arrays (``neigh_indptr``, + ``neigh_indices``) instead of a numba ``List``. +* ``get_edges`` no longer deduplicates edges with a quadratic loop: the arrow + plots on large neighbourhoods are seconds faster. + 1.5.0 (2026-11-01) ------------------ This release contains changes in cell-cell interaction (CCI) mainly diff --git a/docs/conf.py b/docs/conf.py index fb943bf1..df618055 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -21,7 +21,7 @@ exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] TUTORIALS_URL = ( - "https://github.com/GenomicsMachineLearning/stLearn-tutorials/releases/download/V1.5.0/tutorials.zip" + "https://github.com/GenomicsMachineLearning/stLearn-tutorials/releases/download/V1.5.1/tutorials.zip" ) diff --git a/docs/index.rst b/docs/index.rst index 6bade652..42ccfc57 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -38,6 +38,8 @@ undissociated tissue sample. Latest Additions ---------------- +.. include:: release_notes/1.5.1.rst + .. include:: release_notes/1.5.0.rst .. include:: release_notes/1.4.1.rst diff --git a/docs/release_notes/1.5.1.rst b/docs/release_notes/1.5.1.rst new file mode 100644 index 00000000..a95178f6 --- /dev/null +++ b/docs/release_notes/1.5.1.rst @@ -0,0 +1,24 @@ +1.5.1 (2026-11-07) +~~~~~~~~~~~~~~~~~~ + +.. rubric:: Features + +* ``st.tl.cci.run`` now uses a spots-by-spots, sparse graph: + ``adata.obsp["spot_neighbours"]``. It replaces the comma-joined strings + ``adata.obsm["spot_neighbours"]`` and ``adata.obsm["spot_neigh_bcs"]``. +* With ``use_label``, ``run`` counts cell types over the same graph instead of + searching the neighbours a second time (``het.count`` takes it as ``neighbours=``). +* ``st.tl.cache.write_subset_h5ad`` takes ``obsp_keys``, and + ``merge_h5ad_into_adata`` copies ``obsp``, so a cached CCI run keeps its neighbours. +* ``stlearn.tl.cci.neighbours.calc_neighbours`` returns that sparse matrix instead + of a numba ``List``. ``neighbourhood_matrix`` and + ``neighbourhoods_from_barcodes`` are removed. +* ``LRScorer`` takes a sparse matrix. ``get_scores``, ``adjacency`` and + ``flatten_neighbours`` are removed from ``stlearn.tl.cci.scoring``. +* ``st.tl.cci.run_cci`` and ``get_edges`` read the stored graph directly; + ``stlearn.tl.cci.neighbours.get_neighbourhoods`` and the edge helpers in + ``stlearn.tl.cci.edges`` are removed. ``get_interaction_matrix`` and + ``get_interaction_pvals`` take the graph's CSR arrays (``neigh_indptr``, + ``neigh_indices``) instead of a numba ``List``. +* ``get_edges`` no longer deduplicates edges with a quadratic loop: the arrow + plots on large neighbourhoods are seconds faster. \ No newline at end of file diff --git a/docs/release_notes/index.rst b/docs/release_notes/index.rst index 63643a9e..9d3585e0 100644 --- a/docs/release_notes/index.rst +++ b/docs/release_notes/index.rst @@ -1,6 +1,8 @@ Release Notes =================================================== +.. include:: 1.5.1.rst + .. include:: 1.5.0.rst .. include:: 1.4.1.rst diff --git a/pyproject.toml b/pyproject.toml index 4d371fad..f7b17ff4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "stlearn" -version = "1.5.0" +version = "1.5.1" authors = [ { name = "Genomics and Machine Learning Lab", email = "andrew.newman@uq.edu.au" }, ] diff --git a/stlearn/__init__.py b/stlearn/__init__.py index 7c9439b6..d8e882eb 100644 --- a/stlearn/__init__.py +++ b/stlearn/__init__.py @@ -2,7 +2,7 @@ __author__ = """Genomics and Machine Learning Lab""" __email__ = "andrew.newman@uq.edu.au" -__version__ = "1.5.0" +__version__ = "1.5.1" from . import add, datasets, em, io, pl, pp, spatial, tl, types, utils from ._settings import settings diff --git a/stlearn/tl/cache/anndata.py b/stlearn/tl/cache/anndata.py index a208feb8..0d585d4e 100644 --- a/stlearn/tl/cache/anndata.py +++ b/stlearn/tl/cache/anndata.py @@ -3,8 +3,8 @@ import pandas as pd -def write_subset_h5ad(adata, filename, obsm_keys=None, uns_keys=None): - """Write only specific obsm and uns components to H5AD""" +def write_subset_h5ad(adata, filename, obsm_keys=None, uns_keys=None, obsp_keys=None): + """Write only specific obsm, obsp and uns components to H5AD""" # Create a minimal AnnData object with the same structure minimal_adata = ad.AnnData( @@ -32,6 +32,14 @@ def write_subset_h5ad(adata, filename, obsm_keys=None, uns_keys=None): else: print(f"Warning: uns['{key}'] not found") + if obsp_keys: + for key in obsp_keys: + if key in adata.obsp: + minimal_adata.obsp[key] = adata.obsp[key] + print(f"Added obsp['{key}'] with {adata.obsp[key].nnz} entries") + else: + print(f"Warning: obsp['{key}'] not found") + minimal_adata.write_h5ad(filename, compression="gzip", compression_opts=9) print(f"Wrote subset to {filename}") @@ -48,4 +56,8 @@ def merge_h5ad_into_adata(adata_main, h5ad_file): adata_main.uns[key] = value print(f"Added uns['{key}']") + for key, value in adata_subset.obsp.items(): + adata_main.obsp[key] = value + print(f"Added obsp['{key}'] with {value.nnz} entries") + return adata_main diff --git a/stlearn/tl/cci/analysis.py b/stlearn/tl/cci/analysis.py index bb846e44..9a37440b 100644 --- a/stlearn/tl/cci/analysis.py +++ b/stlearn/tl/cci/analysis.py @@ -4,12 +4,12 @@ import warnings +import anndata import numpy as np -import numpy.typing as npt -import pandas as pd -from anndata import AnnData -from scipy import sparse -from tqdm import tqdm +import numpy.typing +import pandas +import scipy +import tqdm from ...utils import _set_n_threads from ._results import ( @@ -28,15 +28,15 @@ label_keys, ) from .ligand_receptors import LRPairs, genes_of, lr_pairs, split_name -from .neighbours import calc_distance, calc_neighbours, get_neighbourhoods +from .neighbours import calc_distance, calc_neighbours, stored_neighbours from .permutation import perform_spot_testing from .scoring import get_lrs_scores def _report_dropped_lrs( - adata: AnnData, + adata: anndata.AnnData, lrs: LRPairs, - n_scored: npt.NDArray[np.int64], + n_scored: numpy.typing.NDArray[np.int64], min_spots: int, min_expr: float, max_named: int = 10, @@ -52,7 +52,8 @@ def _report_dropped_lrs( n_above = np.asarray((min_expr < adata[:, genes].X).sum(axis=0)).ravel() unexpressed = genes[n_above == 0] # sorted, as genes() are by_gene = np.isin(dropped.ligands(), unexpressed) | np.isin( - dropped.receptors(), unexpressed + dropped.receptors(), + unexpressed, ) n_by_gene = int(by_gene.sum()) too_few = len(dropped) - n_by_gene @@ -65,7 +66,7 @@ def _report_dropped_lrs( ) print( f"Dropping {n_by_gene} LR pairs with a gene never expressed above " - f"min_expr={min_expr}: {named}{more}." + f"min_expr={min_expr}: {named}{more}.", ) if too_few: print(f"Dropping {too_few} LR pairs scored in {min_spots} or fewer spots.") @@ -73,7 +74,7 @@ def _report_dropped_lrs( def run( - adata: AnnData, + adata: anndata.AnnData, lrs: LRPairs, min_spots: int = 10, distance: float | None = None, @@ -88,7 +89,7 @@ def run( random_state: int = 0, verbose: bool = True, cell_prop_cutoff: float = 0.2, -) -> AnnData: +) -> anndata.AnnData: """Performs stLearn LR analysis. Parameters @@ -162,8 +163,9 @@ def run( Returns 'lr_species' and 'lrfeatures' With save_bg 'lrs_to_bg' and 'lr_spot_indices'. Relevant information stored in adata: - adata.obsm['spot_neighbours'] and adata.obsm['spot_neigh_bcs'] - The neighbours of each spot. + adata.obsp['spot_neighbours'] + The neighbours of each spot, as a spots-by-spots CSR matrix: row + i's indices are spot i's neighbours. adata.obsm['cci_het'] Only if use_label specified; contains the counts of the cell types found per spot. @@ -180,7 +182,7 @@ def run( raise TypeError( f"lrs must be LRPairs, got {type(lrs).__name__}; use lr_pairs(database), " "read_lr_pairs(path), LRPairs.from_frame(df) or LRPairs.from_names(names, " - "species=...)." + "species=...).", ) if len(lrs) == 0: @@ -202,25 +204,15 @@ def run( # Calculating neighbour & storing # distance = calc_distance(adata, distance) neighbours = calc_neighbours(adata, distance, verbose=verbose) - adata.obsm["spot_neighbours"] = pd.DataFrame( - [",".join(x.astype(str)) for x in neighbours], - index=adata.obs_names, - columns=["neighbour_indices"], - ) - - # Important to store barcodes in-case adata subsetted # - barcodes = adata.obs_names.to_numpy() - adata.obsm["spot_neigh_bcs"] = pd.DataFrame( - list(map(lambda x: ",".join(barcodes[x]), neighbours)), - index=adata.obs_names, - columns=["neighbour_bcs"], - ) + if neighbours.nnz == 0: + raise ValueError( + "All spots have no neighbours at current distance, set distance to " + "higher value, or distance=0 for within-spot mode.", + ) + adata.obsp["spot_neighbours"] = neighbours if verbose: - print( - "Spot neighbour indices stored in adata.obsm['spot_neighbours'] " - "& adata.obsm['spot_neigh_bcs'].", - ) + print("Spot neighbours stored in adata.obsp['spot_neighbours'].") # Weighting by cell heterogeneity if given cell types if use_label is not None: @@ -234,6 +226,7 @@ def run( use_label=use_label, cell_prop_cutoff=cell_prop_cutoff, verbose=verbose, + neighbours=neighbours, ) het_vals = adata.obsm["cci_het"].ravel() else: @@ -242,14 +235,19 @@ def run( # 1. Filter any LRs without stored expression. # Calculating the lr_scores across spots for the inputted lrs # lr_scores, scored = get_lrs_scores( - adata, lrs, neighbours, het_vals, min_expr, sparse_result=True + adata, + lrs, + neighbours, + het_vals, + min_expr, + sparse_result=True, ) # Fail fast. if len(scored) == 0: raise ValueError( f"None of the {len(lrs)} {lrs.species} LR pairs are in adata.var_names; " "check the pairs' species matches the data, e.g. " - "lr_pairs(species='mouse')." + "lr_pairs(species='mouse').", ) n_scored = np.diff(lr_scores.indptr).astype(np.int64) # spots scored per LR lr_bool = n_scored > min_spots @@ -260,7 +258,7 @@ def run( if len(tested) == 0: raise ValueError( f"No LR pairs had scores in more than {min_spots} spots; " - "try lowering min_spots or min_expr." + "try lowering min_spots or min_expr.", ) # Permutation methods generating background per spot, & test lrs in spot. tests = perform_spot_testing( @@ -276,8 +274,9 @@ def run( neg_binom=neg_binom, random_state=random_state, ) - lr_results = AnnData( - obs=pd.DataFrame(index=adata.obs_names), var=pd.DataFrame(index=tested.names()) + lr_results = anndata.AnnData( + obs=pandas.DataFrame(index=adata.obs_names), + var=pandas.DataFrame(index=tested.names()), ) lr_results.var["source"] = tested.to_frame()["source"].to_numpy(dtype=str) @@ -293,13 +292,13 @@ def run( if verbose: print( "\nLR results returned as a table: the summary in .var, the per-spot " - "results in .layers. Pass it to the CCI functions as lr_results=." + "results in .layers. Pass it to the CCI functions as lr_results=.", ) return lr_results def adj_pvals( - lr_results: AnnData, + lr_results: anndata.AnnData, pval_adj_cutoff: float = 0.05, correct_axis: str | None = "spot", adj_method: str = "fdr_bh", @@ -340,10 +339,10 @@ def adj_pvals( def run_lr_go( - lr_results: AnnData, + lr_results: anndata.AnnData, r_path: str, n_top: int = 100, - bg_genes: npt.ArrayLike | None = None, + bg_genes: numpy.typing.ArrayLike | None = None, min_sig_spots: int = 1, species: str | None = None, p_cutoff: float = 0.01, @@ -391,7 +390,7 @@ def run_lr_go( species = pairs_species or "human" elif pairs_species is not None and species != pairs_species: raise ValueError( - f"species is {species!r} but the LR table's pairs are {pairs_species}." + f"species is {species!r} but the LR table's pairs are {pairs_species}.", ) all_species = ["human", "mouse"] if species not in all_species: @@ -425,9 +424,9 @@ def run_lr_go( # Functions for calling Celltype-Celltype interactions def run_cci( - adata: AnnData, + adata: anndata.AnnData, use_label: str, - lr_results: AnnData | None = None, + lr_results: anndata.AnnData | None = None, spot_mixtures: bool = False, min_spots: int = 3, sig_spots: bool = True, @@ -545,9 +544,10 @@ def run_cci( msg = f"Cell type scores misformatted in adata.uns[{uns_key}]:\n" if not cols_present or not rows_present: if not cols_present: + out = [cell for cell in all_set if cell not in adata.uns[uns_key]] msg = ( - msg + f"Cell types missing from adata.uns[{uns_key}] columns:\n" - f"{[cell for cell in all_set if cell not in adata.uns[uns_key]]}\n" + msg + + f"Cell types missing from adata.uns[{uns_key}] columns:{out}\n" ) elif not rows_present: msg = msg + "Rows do not correspond to adata.obs_names.\n" @@ -558,11 +558,9 @@ def run_cci( if len(all_set) < adata.uns[uns_key].shape[1]: all_set = adata.uns[uns_key].columns.to_numpy().astype(str) - # Getting minimum necessary information for edge counting - if verbose: - print("Getting cached neighbourhood information...") - # Getting the neighbourhoods # - _, neighbourhood_indices = get_neighbourhoods(adata) + graph = stored_neighbours(adata) + neigh_indptr = graph.indptr.astype(np.int64) + neigh_indices = graph.indices.astype(np.int64) if verbose: print("Getting information for CCI counting...") @@ -581,10 +579,10 @@ def run_cci( "may need to adjust min_spots, sig_spots parameters, " "or re-run st.tl.cci.run with more relaxed parameters.", ) - lr_expr = sparse.csc_matrix(adata[:, lr_genes].X) # a gene at a time below + lr_expr = scipy.sparse.csc_matrix(adata[:, lr_genes].X) # a gene at a time below gene_column = {gene: i for i, gene in enumerate(lr_genes)} - def expressed(gene: str) -> npt.NDArray[np.bool_]: + def expressed(gene: str) -> numpy.typing.NDArray[np.bool_]: return lr_expr[:, gene_column[gene]].toarray().ravel() > 0 # Sig-CCIs across all LRs # @@ -597,7 +595,7 @@ def expressed(gene: str) -> npt.NDArray[np.bool_]: lr_n_spot_cci = np.zeros(lr_summary.shape[0]) lr_n_spot_cci_sig = np.zeros(lr_summary.shape[0]) lr_n_cci_sig = np.zeros(lr_summary.shape[0]) - with tqdm( + with tqdm.tqdm( total=len(best_lrs), desc="Counting celltype-celltype interactions per LR and permuting " + f"{n_perms} times.", @@ -613,7 +611,8 @@ def expressed(gene: str) -> npt.NDArray[np.bool_]: int_matrix = get_interaction_matrix( cell_data, - neighbourhood_indices, + neigh_indptr, + neigh_indices, all_set, sig_bool, L_bool, @@ -626,7 +625,8 @@ def expressed(gene: str) -> npt.NDArray[np.bool_]: int_matrix, n_perms, cell_data, - neighbourhood_indices, + neigh_indptr, + neigh_indices, all_set, sig_bool, L_bool, @@ -648,9 +648,11 @@ def expressed(gene: str) -> npt.NDArray[np.bool_]: raw_matrix += int_matrix all_matrix += sig_int_matrix - int_df = pd.DataFrame(int_matrix, index=all_set, columns=all_set) - sig_int_df = pd.DataFrame(sig_int_matrix, index=all_set, columns=all_set) - pval_df = pd.DataFrame(int_pvals, index=all_set, columns=all_set) + int_df = pandas.DataFrame(int_matrix, index=all_set, columns=all_set) + sig_int_df = pandas.DataFrame( + sig_int_matrix, index=all_set, columns=all_set + ) + pval_df = pandas.DataFrame(int_pvals, index=all_set, columns=all_set) per_lr_cci[best_lr] = sig_int_df per_lr_cci_pvals[best_lr] = pval_df per_lr_cci_raw[best_lr] = int_df @@ -661,12 +663,12 @@ def expressed(gene: str) -> npt.NDArray[np.bool_]: lr.var[f"n_cci_sig_{use_label}"] = lr_n_cci_sig lr.var[f"n-spot_cci_{use_label}"] = lr_n_spot_cci lr.var[f"n-spot_cci_sig_{use_label}"] = lr_n_spot_cci_sig - lr.uns[f"lr_cci_{use_label}"] = pd.DataFrame( + lr.uns[f"lr_cci_{use_label}"] = pandas.DataFrame( all_matrix, index=all_set, columns=all_set, ) - lr.uns[f"lr_cci_raw_{use_label}"] = pd.DataFrame( + lr.uns[f"lr_cci_raw_{use_label}"] = pandas.DataFrame( raw_matrix, index=all_set, columns=all_set, diff --git a/stlearn/tl/cci/edges.py b/stlearn/tl/cci/edges.py index ec710041..e8ad251c 100644 --- a/stlearn/tl/cci/edges.py +++ b/stlearn/tl/cci/edges.py @@ -2,10 +2,8 @@ import numpy as np from anndata import AnnData -from numba import njit -from numba.typed import List -from .neighbours import get_neighbourhoods +from .neighbours import stored_neighbours def get_edges( @@ -13,13 +11,13 @@ def get_edges( l_bool: np.ndarray, r_bool: np.ndarray, sig_bool: np.ndarray, -): - """Gets a list edges representing significant interactions. +) -> list[list[tuple[str, str]]]: + """Gets the edges representing significant interactions, to draw. Parameters ---------- adata : AnnData - Annotated data object containing spatial transcriptomics data. + The data st.tl.cci.run stored its neighbours in, obsp["spot_neighbours"]. l_bool : np.ndarray of bool, shape (n_spots,) Boolean array indicating spots where the ligand is expressed. r_bool : np.ndarray of bool, shape (n_spots,) @@ -29,123 +27,27 @@ def get_edges( Returns ------- - edge_list_unique: list> Either a list of tuples (directed), or - list of sets (undirected), indicating unique significant - interactions between spots. + [forward_edges, reverse_edges]: forward edges run from each significant spot + expressing the ligand to its neighbours expressing the receptor, reverse + edges from each significant spot expressing the receptor to its neighbours + expressing the ligand. Each is a list of (barcode, barcode), with each + pair of spots once, the way round it was first found. """ - # Getting the neighbourhoods # - neighbourhood_bcs, neighbourhood_indices = get_neighbourhoods(adata) - - # Getting the edges to draw in-between # - l_spot_indices = np.where(np.logical_and(l_bool, sig_bool))[0] - r_spot_indices = np.where(np.logical_and(r_bool, sig_bool))[0] - - # To keep the get_between_spot_edge_array function happy # - cell_data = np.ones((1, len(sig_bool)))[0, :].astype(np.float64) - - # Retrieving the edges # - gene_bools = [r_bool, l_bool] - edge_list = [] - for i, spot_indices in enumerate([l_spot_indices, r_spot_indices]): - # Subsetting to the relevant neighbourhoods # - neigh_bcs = List() - neigh_indices = List() - for j in spot_indices: - neigh_bcs.append(neighbourhood_bcs[j]) - neigh_indices.append(neighbourhood_indices[j]) - - # Getting the edges in this neighbourhood # - edges = init_edge_list(neighbourhood_bcs) # Note this has 1 pseudo edge - if len(spot_indices) != 0: - get_between_spot_edge_array( - edges, - neigh_bcs, - neigh_indices, - gene_bools[i], - cell_data, - ) - edge_list.append(edges[1:]) - - return edge_list - - -@njit -def init_edge_list(neighbourhood_bcs): - """ Initialises the edge-list in a way which ensures consistent typing to \ - prevent errors with Numba. - """ - - # Initialising the edge list # - edge_list = List() - # This ensures consistent typing # - for i in range(len(neighbourhood_bcs)): - if len(neighbourhood_bcs[i][1]) > 0: - edge_list.append((neighbourhood_bcs[i][0], neighbourhood_bcs[i][1][0])) - break - return edge_list - - -@njit -def get_between_spot_edge_array( - edge_list: List, - neighbourhood_bcs: List, - neighbourhood_indices: List, - neigh_bool: np.ndarray, - cell_data: np.ndarray, - cutoff: float = 0, -): - """ Populates edge_list with edges linking spots with a valid neighbour \ - of a given cell type. Validity of neighbour determined by neigh_bool, \ - which can indicate whether the neighbour expresses a certain ligand \ - or receptor. neighbourhood_bcs and neighbourhood_indices hold, per - spot, (barcode, neighbour barcodes) and (index, neighbour indices); - cell_data is spots x cell types, and a neighbour counts as a cell - type when its value is above cutoff. - """ - edge_starts = List() - edge_ends = List() - for i in range(len(neighbourhood_bcs)): - bcs, indices = neighbourhood_bcs[i], neighbourhood_indices[i] - spot_bc, neigh_bcs = bcs - neigh_indices = indices[1] - # Subset the neighbours to only those fitting indicated criteria # - neigh_bcs = neigh_bcs[neigh_bool[neigh_indices]] - neigh_indices = neigh_indices[neigh_bool[neigh_indices]] - - if len(neigh_indices) == 0: # No cases where neighbours meet criteria - continue # Don't add any interactions for this neighbourhood - - # Note that can keep all by inputting cell_data with all 1's # - interact_neigh_bool = cell_data[neigh_indices] > cutoff - - # Retrieving the barcodes of the interacting neighbours # - interact_neigh_bcs = neigh_bcs[interact_neigh_bool] - for interact_neigh_bc in interact_neigh_bcs: - edge_starts.append(spot_bc) - edge_ends.append(interact_neigh_bc) - - # Getting the unique edges # - if len(edge_starts) > 0: - add_unique_edges(edge_list, edge_starts, edge_ends) - - -@njit -def add_unique_edges(edge_list, edge_starts, edge_ends): - """ Adds the unique edges to the given edge list. \ - Complicated in order to satisfy Numba compilation in no-python mode. - """ - n_edges = len(edge_starts) - - # Adding the unique edges # - edge_added = np.zeros((1, len(edge_starts)))[0, :] == 1 - for i in range(n_edges): - if not edge_added[i]: - edge_start, edge_end = edge_starts[i], edge_ends[i] - edge_list.append((edge_start, edge_end)) - for j in range(i, n_edges): - edge_startj, edge_endj = edge_starts[j], edge_ends[j] - # Direction doesn't matter # - if (edge_start == edge_startj and edge_end == edge_endj) or ( - edge_end == edge_startj and edge_start == edge_endj - ): - edge_added[j] = True + graph = stored_neighbours(adata) + names = adata.obs_names.to_numpy().astype(str) + edge_lists = [] + for sources, targets in ((l_bool, r_bool), (r_bool, l_bool)): + spots = np.flatnonzero(np.asarray(sources) & np.asarray(sig_bool)) + rows = graph[spots] + starts = np.repeat(spots, np.diff(rows.indptr)) + ends = rows.indices[np.asarray(targets)[rows.indices]] + starts = starts[np.asarray(targets)[rows.indices]] + # Each pair of spots once, in the order found: direction doesn't matter. + pairs = np.minimum(starts, ends).astype(np.int64) * graph.shape[0] + pairs += np.maximum(starts, ends) + _, first = np.unique(pairs, return_index=True) + first.sort() + edge_lists.append( + list(zip(names[starts[first]], names[ends[first]], strict=True)) + ) + return edge_lists diff --git a/stlearn/tl/cci/het.py b/stlearn/tl/cci/het.py index df537299..1406905e 100644 --- a/stlearn/tl/cci/het.py +++ b/stlearn/tl/cci/het.py @@ -4,10 +4,11 @@ import numpy as np import pandas as pd +import scipy from anndata import AnnData from numba import njit, prange -from .neighbours import calc_distance, neighbourhood_matrix +from .neighbours import calc_distance, calc_neighbours def count( @@ -17,6 +18,7 @@ def count( verbose: bool = True, distance: float | None = None, cell_prop_cutoff: float = 0.2, + neighbours: scipy.sparse.csr_matrix | None = None, ) -> None: """Counts the cell types present in each spot's neighbourhood. @@ -39,6 +41,9 @@ def count( distance=0 means within spot cell_prop_cutoff: float A cell type is present in a spot when its proportion is above this. + neighbours: scipy.sparse.csr_matrix + The spots-by-spots graph calc_neighbours makes, when the caller already + has it (st.tl.cci.run passes its own); otherwise it is found from distance. Returns ------- @@ -53,11 +58,15 @@ def count( """ present = label_proportions(adata, use_label).to_numpy() > cell_prop_cutoff - # between spot if distance != 0: - # The default neighbour distance if not given; distance=0 is within-spot. - neighbourhoods = neighbourhood_matrix(adata, calc_distance(adata, distance)) - # The number of cell types present anywhere in each neighbourhood. + if neighbours is None: + neighbours = calc_neighbours( + adata, calc_distance(adata, distance), verbose=False + ) + spots = scipy.sparse.identity( + neighbours.shape[0], dtype=neighbours.dtype, format="csr" + ) + neighbourhoods = neighbours.maximum(spots) adata.obsm[use_het] = ((neighbourhoods @ present.astype(np.int32)) > 0).sum(1) # within spot @@ -129,7 +138,8 @@ def get_interaction_pvals( int_matrix, n_perms, cell_data, - neighbourhood_indices, + neigh_indptr, + neigh_indices, all_set, sig_bool, L_bool, @@ -176,7 +186,8 @@ def get_interaction_pvals( # Calculating interactions for permuted labels # perm_matrix = get_interaction_matrix( perm_data, - neighbourhood_indices, + neigh_indptr, + neigh_indices, all_set, sig_bool, L_bool, @@ -204,7 +215,8 @@ def get_interaction_pvals( @njit def get_interaction_matrix( cell_data, - neighbourhood_indices, + neigh_indptr, + neigh_indices, all_set, sig_bool, L_bool, @@ -228,10 +240,8 @@ def get_interaction_matrix( if cell_data[i, t1] <= cell_prop_cutoff: continue - neighs = neighbourhood_indices[i][1] - - for k in range(len(neighs)): - n_idx = neighs[k] + for k in range(neigh_indptr[i], neigh_indptr[i + 1]): + n_idx = neigh_indices[k] if cell_data[n_idx, t2] > cell_prop_cutoff: valid = False if L_bool[i] and R_bool[n_idx]: diff --git a/stlearn/tl/cci/neighbours.py b/stlearn/tl/cci/neighbours.py index 1169b801..6626ef96 100644 --- a/stlearn/tl/cci/neighbours.py +++ b/stlearn/tl/cci/neighbours.py @@ -1,12 +1,7 @@ -"""Spot neighbourhoods: computing them, and reading back the ones st.tl.cci.run -stored in adata.obsm. -""" - import numpy as np import scipy.sparse import scipy.spatial as spatial from anndata import AnnData -from numba.typed import List from ...utils import _image_scalef @@ -37,9 +32,8 @@ def calc_distance(adata: AnnData, distance: float | None) -> float: def calc_neighbours( adata: AnnData, distance: float | None = None, - index: bool = True, verbose: bool = True, -) -> List: +) -> scipy.sparse.csr_matrix: """Each spot's neighbours: the other spots within distance, or the spot itself when distance=0 (within-spot mode). @@ -50,140 +44,61 @@ def calc_neighbours( distance (float): Distance to determine the neighbours (default: closest), distance=0 means within spot - index (bool): - Indicates whether to return neighbours as indices to other spots or names of - other spots. verbose (bool): Display debugging information Returns ------- - neighbours (numba.typed.List): - List of np.array's indicating neighbours by indices for each spot. + neighbours (scipy.sparse.csr_matrix): + Spots by spots: row i's indices are spot i's neighbours, in the order + the tree finds them, as st.tl.cci.run stores in obsp["spot_neighbours"]. - Raises - ------ - ValueError - If no spot has a neighbour at this distance. """ if verbose: print("Calculating neighbours...") - # get neighbour spots for each spot according to the specified distance - coor = adata.obs[["imagerow", "imagecol"]].to_numpy() - point_tree = spatial.cKDTree(coor) - # One query for all spots, across all CPUs. - found = point_tree.query_ball_point( - coor, - distance, - workers=-1, - return_sorted=False, - ) - neighbours = [] - for i, spot in enumerate(adata.obs_names): - if distance == 0: - neighbours.append(np.array([i if index else spot])) - else: - n_index = found[i] - if index: - n_index = np.array(n_index, dtype=np.int32) - neighbours.append(n_index[n_index != i]) - else: - n_spots = adata.obs_names[n_index] - neighbours.append(n_spots[n_spots != spot]) - typed_neighs = List() - [typed_neighs.append(neigh) for neigh in neighbours] + n_spots = adata.n_obs + if distance == 0: + # Within-spot mode: each spot is its only neighbour. + indices = np.arange(n_spots, dtype=np.int32) + indptr = np.arange(n_spots + 1, dtype=np.int64) + else: + coor = adata.obs[["imagerow", "imagecol"]].to_numpy() + # One query for all spots, across all CPUs. + found = spatial.cKDTree(coor).query_ball_point( + coor, + distance, + workers=-1, + return_sorted=False, + ) + # Each row's spots as found, the spot itself left out. + lengths = np.fromiter(map(len, found), np.int64, count=n_spots) + rows = np.repeat(np.arange(n_spots), lengths) + flat = np.concatenate(found).astype(np.int32) + other = flat != rows + indices = flat[other] + indptr = np.zeros(n_spots + 1, np.int64) + np.cumsum(np.bincount(rows[other], minlength=n_spots), out=indptr[1:]) + + graph = scipy.sparse.csr_matrix( + (np.ones(len(indices), np.int8), indices, indptr), + shape=(n_spots, n_spots), + ) + n_neighs = np.diff(graph.indptr) - n_neighs = np.array([len(neigh) for neigh in neighbours]) if verbose: print( - f"{len(np.where(n_neighs == 0)[0])} spots with no neighbours, " + f"{np.count_nonzero(n_neighs == 0)} spots with no neighbours, " f"{int(np.median(n_neighs))} median spot neighbours.", ) + return graph + - if np.all(n_neighs == 0): +def stored_neighbours(adata) -> scipy.sparse.csr_matrix: + """The graph st.tl.cci.run stored in obsp["spot_neighbours"], as CSR: row i's + indices are spot i's neighbours. Raises ValueError when there is none.""" + if "spot_neighbours" not in adata.obsp: raise ValueError( - "All spots have no neighbours at current distance," - " set distance to higher value, or distance=0 for " - "within-spot mode.", + "No neighbours stored in adata.obsp['spot_neighbours']: call st.tl.cci.run." ) - return typed_neighs - - -def neighbourhood_matrix(adata: AnnData, distance: float) -> scipy.sparse.csr_matrix: - """Spots by spots; row i marks spot i's neighbourhood within distance, the - spot itself included.""" - coor = adata.obs[["imagerow", "imagecol"]].to_numpy() - found = spatial.cKDTree(coor).query_ball_point(coor, distance, workers=-1) - rows = np.repeat(np.arange(len(found)), list(map(len, found))) - return scipy.sparse.csr_matrix( - (np.ones(len(rows), np.int32), (rows, np.concatenate(found))), - shape=(len(found), len(found)), - ) - - -def neighbourhoods_from_barcodes( - spot_bcs: np.ndarray, - spot_neigh_bcs: np.ndarray, - str_dtype: str, -): - """(barcode, neighbour barcodes) and (index, neighbour indices) per spot, from - the comma-joined neighbour barcodes st.tl.cci.run stores. Neighbours no longer - in the data, after subsetting, are skipped. - """ - neighbourhood_bcs, neighbourhood_indices = [], [] - spot_bc_to_index = {bc: i for i, bc in enumerate(spot_bcs)} - - for i in range(spot_neigh_bcs.shape[0]): - neigh_bcs = np.array(spot_neigh_bcs[i, :][0].split(",")) - neigh_bcs = neigh_bcs[neigh_bcs != ""] - - neigh_bcs_array, neigh_indices = [], [] - for neigh_bc in neigh_bcs: - if neigh_bc in spot_bc_to_index: - neigh_bcs_array.append(neigh_bc) - neigh_indices.append(spot_bc_to_index[neigh_bc]) - - neigh_bcs_array = np.array(neigh_bcs_array, dtype=str_dtype) - neigh_indices = np.array(neigh_indices, dtype=np.int64) - - neighbourhood_indices.append((i, neigh_indices)) - neighbourhood_bcs.append((spot_bcs[i], neigh_bcs_array)) - - return List(neighbourhood_bcs), List(neighbourhood_indices) - - -def get_neighbourhoods(adata): - """Gets the neighbourhood information stored by st.tl.cci.run.""" - - # Old stlearn version where didn't store neighbourhood barcodes, not good - # for anndata subsetting!! - if "spot_neigh_bcs" not in adata.obsm: - # Determining the neighbour spots used for significance testing # - neighbours = List() - for i in range(adata.obsm["spot_neighbours"].shape[0]): - neighs = np.array(adata.obsm["spot_neighbours"].values[i, :][0].split(",")) - neighs = neighs[neighs != ""].astype(int) - neighbours.append(neighs) - - # Getting the neighbourhood information # - neighbourhood_bcs = List() - neighbourhood_indices = List() - spot_bcs = adata.obs_names.to_numpy().astype(str) - for spot_i in range(len(spot_bcs)): - neighbourhood_indices.append((spot_i, neighbours[spot_i])) - neighbourhood_bcs.append((spot_bcs[spot_i], spot_bcs[neighbours[spot_i]])) - - return neighbourhood_bcs, neighbourhood_indices - else: # Newer version - spot_bcs = adata.obs_names.to_numpy().astype(str) - spot_neigh_bcs = adata.obsm["spot_neigh_bcs"].to_numpy().astype(str) - - max_len = 0 - for bc in spot_bcs: - if len(bc) > max_len: - max_len = len(bc) - - str_dtype = f" Expression: if self.columns is not None: return self.columns[:, genes].tocsr() expr = self.adata[:, genes].X - if sparse.issparse(expr): - return sparse.csr_matrix(expr) + if scipy.sparse.issparse(expr): + return scipy.sparse.csr_matrix(expr) return np.asfortranarray(expr) def score( @@ -70,7 +74,7 @@ def score( filter_pairs: bool = True, spot_indices: npt.NDArray[np.int32] | None = None, sparse_result: bool = False, - ) -> tuple[npt.NDArray[np.float64] | sparse.csc_matrix, LRPairs]: + ) -> tuple[npt.NDArray[np.float64] | scipy.sparse.csc_matrix, LRPairs]: """The scores of the LR pairs, and the pairs scored; see get_lrs_scores.""" genes = np.column_stack((lrs.ligands(), lrs.receptors())) pairs = self.adata.var_names.get_indexer(genes.ravel()).reshape(genes.shape) @@ -93,7 +97,7 @@ def score_indices( min_expr: float, spot_indices: npt.NDArray[np.int32] | None = None, sparse_result: bool = False, - ) -> npt.NDArray[np.float64] | sparse.csc_matrix: + ) -> npt.NDArray[np.float64] | scipy.sparse.csc_matrix: """The scores of pairs given as positions in adata.var_names, so the thousands of random pairs of a background are never named. @@ -127,13 +131,13 @@ def score_indices( def get_lrs_scores( adata: AnnData, lrs: LRPairs, - neighbours: np.ndarray, + neighbours: scipy.sparse.csr_matrix, het_vals: np.ndarray, min_expr: float, filter_pairs: bool = True, spot_indices: npt.NDArray[np.int32] | None = None, sparse_result: bool = False, -) -> tuple[npt.NDArray[np.float64] | sparse.csc_matrix, LRPairs]: +) -> tuple[npt.NDArray[np.float64] | scipy.sparse.csc_matrix, LRPairs]: """Gets the scores for the indicated set of LR pairs & the heterogeneity values. For repeated calls on the same sample, make an LRScorer once and call its @@ -145,13 +149,13 @@ def get_lrs_scores( See run() doc-string. lrs: LRPairs See run() doc-string. - neighbours: np.ndarray - Array of arrays with indices specifying neighbours of each spot. + neighbours: scipy.sparse.csr_matrix + CSR Matrix with indices specifying neighbours of each spot. het_vals: np.ndarray Cell heterogeneity counts per spot. min_expr: float Minimum gene expression of either L or R for spot to be considered to - have reasonable score. + have a reasonable score. filter_pairs: bool Whether to drop the pairs with a gene not in adata.var_names; if not, such a pair raises KeyError. @@ -198,26 +202,6 @@ def direction_score(at_spot: float, mean_partner: float, min_expr: float) -> flo return at_spot * (mean_partner > min_expr) + (at_spot > min_expr) * mean_partner -@njit(cache=True) -def flatten_neighbours(neighbours: List) -> tuple[np.ndarray, np.ndarray]: - """The neighbours as one array, spot i's at flat[starts[i] : starts[i + 1]].""" - starts = np.zeros(len(neighbours) + 1, np.int64) - for i in range(len(neighbours)): - starts[i + 1] = starts[i] + len(neighbours[i]) - flat = np.empty(starts[-1], np.int64) - for i in range(len(neighbours)): - flat[starts[i] : starts[i + 1]] = neighbours[i] - return starts, flat - - -def adjacency(neighbours: List) -> sparse.csr_matrix: - """The neighbours as a spots-by-spots CSR matrix: row i's indices are spot - i's neighbours, in their order, so rows of spots slice out cheaply.""" - starts, flat = flatten_neighbours(neighbours) - n = len(starts) - 1 - return sparse.csr_matrix((np.ones(len(flat), np.int8), flat, starts), shape=(n, n)) - - @njit(parallel=True) def neighbourhood_means( expr: np.ndarray, @@ -233,7 +217,8 @@ def neighbourhood_means( expr: np.ndarray Spots*Genes expression, column-major. starts, flat: np.ndarray - The neighbours from flatten_neighbours, as rows of expr. + A CSR matrix's indptr and indices: spot i's neighbours, as rows of expr, + are flat[starts[i] : starts[i + 1]]. spot_indices: npt.NDArray[np.int32] The spots to compute the means at, as positions in starts. Returns @@ -306,40 +291,47 @@ def _dense_rows(expr: Expression, rows: np.ndarray) -> np.ndarray: """These rows of expr, dense and column-major.""" taken = expr[rows] return ( - taken.toarray(order="F") if sparse.issparse(taken) else np.asfortranarray(taken) + taken.toarray(order="F") + if scipy.sparse.issparse(taken) + else np.asfortranarray(taken) ) def score_spots( expr: Expression, pairs: np.ndarray, - adjacency: sparse.csr_matrix, + adjacency: scipy.sparse.csr_matrix, het_vals: np.ndarray, min_expr: float, spot_indices: npt.NDArray[np.int32], sparse_result: bool = False, -) -> npt.NDArray[np.float64] | sparse.csc_matrix: +) -> npt.NDArray[np.float64] | scipy.sparse.csc_matrix: """Scores the pairs at the spots, a chunk of spots at a time, each chunk's rows of expr dense and its neighbourhood means held only while it is scored. - See get_scores for the arguments and result; with sparse_result the scores - come back as get_lrs_scores describes. """ het_vals = np.asarray(het_vals) n_pairs, n_spots = pairs.shape[0], len(spot_indices) chunk = spots_per_chunk(expr.shape[1], n_pairs) - pieces: list[sparse.csr_matrix] = [] + pieces: list[scipy.sparse.csr_matrix] = [] spot_scores = np.empty((n_pairs, 0 if sparse_result else n_spots), np.float64) + + n_rows = expr.shape[0] + mask = np.zeros(n_rows, np.bool_) # the chunk's rows; reset after each chunk + position = np.empty(n_rows, np.int64) # a row's index in the chunk's block + for start in range(0, n_spots, chunk): spots = spot_indices[start : start + chunk] neighbours = adjacency[spots] - rows, local = np.unique( - np.concatenate([spots, neighbours.indices]), return_inverse=True - ) + mask[spots] = True + mask[neighbours.indices] = True + rows = np.flatnonzero(mask) + mask[rows] = False + position[rows] = np.arange(len(rows)) block = _dense_rows(expr, rows) means = neighbourhood_means( block, neighbours.indptr.astype(np.int64), - local[len(spots) :].astype(np.int64), + position[neighbours.indices], np.arange(len(spots), dtype=np.int32), ) scores = ( @@ -347,50 +339,19 @@ def score_spots( if sparse_result else spot_scores[:, start : start + chunk] ) - block_spots = local[: len(spots)].astype(np.int32) + block_spots = position[spots].astype(np.int32) score_pairs(block, means, pairs, het_vals[rows], min_expr, block_spots, scores) if sparse_result: by_spot = scores.T above = by_spot > 0 pieces.append( - sparse.csr_matrix((by_spot[above], np.nonzero(above)), by_spot.shape) + scipy.sparse.csr_matrix( + (by_spot[above], np.nonzero(above)), by_spot.shape + ) ) del scores, by_spot, above # before the next chunk's are made if not sparse_result: return spot_scores.T if not pieces: - return sparse.csc_matrix((n_spots, n_pairs), dtype=np.float64) - return sparse.vstack(pieces, format="csc") - - -def get_scores( - expr: np.ndarray, - pairs: np.ndarray, - neighbours: List, - het_vals: np.ndarray, - min_expr: float, - spot_indices: npt.NDArray[np.int32], -) -> npt.NDArray[np.float64]: - """Calculates the scores. - Parameters - ---------- - expr: np.ndarray - Spots*Genes expression of the genes in the pairs. - pairs: np.ndarray - Pairs*2 column indices into expr, as (ligand, receptor). - neighbours: List - List of np.array's indicating neighbours by indices for each spot. - het_vals: np.ndarray - Spots*Het counts - min_expr: float - Minimum expression for gene to be considered expressed. - spot_indices: npt.NDArray[np.int32] - Subset of spots to score, given as their integer row positions. - Returns - ------- - spot_scores: np.ndarray - Spots*LR pair of the LR scores per spot. - """ - return score_spots( - expr, pairs, adjacency(neighbours), het_vals, min_expr, spot_indices - ) + return scipy.sparse.csc_matrix((n_spots, n_pairs), dtype=np.float64) + return scipy.sparse.vstack(pieces, format="csc") diff --git a/tests/acceptance/run_cci.py b/tests/acceptance/run_cci.py index 5efa611d..a89ba28b 100644 --- a/tests/acceptance/run_cci.py +++ b/tests/acceptance/run_cci.py @@ -5,6 +5,10 @@ python -m tests.acceptance.run_cci SAMPLE.h5ad RESULT.h5ad [--use-label] +The acceptance test runs it for every scenario: + + pytest -m acceptance --acceptance-bins 500000 --acceptance-baseline=write + numba's functions are compiled on a small slice first, so the timed run is steady state. Results are stored sparse, as their difference from the value most spots have (score 0, p-value 1): the entries of the LR table's sparse @@ -40,7 +44,8 @@ "neg_log10_p_adjs": 0.0, "lr_sig_scores": 0.0, } -NEIGHBOURS = ["spot_neighbours", "spot_neigh_bcs"] +# The neighbours, obsp["spot_neighbours"], kept in the result as its CSR arrays. +NEIGHBOURS = ["neighbours_indptr", "neighbours_indices"] # Phases of st.tl.cci.run, timed by wrapping the functions it calls. PHASES = { "neighbours": ("stlearn.tl.cci.analysis", "calc_neighbours"), @@ -185,8 +190,9 @@ def main() -> None: result = anndata.AnnData(obs=pandas.DataFrame(index=adata.obs_names)) for key, usual in RESULTS.items(): result.obsm[key] = shifted(lr_results.layers[key], usual) - for key in NEIGHBOURS: - result.obs[key] = adata.obsm[key].iloc[:, 0].to_numpy() + graph = scipy.sparse.csr_matrix(adata.obsp["spot_neighbours"]) + result.uns["neighbours_indptr"] = graph.indptr + result.uns["neighbours_indices"] = graph.indices if args.use_label: result.obs["cci_het"] = numpy.asarray(adata.obsm["cci_het"]).ravel() # The layout of baselines written before the LR table. diff --git a/tests/acceptance/test_cci_acceptance.py b/tests/acceptance/test_cci_acceptance.py index cc6ef935..eaf7a6ee 100644 --- a/tests/acceptance/test_cci_acceptance.py +++ b/tests/acceptance/test_cci_acceptance.py @@ -231,7 +231,13 @@ def differences(baseline: anndata.AnnData, current: anndata.AnnData, atol: float current_features = current.uns["lrfeatures"].loc[features.index, features.columns] if not np.allclose(features, current_features, rtol=0, atol=atol, equal_nan=True): problems.append("lrfeatures differ") - for key in [*NEIGHBOURS, "cci_het"]: + for key in NEIGHBOURS: + before, now = baseline.uns.get(key), current.uns.get(key) + if (before is None) != (now is None): + problems.append(f"{key}: in only one run") + elif before is not None and not np.array_equal(before, now): + problems.append(f"{key}: the neighbourhoods differ") + for key in ["cci_het"]: in_baseline, in_current = key in baseline.obs, key in current.obs if in_baseline != in_current: problems.append(f"{key}: in only one run") diff --git a/tests/io/test_cci_round_trip.py b/tests/io/test_cci_round_trip.py index d8f4ee40..16f4f316 100644 --- a/tests/io/test_cci_round_trip.py +++ b/tests/io/test_cci_round_trip.py @@ -109,7 +109,7 @@ def test_every_key(self): pd.testing.assert_frame_equal( getattr(read[name], slot), getattr(table, slot) ) - for slot in ("obsm", "layers", "uns"): + for slot in ("obsm", "obsp", "layers", "uns"): written, back = getattr(table, slot), getattr(read[name], slot) self.assertEqual(sorted(back), sorted(written), msg=f"{name}.{slot}") for key in written: diff --git a/tests/test_utils.py b/tests/test_utils.py index 923fe65e..bd873787 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -61,7 +61,7 @@ def test_cci_neighbours(self): ): distance = st.utils.microns_to_pixels(adata, distance_um) neighbours = calc_neighbours(adata, distance, verbose=False) - self.assertEqual(len(neighbours[centre]), n_expected) + self.assertEqual(neighbours[centre].nnz, n_expected) def test_arrays_and_round_trip(self): adata = _spots(microns_per_pixel=0.2125, scalef=0.1) diff --git a/tests/tl/test_cci.py b/tests/tl/test_cci.py index 9b6fded5..6e12855c 100644 --- a/tests/tl/test_cci.py +++ b/tests/tl/test_cci.py @@ -1,24 +1,17 @@ #!/usr/bin/env python -"""Tests for the ligand-receptor and cell-cell interaction functions in `stlearn`. - -TestCCI holds unit tests of the neighbourhood and counting helpers on a toy -seven-spot neighbourhood. TestRunCCILabels runs st.tl.cci.run_cci on the test data -with cell-type labels stored as a categorical and as plain strings. -""" +"""Neighbourhood and counting helpers on a toy seven-spot neighbourhood.""" import unittest -from typing import ClassVar +from typing import Any import numpy as np import pandas as pd +import scipy.sparse from anndata import AnnData -from numba.typed import List -import stlearn as st -import stlearn.tl.cci.edges as edges import stlearn.tl.cci.het as het -from tests.utils import read_test_data +from stlearn.tl.cci import get_edges # Per line - which cells to annotate CELL_TYPE_ANNOTATIONS = ["CT1", "CT2", "CT3", "CT2", "CT1", "CT3", "CT2"] @@ -27,35 +20,9 @@ CELL_TYPE_LABELS = np.array(["CT1", "CT2", "CT3"]) -def read_test_data_with_lr_pairs(n_pairs: int = 20): - """Test data that ligand-receptor analysis can run on. - - The 100 genes in the test data contain no pair from the ligand-receptor - database and the matrix is scaled, so `st.tl.cci.run` has nothing to test. This - renames the first `2 * n_pairs` genes to database pairs and clips negative values - to zero. It exercises the code path; the scores carry no biological meaning. - """ - adata = read_test_data() - lrs = st.tl.cci.lr_pairs("connectomeDB2020") - - genes: list[str] = [] - for ligand, receptor in zip( - lrs.ligands().tolist(), lrs.receptors().tolist(), strict=True - ): - if ligand != receptor and ligand not in genes and receptor not in genes: - genes += [ligand, receptor] - if len(genes) == 2 * n_pairs: - break - - names = list(adata.var_names) - names[: len(genes)] = genes - adata.var_names = names - adata.X = np.clip(adata.X, 0, None) - return adata, lrs - - class TestCCI(unittest.TestCase): """Tests for `stlearn` CCI capability.""" + all_edges: list[Any] def setUp(self) -> None: """Setup some basic test-cases as sanity checks.""" @@ -75,141 +42,78 @@ def setUp(self) -> None: "F": ["G", "A", "E"], "G": ["F", "A", "B"], } - neighbourhood_bcs = List() - neighbourhood_indices = List() - for i, bc in enumerate(spot_bcs): - neigh_bcs = np.array(neigh_dict[bc]) - neigh_indices = np.array( - [np.where(spot_bcs == n_bc)[0][0] for n_bc in neigh_bcs] - ) - neighbourhood_bcs.append((bc, neigh_bcs)) - neighbourhood_indices.append((i, neigh_indices)) - - self.neighbourhood_bcs = neighbourhood_bcs - self.neighbourhood_indices = neighbourhood_indices - self.neigh_dict = neigh_dict - - # Important, granular tests related to LR scoring - - # Important, granular tests related to CCI counting - def test_edge_retrieval_basic(self): - """ Basic test of functionality to retrieve edges via \ - get_between_spot_edge_array. - """ - neighbourhood_bcs = self.neighbourhood_bcs - neighbourhood_indices = self.neighbourhood_indices - - # Initialising the edge list # - edge_list = edges.init_edge_list(neighbourhood_bcs) - - # Basic case, should populate with all edges - neigh_bool = np.array([True] * len(neighbourhood_bcs)) - cell_data = np.array([1] * len(neighbourhood_bcs), dtype=np.float64) - edges.get_between_spot_edge_array( - edge_list, neighbourhood_bcs, neighbourhood_indices, neigh_bool, cell_data - ) - edge_list = edge_list[1:] # Remove initialisation edge - - # Getting all the edges in the neighbourhood # - all_edges = [] - for bc in self.neigh_dict: - neigh_bcs = self.neigh_dict[bc] - for neigh_bc in neigh_bcs: - if (bc, neigh_bc) not in all_edges and (neigh_bc, bc) not in all_edges: - all_edges.append((bc, neigh_bc)) - n_edges = len(all_edges) - - self.assertEqual(len(edge_list), n_edges) - self.assertTrue( - np.all([edge in all_edges or edge[::-1] in all_edges for edge in edge_list]) - ) - - # Some neighbours not valid but no effect on edge list - # No effect since though not a valid neighbour, still a valid spot # - edge_list = edges.init_edge_list(neighbourhood_bcs) - invalid_neighs = ["B", "E"] - neigh_bool = np.array([bc not in invalid_neighs for bc in self.neigh_dict]) - edges.get_between_spot_edge_array( - edge_list, neighbourhood_bcs, neighbourhood_indices, neigh_bool, cell_data + position = {bc: i for i, bc in enumerate(spot_bcs)} + lengths = [len(neigh_dict[bc]) for bc in spot_bcs] + self.indptr = np.concatenate([[0], np.cumsum(lengths)]).astype(np.int64) + self.indices = np.array( + [position[n] for bc in spot_bcs for n in neigh_dict[bc]], dtype=np.int64 ) - edge_list = edge_list[1:] # Remove initialisation edge - - self.assertEqual(len(edge_list), n_edges) - self.assertTrue( - np.all([edge in all_edges or edge[::-1] in all_edges for edge in edge_list]) + self.adata = AnnData( + X=np.zeros((len(spot_bcs), 1)), obs=pd.DataFrame(index=spot_bcs) ) - - # Some neighbours not valid, effects the edge list - # Two neighbouring spots no longer valid neighbours # - edge_list = edges.init_edge_list(neighbourhood_bcs) - invalid_neighs = ["B", "C"] - neigh_bool = np.array([bc not in invalid_neighs for bc in self.neigh_dict]) - edges.get_between_spot_edge_array( - edge_list, neighbourhood_bcs, neighbourhood_indices, neigh_bool, cell_data + self.adata.obsp["spot_neighbours"] = scipy.sparse.csr_matrix( + (np.ones(len(self.indices), np.int8), self.indices, self.indptr), + shape=(len(spot_bcs), len(spot_bcs)), ) - edge_list = edge_list[1:] # Remove initialisation edge - - sub_edges1 = [ - edge for edge in all_edges if edge not in [("B", "C"), ("C", "B")] - ] - n_subedges = len(sub_edges1) - - self.assertEqual(len(edge_list), n_subedges) - self.assertTrue( - np.all( - [edge in sub_edges1 or edge[::-1] in sub_edges1 for edge in edge_list] - ) - ) - - # Middle spot not neighbour, cell type, or spot of interest - # Removing the centre-spot as being the cell type of interest # - neigh_bool = np.array([True] * len(neighbourhood_bcs)) - neigh_bool[0] = False - cell_data = np.array([1] * len(neighbourhood_bcs), dtype=np.float64) - cell_data[0] = 0 - neigh_bcs = neighbourhood_bcs[1:] - neigh_indices = neighbourhood_indices[1:] - - edge_list = edges.init_edge_list(neighbourhood_bcs) - edges.get_between_spot_edge_array( - edge_list, neigh_bcs, neigh_indices, neigh_bool, cell_data - ) - edge_list = edge_list[1:] # Remove initialisation edge + self.spot_bcs = spot_bcs + self.neigh_dict = neigh_dict + # Every edge of the neighbourhood, each pair of spots once. + self.all_edges = [] + for bc, neigh_bcs in neigh_dict.items(): + for neigh_bc in neigh_bcs: + if (neigh_bc, bc) not in self.all_edges: + self.all_edges.append((bc, neigh_bc)) - sub_edges2 = [edge for edge in all_edges if "A" not in edge] - n_subedges2 = len(sub_edges2) + def _bools(self, excluded=()): + return np.array([bc not in excluded for bc in self.spot_bcs]) - self.assertEqual(len(edge_list), n_subedges2) - self.assertTrue( - np.all( - [edge in sub_edges2 or edge[::-1] in sub_edges2 for edge in edge_list] - ) - ) + def assertSameEdges(self, edges, expected): + """The same pairs of spots, either way round.""" + self.assertEqual(len(edges), len(expected)) + self.assertTrue(all(e in expected or e[::-1] in expected for e in edges)) - # Corner spot valid neighbour, not cell type, not spot of interest - neigh_bool = np.array([True] * len(neighbourhood_bcs)) - cell_data = np.array([1] * len(neighbourhood_bcs), dtype=np.float64) - cell_data[1] = 0 - neigh_bcs, neigh_indices = List(), List() - for i in range(len(cell_data)): - if i != 1: - neigh_bcs.append(neighbourhood_bcs[i]) - neigh_indices.append(neighbourhood_indices[i]) + # Important, granular tests related to LR scoring - edge_list = edges.init_edge_list(neighbourhood_bcs) - edges.get_between_spot_edge_array( - edge_list, neigh_bcs, neigh_indices, neigh_bool, cell_data + # Important, granular tests related to CCI counting + def test_edges(self): + """get_edges draws, forward, from each significant spot with the ligand + to its neighbours with the receptor, and reverse the other way round, + each pair of spots once.""" + everywhere = self._bools() + forward, reverse = get_edges(self.adata, everywhere, everywhere, everywhere) + self.assertSameEdges(forward, self.all_edges) + self.assertSameEdges(reverse, self.all_edges) + self.assertEqual(len(self.all_edges), 12) + # Each edge once, the way round it was first found: A first, then B, ... + self.assertEqual(forward[:6], [("A", n) for n in self.neigh_dict["A"]]) + self.assertEqual(forward[6], ("B", "G")) + + # B and E without the receptor: their edges still come from their own + # neighbourhoods, where they are the significant ligand spot. + forward, _ = get_edges(self.adata, everywhere, self._bools("BE"), everywhere) + self.assertSameEdges(forward, self.all_edges) + + # Neighbouring B and C without the receptor: only their edge is lost. + forward, _ = get_edges(self.adata, everywhere, self._bools("BC"), everywhere) + self.assertSameEdges( + forward, [e for e in self.all_edges if set(e) != {"B", "C"}] ) - edge_list = edge_list[1:] # Remove initialisation edge - - sub_edges3 = [edge for edge in all_edges if "B" not in edge] - n_subedges3 = len(sub_edges3) - self.assertEqual(len(edge_list), n_subedges3) - self.assertTrue( - np.all( - [edge in sub_edges3 or edge[::-1] in sub_edges3 for edge in edge_list] - ) + # The middle spot neither significant nor with the receptor: the ring. + no_a = self._bools("A") + forward, _ = get_edges(self.adata, everywhere, no_a, no_a) + self.assertSameEdges(forward, [e for e in self.all_edges if "A" not in e]) + + # Forward and reverse differ when the ligand and receptor spots do. + only_a, not_a = ~no_a, no_a + forward, reverse = get_edges(self.adata, only_a, not_a, everywhere) + self.assertEqual(forward, [("A", n) for n in self.neigh_dict["A"]]) + self.assertSameEdges(reverse, [e for e in self.all_edges if "A" in e]) + self.assertEqual(reverse[0], ("B", "A")) + + # Nothing significant: no edges, and lists all the same. + self.assertEqual( + get_edges(self.adata, everywhere, everywhere, ~everywhere), [[], []] ) def test_get_interaction_matrix(self): @@ -230,7 +134,8 @@ def test_get_interaction_matrix(self): # Get interaction matrix int_matrix = het.get_interaction_matrix( cell_data, - self.neighbourhood_indices, + self.indptr, + self.indices, CELL_TYPE_LABELS, sig_bool, ligand_bool, @@ -254,52 +159,3 @@ def create_cci(cell_annotations: list[str], unique_cell_type_labels): ct_index = np.where(unique_cell_type_labels == annot)[0][0] cell_data[i, ct_index] = 1 return cell_data - - # TODO next things to test: - # 1. Getting the LR scores. - - -class TestRunCCILabels(unittest.TestCase): - """run_cci with cell-type labels stored as a categorical and as plain strings. - - Under pandas 3, `Series.values.astype(str)` returns an Arrow-backed string array, - which numba rejects, so run_cci failed for both dtypes. - """ - - adata: ClassVar[AnnData] - lr_results: ClassVar[AnnData] - - @classmethod - def setUpClass(cls): - adata, lrs = read_test_data_with_lr_pairs() - cls.lr_results = st.tl.cci.run( - adata, lrs, min_spots=5, distance=None, n_pairs=200, n_cpus=1, verbose=False - ) - # Three bands of cell types across the tissue. - bands = pd.qcut(adata.obs["imagecol"], 3, labels=CELL_TYPE_LABELS) - adata.obs["cell_type"] = bands - adata.obs["cell_type_str"] = bands.to_numpy().astype(str) - cls.adata = adata - - def assert_run_cci_runs(self, use_label: str): - adata = self.adata.copy() - lr_results = self.lr_results.copy() - st.tl.cci.run_cci( - adata, - use_label, - lr_results=lr_results, - min_spots=3, - n_perms=20, - n_cpus=1, - verbose=False, - ) - - matrix = lr_results.uns[f"lr_cci_{use_label}"] - self.assertEqual(list(matrix.index), CELL_TYPE_LABELS.tolist()) - self.assertEqual(list(matrix.columns), CELL_TYPE_LABELS.tolist()) - - def test_categorical_labels(self): - self.assert_run_cci_runs("cell_type") - - def test_string_labels(self): - self.assert_run_cci_runs("cell_type_str") diff --git a/tests/tl/test_cci_labels.py b/tests/tl/test_cci_labels.py new file mode 100644 index 00000000..f646cfc8 --- /dev/null +++ b/tests/tl/test_cci_labels.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python + +"""Runs st.tl.cci.run_cci on the test data with cell-type labels stored as a +categorical and as plain strings. +""" + +import unittest +from typing import ClassVar + +import numpy as np +import pandas as pd +from anndata import AnnData + +import stlearn as st +from tests.utils import read_test_data + +CELL_TYPE_LABELS = np.array(["CT1", "CT2", "CT3"]) + + +def read_test_data_with_lr_pairs(n_pairs: int = 20): + """Test data that ligand-receptor analysis can run on. + + The 100 genes in the test data contain no pair from the ligand-receptor + database and the matrix is scaled, so `st.tl.cci.run` has nothing to test. This + renames the first `2 * n_pairs` genes to database pairs and clips negative values + to zero. It exercises the code path; the scores carry no biological meaning. + """ + adata = read_test_data() + lrs = st.tl.cci.lr_pairs("connectomeDB2020") + + genes: list[str] = [] + for ligand, receptor in zip( + lrs.ligands().tolist(), lrs.receptors().tolist(), strict=True + ): + if ligand != receptor and ligand not in genes and receptor not in genes: + genes += [ligand, receptor] + if len(genes) == 2 * n_pairs: + break + + names = list(adata.var_names) + names[: len(genes)] = genes + adata.var_names = names + adata.X = np.clip(adata.X, 0, None) + return adata, lrs + + +class TestCCILabels(unittest.TestCase): + """run_cci with cell-type labels stored as a categorical and as plain strings. + + Under pandas 3, `Series.values.astype(str)` returns an Arrow-backed string array, + which numba rejects, so run_cci failed for both dtypes. + """ + + adata: ClassVar[AnnData] + lr_results: ClassVar[AnnData] + + @classmethod + def setUpClass(cls): + adata, lrs = read_test_data_with_lr_pairs() + cls.lr_results = st.tl.cci.run( + adata, lrs, min_spots=5, distance=None, n_pairs=200, n_cpus=1, verbose=False + ) + # Three bands of cell types across the tissue. + bands = pd.qcut(adata.obs["imagecol"], 3, labels=CELL_TYPE_LABELS) + adata.obs["cell_type"] = bands + adata.obs["cell_type_str"] = bands.to_numpy().astype(str) + cls.adata = adata + + def assert_run_cci_runs(self, use_label: str): + adata = self.adata.copy() + lr_results = self.lr_results.copy() + st.tl.cci.run_cci( + adata, + use_label, + lr_results=lr_results, + min_spots=3, + n_perms=20, + n_cpus=1, + verbose=False, + ) + + matrix = lr_results.uns[f"lr_cci_{use_label}"] + self.assertEqual(list(matrix.index), CELL_TYPE_LABELS.tolist()) + self.assertEqual(list(matrix.columns), CELL_TYPE_LABELS.tolist()) + + def test_categorical_labels(self): + self.assert_run_cci_runs("cell_type") + + def test_string_labels(self): + self.assert_run_cci_runs("cell_type_str") diff --git a/tests/tl/test_het.py b/tests/tl/test_het.py index 6e3301ad..feaf8558 100644 --- a/tests/tl/test_het.py +++ b/tests/tl/test_het.py @@ -12,6 +12,7 @@ import stlearn as st from stlearn.tl.cci import lr_result from stlearn.tl.cci.het import count, get_data_for_counting +from stlearn.tl.cci.neighbours import calc_neighbours from tests.utils import CCI_LRS, cci_table LRS = st.tl.cci.LRPairs.from_names( @@ -196,6 +197,31 @@ def test_spot_mixtures_counts_cell_types_never_dominant(self): self.assertIn("CT4", counts[True].index) self.assertGreater(counts[True].loc["CT4"].sum(), 0) + def test_counts_over_given_neighbours(self): + """Given the neighbour graph st.tl.cci.run has already built, count uses + it instead of searching again, the spot itself still in its neighbourhood, + and leaves the graph as it was.""" + searched, given, narrow = _make_adata(), _make_adata(), _make_adata() + count(searched, use_label="cell_type", verbose=False, distance=1.5) + count(narrow, use_label="cell_type", verbose=False, distance=1.0) + # 8 neighbours against 4: the two distances must count differently. + self.assertFalse( + np.array_equal(searched.obsm["cci_het"], narrow.obsm["cci_het"]) + ) + + graph = calc_neighbours(given, 1.5, verbose=False) + nnz = graph.nnz + # The graph decides the neighbourhood, whatever distance says. + count( + given, use_label="cell_type", verbose=False, distance=1.0, neighbours=graph + ) + np.testing.assert_array_equal(given.obsm["cci_het"], searched.obsm["cci_het"]) + # The spot was added to its own neighbourhood in a copy: run stores this + # graph as obsp["spot_neighbours"], where a diagonal would make every spot + # its own interaction partner in run_cci. + self.assertEqual(graph.nnz, nnz) + self.assertEqual(graph.diagonal().sum(), 0) + if __name__ == "__main__": unittest.main() diff --git a/tests/tl/test_lr_table.py b/tests/tl/test_lr_table.py index dc0ea841..8c73ef52 100644 --- a/tests/tl/test_lr_table.py +++ b/tests/tl/test_lr_table.py @@ -109,10 +109,8 @@ def test_layers_are_sparse_with_the_tested_spots(self): np.testing.assert_array_equal(lr_result(lr, key, j), dense[:, j]) def test_expression_table_keeps_only_what_is_per_spot(self): - self.assertEqual( - sorted(self.adata.obsm), - ["cci_het", "spatial", "spot_neigh_bcs", "spot_neighbours"], - ) + self.assertEqual(sorted(self.adata.obsm), ["cci_het", "spatial"]) + self.assertEqual(list(self.adata.obsp), ["spot_neighbours"]) self.assertFalse(set(self.adata.uns) & set(OLD_UNS)) def test_without_lr_results(self): diff --git a/tests/tl/test_neighbourhoods.py b/tests/tl/test_neighbourhoods.py new file mode 100644 index 00000000..a98da88c --- /dev/null +++ b/tests/tl/test_neighbourhoods.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python + +import unittest + +import numpy as np + +import stlearn as st +from stlearn.tl.cci.neighbours import calc_neighbours, stored_neighbours +from tests.utils import CCI_LRS, cci_table + + +class TestStoredNeighbours(unittest.TestCase): + """stored_neighbours reads the graph st.tl.cci.run stores in obsp.""" + + @classmethod + def setUpClass(cls): + cls.distance = 1.5 + adata = cci_table() + st.tl.cci.run( + adata, + CCI_LRS, + min_spots=2, + distance=cls.distance, + n_pairs=100, + verbose=False, + ) + cls.adata = adata + graph = calc_neighbours(adata, cls.distance, verbose=False) + cls.found = np.split(graph.indices, graph.indptr[1:-1]) + + def test_matches_calc_neighbours(self): + """The stored graph is the one calc_neighbours finds, as CSR, each + spot's neighbours in its order.""" + graph = stored_neighbours(self.adata) + self.assertEqual(graph.format, "csr") + self.assertEqual(graph.shape, (self.adata.n_obs, self.adata.n_obs)) + for i, found in enumerate(self.found): + np.testing.assert_array_equal( + graph.indices[graph.indptr[i] : graph.indptr[i + 1]], found + ) + + def test_subset_keeps_neighbours_still_in_it(self): + """Subsetting slices obsp on both axes, so a spot keeps only the + neighbours still in the subset, at their positions in it.""" + keep = np.arange(self.adata.n_obs) % 3 != 0 + subset = self.adata[keep].copy() + graph = stored_neighbours(subset) + names = self.adata.obs_names.to_numpy() + self.assertEqual(graph.shape, (subset.n_obs, subset.n_obs)) + for row, i in enumerate(np.flatnonzero(keep)): + expected = names[[n for n in self.found[i] if keep[n]]] + neighbours = graph.indices[graph.indptr[row] : graph.indptr[row + 1]] + self.assertEqual(sorted(subset.obs_names[neighbours]), sorted(expected)) + + def test_without_neighbours(self): + """A table without obsp["spot_neighbours"] says to run st.tl.cci.run.""" + adata = self.adata.copy() + del adata.obsp["spot_neighbours"] + with self.assertRaisesRegex(ValueError, r"st\.tl\.cci\.run"): + stored_neighbours(adata) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/tl/test_scoring.py b/tests/tl/test_scoring.py index a570a738..1639d102 100644 --- a/tests/tl/test_scoring.py +++ b/tests/tl/test_scoring.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python + """Tests for stlearn.tl.cci.scoring: LR scores against a plain-numpy reference.""" import unittest @@ -7,22 +9,20 @@ import pandas as pd import scipy.sparse from anndata import AnnData -from numba.typed import List from stlearn.tl.cci import scoring from stlearn.tl.cci.neighbours import calc_neighbours from stlearn.tl.cci.scoring import ( LRPairs, LRScorer, - flatten_neighbours, get_lrs_scores, - get_scores, neighbourhood_means, + score_spots, ) def _reference_scores(expr, pairs, neighbour_lists, het_vals, min_expr, spot_indices): - """Plain-numpy oracle for get_scores (slow but obvious). + """Plain-numpy oracle for score_spots (slow but obvious). For each pair: the (ligand, receptor) expression at the spot, against the (receptor, ligand) expression averaged over the spot's neighbours (0 with no @@ -44,6 +44,18 @@ def _reference_scores(expr, pairs, neighbour_lists, het_vals, min_expr, spot_ind return scores +def _graph(neighbour_lists): + """Each spot's neighbours as the spots-by-spots CSR matrix calc_neighbours + makes: row i's indices are neighbour_lists[i], in order.""" + n_spots = len(neighbour_lists) + lengths = [len(neighbours) for neighbours in neighbour_lists] + indptr = np.concatenate([[0], np.cumsum(lengths)]).astype(np.int64) + indices = np.concatenate([np.asarray(n, np.int32) for n in neighbour_lists]) + return scipy.sparse.csr_matrix( + (np.ones(len(indices), np.int8), indices, indptr), shape=(n_spots, n_spots) + ) + + def _random_neighbours(rng, n_spots, max_neighbours=4): """Random neighbour sets of up to max_neighbours - 1 other spots, some empty.""" neighbour_lists = [] @@ -56,33 +68,27 @@ def _random_neighbours(rng, n_spots, max_neighbours=4): return neighbour_lists -class TestGetScores(unittest.TestCase): +class TestScoreSpots(unittest.TestCase): def test_worked_example(self): """Spot 0: ligand 1 with receptor mean (3 + 7) / 2 = 5 gives 1 + 5; its receptor is 0, so (6 + 0) / 2 = 3. Spot 1: ligand 5 with receptor mean 0 gives 0; receptor 3 with ligand mean 1 gives 3 + 1; (0 + 4) / 2 = 2. Spot 2 has no neighbours, so 0. Then each is multiplied by its het value.""" expr = np.array([[1.0, 0.0], [5.0, 3.0], [0.0, 7.0]]) # ligand, receptor - neighbours = List( - [ - np.array([1, 2], dtype=np.int32), - np.array([0], dtype=np.int32), - np.array([], dtype=np.int32), - ], - ) + neighbours = _graph([[1, 2], [0], []]) pairs = np.array([[0, 1]]) spots = np.arange(3, dtype=np.int32) - scores = get_scores(expr, pairs, neighbours, np.ones(3), 0.0, spots) + scores = score_spots(expr, pairs, neighbours, np.ones(3), 0.0, spots) self.assertEqual(scores.dtype, np.float64) np.testing.assert_array_equal(scores[:, 0], [3.0, 2.0, 0.0]) het_vals = np.array([1, 2, 3]) - scores = get_scores(expr, pairs, neighbours, het_vals, 0.0, spots) + scores = score_spots(expr, pairs, neighbours, het_vals, 0.0, spots) np.testing.assert_array_equal(scores[:, 0], [3.0, 4.0, 0.0]) def test_matches_reference_on_random_inputs(self): - """get_scores agrees with the numpy oracle across shapes, thresholds, het + """score_spots agrees with the numpy oracle across shapes, thresholds, het values and subsets of spots. Integer-valued expression puts values exactly on min_expr, which is not above it.""" rng = np.random.default_rng(0) @@ -108,10 +114,10 @@ def test_matches_reference_on_random_inputs(self): some_spots = np.sort(rng.choice(n_spots, size=n_spots // 2, replace=False)) for spots in (all_spots, some_spots): - got = get_scores( + got = score_spots( np.asfortranarray(expr), pairs, - List(neighbour_lists), + _graph(neighbour_lists), het_vals, min_expr, spots, @@ -142,7 +148,8 @@ def test_means_once_per_gene(self): neighbour_lists = _random_neighbours(rng, n_spots) spots = np.sort(rng.choice(n_spots, size=12, replace=False)).astype(np.int32) - starts, flat = flatten_neighbours(List(neighbour_lists)) + graph = _graph(neighbour_lists) + starts, flat = graph.indptr.astype(np.int64), graph.indices.astype(np.int64) means = neighbourhood_means(expr, starts, flat, spots) self.assertEqual(means.shape, (n_genes, len(spots))) @@ -160,16 +167,16 @@ def test_spots_scored_in_chunks(self): n_spots, n_genes = 23, 4 expr = np.asfortranarray(rng.random((n_spots, n_genes))) pairs = np.array([[0, 1], [2, 3], [1, 2]]) - neighbours = List(_random_neighbours(rng, n_spots)) + neighbours = _graph(_random_neighbours(rng, n_spots)) het_vals = np.ones(n_spots) spots = np.arange(n_spots, dtype=np.int32) - whole = get_scores(expr, pairs, neighbours, het_vals, 0.0, spots) + whole = score_spots(expr, pairs, neighbours, het_vals, 0.0, spots) chunk_bytes = scoring.MEANS_CHUNK_BYTES try: # 8 bytes per mean, n_genes genes: 5 spots of means per chunk. scoring.MEANS_CHUNK_BYTES = 8 * n_genes * 5 - chunked = get_scores(expr, pairs, neighbours, het_vals, 0.0, spots) + chunked = score_spots(expr, pairs, neighbours, het_vals, 0.0, spots) finally: scoring.MEANS_CHUNK_BYTES = chunk_bytes np.testing.assert_array_equal(chunked, whole) @@ -216,7 +223,7 @@ def test_matches_reference_through_adata(self): expected = _reference_scores( dense.to_numpy(dtype=np.float64), pairs, - list(neighbours), + np.split(neighbours.indices, neighbours.indptr[1:-1]), het_vals, 0.5, np.arange(len(adata)), @@ -262,17 +269,18 @@ def test_homophilic_pair_scores_as_one_direction(self): self.assertEqual(list(scored.names()), ["Lig0_Lig0", "Lig0_Rec0"]) dense = adata.to_df().to_numpy(dtype=np.float64) + rows = np.split(neighbours.indices, neighbours.indptr[1:-1]) expected = _reference_scores( dense, [[0, 0]], - list(neighbours), + rows, het_vals, 0.5, np.arange(len(adata)), ) np.testing.assert_allclose(scores[:, [0]], expected, rtol=1e-6) at_spot = dense[:, 0] - own_mean = np.array([dense[n, 0].mean() if len(n) else 0.0 for n in neighbours]) + own_mean = np.array([dense[n, 0].mean() if len(n) else 0.0 for n in rows]) one_direction = at_spot * (own_mean > 0.5) + (at_spot > 0.5) * own_mean np.testing.assert_allclose(scores[:, 0], one_direction, rtol=1e-6) @@ -488,7 +496,7 @@ def test_only_a_chunks_rows_are_made_dense(self): neighbours = calc_neighbours(adata, distance=1.5, verbose=False) scorer = LRScorer(adata, neighbours) spots = np.array([0, 1, 35], dtype=np.int32) - needed = set(spots) | {n for s in spots for n in neighbours[s]} + needed = set(spots) | {n for s in spots for n in neighbours[s].indices} taken = [] original = scoring._dense_rows diff --git a/uv.lock b/uv.lock index 10b92c5c..6110250c 100644 --- a/uv.lock +++ b/uv.lock @@ -5610,7 +5610,7 @@ wheels = [ [[package]] name = "stlearn" -version = "1.5.0" +version = "1.5.1" source = { editable = "." } dependencies = [ { name = "anndata" },