Feat/pypangraph junctions

[mmolari]

The aim of this PR is to add functionalities to PyPangraph to analyze core-genome junctions.

Core-genome junctions are defined as regions in the graph separated by two adjacent flanking core blocks.

The main additions of the PR are:

• a class BackboneJunctions that serves as the interface to functions for junction analysis
• the method stats of this class, returning summary statistics for each junction.
• the method positions that computed junction coordinates in genomes frames of reference.
• the method sequences, that extracts from the graph unaligned sequences for each junction.
• updates to the documentation

[ivan-aksamentov]

Will check later, but I am somewhat weak on both the junctions and the snake language, so not much hope. Pasting AI review in the meantime

---

⚠️ AI-generated content below. Verify all claims.

Overview

Click to expand

This branch adds a junctions sub-package to pypangraph for analyzing accessory genome regions flanked by shared core blocks. It introduces BackboneJunctions as the main entry point, with methods for statistics, genomic positions, and sequence extraction per junction. Supporting this, a new topology_utils module provides OrientedBlock, Walk, and Edge primitives shared by both junctions and the refactored MSU code. Three tutorial pages document the analysis workflow on a real S. aureus dataset.

Changes span the pypangraph Python package (13 production files, 7 test files), documentation (3 new tutorials, image assets, changelog), and a minor README update.

Background

Pangenome junctions and structural variation [click to expand]

In a pangenome graph, genomes are decomposed into blocks (pancontigs) of homologous sequence. A junction is the accessory region between two consecutive core blocks shared across all (or most) genomes. By comparing junctions at the same pair of core flanks across isolates, one can identify insertions, deletions, rearrangements, and mobile element activity. This is the analytical framework described in Molari, Shaw, and Neher (2025) [1].

The branch implements this framework as a Python library, building on pangraph's block/path/node data model. Each genome's path is split at backbone boundaries (core blocks above a length threshold), producing per-isolate junctions keyed by their flanking edge. Reverse-complement symmetry is handled by canonicalizing junction orientation before comparison.

Blocking issues

Correctness concerns worth addressing before merge.

🔴 H1. Walk.__eq__ and __hash__ ignore the circular field [click to expand]

Walk.__eq__ compares only oriented_blocks, ignoring the circular attribute. __hash__ also excludes it. A circular walk [A+, B+] and a linear walk [A+, B+] compare equal and hash identically, despite representing different biological structures (a cycle vs a path). [src]

Effect: Counter(center_paths.values()) in _edge_stats groups walks by equality. Center walks currently have circular=None, so the bug does not fire today. But walk_categories() in the MSU workflow operates on walks where circularity differs between isolates, and any future use of Walk equality in a mixed-circularity context will silently conflate distinct structures.

Fix: include circular in both __eq__ and __hash__:

def __eq__(self, o: object) -> bool:
    if not isinstance(o, Walk):
        return NotImplemented
    return self.oriented_blocks == o.oriented_blocks and self.circular == o.circular

def __hash__(self) -> int:
    return hash((tuple(self.oriented_blocks), self.circular))

Verify test fixtures that construct walks with unset circularity still compare as expected.

🔴 H2. OrientedBlock.from_str_id cannot round-trip underscore-containing IDs [click to expand]

OrientedBlock.from_str_id() splits on _ to separate the block ID from the strand suffix. [src] But minimal_synteny_units() creates block IDs like MSU_0, MSU_1 [src]. Calling to_str_id() on such a block produces MSU_0_f, which from_str_id() cannot unpack:

>>> OrientedBlock.from_str_id("MSU_0_f")
ValueError: too many values to unpack

The same defect applies to Edge.from_str_id() which delegates to OrientedBlock.from_str_id(). [src]

Effect: any downstream code that serializes MSU walks via to_str_id() and attempts to parse them back crashes.

Fix: split from the right:

bid_str, strand_str = t.rsplit("_", 1)

🔴 H3. Junction.__init__ type hints contradict runtime behavior [click to expand]

Junction.__init__ declares left: OrientedBlock and right: OrientedBlock [src], but terminal junctions on linear paths have None flanks. The constructor is called with None at lines 155 and 170 [src], and invert() explicitly handles None [src].

Effect: type checkers in strict mode will not flag unsafe attribute access on left/right without a guard. Library consumers writing typed code will not know None is possible without reading the implementation.

Fix: change the signature to left: OrientedBlock | None, center: Walk, right: OrientedBlock | None. Since from __future__ import annotations is already imported, the X | Y syntax works on Python 3.9+.

Non-blocking issues

Test gaps, API roughness, performance, and documentation. Fix if time allows.

🟡 M1. _ensure_split crashes entire analysis if any genome has < 2 backbone blocks [click to expand]

BackboneJunctions._ensure_split() iterates all genomes and calls path_junction_split() for each [src]. If any single genome has fewer than 2 backbone blocks (plausible with high L_thr), path_junction_split raises ValueError and the entire analysis aborts.

Effect: a single short or divergent genome prevents analysis of all other genomes. The error message does not identify which genome failed.

Fix: catch ValueError per genome with a warning, or add the genome name to the error message so users can exclude it.

🟡 M2. Tutorial exports the wrong JSON filename [click to expand]

The tutorial builds junction.json (singular) on line 146, then runs pangraph export gfa junctions.json (plural) on line 147. [src]

Effect: following the tutorial as written fails because junctions.json was never created.

Fix: change junctions.json to junction.json in the export command.

🟡 M3. Plot function accesses private bj._bdf [click to expand]

linear_junction_plot reads bj._bdf["len"] directly [src], coupling the plot module to a private cache layout in BackboneJunctions.

Effect: if _bdf is renamed or restructured, the plot function breaks silently.

Fix: expose a public accessor on BackboneJunctions for block lengths, or pass block stats explicitly.

🟡 M4. Pandas scalar lookups in backbone splitting and stats hot paths [click to expand]

_is_backbone() performs two pandas .loc lookups per call [src], and path_junction_split calls it twice per block (once in the count pre-pass, once in the splitting loop) [src]. Similarly, _edge_stats() reads block lengths via .loc for each edge and unique accessory block [src].

Effect: the constant factor scales with node count on the documented S. aureus dataset (6,817 nodes, 151 edges). Not a correctness issue, but avoidable overhead.

Fix: precompute a backbone ID set in __init__:

is_bb = self._bdf["core"] & (self._bdf["len"] >= self.L_thr)
self._backbone_ids = set(self._bdf.index[is_bb])

For stats, convert bdf["len"] to a dict once before the edge loop.

🟡 M5. find_mergers has quadratic component merging [click to expand]

When a qualifying edge joins two merger groups, find_mergers scans all entries to relabel one component [src]. For a chain of B core blocks, this gives O(E * B) relabeling work.

Effect: MSU extraction slows on chromosomal graphs with many conserved core adjacencies.

Fix: use union-find for merger groups.

🟡 M6. Missing scientific citations in tutorial t08 [click to expand]

The tutorial on interesting junctions makes specific biological claims about IS256, clfA, staphylocoagulase, and SCCmec without citing primary literature [src]. The SCCmec claim links to Wikipedia.

Effect: scientific claims in user-facing documentation lack traceability.

Fix: add DOI-backed citations. Suggested references:

• IS256 in S. aureus: Lyon, Gillespie, and Skurray 1987. "Detection and Characterization of IS256, an Insertion Sequence in Staphylococcus aureus." Journal of General Microbiology 133 (11): 3031-38. https://doi.org/10.1099/00221287-133-11-3031
• ClfA as adhesion virulence factor: Foster et al. 2014. "Adhesion, Invasion and Evasion: The Many Functions of the Surface Proteins of Staphylococcus aureus." Nature Reviews Microbiology 12 (1): 49-62. https://doi.org/10.1038/nrmicro3161
• Staphylocoagulase as virulence factor: Cheng et al. 2010. "Contribution of Coagulases towards Staphylococcus aureus Disease and Protective Immunity." PLoS Pathogens 6 (8): e1001036. https://doi.org/10.1371/journal.ppat.1001036
• SCCmec (replacing Wikipedia link): IWG-SCC 2009. "Classification of Staphylococcal Cassette Chromosome Mec (SCCmec): Guidelines for Reporting Novel SCCmec Elements." Antimicrobial Agents and Chemotherapy 53 (12): 4961-67. https://doi.org/10.1128/AAC.00579-09

🟡 M7. Documentation prose style [click to expand]

The new tutorials and CHANGELOG use em dashes (U+2014) throughout (~20 instances across 4 files), and t07-junction-stats.md uses bold definition labels. Em dashes are a common marker of AI-generated prose; human-written technical documentation typically uses commas, periods, or parentheses instead. Bold on every label in a definition list creates uniform visual weight, so nothing stands out when scanning for a specific column name.

Fix: replace em dashes with -- or restructure sentences. Use plain labels in definition lists.

🟡 M8. TODO placeholder in production stats [click to expand]

junction_stats contains a speculative TODO at line 102 [src].

Fix: remove if the current stat set is intentional, or convert to a tracked issue.

Low

Cosmetic or minor issues.

🔵 L1. Parameter named node_id receives a block ID [click to expand]

The inner function is_core(node_id) in core_paths actually receives a block_id from filter_walks calling keep_f(ob.id) [src].

Fix: rename to bid or block_id.

🔵 L2. Implicit exception chaining in IndexedCollection.__getitem__ [click to expand]

except KeyError: raise KeyError(...) produces a confusing double traceback [src].

Fix: use explicit chaining (raise KeyError(...) from ex) to preserve the original key for debugging while producing a cleaner "direct cause" traceback instead of "during handling of."

Dismissed

Reported by automated reviewers, investigated and found non-issues.

⚪ X1. Core-genome alignment uses guide strand for all strains [click to expand]

Reported: core_genome_alignment() orients every strain's sequence using only the guide strain's strand, ignoring per-strain node orientation. Claimed this produces wrong alignments for inverted core blocks.

Dismissed: Block.to_alignment() returns sequences in consensus orientation regardless of the node's genomic strand. The edits stored per node are in consensus coordinates. The guide_strand flag determines whether to present the entire block in consensus orientation (forward) or reverse-complement, so the alignment follows the guide strain's genomic order. Applying per-node strand corrections would double-flip sequences and produce wrong results. The current behavior is correct for core-genome multiple sequence alignment.

⚪ X2. Hash over mutable fields on topology value objects [click to expand]

Reported: OrientedBlock, Walk, Edge, and Junction define __hash__ over public mutable attributes. Mutating after insertion into hashed containers breaks lookup.

Dismissed: these are value objects used as dictionary keys and Counter entries. The library never mutates them after hashing. Walk.add_left()/add_right() are construction helpers used before the walk is inserted into any container. Mutation after hashing would be a caller error, not a library defect.

⚪ X3. is_singleton definition at n=2 [click to expand]

Reported: with exactly 2 isolates and 2 different variants, is_singleton yields True even though neither variant has a meaningful "majority."

Dismissed: the formula n_majority_category == n_isolates - 1 is mathematically correct and consistent with its docstring. At n=2, the flag correctly identifies that one variant is unique. This is a useful signal for downstream filtering and users can add their own n_isolates >= 3 guard if needed.

⚪ X4. Palindromic edge hash collision [click to expand]

Reported: Edge.__hash__ and Junction.__hash__ use XOR of forward and inverted side hashes. Palindromic structures (where both sides hash identically) all hash to 0.

Dismissed: RC-palindromic edges are rare in bacterial pangenomes. Hash collision degrades performance (O(n) probe chains) but does not affect correctness. The collision rate in practice is negligible.

⚪ X5. pypangraph.msu module renamed without compatibility shim [click to expand]

Reported: deleting pypangraph.msu and moving to pypangraph.minimal_synteny_units breaks downstream users.

Dismissed: the CHANGELOG documents this as an intentional change in version 1.1.0, and pypangraph.__init__ exports minimal_synteny_units at the top level. The package uses semver minor bumps for API additions. A smoke test verifies the new import paths.

Notes

Positive observations about the change.

Click to expand

• __eq__ methods guard against foreign types returning NotImplemented throughout: OrientedBlock, Walk, Edge, Junction [src]
• Junction.invert() handles None flanks for terminal junctions on linear paths [src]
• The is_canonical()/to_canonical() pattern with explicit ValueError on terminal junctions prevents silent misuse [src]
• PangraphLoadError with chained exceptions replaces bare print(ex) and naked assert [src]
• Test fixtures in conftest.py are well-designed: the junction, linear, sequence, and inversion pangraphs each exercise distinct edge cases with hand-verifiable expected values
• Sequence extraction correctly applies reverse-complement based on post-canonicalization strand [src]
• Smoke tests on real data (plasmids.json, staph.json.gz) complement synthetic fixtures

Glossary

Click to expand

• backbone block: a core block whose consensus length meets a minimum threshold (L_thr)
• canonical orientation: a deterministic direction chosen for an edge or junction to ensure reverse-complement equivalents are recognized as identical
• center path: the walk of accessory blocks between two flanking core blocks in a junction
• core block: a block present in every genome in the pangenome
• edge: the pair of oriented core blocks flanking a junction
• junction: a segment of accessory genome flanked by two backbone blocks, with RC symmetry
• MSU (minimal synteny unit): a maximal contiguous set of core blocks always adjacent across all genomes
• walk: an ordered list of oriented blocks representing a genome's traversal through the block graph

References

Click to expand

[1] Molari, Marco, Liam P. Shaw, and Richard A. Neher. 2025. "Quantifying the Evolutionary Dynamics of Structure and Content in Closely Related E. coli Genomes." Molecular Biology and Evolution 42 (1): msae272. https://doi.org/10.1093/molbev/msae272

[mmolari]

Thank you! Still very appreciated! I'll go through these and after this polishing I'll release next week!