Skip to content

Commit d67f17f

Browse files
committed
More data type fixes and improve speed of generating random pairs.
1 parent 718be6e commit d67f17f

4 files changed

Lines changed: 57 additions & 17 deletions

File tree

stlearn/tl/cci/analysis.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -259,7 +259,7 @@ def run(
259259
Whether print dialogue to user during run-time.
260260
Returns
261261
--------
262-
adata: None
262+
None
263263
Relevant information stored:
264264
adata.uns['lr_summary']
265265
Summary of significant spots detected per LR,
@@ -398,7 +398,7 @@ def adj_pvals(
398398
-------
399399
adata: AnnData
400400
Adjusts all of the LR results; warning, does not adjust
401-
celltype-celltype results from running ran st.tl.run_cci downstream.
401+
celltype-celltype results from running st.tl.run_cci downstream.
402402
"""
403403
if "lr_summary" not in adata.uns:
404404
raise Exception("Need to run st.tl.cci.run first.")

stlearn/tl/cci/base.py

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,13 +42,13 @@ def calc_distance(adata: AnnData, distance: float | None) -> float:
4242

4343
def get_lrs_scores(
4444
adata: AnnData,
45-
lrs: np.ndarray,
45+
lrs: npt.NDArray[np.str_],
4646
neighbours: np.ndarray,
4747
het_vals: np.ndarray,
4848
min_expr: float,
4949
filter_pairs: bool = True,
5050
spot_indices: npt.NDArray[np.int32] | None = None,
51-
):
51+
) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.str_]]:
5252
"""Gets the scores for the indicated set of LR pairs & the heterogeneity values.
5353
Parameters
5454
----------
@@ -69,7 +69,12 @@ def get_lrs_scores(
6969
Subset of spots to score, given as their integer row positions.
7070
Returns
7171
-------
72-
lrs: np.ndarray lr pairs from the database in format ['L1_R1', 'LN_RN']
72+
lr_scores: npt.NDArray[np.float64]
73+
Shape (n_scored_spots, n_pairs). LR score for each scored spot (rows in
74+
the order of spot_indices) and each ligand-receptor pair (columns).
75+
new_lrs: npt.NDArray[np.str_]
76+
Shape (n_pairs,). Ligand-receptor pair labels in 'L_R' format
77+
(e.g. 'L1_R1'), column-aligned with lr_scores.
7378
"""
7479
if spot_indices is None:
7580
spot_indices = np.arange(len(adata), dtype=np.int32)

stlearn/tl/cci/perm_utils.py

Lines changed: 40 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,9 @@
44
from numba.typed import List
55
from scipy.spatial.distance import canberra
66
from sklearn.preprocessing import MinMaxScaler
7+
from numba.core import types
8+
import numpy.typing as npt
9+
from numba.typed import Dict, List
710

811
from .base import get_lrs_scores
912

@@ -264,23 +267,54 @@ def get_similar_genes_fast(
264267
return similar_genes
265268

266269

267-
@njit
268-
def gen_rand_pairs(genes1: np.ndarray, genes2: np.ndarray, n_pairs: int, seed: int):
269-
"""Generates random pairs of genes."""
270+
def gen_rand_pairs(
271+
genes1: npt.NDArray[np.str_],
272+
genes2: npt.NDArray[np.str_],
273+
n_pairs: int,
274+
seed: int,
275+
) -> npt.NDArray[np.str_]:
276+
"""Generate unique random gene pairs for building background LR scores.
277+
278+
Each pair is formed by drawing one gene from genes1 and one from genes2,
279+
formatted as 'gene1_gene2'. Self-pairs (the same gene on both sides) and
280+
duplicate pairs are rejected, so every returned pair is unique and its two
281+
genes are distinct.
282+
283+
Parameters
284+
----------
285+
genes1: npt.NDArray[np.str_]
286+
Candidate genes for the first (ligand) position of each pair.
287+
genes2: npt.NDArray[np.str_]
288+
Candidate genes for the second (receptor) position of each pair.
289+
n_pairs: int
290+
Number of unique pairs to generate.
291+
seed: int
292+
Seed for the random generator, for reproducible pair selection.
293+
"""
294+
n_possible = len(genes1) * len(genes2) - np.intersect1d(genes1, genes2).size
295+
if n_pairs > n_possible:
296+
raise ValueError(
297+
f"Requested {n_pairs} unique pairs but only {n_possible} are "
298+
f"possible from {len(genes1)}×{len(genes2)} genes."
299+
)
300+
return np.array(list(_gen_rand_pairs(genes1, genes2, n_pairs, seed)))
301+
270302

303+
@njit
304+
def _gen_rand_pairs(genes1, genes2, n_pairs, seed):
271305
np.random.seed(seed) # noqa: NPY002 (numba requires legacy API)
272306
rand_pairs = List()
307+
seen = Dict.empty(types.unicode_type, types.boolean) # O(1) membership
273308
for _j in range(0, n_pairs):
274309
l_rand = np.random.choice(genes1, 1)[0] # noqa: NPY002
275310
r_rand = np.random.choice(genes2, 1)[0] # noqa: NPY002
276311
rand_pair = "_".join([l_rand, r_rand])
277-
while rand_pair in rand_pairs or l_rand == r_rand:
312+
while rand_pair in seen or l_rand == r_rand: # was: in rand_pairs
278313
l_rand = np.random.choice(genes1, 1)[0] # noqa: NPY002
279314
r_rand = np.random.choice(genes2, 1)[0] # noqa: NPY002
280315
rand_pair = "_".join([l_rand, r_rand])
281-
282316
rand_pairs.append(rand_pair)
283-
317+
seen[rand_pair] = True
284318
return rand_pairs
285319

286320

stlearn/tl/cci/permutation.py

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@
1313

1414
def perform_spot_testing(
1515
adata: AnnData,
16-
lr_scores: np.ndarray,
16+
lr_scores: npt.NDArray[np.float64],
1717
lrs: npt.NDArray[np.str_],
1818
n_pairs: int,
1919
neighbours: List,
@@ -145,14 +145,15 @@ def perform_spot_testing(
145145
else: # Fitting NB per LR
146146
lr_j_scores = lr_scores[spot_indices, lr_j]
147147
bg_ = background.ravel()
148-
bg_wScore = np.array(list(lr_j_scores) + list(bg_))
148+
bg_w_score = np.array(list(lr_j_scores) + list(bg_))
149149

150150
# 1) rounding discretisation
151151
# First multiple to get minimum value to be one before rounding #
152-
bg_1 = bg_wScore * (1 / min(bg_wScore[bg_wScore != 0]))
153-
bg_1 = np.round(bg_1)
154-
lr_j_scores_1 = bg_1[0 : len(lr_j_scores)]
155-
bg_1 = bg_1[len(lr_j_scores) : len(bg_1)]
152+
scale = 1.0 / min(bg_w_score[bg_w_score != 0])
153+
scaled = np.round(bg_w_score * scale)
154+
n_obs = len(lr_j_scores)
155+
lr_j_scores_1 = scaled[:n_obs]
156+
bg_1 = scaled[n_obs:]
156157

157158
# Getting the pvalue from negative binomial approach
158159
round_pvals, _, _, _ = get_stats(

0 commit comments

Comments
 (0)