persistence plot

Author: chengmingbo

README.md

@@ -75,6 +75,7 @@ We have tested python version `3.9.0`, `3.10.8`, `3.10.14`, `3.11.0`, `3.11.5`,
   - scanpy (>=1.9.3)
   - adjustText (>=0.8)
   - pygraphviz (>=1.11)
+  - gudhi (>=3.10.1)
   - magic-impute (>=3.0.0)
   - anndata (>=0.9.2)
 

phlower/plot/persistence.py

@@ -11,12 +11,13 @@
 import gudhi as gd
 
 def persisitence_barcodes(adata: AnnData,
-                          include_holes=True,
-                          persistence='persistence',
-                          barcodes_dim=[0],
-                          min_persistence=0.01,
-                          show_threshold=True,
-                          show_legend=True,
+                          include_holes:bool=True,
+                          persistence:str='persistence',
+                          barcodes_dim:List[int]=[0],
+                          min_persistence:float=0.01,
+                          show_threshold:bool=True,
+                          show_legend:bool=True,
+                          manual_threshold:Optional[List[float]] = None,
                           ax=None,
                           **args):
     """
@@ -34,6 +35,7 @@ def persisitence_barcodes(adata: AnnData,
     if not set(barcodes_dim) <= set([0, 1, 2, 3]):
         raise ValueError(f"barcodes_dim: {barcodes_dim} should be in [0, 1, 2, 3]")
 
+
     if include_holes:
         simplex_tree = adata.uns[persistence]['simplextree_tri_ae']
     else:
@@ -45,21 +47,26 @@ def persisitence_barcodes(adata: AnnData,
     limit = max(adata.uns[persistence]['filter_num']+10, ax.get_xlim()[1])
     ax.set_xlim(0, limit)
     if show_threshold:
-        ax.axvline(x=adata.uns[persistence]['filter_num'], ls='--', color='blue')
+        if not manual_threshold:
+            ax.axvline(x=adata.uns[persistence]['filter_num'], ls='--', color='blue')
+        else:
+            for th in manual_threshold:
+                ax.axvline(x=th, ls='--', color='blue')
     if not show_legend:
         ax.get_legend().remove()
 
     return ax
 
 
 def persisitence_birth_death(adata: AnnData,
-                          include_holes=True,
-                          persistence='persistence',
-                          min_persistence=0.01,
-                          show_threshold=True,
-                          show_legend=True,
-                          ax=None,
-                          **args):
+                             include_holes:bool=True,
+                             persistence:str='persistence',
+                             min_persistence:float=0.01,
+                             show_threshold:bool=True,
+                             show_legend:bool=True,
+                             manual_threshold:Optional[List[float]] = None,
+                             ax=None,
+                             **args):
     """
     Plot the persistence birth death diagram
 
@@ -81,7 +88,11 @@ def persisitence_birth_death(adata: AnnData,
     pers = simplex_tree.persistence(min_persistence=min_persistence)
     ax = gd.plot_persistence_diagram(pers, axes=ax, **args)
     if show_threshold:
-        ax.axhline(y=adata.uns[persistence]['filter_num'], ls='--', color='blue')
+        if not manual_threshold:
+            ax.axhline(y=adata.uns[persistence]['filter_num'], ls='--', color='blue')
+        else:
+            for th in manual_threshold:
+                ax.axhline(y=th, ls='--', color='blue')
     if not show_legend:
         ax.get_legend().remove()
     return ax

phlower/tools/triangulation.py

@@ -1069,6 +1069,6 @@ def barcode0_max(barcode0):
 
 
 def persistence_tree_edges(g, filtration_list, edge_list, filter_num):
-    index = np.where(np.array(filtration_list) < filter_num)[0]
+    index = np.where(np.array(filtration_list) <= filter_num)[0]
     keep_edges = np.array(edge_list)[index]
     return keep_edges