Add make_factor_obs and bugfixes

Author: pedrofale

src/scdef/models/_scdef.py

@@ -1556,7 +1556,8 @@ def filter_factors(
         self,
         thres: Optional[float] = 0.0,
         iqr_mult: Optional[float] = 0.0,
-        min_cells: Optional[float] = 0.001,
+        min_cells_upper: Optional[float] = 0.001,
+        min_cells_lower: Optional[float] = 0.0,
         filter_up: Optional[bool] = True,
         normalized: Optional[bool] = False,
     ):
@@ -1568,17 +1569,20 @@ def filter_factors(
             min_cells: minimum number of cells that each factor must have attached to it for it to be kept. If between 0 and 1, fraction. Otherwise, absolute value
             filter_up: whether to remove factors in upper layers via inter-layer attachments
         """
-        if min_cells != 0:
-            if min_cells < 1.0:
-                min_cells = max(min_cells * self.adata.shape[0], 10)
+        if min_cells_upper != 0:
+            if min_cells_upper < 1.0:
+                min_cells_upper = max(min_cells_upper * self.adata.shape[0], 10)
+        if min_cells_lower != 0:
+            if min_cells_lower < 1.0:
+                min_cells_lower = max(min_cells_lower * self.adata.shape[0], 10)
 
         self.factor_lists = []
         for i, layer_name in enumerate(self.layer_names):
             if i == 0:
                 keep = self.get_effective_factors(
                     thres=thres,
                     iqr_mult=iqr_mult,
-                    min_cells=min_cells,
+                    min_cells=min_cells_lower,
                     normalized=normalized,
                 )
             else:
@@ -1590,7 +1594,7 @@ def filter_factors(
                     ]
                 )
                 keep = np.array(range(self.layer_sizes[i]))[
-                    np.where(counts >= min_cells)[0]
+                    np.where(counts >= min_cells_upper)[0]
                 ]
                 if filter_up:
                     mat = self.pmeans[f"{layer_name}W"][keep]

src/scdef/plotting/graph.py

@@ -939,22 +939,22 @@ def plot_biological_hierarchy(model, **kwargs):
     return g
 
 
-def plot_technical_hierarchy(model, **kwargs):
+def plot_technical_hierarchy(model, show_signatures=True, **kwargs):
     technical_signature = None
     technical_scores = None
-    if "show_signatures" in kwargs:
-        if kwargs["show_signatures"]:
-            technical_signature, technical_scores = get_technical_signature(
-                model,
-                return_scores=True,
-                top_genes=None,
-            )
+    if show_signatures:
+        technical_signature, technical_scores = get_technical_signature(
+            model,
+            return_scores=True,
+            top_genes=None,
+        )
     g = make_technical_hierarchy_graph(
         model,
         hierarchy=model.adata.uns["technical_hierarchy"],
         root_gene_rankings=[technical_signature],
         root_gene_scores=[technical_scores],
         root_name="tech_top",
+        show_signatures=show_signatures,
         **kwargs,
     )
     return g

src/scdef/tools/factor.py

@@ -1,5 +1,23 @@
 import numpy as np
 import pandas as pd
+from .hierarchy import get_hierarchy, compute_hierarchy_scores
+
+
+def make_factor_obs(model):
+    res = compute_hierarchy_scores(model)
+    model.adata.uns["factor_obs"] = res["per_factor"].set_index("child_factor")
+    model.adata.uns["factor_obs"]["ARD"] = np.array(
+        [np.nan] * len(model.adata.uns["factor_obs"])
+    )
+    model.adata.uns["factor_obs"]["BRD"] = np.array(
+        [np.nan] * len(model.adata.uns["factor_obs"])
+    )
+    model.adata.uns["factor_obs"].loc[model.factor_names[0], "ARD"] = np.asarray(
+        model.pmeans["factor_means"]
+    )[model.factor_lists[0]].ravel()
+    model.adata.uns["factor_obs"].loc[model.factor_names[0], "BRD"] = np.asarray(
+        model.pmeans["factor_concentrations"]
+    )[model.factor_lists[0]].ravel()
 
 
 def set_factor_signatures(model, signatures=None, top_genes=10):
@@ -10,16 +28,30 @@ def set_factor_signatures(model, signatures=None, top_genes=10):
 
 
 def set_technical_factors(model, factors=None):
-    """Set the technical factors of the model."""
+    """Set the technical factors of the model. They must be layer 0 factors."""
     # in model.adata.uns["factor_obs"], annotate as technical or not.
+    all_factor_names = [name for names in model.factor_names for name in names][
+        :-1
+    ]  # do not use root
     if "factor_obs" not in model.adata.uns:
         # Collect all factor names from all layers into a flat list
-        all_factor_names = [name for names in model.factor_names for name in names]
-        model.adata.uns["factor_obs"] = pd.DataFrame(index=all_factor_names)
+        make_factor_obs(model)
     model.adata.uns["factor_obs"]["technical"] = np.array(
         [False] * len(all_factor_names)
     )
-    model.adata.uns["factor_obs"]["technical"].loc[factors] = True
+    model.adata.uns["factor_obs"].loc[factors, "technical"] = True
+
+    # Get complete hierarchy
+    complete_hierarchy = get_hierarchy(model, simplified=False)
+    # Traverse hierarchy. If all the children of a factor are technical, set the factor as technical.
+    for factor, children in complete_hierarchy.items():
+        if all(
+            [
+                model.adata.uns["factor_obs"].loc[child, "technical"]
+                for child in children
+            ]
+        ):
+            model.adata.uns["factor_obs"].loc[factor, "technical"] = True
 
 
 def __build_consensus_signature(var_names, gene_scores_array, sizes_array):

src/scdef/tools/hierarchy.py

@@ -159,9 +159,6 @@ def _usage_weights(child_layer_idx, child_factor_indices):
             np.sum(scores_concat * weights_concat) / (np.sum(weights_concat) + eps)
         )
 
-    if "factor_obs" not in model.adata.uns:
-        model.adata.uns["factor_obs"] = per_factor.set_index("child_factor")
-
     return {
         "per_factor": per_factor,
         "per_transition": per_transition,
@@ -185,7 +182,9 @@ def make_technical_hierarchy(model):
     """Make the technical hierarchy of the model."""
     technical_factors = model.adata.uns["factor_obs"][
         model.adata.uns["factor_obs"]["technical"]
-    ].index.tolist()  # layer 0
+    ][
+        model.adata.uns["factor_obs"]["child_layer"] == "L0"
+    ].index.tolist()  # only layer 0 factors
     # technical hierarchy is a root with all technical factors as direct children.
     # connection weights are proportional to the usage of each factor
     technical_hierarchy = dict()

tests/test_integration.py

@@ -97,7 +97,7 @@ def test_scdef():
 
     model.fit(n_epoch=3)
 
-    model.filter_factors(thres=0.0, min_cells=0)  # make sure we keep factors
+    model.filter_factors(thres=0.0, min_cells_lower=0)  # make sure we keep factors
 
     model.logger.info(model.factor_lists)