Add parsimony and fitch_count; consolidate linkage tests (v0.3.0)

CHANGELOG.md

@@ -11,6 +11,15 @@ and this project adheres to [Semantic Versioning][].
 ## Unreleased
 
 ### Added
+
+### Changed
+
+### Fixed
+
+## [0.3.0] - 2026-07-08
+
+### Added
+- `pycea.tl.parsimony` and `pycea.tl.fitch_count` for small-parsimony scoring and the FitchCount transition-count algorithm (migrated from Cassiopeia) (#59)
 - `pycea.pl.ancestral_linkage` for plotting the pairwise linkage matrix as a clustered heatmap (#58)
 - Added fast `"sum"` method to `pycea.tl.ancestral_states` (#54, #56)
 - Added additional node ploting customization - `outline_width` and option to directly specify color with hex code (#53)

docs/api.md

@@ -22,6 +22,7 @@
     :toctree: generated
 
     tl.ancestral_states
+    tl.ancestral_linkage
     tl.autocorr
     tl.clades
     tl.compare_distance
@@ -34,6 +35,8 @@
     tl.fitness
     tl.partition_test
     tl.expansion_test
+    tl.parsimony
+    tl.fitch_count
 ```
 
 ## Plotting

docs/conf.py

@@ -87,6 +87,11 @@
 nb_execution_mode = "off"
 nb_merge_streams = True
 typehints_defaults = "braces"
+# Render @typing.overload'd functions as a single entry (the implementation).
+# Otherwise sphinx-autodoc-typehints emits every overload signature, which both
+# drops the short description from the API summary table and clutters each
+# function's page with the (copy-only) overload stubs.
+typehints_document_overloads = False
 
 source_suffix = {
     ".rst": "restructuredtext",

docs/references.bib

@@ -102,3 +102,18 @@ @article{Griffiths_1998
 	keywords = {Age of a mutation, Coalescent process, Population Genetics, Samples of DNA},
 	pages = {273--295},
 }
+
+@article{Quinn_2021,
+	title = {Single-cell lineages reveal the rates, routes, and drivers of metastasis in cancer xenografts},
+	doi = {10.1126/science.abc1944},
+	abstract = {Detailed phylogenies of tumor populations can recount the history and chronology of critical events during cancer progression, such as metastatic dissemination. We applied a Cas9-based, single-cell lineage tracer to study the rates, routes, and drivers of metastasis in a lung cancer xenograft mouse model. We report deeply resolved phylogenies for tens of thousands of cancer cells traced over months of growth and dissemination. This revealed stark heterogeneity in metastatic capacity, arising from preexisting and heritable differences in gene expression. We demonstrate that these identified genes can drive invasiveness and uncovered an unanticipated suppressive role for KRT17. We also show that metastases disseminated via multidirectional tissue routes and complex seeding topologies. Overall, we demonstrate the power of tracing cancer progression at subclonal resolution and vast scale.},
+	language = {en},
+	journal = {Science},
+	volume = {371},
+	number = {6532},
+	pages = {eabc1944},
+	year = {2021},
+	month = feb,
+	author = {Quinn, Jeffrey J. and Jones, Matthew G. and Okimoto, Ross A. and Nanjo, Shigeki and Chan, Michelle M. and Yosef, Nir and Bivona, Trever G. and Weissman, Jonathan S.},
+
+}

pyproject.toml

@@ -4,7 +4,7 @@ requires = [ "hatchling" ]
 
 [project]
 name = "pycea-lineage"
-version = "0.2.0"
+version = "0.3.0"
 description = "Scverse lineage tracing toolkit"
 readme = "README.md"
 license = { file = "LICENSE" }

src/pycea/pl/plot_ancestral_linkage.py

@@ -142,8 +142,9 @@ def ancestral_linkage(
             groupby = candidates[0]
         stats = tdata.uns.get(f"{groupby}{_STATS_SUFFIX}")
         if stats is None:
+            key_added = groupby if groupby is not None else "<groupby>"
             raise KeyError(
-                f"{groupby + _STATS_SUFFIX!r} not found in tdata.uns. "
+                f"{key_added + _STATS_SUFFIX!r} not found in tdata.uns. "
                 f"Run pycea.tl.ancestral_linkage(..., key_added={groupby!r}) first."
             )
         stats = stats.copy()
@@ -256,6 +257,9 @@ def ancestral_linkage(
     ax.set_ylabel("")
 
     if cbar:
-        ax.get_figure().colorbar(mesh, ax=ax)
+        fig = ax.get_figure()
+        if fig is None:
+            raise RuntimeError("no figure associated with ax; cannot add colorbar")
+        fig.colorbar(mesh, ax=ax)
 
     return ax

src/pycea/tl/__init__.py

@@ -6,6 +6,7 @@
 from .fitness import fitness
 from .n_extant import n_extant
 from .neighbor_distance import neighbor_distance
+from .parsimony import fitch_count, parsimony
 from .partition_test import partition_test
 from .sort import sort
 from .topology import expansion_test

src/pycea/tl/_metrics.py

@@ -40,6 +40,19 @@ def __call__(self, a, b):
     def pairwise(self, X, Y):
         return pairwise_distances(X, Y, metric=self.__call__)
 
+
+class _CallableMetric:
+    """Adapts a plain ``metric(a, b) -> float`` callable to the ``.pairwise`` API."""
+
+    def __init__(self, func: _MetricFn):
+        self.func = func
+
+    def __call__(self, a, b):
+        return self.func(a, b)
+
+    def pairwise(self, X, Y):
+        return pairwise_distances(X, Y, metric=self.func)
+
 def _lca_distance(tree, depth_key, node1, node2, lca):
     """Compute the lca distance between two nodes in a tree."""
     if node1 == node2:

src/pycea/tl/ancestral_linkage.py

@@ -4,7 +4,7 @@
 import sys
 import warnings
 from collections import defaultdict
-from collections.abc import Callable, Sequence
+from collections.abc import Callable, Mapping, Sequence
 from typing import Literal, overload
 
 import networkx as nx
@@ -123,7 +123,7 @@ def _max_lca_depth_scores(
 
 def _all_pairs_scores(
     tdata: td.TreeData,
-    trees: dict,
+    trees: dict | Mapping,
     source_leaves_by_tree: dict,
     target_cats: list,
     cat_to_leaves_by_tree: dict,
@@ -135,7 +135,7 @@ def _all_pairs_scores(
     from pycea.tl.tree_distance import tree_distance as _td
 
     tree_keys = list(trees.keys())
-    result = _td(tdata, depth_key=depth_key, metric=metric, tree=tree_keys, copy=True)
+    result = _td(tdata, depth_key=depth_key, metric=metric, tree=tree_keys, copy=True)  # type: ignore
     precomputed = result.toarray() if sp.sparse.issparse(result) else result
 
     scores: dict = {}
@@ -160,7 +160,7 @@ def _all_pairs_scores(
 
 def _compute_scores(
     tdata: td.TreeData,
-    trees: dict,
+    trees: dict | Mapping,
     leaf_to_cat: dict,
     target_cats: list,
     aggregate: str | Callable,
@@ -416,7 +416,7 @@ def _compute_p_values(
 
 def _run_permutation_test(
     tdata: td.TreeData,
-    trees: dict,
+    trees: dict | Mapping,
     leaf_to_cat: dict,
     all_cats: list,
     target_cats: list,
@@ -470,7 +470,7 @@ def _run_permutation_test(
 
 def _run_permutation_test_non_target(
     tdata: td.TreeData,
-    trees: dict,
+    trees: dict | Mapping,
     leaf_to_cat: dict,
     all_cats: list,
     observed_df: pd.DataFrame,
@@ -954,9 +954,9 @@ def _run_single_perm(single_tree, tree_lc, tree_sm, tree_cl, extra_row_fields=No
                         perm_val = cat_null_mean.get(cat, np.nan) if cat is not None else np.nan
                         merged_norm_map[leaf] = (score - perm_val) if not np.isnan(score) else np.nan
 
-            tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(merged_score_map, dtype=float))
+            tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(merged_score_map, dtype=float).to_dict())
             if normalize:
-                tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(merged_norm_map, dtype=float))
+                tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(merged_norm_map, dtype=float).to_dict())
             # Return per-category means from whichever map was written to obs.
             linkage_map = merged_norm_map if normalize else merged_score_map
             if test == "permutation":
@@ -977,9 +977,9 @@ def _run_single_perm(single_tree, tree_lc, tree_sm, tree_cl, extra_row_fields=No
 
         else:
             # Global (non-by_tree) path
-            all_scores = _compute_scores(tdata, trees, leaf_to_cat, [target], single_agg, metric, depth_key)
+            all_scores = _compute_scores(tdata, trees, leaf_to_cat, [target], single_agg, metric, depth_key)  # type: ignore
             score_map = {leaf: scores.get(target, np.nan) for leaf, scores in all_scores.items()}
-            tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(score_map, dtype=float))
+            tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(score_map, dtype=float).to_dict())
 
             if run_perm:
                 rows, cat_null_mean = _run_single_perm(trees, leaf_to_cat, score_map, cat_to_leaves)
@@ -992,7 +992,7 @@ def _run_single_perm(single_tree, tree_lc, tree_sm, tree_cl, extra_row_fields=No
                         else np.nan
                         for leaf, score in score_map.items()
                     }
-                    tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(score_map, dtype=float))
+                    tdata.obs[f"{target}_linkage"] = tdata.obs.index.map(pd.Series(score_map, dtype=float).to_dict())
                 if test == "permutation":
                     test_df = pd.DataFrame(rows)
                     tdata.uns[f"{key_added}_test"] = test_df

src/pycea/tl/ancestral_states.py

@@ -44,39 +44,59 @@ def _remove_node_attributes(tree: nx.DiGraph, key: str) -> None:
             del tree.nodes[node][key]
 
 
-def _reconstruct_fitch_hartigan(
-    tree: nx.DiGraph, key: str, missing: str | None = None, index: int | None = None, fixed_nodes: set | None = None
+def _fitch_hartigan_downpass(
+    tree: nx.DiGraph,
+    key: str,
+    missing: str | None = None,
+    index: int | None = None,
+    fixed_nodes: set | None = None,
+    set_key: str = "value_set",
 ) -> None:
-    """Reconstructs ancestral states using the Fitch-Hartigan algorithm."""
+    """Populate ``tree.nodes[node][set_key]`` with Fitch-Hartigan state sets (down-pass).
+
+    Each node is annotated with the set of states it may take on in an
+    equally-parsimonious solution. Nodes for which all descendants are missing are
+    annotated with the ``missing`` sentinel instead of a set.
+    """
 
     def _is_fixed(node, value):
         return fixed_nodes is not None and node in fixed_nodes and value is not None and value != missing
 
-    # Recursive function to calculate the downpass
     def downpass(node):
         value = _get_node_value(tree, node, key, index)
         # Base case: leaf or fixed observed internal node
         if tree.out_degree(node) == 0 or _is_fixed(node, value):
             if value == missing:
-                tree.nodes[node]["value_set"] = missing
+                tree.nodes[node][set_key] = missing
             else:
-                tree.nodes[node]["value_set"] = {value}
+                tree.nodes[node][set_key] = {value}
         # Recursive case: internal node
         else:
             value_sets = []
             for child in tree.successors(node):
                 downpass(child)
-                value_set = tree.nodes[child]["value_set"]
+                value_set = tree.nodes[child][set_key]
                 if value_set != missing:
                     value_sets.append(value_set)
             if len(value_sets) > 0:
                 intersection = set.intersection(*value_sets)
                 if intersection:
-                    tree.nodes[node]["value_set"] = intersection
+                    tree.nodes[node][set_key] = intersection
                 else:
-                    tree.nodes[node]["value_set"] = set.union(*value_sets)
+                    tree.nodes[node][set_key] = set.union(*value_sets)
             else:
-                tree.nodes[node]["value_set"] = missing
+                tree.nodes[node][set_key] = missing
+
+    downpass(get_root(tree))
+
+
+def _reconstruct_fitch_hartigan(
+    tree: nx.DiGraph, key: str, missing: str | None = None, index: int | None = None, fixed_nodes: set | None = None
+) -> None:
+    """Reconstructs ancestral states using the Fitch-Hartigan algorithm."""
+
+    def _is_fixed(node, value):
+        return fixed_nodes is not None and node in fixed_nodes and value is not None and value != missing
 
     # Recursive function to calculate the uppass
     def uppass(node, parent_state=None):
@@ -100,7 +120,7 @@ def uppass(node, parent_state=None):
 
     # Run the algorithm
     root = get_root(tree)
-    downpass(root)
+    _fitch_hartigan_downpass(tree, key, missing, index, fixed_nodes)
     uppass(root)
     # Clean up
     for node in tree.nodes: