v0.3.1: Mini update

Author: xiachenrui

decipher/cfg.yaml

@@ -73,7 +73,7 @@ sp_model:
     warmup_steps: 1000
   loader:
     batch_size: 256
-    num_workers: 4
+    num_workers: 8
     pin_memory: true
     persistent_workers: true
   num_neighbors: [-1]

decipher/cls.py

@@ -58,6 +58,7 @@ def register_data(
         split_by: str = None,
         preprocess: bool = True,
         edge_index: np.ndarray = None,
+        save_data: bool = True,
     ) -> None:
         r"""
         Register spatial omics data
@@ -73,10 +74,11 @@ def register_data(
         edge_index
             use self-defined spatial neighbor edges (PyG format), only for advanced users
         """
-        # preprocess single cell data
-        adata, self.batch_idx = omics_data_process(adata, split_by, preprocess)
-        self.x = adata.X.astype(np.float32)
-        np.save(self.work_dir / "x.npy", self.x)
+        # process spatial omics data
+        self.x, coords, self.batch_idx = omics_data_process(adata, split_by, preprocess)
+        del adata
+        if save_data:
+            np.save(self.work_dir / "x.npy", self.x)
 
         # mnn correction
         if self.batch_idx is not None:
@@ -87,9 +89,10 @@ def register_data(
 
         # build spatial graph
         if edge_index is None:
-            self.edge_index = build_graph(adata.obsm["spatial"], self.batch_idx, **CFG.graph)
+            self.edge_index = build_graph(coords, self.batch_idx, **CFG.graph)
         else:
             logger.info("Use self-defined edge index.")
+            assert edge_index.max() < self.x.shape[0], "Edge index out of range."
             self.edge_index = edge_index
         np.save(self.work_dir / "edge_index.npy", self.edge_index.numpy())
 

decipher/data/process.py

@@ -13,7 +13,7 @@
 
 
 def omics_data_process(
-    adata: list[AnnData] | AnnData,
+    adata_raw: list[AnnData] | AnnData,
     split_by: str = None,
     preprocess: bool = True,
 ) -> tuple[AnnData, np.ndarray | None]:
@@ -34,20 +34,22 @@ def omics_data_process(
     Preprocessed AnnData
     """
     batch_idx = None
-    if isinstance(adata, AnnData):
-        if split_by is not None and split_by in adata.obs.columns:
-            batch_idx = LabelEncoder().fit_transform(adata.obs[split_by])
+    if isinstance(adata_raw, AnnData):
+        if split_by is not None and split_by in adata_raw.obs.columns:
+            batch_idx = LabelEncoder().fit_transform(adata_raw.obs[split_by])
+        adata = adata_raw.copy()
     else:
-        adata = adata[0].concatenate(adata[1:], batch_key="batch", uns_merge="same")
+        adata = adata_raw[0].concatenate(adata_raw[1:], batch_key="batch", uns_merge="same")
         batch_idx = LabelEncoder().fit_transform(adata.obs["batch"])
+    del adata_raw
 
     if batch_idx is not None and not CFG.pp.ignore_batch:
         logger.info(f"Detected {np.unique(batch_idx).shape[0]} batches.")
         adata.obs["_batch"] = batch_idx
 
     if not preprocess:
         logger.warning("Skip preprocessing steps, only for advanced users.")
-        return adata.copy(), batch_idx
+        return adata.X.astype(np.float32), adata.obsm["spatial"], batch_idx
 
     logger.info(f"Preprocessing {adata.n_obs} cells.")
     start_time = time.time()
@@ -79,4 +81,4 @@ def omics_data_process(
     sc.pp.log1p(adata)
     sc.pp.scale(adata, max_value=10)
     logger.success(f"Preprocessing finished in {time.time() - start_time:.2f} seconds.")
-    return adata.copy(), batch_idx
+    return adata.X.astype(np.float32), adata.obsm["spatial"], batch_idx

decipher/explain/gene/lr.py

@@ -14,10 +14,11 @@
 
 def get_lr_expr(
     adata: sc.AnnData,
-    lr: pd.DataFrame,
+    lr_df: pd.DataFrame,
     radius: float,
     aggr: str = "mean",
     binary: bool = False,
+    threshold: float = 0.01,
 ) -> tuple[torch.Tensor, pd.DataFrame]:
     r"""
     Get ligand-receptor activity from spatial data.
@@ -26,23 +27,24 @@ def get_lr_expr(
 
     Args:
         adata: AnnData object
-        lr: Ligand receptor pairs, must contain columns `ligand.symbol` and `receptor.symbol`
+        lr_df: Ligand receptor pairs, must contain columns `ligand.symbol` and `receptor.symbol`
         radius: float radius for valid ligand
         aggr: str aggregation method of ligand expression, default is `mean`
         binary: bool if convert the expression to binary, default is `False`
+        threshold: float threshold to filter low-expressed LR pairs, default is `0.01`
     Returns:
         LR activity: Ligand-receptor pair activity
-        lr: filtered ligand receptor pairs
+        lr_df: filtered ligand receptor pairs
     Note:
         Not support batched data yet.
     """
     # check
-    assert "ligand.symbol" in lr.columns
-    assert "receptor.symbol" in lr.columns
+    assert "ligand.symbol" in lr_df.columns
+    assert "receptor.symbol" in lr_df.columns
 
     # filter non-LR genes (reduce memory usage)
     lr_genes = []
-    for _, row in lr.iterrows():
+    for _, row in lr_df.iterrows():
         ligand = row["ligand.symbol"]
         if "," in ligand:
             lr_genes += ligand.split(",")
@@ -57,20 +59,21 @@ def get_lr_expr(
     lr_genes = list(set(lr_genes))
     lr_idx = [x in lr_genes for x in adata.var.index]
     adata = adata[:, lr_idx]
-    sc.pp.filter_genes(adata, min_cells=30)
+
     logger.info(f"Find {adata.n_vars} LR-related genes.")
 
     # merge neighbor expr
     edge_index = build_graph(adata.obsm["spatial"], radius=radius, mode="radius")
     aggr_net = SimpleConv(aggr=aggr, combine_root=None)
     expr = torch.from_numpy(adata.X.toarray()) if isspmatrix(adata.X) else torch.from_numpy(adata.X)
+    # expr = expr.to(torch.float32)
     expr_neighbors = aggr_net(expr, edge_index)
 
     # get the ligand and receptor gene index in adata.var
     lr_activity_list = []
     gene_set = set(adata.var.index)
-    lr_filter = np.zeros(len(lr), dtype=bool)
-    for i, lr_pair in lr.iterrows():
+    lr_filter = np.zeros(len(lr_df), dtype=bool)
+    for i, lr_pair in lr_df.iterrows():
         ligand = lr_pair["ligand.symbol"]
         receptor = lr_pair["receptor.symbol"]
         ligand = [ligand] if "," not in ligand else ligand.split(",")
@@ -82,10 +85,16 @@ def get_lr_expr(
         receptor_idx = [adata.var.index.get_loc(x) for x in receptor if x in gene_set]
         ligand_expr = expr_neighbors[:, ligand_idx].mean(dim=1)  # from neighbor
         receptor_expr = expr[:, receptor_idx].mean(dim=1)  # from cell
+        lr_activity = ligand_expr * receptor_expr
+        if (lr_activity > 0).float().mean().item() < threshold:
+            continue
         lr_activity_list.append(ligand_expr * receptor_expr)
         lr_filter[i] = True
+    if len(lr_activity_list) == 0:
+        logger.warning(f"No ligand-receptor pairs (threshold: {threshold}).")
+        return None, None
     lr_activity = torch.stack(lr_activity_list, dim=1)
     if binary:
         lr_activity = (lr_activity > 0).float()
-    lr = lr[lr_filter]
-    return lr_activity, lr
+    lr_df = lr_df[lr_filter]
+    return lr_activity, lr_df

decipher/graphic/build.py

@@ -7,7 +7,7 @@
 from loguru import logger
 from torch import Tensor
 from torch_geometric.nn import knn_graph, radius_graph
-from torch_geometric.utils import to_undirected
+from torch_geometric.utils import sort_edge_index, to_undirected
 
 from ..utils import estimate_spot_distance
 from .knn import knn
@@ -108,13 +108,13 @@ def build_graph(
     )
     if mode == "radius":
         num_neighbors = edge_index[0].bincount().float()
-        logger.info(f"Mean number of neighbors: {num_neighbors.mean().item():.2f}")
-        logger.info(f"Max number of neighbors: {num_neighbors.max().item()}")
-        logger.info(f"Min number of neighbors: {num_neighbors.min().item()}")
-        logger.info(f"Median number of neighbors: {num_neighbors.median().item()}")
+        logger.debug(f"Mean number of neighbors: {num_neighbors.mean().item():.2f}")
+        logger.debug(f"Max number of neighbors: {num_neighbors.max().item()}")
+        logger.debug(f"Min number of neighbors: {num_neighbors.min().item()}")
+        logger.debug(f"Median number of neighbors: {num_neighbors.median().item()}")
         for percentile in [0.05, 0.25, 0.75, 0.95]:
             tile = num_neighbors.kthvalue(int(num_neighbors.numel() * percentile)).values.item()
-            logger.info(f"{percentile * 100}th percentile of number of neighbors: {tile}")
+            logger.debug(f"{percentile * 100}th percentile of number of neighbors: {tile}")
     return edge_index
 
 
@@ -135,4 +135,4 @@ def knn_to_edge_index(knn_result: np.ndarray) -> Tensor:
 
     edge_index = torch.stack([src_nodes, dst_nodes], dim=0)
 
-    return edge_index
+    return sort_edge_index(edge_index)

pyproject.toml

@@ -6,7 +6,7 @@ build-backend = "poetry.core.masonry.api"
 [tool.poetry]
 name = "cell-decipher"
 packages = [{include = "decipher"}]
-version = "0.3.0"
+version = "0.3.1"
 description = "DECIPHER for learning disentangled cellular embeddings in large-scale heterogeneous spatial omics data"
 readme = "README.md"
 license = "MIT"