Use more modern approach to loading included databases and fix mouse …

stlearn/tl/cci/analysis.py

@@ -3,6 +3,7 @@
 """
 
 import os
+from importlib.resources import files
 
 import numba
 import numpy as np
@@ -27,7 +28,7 @@
 
 # Functions related to Ligand-Receptor interactions
 def load_lrs(
-    names: str | list | None = None, species: str = "human"
+        names: str | list[str] | None = None, species: str = "human",
 ) -> npt.NDArray[np.str_]:
     """Loads inputted LR database, & concatenates into consistent database set of
     pairs without duplicates. If None loads 'connectomeDB2020_lit'.
@@ -50,38 +51,30 @@ def load_lrs(
     if isinstance(names, str):
         names = [names]
 
-    path = os.path.dirname(os.path.realpath(__file__))
-    dbs = [pd.read_csv(f"{path}/databases/{name}.txt", sep="\t") for name in names]
-    lrs_full = []
-    for db in dbs:
-        lrs = [f"{db.values[i, 0]}_{db.values[i, 1]}" for i in range(db.shape[0])]
-        lrs_full.extend(lrs)
-    lrs_full_arr = np.unique(np.array(lrs_full))
+    db_dir = files("stlearn.tl.cci") / "databases"
+    lrs: set[str] = set()
+    for name in names:
+        with (db_dir / f"{name}.txt").open("rb") as fh:
+            db = pd.read_csv(fh, sep="\t")
+        lrs.update(
+            f"{ligand}_{receptor}" for ligand, receptor in db.iloc[:, :2].values
+        )
     # If dealing with mouse, need to reformat #
     if species == "mouse":
-        genes1 = [lr_.split("_")[0] for lr_ in lrs_full]
-        genes2 = [lr_.split("_")[1] for lr_ in lrs_full]
-        lrs_full_arr = np.array(
-            [
-                genes1[i][0]
-                + genes1[i][1:].lower()
-                + "_"
-                + genes2[i][0]
-                + genes2[i][1:].lower()
-                for i in range(len(lrs_full))
-            ],
-        )
-
-    return lrs_full_arr
+        lrs = {
+            f"{ligand[0]}{ligand[1:].lower()}_{receptor[0]}{receptor[1:].lower()}"
+            for ligand, receptor in (lr_.split("_") for lr_ in lrs)
+        }
+    return np.array(sorted(lrs), dtype=np.str_)
 
 
 def grid(
-    adata: AnnData,
-    n_row: int = 10,
-    n_col: int = 10,
-    use_label: str | None = None,
-    n_cpus: int | None = None,
-    verbose: bool = True,
+        adata: AnnData,
+        n_row: int = 10,
+        n_col: int = 10,
+        use_label: str | None = None,
+        n_cpus: int | None = None,
+        verbose: bool = True,
 ) -> AnnData:
     """Creates a new anndata representing a gridded version of the data; can be
         used upstream of CCI pipeline. NOTE: intended use is for single cell
@@ -199,20 +192,20 @@ def grid(
 
 
 def run(
-    adata: AnnData,
-    lrs: npt.NDArray[np.str_],
-    min_spots: int = 10,
-    distance: float | None = None,
-    n_pairs: int = 1000,
-    n_cpus: int | None = None,
-    use_label: str | None = None,
-    adj_method: str = "fdr_bh",
-    pval_adj_cutoff: float = 0.05,
-    min_expr: float = 0.0,
-    save_bg: bool = False,
-    neg_binom: bool = False,
-    random_state: int = 0,
-    verbose: bool = True,
+        adata: AnnData,
+        lrs: npt.NDArray[np.str_],
+        min_spots: int = 10,
+        distance: float | None = None,
+        n_pairs: int = 1000,
+        n_cpus: int | None = None,
+        use_label: str | None = None,
+        adj_method: str = "fdr_bh",
+        pval_adj_cutoff: float = 0.05,
+        min_expr: float = 0.0,
+        save_bg: bool = False,
+        neg_binom: bool = False,
+        random_state: int = 0,
+        verbose: bool = True,
 ) -> None:
     """Performs stLearn LR analysis.
 
@@ -373,10 +366,10 @@ def run(
 
 
 def adj_pvals(
-    adata,
-    pval_adj_cutoff: float = 0.05,
-    correct_axis: str = "spot",
-    adj_method: str = "fdr_bh",
+        adata,
+        pval_adj_cutoff: float = 0.05,
+        correct_axis: str = "spot",
+        adj_method: str = "fdr_bh",
 ):
     """Performs p-value adjustment and determination of significant spots.
         Default settings of this function are already run in st.tl.cci.run.
@@ -455,16 +448,16 @@ def adj_pvals(
 
 
 def run_lr_go(
-    adata: AnnData,
-    r_path: str,
-    n_top: int = 100,
-    bg_genes: np.ndarray | None = None,
-    min_sig_spots: int = 1,
-    species: str = "human",
-    p_cutoff: float = 0.01,
-    q_cutoff: float = 0.5,
-    onts: str = "BP",
-    verbose: bool = True,
+        adata: AnnData,
+        r_path: str,
+        n_top: int = 100,
+        bg_genes: np.ndarray | None = None,
+        min_sig_spots: int = 1,
+        species: str = "human",
+        p_cutoff: float = 0.01,
+        q_cutoff: float = 0.5,
+        onts: str = "BP",
+        verbose: bool = True,
 ):
     """Runs a basic GO analysis on the genes in the top ranked LR pairs.
         Only supported for human and mouse species.
@@ -533,17 +526,17 @@ def run_lr_go(
 
 # Functions for calling Celltype-Celltype interactions
 def run_cci(
-    adata: AnnData,
-    use_label: str,
-    spot_mixtures: bool = False,
-    min_spots: int = 3,
-    sig_spots: bool = True,
-    cell_prop_cutoff: float = 0.2,
-    p_cutoff: float = 0.05,
-    n_perms: int = 100,
-    n_cpus: int | None = None,
-    random_state: int = 0,
-    verbose: bool = True,
+        adata: AnnData,
+        use_label: str,
+        spot_mixtures: bool = False,
+        min_spots: int = 3,
+        sig_spots: bool = True,
+        cell_prop_cutoff: float = 0.2,
+        p_cutoff: float = 0.05,
+        n_perms: int = 100,
+        n_cpus: int | None = None,
+        random_state: int = 0,
+        verbose: bool = True,
 ):
     """Calls significant celltype-celltype interactions based on cell-type data
     randomisation.
@@ -658,8 +651,8 @@ def run_cci(
         if not cols_present or not rows_present:
             if not cols_present:
                 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:\n"
+                              f"{[cell for cell in all_set if cell not in adata.uns[uns_key]]}\n"
                 )
             elif not rows_present:
                 msg = msg + "Rows do not correspond to adata.obs_names.\n"
@@ -706,11 +699,11 @@ def run_cci(
     lr_n_spot_cci_sig = np.zeros(lr_summary.shape[0])
     lr_n_cci_sig = np.zeros(lr_summary.shape[0])
     with tqdm(
-        total=len(best_lrs),
-        desc="Counting celltype-celltype interactions per LR and permuting "
-        + f"{n_perms} times.",
-        bar_format="{l_bar}{bar} [ time left: {remaining} ]",
-        disable=not verbose,
+            total=len(best_lrs),
+            desc="Counting celltype-celltype interactions per LR and permuting "
+                 + f"{n_perms} times.",
+            bar_format="{l_bar}{bar} [ time left: {remaining} ]",
+            disable=not verbose,
     ) as pbar:
         for i, best_lr in enumerate(best_lrs):
             ligand, receptor = best_lr.split("_")

tests/tl/test_cci.py

@@ -81,6 +81,10 @@ def test_load_lrs(self):
         self.assertTrue(np.all([gene[0].isupper() for gene in genes2]))
         self.assertTrue(np.all([gene[1:] == gene[1:].lower() for gene in genes2]))
 
+        # Should not have duplicates.
+        self.assertEqual(len(lrs), len(set(lrs)))
+        self.assertLessEqual(len(lrs), sizes[1])
+
     # Important, granular tests related to LR scoring
 
     # Important, granular tests related to CCI counting