Upgrade to pyham 2 and ete4, make --species_tree optional

.github/workflows/run-tests.yml

@@ -10,9 +10,15 @@ jobs:
   build:
 
     runs-on: ubuntu-latest
+    continue-on-error: ${{ matrix.experimental }}
     strategy:
+      fail-fast: false
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13", "3.14"]
+        experimental: [false]
+        include:
+          - python-version: "3.15"
+            experimental: true
 
     steps:
     - uses: actions/checkout@v4
@@ -21,7 +27,19 @@ jobs:
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies and package
+      if: ${{ !matrix.experimental }}
       run: uv sync --locked
+    # Python versions not yet allowed by requires-python: widen the bound in
+    # this checkout only and resolve dependencies afresh instead of using uv.lock
+    - name: Install dependencies and package (unreleased Python version)
+      if: ${{ matrix.experimental }}
+      env:
+        PY_VERSION: ${{ matrix.python-version }}
+      run: |
+        upper="${PY_VERSION%%.*}.$(( ${PY_VERSION#*.} + 1 ))"
+        sed -i -E "s/^(requires-python = \">=[0-9.]+,<)[0-9.]+\"/\1${upper}\"/" pyproject.toml
+        grep -q "^requires-python = \".*<${upper}\"" pyproject.toml
+        uv sync
     - name: Check edgehog command-line entry point
       run: uv run edgehog --version
     - name: Test with pytest

CLAUDE.md

@@ -0,0 +1,55 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Project
+
+edgeHOG infers ancestral gene order (synteny graphs) at every internal node of a species tree, using HOGs (Hierarchical Orthologous Groups) from an orthoXML file plus extant genome annotations (GFF directory or an OMA HDF5 database). Python >=3.9,<3.15, managed with uv (build backend `uv_build`, flat layout via `module-root = ""`); CLI entry point `edgehog = edgehog.edgehog:main`.
+
+## Commands
+
+```bash
+uv sync --extra oma             # dev install incl. the default `dev` group; the `oma` extra (pyoma + pytables) is needed for --hdf5 input/output
+uv run edgehog --help
+uv lock                         # after changing dependencies; CI uses `uv sync --locked`, so uv.lock must be committed and current
+
+# Reference end-to-end run (the same run tests/test_pipeline.py does)
+uv run edgehog --hog test_data/FastOMA_HOGs.orthoxml \
+    --species_tree test_data/species_tree.nwk \
+    --gff_directory test_data/gff3/ \
+    --date_edges --orient_edges \
+    --output_directory test-results
+
+# Check a known ancestral edge by hand: exactly one matching line must exist
+gunzip -c test-results/12_linearized_synteny_graph_edges.tsv.gz \
+  | awk '/(HOG_HOG:0000285_13\tHOG_HOG:0000284_13)|(HOG_HOG:0000284_13\tHOG_HOG:0000285_13)/{c++} END{exit c==1?0:1}'
+```
+
+```bash
+uv run pytest                   # ~40s; tests/test_pipeline.py runs the test_data pipeline above in-process
+uv run pytest tests/test_pipeline.py::test_known_ancestral_adjacency
+uv run --python 3.9 pytest      # check the oldest supported Python
+```
+
+`tests/test_pipeline.py` runs the test_data pipeline above once per module and checks the same edge. Unit tests per phase (`test_read_gff.py`, `test_extant_graphs.py`, `test_bottom_up.py`, `test_top_down.py`, `test_linearize.py`) use stubs from `tests/synteny_helpers.py` (genes with `.parent`, an ete4 tree with fake pyham genomes patched in via `patch_pyham`) instead of real orthoXML; `test_synthetic_pipeline.py` runs the phases with real pyham on a hand-made 4-species dataset with known ancestral orders. Known deviations from the paper's Supplementary Information and other bugs are documented as `xfail(strict=True)` tests, so a fix makes them fail until the marker is removed. CI (`.github/workflows/run-tests.yml`, Python 3.9–3.14) runs them after `uv run edgehog --version`, which checks the console-script entry point. An `experimental` matrix entry (currently 3.15) tests the next Python version with `continue-on-error`: it widens `requires-python` in its checkout to the matrix version plus one minor and runs `uv sync` without `--locked`; once its dependencies have wheels, move it into the main version list and bump `requires-python`. The `oma` extra is limited to Python <3.13 because pyoma requires numpy<2, which has no wheels for 3.13+, so `--hdf5` only works there. `test_data/oma-hdf5-input/` contains a small OMA HDF5 dataset (`OmaServer.h5`, `edgehog-hogs.orthoxml.gz`, `simplified_speciestree.nwk`) for exercising the `--hdf5` code path. `utils/extract_age_timetree.py` is a standalone post-processing script (dates extant edges in million years using a TimeTree tree).
+
+## Architecture
+
+The pipeline is orchestrated linearly in `edgehog/edgehog.py:main()`. All state lives on a single **pyham `Ham` object**: each stage stores its networkx graphs as ete4 props on the nodes of `ham.taxonomy.tree` (pyham 2 / ete4: write with `tree_node.add_prop(name, value)`, read values edgehog always sets with `tree_node.props[name]`, and the genome with `tree_node.props.get('genome')`, which is None for nodes without genome; attribute access and setting arbitrary attributes on ete4 nodes do not work, `is_leaf` is a property, and missing branch lengths are `None`, which `characterize_species_tree.branch_length` counts as 1.0). Graph nodes are pyham gene/HOG objects (extant genes at leaves, ancestral HOGs at internal nodes), so `gene.parent` maps a node to its HOG in the parent level.
+
+Stages and the tree-node features they produce:
+
+1. `process_hogs.map_hogs_onto_tree` — builds the `Ham` from tree + orthoXML. Whether pyham uses internal node names depends on whether the newick root has a non-numeric name; with `--hdf5`, OMA species resolve mode is used. Without `--species_tree`, pyham builds the tree from the `<taxonomy>` section of the orthoXML (if both exist, the newick wins).
+2. `characterize_species_tree.characterize_tree` — computes per-node stats (RED score, level from root) used in output/genome_dict.tsv.
+3. Extant graphs → `bottom_up_synteny` and `contiguity_dict` on leaves:
+   - GFF path: `init_extant_synteny_graphs.init_extant_graphs` with `read_gff.py`; orthoXML `protId` (or its prefix before a space) must match the GFF CDS `protein_id`, and GFF files must be named after tree leaves.
+   - HDF5 path: `hdf5.init_extant_graphs_from_hdf5` reads OMA main isoforms sorted by chromosome/locus. Genes absent from the orthoXML are skipped (including the first one). Also sets `genome_code` and `taxid`.
+   - Adjacent genes are only linked when on the same contig. Edge attrs: `weight`, `unidirectional`/`convergent`/`divergent` (orientation counts, only filled with `--orient_edges`), `children`, `extant_descendants`.
+4. `infer_ancestral_synteny_graphs.leaves_to_root_synteny_propagation` (postorder) — projects each child's graph onto parent HOGs, summing weights; `build_transitive_graph` bridges over up to `--max_gaps` child genes that have no parent HOG. Sets `bottom_up_synteny` on ancestors and `transitive_synteny` / `real_edge_to_transitive_edge` on children.
+5. `root_to_leaves_edge_trimming` (preorder) — removes edges not supported parsimoniously (seen in only one child lineage) → `top_down_synteny`.
+6. `add_ons.date_edges` (optional, `--date_edges`), then `linearize_graphs` → `linear_synteny`, then `add_ons.phylostratify` (optional) → `stratigraphy`.
+7. `write_output.write_output` writes per-node gzipped TSVs (`<id>_bottom-up|top-down|linearized|extant_synteny_graph_edges.tsv.gz`) plus `genome_dict.tsv` mapping numeric prefixes to tree nodes; `write_as_hdf5` (`--out-format HDF5`) uses `hdf5.HDF5Writer`, which requires the OMA HDF5 DB to map level names to NCBI taxids and HOG IDs to row numbers.
+
+`hdf5.py` imports pyoma/tables at module level and exits if missing, so import it lazily (as `init_extant_synteny_graphs.py` does) to keep the GFF path working without the `oma` extra.
+
+Version is read from installed package metadata (`edgehog/__init__.py`), so bump it in `pyproject.toml`.

Dockerfile

@@ -7,6 +7,9 @@ ENV UV_COMPILE_BYTECODE=1 \
 
 WORKDIR /app
 
+# ete4 is only distributed as source and needs a C compiler to build
+RUN apt-get update && apt-get install -y --no-install-recommends gcc libc6-dev && rm -rf /var/lib/apt/lists/*
+
 # install dependencies first so they are cached independently of source changes
 COPY pyproject.toml uv.lock README.md LICENSE.txt ./
 RUN --mount=type=cache,target=/root/.cache/uv uv sync --locked --no-dev --no-install-project

README.md

@@ -90,7 +90,7 @@ pip install ".[oma]"
 
 ```
 usage: edgehog [-h] [--version] [--output_directory OUTPUT_DIRECTORY] 
-               --species_tree SPECIES_TREE  --hogs HOGS 
+               [--species_tree SPECIES_TREE] --hogs HOGS 
                [--gff_directory GFF_DIRECTORY] [--hdf5 HDF5] [--orient_edges] 
                [--date_edges] [--phylostratify] [--max_gaps MAX_GAPS] 
                [--include_extant_genes] [--out-format {TSV,HDF5}]
@@ -104,7 +104,9 @@ optional arguments:
   --output_directory OUTPUT_DIRECTORY
                         path to output directory (default is ./edgehog_output)
   --species_tree SPECIES_TREE
-                        path to species/genomes phylogenetic tree (newick format)
+                        path to species/genomes phylogenetic tree (newick format). Can be omitted if the HOGs
+                        file contains the species tree as <taxonomy> section; if both are given, the species
+                        tree is used
   --hogs HOGS           path to the HierarchicalGroups.orthoxml file in which HOGs are stored
   --gff_directory GFF_DIRECTORY
                         path to directory with the gffs of extant genomes (each gff file must be named according 
@@ -131,7 +133,7 @@ optional arguments:
 ## Input data
 
 Three types of input data are needed for ```edgeHOG``` to run:
-* a phylogenetic tree of species/genomes of interest (in newick format)
+* a phylogenetic tree of species/genomes of interest (in newick format, or as `<taxonomy>` section of the HierarchicalGroups.orthoxml file)
 * the annotation of each of these genomes (e.g. in the form of a directory of gff files)
 * an HierarchicalGroups.orthoxml file corresponding to the extant genomes
 
@@ -144,11 +146,13 @@ Since these input data are intersected, they must comply with the following requ
 
 A phylogenetic tree of the input genomes/species must be provided in the newick format. If internal nodes are not named, they will be named based on the concatenation of the names of their descendant leaves
 
+Alternatively, the species tree can be taken from the HierarchicalGroups.orthoxml file, if it contains the tree as `<taxonomy>` section (e.g. orthoXML files exported from the OMA browser). In this case, omit ```--species_tree```. If a species tree is given as well, it is used instead of the `<taxonomy>` section. Missing branch lengths count as 1 for the relative evolutionary divergence (RED) scores. The inferred synteny graphs are the same with either input, but the `<taxonomy>` section provides a taxon id for each level, so the HOGs that are not explicitly listed in the orthoXML get ids ending with the taxon id of their own level (e.g. `HOG:0000118_11` at Hominidae instead of `HOG:0000118_13`, the id of its descendant at Homininae); the labels in the output files then always refer to the level they belong to.
+
 * The species phylogenetic tree used in the OMA database can be downloaded [here](https://omabrowser.org/All/speciestree.nwk)
 * The high-quality GTDB archaeal and bacterial species trees, along with metadata can be found [here](https://data.gtdb.ecogenomic.org/releases/latest/)
-* To use the tree or a subtree of the NCBI taxonomy database, the ete3 python package has some [useful build-in functions](http://etetoolkit.org/docs/latest/tutorial/tutorial_ncbitaxonomy.html)
+* To use the tree or a subtree of the NCBI taxonomy database, the ete4 python package has some [useful build-in functions](https://etetoolkit.github.io/ete/tutorial/tutorial_taxonomy.html)
 
-To prune a tree in order to obtain only the phylogeny of your genomes of interest, please refer to the corresponding [ete3 tutorial](http://etetoolkit.org/docs/latest/tutorial/tutorial_trees.html#pruning-trees)
+To prune a tree in order to obtain only the phylogeny of your genomes of interest, please refer to the corresponding [ete4 tutorial](https://etetoolkit.github.io/ete/tutorial/tutorial_trees.html#pruning-trees)
 
 If you don't have a species tree available for your genomes, you can follow [this tutorial](https://github.com/DessimozLab/f1000_PhylogeneticTree) on how to use OMA Browser and OMA standalone for species tree inference.
 

edgehog/add_ons.py

@@ -7,17 +7,16 @@ def date_edges(ham):
     print('###################################')
     print('Predicting age of extant edges ...')
     for tree_node in ham.taxonomy.tree.traverse('preorder'):
-        try:
-            genome = tree_node.genome
-        except AttributeError:
+        genome = tree_node.props.get('genome')
+        if genome is None:
             print('No genome stored for {}'.format(tree_node.name))
             continue
         if isinstance(genome, pyham.AncestralGenome):
-            graph = tree_node.top_down_synteny
+            graph = tree_node.props['top_down_synteny']
         elif isinstance(genome, pyham.ExtantGenome):
-            graph = tree_node.bottom_up_synteny
+            graph = tree_node.props['bottom_up_synteny']
         if tree_node.up:
-            parent_graph = tree_node.up.top_down_synteny
+            parent_graph = tree_node.up.props['top_down_synteny']
         for gene, adjacent_gene, edge_attributes in graph.edges.data():
             if not tree_node.up:
                 graph[gene][adjacent_gene]["lca_node"] = tree_node
@@ -30,25 +29,24 @@ def date_edges(ham):
                         graph[gene][adjacent_gene]["lca_node"] = parent_graph[gene.parent][adjacent_gene.parent]["lca_node"]
                     else:
                         graph[gene][adjacent_gene]["lca_node"] = tree_node
-            graph[gene][adjacent_gene]["lca"] = graph[gene][adjacent_gene]["lca_node"].genome.name
-            graph[gene][adjacent_gene]["age"] = 1 - graph[gene][adjacent_gene]["lca_node"].red
+            graph[gene][adjacent_gene]["lca"] = graph[gene][adjacent_gene]["lca_node"].props['genome'].name
+            graph[gene][adjacent_gene]["age"] = 1 - graph[gene][adjacent_gene]["lca_node"].props['red']
 
 
 def phylostratify(ham):
     print('###################################')
     print('Phylostratigraphy ...')
     for tree_node in ham.taxonomy.tree.traverse('preorder'):
-        try:
-            genome = tree_node.genome
-        except AttributeError:
+        genome = tree_node.props.get('genome')
+        if genome is None:
             print('No genome stored for {}'.format(tree_node.name))
             continue
         if tree_node.up:
             if isinstance(genome, pyham.AncestralGenome):
-                graph = tree_node.top_down_synteny
+                graph = tree_node.props['top_down_synteny']
             elif isinstance(genome, pyham.ExtantGenome):
-                graph = tree_node.bottom_up_synteny
-            parent_graph = tree_node.up.top_down_synteny
+                graph = tree_node.props['bottom_up_synteny']
+            parent_graph = tree_node.up.props['top_down_synteny']
 
             parent_occurrence = dict()
             for gene in graph:
@@ -89,4 +87,4 @@ def phylostratify(ham):
                     except:
                         stratigraphy["lost_edges_due_to_gene_loss"] += 1
 
-            tree_node.add_feature('stratigraphy', stratigraphy)      
+            tree_node.add_prop('stratigraphy', stratigraphy)      

edgehog/characterize_species_tree.py

@@ -6,14 +6,20 @@
 @author: charles
 """
 
+# length of the branch to the parent node; missing branch lengths count as 1.0 (as in ete3),
+# e.g. for newick trees without branch lengths or taxonomies read from the orthoxml
+def branch_length(node):
+    return node.dist if node.dist is not None else 1.0
+
+
 def compute_sum_dist_to_desc_leaves(node, sum_d = 0, n_leaves = 0, n_internal_nodes = 0):
     n_internal_nodes = 0
-    if node.is_leaf():
+    if node.is_leaf:
         n_leaves = 1
-        sum_d = node.get_distance(node.up)
-        node.n_desc_leaves = 0
-        node.sum_dist_to_desc_leaves = 0
-        node.n_internal_nodes = 0
+        sum_d = branch_length(node)
+        node.add_prop('n_desc_leaves', 0)
+        node.add_prop('sum_dist_to_desc_leaves', 0)
+        node.add_prop('n_internal_nodes', 0)
     else:
         n_leaves = 0
         sum_d = 0
@@ -23,35 +29,35 @@ def compute_sum_dist_to_desc_leaves(node, sum_d = 0, n_leaves = 0, n_internal_no
             n_leaves += res[1]
             n_internal_nodes += res[2]
         n_internal_nodes += 1
-        node.sum_dist_to_desc_leaves = sum_d
-        node.n_desc_leaves = n_leaves
-        node.n_internal_nodes = n_internal_nodes
+        node.add_prop('sum_dist_to_desc_leaves', sum_d)
+        node.add_prop('n_desc_leaves', n_leaves)
+        node.add_prop('n_internal_nodes', n_internal_nodes)
         if node.up:
-            sum_d += node.get_distance(node.up) * n_leaves
+            sum_d += branch_length(node) * n_leaves
     if node.up:
         return sum_d, n_leaves, n_internal_nodes
 
 
 def compute_red_score(node, red = 0, level_from_root = 0):
-    if node.is_leaf():
+    if node.is_leaf:
         red = 1
     elif not node.up:
         red = 0
     else:
         p = red
-        d = node.get_distance(node.up)
-        u = (node.sum_dist_to_desc_leaves + (d * node.n_desc_leaves)) / node.n_desc_leaves
+        d = branch_length(node)
+        u = (node.props['sum_dist_to_desc_leaves'] + (d * node.props['n_desc_leaves'])) / node.props['n_desc_leaves']
         red = p + (d/u) * (1-p)
-    node.level = level_from_root
-    node.red = red
+    node.add_prop('level', level_from_root)
+    node.add_prop('red', red)
     for child in node.children:
         compute_red_score(child, red, level_from_root + 1)
 
 
 # the list of nodes per level will be useful for subsequent parallelization of the code
 # def get_nodes_per_level_from_root(node, level_from_root, nodes_per_level):
 #     node.level = level_from_root
-#     if not node.is_leaf():
+#     if not node.is_leaf:
 #         if not level_from_root in nodes_per_level:
 #             nodes_per_level[level_from_root] = list()
 #         nodes_per_level[level_from_root].append(node)

edgehog/check_args.py

@@ -30,7 +30,7 @@ def check_file_extension(file, file_type):
 
 def print_args(args, out_dir):
     print('output directory: \'%s\'' % out_dir)
-    print('species tree:     \'%s\'' % args.species_tree)
+    print('species tree:     \'%s\'' % (args.species_tree if args.species_tree else '<taxonomy> of the HOGs file'))
     print('HOGs.orthoxml:    \'%s\'' % args.hogs)
     print('gff directory:    \'%s\'' % args.gff_directory)
 
@@ -48,7 +48,8 @@ def check_args(args):
     # check output directory exists
     out_dir = check_out_dir(os.path.abspath(args.output_directory))
     # check mandatory input files exist and have the correct extension
-    check_path_exists(args.species_tree)
+    if args.species_tree:
+        check_path_exists(args.species_tree)
     check_path_exists(args.hogs)
     check_file_extension(args.hogs, 'orthoxml')
     if args.gff_directory: