Merge pull request #21 from Image-Analysis-Hub/feature/more-mandatory-cycle-feats

Feature/more mandatory cycle feats

- Add cycle_duration as mandatory cycle feature
- Make DivisionTime and DivisionRate computation more robust
- Rename cell_cycle_completeness feature into cycle_completeness

Author: lxenard

notebooks/Managing features.ipynb

@@ -308,7 +308,7 @@ def cells_Feature(provenance: str = "Pycellin") -> Feature:
 def cycle_length_Feature(provenance: str = "Pycellin") -> Feature:
     feat = Feature(
         name="cycle_length",
-        description="Number of cells in the cell cycle",
+        description="Number of cells in the cell cycle, minding gaps",
         provenance=provenance,
         feat_type="node",
         lin_type="CycleLineage",
@@ -317,7 +317,17 @@ def cycle_length_Feature(provenance: str = "Pycellin") -> Feature:
     return feat
 
 
-# TODO: define_cycle_duration_Feature
+def cycle_duration_Feature(provenance: str = "Pycellin") -> Feature:
+    feat = Feature(
+        name="cycle_duration",
+        description="Number of frames in the cell cycle, regardless of gaps",
+        provenance=provenance,
+        feat_type="node",
+        lin_type="CycleLineage",
+        data_type="int",
+        unit="frame",
+    )
+    return feat
 
 
 def level_Feature(provenance: str = "Pycellin") -> Feature:
@@ -574,8 +584,9 @@ def _add_cycle_lineage_features(self) -> None:
         feat_ID = cycle_ID_Feature()
         feat_cells = cells_Feature()
         feat_length = cycle_length_Feature()
+        feat_duration = cycle_duration_Feature()
         feat_level = level_Feature()
-        for feat in [feat_ID, feat_cells, feat_length, feat_level]:
+        for feat in [feat_ID, feat_cells, feat_length, feat_duration, feat_level]:
             if feat.name not in self.feats_dict:
                 self._add_feature(feat)
                 self._protect_feature(feat.name)

pycellin/classes/feature.py

@@ -1218,19 +1218,24 @@ def __init__(self, cell_lineage: CellLineage | None = None) -> None:
             # Adding node and graph features.
             self.graph["lineage_ID"] = cell_lineage.graph["lineage_ID"]
             for n in divs + leaves:
+                cells_in_cycle = cell_lineage.get_cell_cycle(n)
+                first = cells_in_cycle[0]
+                last = cells_in_cycle[-1]
                 self.nodes[n]["cycle_ID"] = n
-                self.nodes[n]["cells"] = cell_lineage.get_cell_cycle(n)
-                self.nodes[n]["cycle_length"] = len(self.nodes[n]["cells"])
+                self.nodes[n]["cells"] = cells_in_cycle
+                # How many cells in the cycle?
+                self.nodes[n]["cycle_length"] = len(cells_in_cycle)
+                # How many frames in the cycle?
+                self.nodes[n]["cycle_duration"] = (
+                    cell_lineage.nodes[last]["frame"]
+                    - cell_lineage.nodes[first]["frame"]
+                ) + 1
                 root = self.get_root()
                 if isinstance(root, list):
                     raise LineageStructureError(
                         "A cycle lineage cannot have multiple roots."
                     )
                 self.nodes[n]["level"] = nx.shortest_path_length(self, root, n)
-            # cell_cycle completeness?
-            # div_time?
-            # cell cycle duration?
-            # Or I add it later with add_custom_feature()?
 
     def __str__(self) -> str:
         name_txt = f" named {self.graph['name']}" if "name" in self.graph else ""

pycellin/classes/lineage.py

@@ -1063,7 +1063,7 @@ def add_branch_total_displacement(
         )
         self.add_custom_feature(motion.BranchTotalDisplacement(feat))
 
-    def add_cell_cycle_completeness(
+    def add_cycle_completeness(
         self,
         rename: str | None = None,
     ) -> None:
@@ -1083,15 +1083,15 @@ def add_cell_cycle_completeness(
             New name for the feature (default is None).
         """
         feat = Feature(
-            name=rename if rename else "cell_cycle_completeness",
+            name=rename if rename else "cycle_completeness",
             description="Completeness of the cell cycle",
             provenance="Pycellin",
             feat_type="node",
             lin_type="CycleLineage",
             data_type="bool",
             unit="none",
         )
-        self.add_custom_feature(tracking.CellCycleCompleteness(feat))
+        self.add_custom_feature(tracking.CycleCompleteness(feat))
 
     def add_cell_displacement(
         self,

pycellin/classes/model.py

@@ -6,7 +6,7 @@
 from .tracking import (
     AbsoluteAge,
     RelativeAge,
-    CellCycleCompleteness,
+    CycleCompleteness,
     DivisionTime,
     DivisionRate,
 )

pycellin/graph/features/__init__.py

@@ -30,9 +30,8 @@
 """
 
 
-from pycellin.classes import CellLineage, CycleLineage
-from pycellin.classes import Feature
-from pycellin.classes import Data
+from pycellin.classes import CellLineage, CycleLineage, Data, Feature
+from pycellin.classes.exceptions import FusionError
 from pycellin.classes.feature_calculator import NodeGlobalFeatureCalculator
 
 # TODO: should I add the word Calc or Calculator to the class names?
@@ -148,7 +147,7 @@ def compute(  # type: ignore[override]
         return age_in_frame * self.time_step
 
 
-class CellCycleCompleteness(NodeGlobalFeatureCalculator):
+class CycleCompleteness(NodeGlobalFeatureCalculator):
     """
     Calculator to compute the cell cycle completeness.
 
@@ -202,14 +201,97 @@ def compute(  # type: ignore[override]
                 return True
 
 
+def _get_cell_lin_frames(lineage: CellLineage, noi: int) -> tuple[int, int]:
+    """
+    Get the frames of the divisions defining the cell cycle.
+
+    This function is used by the DivisionTime and DivisionRate calculators.
+
+    Parameters
+    ----------
+    lineage : CellLineage
+        Lineage graph containing the node of interest.
+    noi : int
+        Node ID (cell_ID) of the cell of interest.
+
+    Returns
+    -------
+    tuple[int, int]
+        Frames of the current and previous division.
+
+    Raises
+    ------
+    KeyError
+        If the cell is not in the lineage.
+    FusionError
+        If the cell has more than one ancestor.
+    """
+    if noi not in lineage.nodes:
+        raise KeyError(f"Cell {noi} not in the lineage.")
+    cells = lineage.get_cell_cycle(noi)
+    frame_current_div = lineage.nodes[cells[-1]]["frame"]
+    ancestors = list(lineage.predecessors(cells[0]))
+    if len(ancestors) > 1:
+        raise FusionError(noi, lineage.graph["lineage_ID"])
+    elif len(ancestors) == 0:
+        frame_prev_div = lineage.nodes[cells[0]]["frame"]
+    else:
+        frame_prev_div = lineage.nodes[ancestors[0]]["frame"]
+    return frame_current_div, frame_prev_div
+
+
+def _get_cycle_lin_frames(
+    data: Data, lineage: CycleLineage, noi: int
+) -> tuple[int, int]:
+    """
+    Get the frames of the divisions defining the cell cycle.
+
+    This function is used by the DivisionTime and DivisionRate calculators.
+
+    Parameters
+    ----------
+    data : Data
+        Data object containing the lineage.
+    lineage : CycleLineage
+        Lineage graph containing the node of interest.
+    noi : int
+        Node ID (cell_ID) of the cell of interest.
+
+    Returns
+    -------
+    tuple[int, int]
+        Frames of the current and previous division.
+
+    Raises
+    ------
+    KeyError
+        If the cycle is not in the lineage.
+    FusionError
+        If the cycle has more than one ancestor.
+    """
+    if noi not in lineage.nodes:
+        raise KeyError(f"Cycle {noi} not in the lineage.")
+    cells = lineage.nodes[noi]["cells"]
+    cell_lin = data.cell_data[lineage.graph["lineage_ID"]]
+    frame_current_div = cell_lin.nodes[cells[-1]]["frame"]
+    ancestors = list(lineage.predecessors(noi))
+    if len(ancestors) > 1:
+        raise FusionError(noi, lineage.graph["lineage_ID"])
+    elif len(ancestors) == 0:
+        frame_prev_div = cell_lin.nodes[cells[0]]["frame"]
+    else:
+        prev_cells = lineage.nodes[ancestors[0]]["cells"]
+        frame_prev_div = cell_lin.nodes[prev_cells[-1]]["frame"]
+    return frame_current_div, frame_prev_div
+
+
 class DivisionTime(NodeGlobalFeatureCalculator):
     """
     Calculator to compute the division time of cells.
 
-    Division time is defined as the time between 2 divisions.
-    It is also the length of the cell cycle of the cell of interest.
-    It is given in frames by default, but can be converted
-    to the time unit of the model if specified.
+    Division time is defined as the time elapsed between the 2 divisions
+    that define the cell cycle. It is given in frames by default, but can
+    be converted to the time unit of the model if specified.
     """
 
     def __init__(self, feature: Feature, time_step: int | float = 1):
@@ -250,14 +332,16 @@ def compute(  # type: ignore[override]
             If the cell or cycle is not in the lineage.
         """
         if isinstance(lineage, CellLineage):
-            if noi not in lineage.nodes:
-                raise KeyError(f"Cell {noi} not in the lineage.")
-            cell_cycle = lineage.get_cell_cycle(noi)
-            return len(cell_cycle) * self.time_step
+            frame_curr_div, frame_prev_div = _get_cell_lin_frames(lineage, noi)
         elif isinstance(lineage, CycleLineage):
-            if noi not in lineage.nodes:
-                raise KeyError(f"Cycle {noi} not in the lineage.")
-            return lineage.nodes[noi]["cycle_length"] * self.time_step
+            frame_curr_div, frame_prev_div = _get_cycle_lin_frames(data, lineage, noi)
+        else:
+            raise TypeError(
+                f"Lineage must be of type CellLineage or CycleLineage, "
+                f"not {type(lineage)}."
+            )
+
+        return (frame_curr_div - frame_prev_div) * self.time_step
 
 
 class DivisionRate(NodeGlobalFeatureCalculator):
@@ -270,17 +354,29 @@ class DivisionRate(NodeGlobalFeatureCalculator):
     to divisions per time unit of the model if specified.
     """
 
-    def __init__(self, feature: Feature, time_step: int | float = 1):
+    def __init__(
+        self, feature: Feature, time_step: int | float = 1, use_div_time: bool = False
+    ):
         """
         Parameters
         ----------
         feature : Feature
             Feature object to which the calculator is associated.
         time_step : int | float, optional
             Time step between 2 frames, in time unit. Default is 1.
+        use_div_time : bool, optional
+            If True, use the division time already computed in the lineage.
+            If False, compute the division time from the lineage. Default is False.
+            The first option is faster but you need to ensure that the division time
+            is computed and updated BEFORE division rate. This can be ensured
+            by adding to the model the division time feature before the division
+            rate feature. Moreover, if `use_div_time` is True, `time_step` will be
+            ignored: division rate will use the division time unit (e.g. if division
+            time is in frames, division rate will be in divisions per frame).
         """
         super().__init__(feature)
         self.time_step = time_step
+        self.use_div_time = use_div_time
 
     def compute(  # type: ignore[override]
         self, data: Data, lineage: CellLineage | CycleLineage, noi: int
@@ -307,15 +403,40 @@ def compute(  # type: ignore[override]
         KeyError
             If the cell or cycle is not in the lineage.
         """
-        if isinstance(lineage, CellLineage):
+        if self.use_div_time:
             if noi not in lineage.nodes:
-                raise KeyError(f"Cell {noi} not in the lineage.")
-            cell_cycle = lineage.get_cell_cycle(noi)
-            return 1 / (len(cell_cycle) * self.time_step)
+                if isinstance(lineage, CellLineage):
+                    lin_txt = "Cell"
+                elif isinstance(lineage, CycleLineage):
+                    lin_txt = "Cycle"
+                else:
+                    raise TypeError(
+                        f"Lineage must be of type CellLineage or CycleLineage, "
+                        f"not {type(lineage)}."
+                    )
+                raise KeyError(f"{lin_txt} {noi} not in the lineage.")
+            try:
+                div_time = lineage.nodes[noi]["division_time"]
+            except KeyError:
+                raise KeyError(
+                    f"Division time not present for cell {noi} in lineage "
+                    f"{lineage.graph['lineage_ID']}."
+                )
+            return 1 / div_time
+
+        if isinstance(lineage, CellLineage):
+            frame_curr_div, frame_prev_div = _get_cell_lin_frames(lineage, noi)
         elif isinstance(lineage, CycleLineage):
-            if noi not in lineage.nodes:
-                raise KeyError(f"Cycle {noi} not in the lineage.")
-            return 1 / (lineage.nodes[noi]["cycle_length"] * self.time_step)
+            frame_curr_div, frame_prev_div = _get_cycle_lin_frames(data, lineage, noi)
+        else:
+            raise TypeError(
+                f"Lineage must be of type CellLineage or CycleLineage, "
+                f"not {type(lineage)}."
+            )
+
+        div_time = (frame_curr_div - frame_prev_div) * self.time_step
+        div_rate = 1 / div_time
+        return div_rate
 
 
 # class CellPhase(NodeGlobalFeatureCalculator):

pycellin/graph/features/tracking.py

@@ -42,7 +42,7 @@ def get_pycellin_cycle_lineage_features() -> dict[str, str]:
             "Mean displacement of the cell during the cell cycle"
         ),
         "branch_mean_speed": "Mean speed of the cell during the cell cycle",
-        "cell_cycle_completeness": (
+        "cycle_completeness": (
             "Completeness of the cell cycle, i.e. does it start and end with a division"
         ),
         "division_time": "Time elapsed between the birth of a cell and its division",