diff --git a/README.md b/README.md
index 1789b98a..06d74272 100644
--- a/README.md
+++ b/README.md
@@ -74,7 +74,7 @@ PyPangraph is a python package with convenient utilities to load and explore the
```python
import pypangraph as pp
-graph = pp.Pangraph.load_graph("graph.json")
+graph = pp.Pangraph.from_json("graph.json")
print(graph)
# pangraph object with 15 paths, 137 blocks and 1042 nodes
```
diff --git a/docs/docs/assets/pp_t6_core_edge_scheme.png b/docs/docs/assets/pp_t6_core_edge_scheme.png
new file mode 100644
index 00000000..f814b40c
Binary files /dev/null and b/docs/docs/assets/pp_t6_core_edge_scheme.png differ
diff --git a/docs/docs/assets/pp_t6_junction_scheme.png b/docs/docs/assets/pp_t6_junction_scheme.png
new file mode 100644
index 00000000..d368c16b
Binary files /dev/null and b/docs/docs/assets/pp_t6_junction_scheme.png differ
diff --git a/docs/docs/assets/pp_t7_linear_junction_plot.png b/docs/docs/assets/pp_t7_linear_junction_plot.png
new file mode 100644
index 00000000..137b177c
Binary files /dev/null and b/docs/docs/assets/pp_t7_linear_junction_plot.png differ
diff --git a/docs/docs/assets/pp_t7_linear_junction_plot_hotspot.png b/docs/docs/assets/pp_t7_linear_junction_plot_hotspot.png
new file mode 100644
index 00000000..8fdaee08
Binary files /dev/null and b/docs/docs/assets/pp_t7_linear_junction_plot_hotspot.png differ
diff --git a/docs/docs/assets/pp_t7_stats_scheme.png b/docs/docs/assets/pp_t7_stats_scheme.png
new file mode 100644
index 00000000..0e3ceb10
Binary files /dev/null and b/docs/docs/assets/pp_t7_stats_scheme.png differ
diff --git a/docs/docs/assets/pp_t7_stripplot.png b/docs/docs/assets/pp_t7_stripplot.png
new file mode 100644
index 00000000..acb76175
Binary files /dev/null and b/docs/docs/assets/pp_t7_stripplot.png differ
diff --git a/docs/docs/assets/pp_t7_transitive_scheme.png b/docs/docs/assets/pp_t7_transitive_scheme.png
new file mode 100644
index 00000000..78637ea7
Binary files /dev/null and b/docs/docs/assets/pp_t7_transitive_scheme.png differ
diff --git a/docs/docs/assets/pp_t8_bandage.png b/docs/docs/assets/pp_t8_bandage.png
new file mode 100644
index 00000000..2e0e501c
Binary files /dev/null and b/docs/docs/assets/pp_t8_bandage.png differ
diff --git a/docs/docs/assets/pp_t8_mafft.png b/docs/docs/assets/pp_t8_mafft.png
new file mode 100644
index 00000000..e04e776a
Binary files /dev/null and b/docs/docs/assets/pp_t8_mafft.png differ
diff --git a/docs/docs/assets/pp_t8_ncbi_viewer_IS_screenshot.png b/docs/docs/assets/pp_t8_ncbi_viewer_IS_screenshot.png
new file mode 100644
index 00000000..f527d5c9
Binary files /dev/null and b/docs/docs/assets/pp_t8_ncbi_viewer_IS_screenshot.png differ
diff --git a/docs/docs/assets/pp_t8_ncbi_viewer_interrupted_gene.png b/docs/docs/assets/pp_t8_ncbi_viewer_interrupted_gene.png
new file mode 100644
index 00000000..925a86a6
Binary files /dev/null and b/docs/docs/assets/pp_t8_ncbi_viewer_interrupted_gene.png differ
diff --git a/docs/docs/assets/pp_t8_positions_scheme.png b/docs/docs/assets/pp_t8_positions_scheme.png
new file mode 100644
index 00000000..aa72aaac
Binary files /dev/null and b/docs/docs/assets/pp_t8_positions_scheme.png differ
diff --git a/docs/docs/assets/pp_t9_large_pangraph.png b/docs/docs/assets/pp_t9_large_pangraph.png
new file mode 100644
index 00000000..f03fab4a
Binary files /dev/null and b/docs/docs/assets/pp_t9_large_pangraph.png differ
diff --git a/docs/docs/assets/pp_t9_scaffold_all.png b/docs/docs/assets/pp_t9_scaffold_all.png
new file mode 100644
index 00000000..939a6864
Binary files /dev/null and b/docs/docs/assets/pp_t9_scaffold_all.png differ
diff --git a/docs/docs/assets/pp_t9_scheme.png b/docs/docs/assets/pp_t9_scheme.png
new file mode 100644
index 00000000..edf63b8c
Binary files /dev/null and b/docs/docs/assets/pp_t9_scheme.png differ
diff --git a/docs/docs/assets/pp_t9_tangled.png b/docs/docs/assets/pp_t9_tangled.png
new file mode 100644
index 00000000..9685bf11
Binary files /dev/null and b/docs/docs/assets/pp_t9_tangled.png differ
diff --git a/docs/docs/assets/pp_t9_untangled.png b/docs/docs/assets/pp_t9_untangled.png
new file mode 100644
index 00000000..a1f561dc
Binary files /dev/null and b/docs/docs/assets/pp_t9_untangled.png differ
diff --git a/docs/docs/assets/svgs/junctions.svg b/docs/docs/assets/svgs/junctions.svg
new file mode 100644
index 00000000..fb4512f5
--- /dev/null
+++ b/docs/docs/assets/svgs/junctions.svg
@@ -0,0 +1,4545 @@
+
+
+
+
diff --git a/docs/docs/assets/t1_gfa_untangled.png b/docs/docs/assets/t1_gfa_untangled.png
new file mode 100644
index 00000000..bf85f704
Binary files /dev/null and b/docs/docs/assets/t1_gfa_untangled.png differ
diff --git a/docs/docs/pypangraph/t01-load-graph.md b/docs/docs/pypangraph/t01-load-graph.md
index 672f2704..192c29b0 100644
--- a/docs/docs/pypangraph/t01-load-graph.md
+++ b/docs/docs/pypangraph/t01-load-graph.md
@@ -4,7 +4,7 @@ sidebar_position: 3
# Loading and exploring a graph
-We start the tutorial by loading a pangenome graph object and exploring its properties. For this tutorial we will use the `plasmids.json` file that you can find in pypangraph repository under `tests/data/plasmids.json`
+We start the tutorial by loading a pangenome graph object and exploring its properties. For this tutorial we will use the `plasmids.json` file, a pangraph of 15 plasmid sequences. You can download it from the pypangraph repository at [`packages/pypangraph/tests/data/plasmids.json`](https://raw.githubusercontent.com/neherlab/pangraph/master/packages/pypangraph/tests/data/plasmids.json).
We can load the graph object with:
diff --git a/docs/docs/pypangraph/t04-core-synteny.md b/docs/docs/pypangraph/t04-core-synteny.md
index ca72a806..c3e8c043 100644
--- a/docs/docs/pypangraph/t04-core-synteny.md
+++ b/docs/docs/pypangraph/t04-core-synteny.md
@@ -71,10 +71,12 @@ From this plot we observe a strong conservation in the order of core blocks. Thi
pangraph export gfa --no-duplicated --minimum-depth 15 plasmids.json -o plasmids_core.gfa
```
-Moreover, we can save the block colors that we used in the previous plot in a csv file, that can be loaded by Bandage to color the blocks.
+Moreover, we can save the block colors that we used in the previous plot in a csv file, that can be loaded by Bandage to color the blocks. We quote every field (`quoting=csv.QUOTE_ALL`) so that Bandage matches the long numeric block ids as strings instead of misreading them as numbers.
```python
-pd.Series(block_color, name="Colour").to_csv("block_colors.csv")
+import csv
+
+pd.Series(block_color, name="Colour").to_csv("block_colors.csv", quoting=csv.QUOTE_ALL)
```
After loading the graph and coloring it we obtain the following picture:
@@ -92,9 +94,11 @@ For these cases, pypangraph provides a method to quickly survey all changes in c
-For this part of the tutorial we will analyze the `graph.json` file created [in the first tutorial](../tutorial/t01-building-pangraph.md#building-the-pangraph), containing 10 _E. coli_ chromosomes. The minimal sinteny units for this graph can be extracted with the function:
+For this part of the tutorial we will analyze the `graph.json` file created [in the first tutorial](../tutorial/t01-building-pangraph.md#building-the-pangraph), containing 10 _E. coli_ chromosomes. The minimal synteny units for this graph can be extracted with the function:
```python
+import pypangraph as pp
+
graph = pp.Pangraph.from_json("graph.json")
# find MSUs
@@ -105,7 +109,7 @@ MSU_mergers, MSU_paths, MSU_len = pp.minimal_synteny_units(graph, threshold_len)
This returns three objects:
- `MSU_mergers`: a dictionary where keys are core block ids and the values are the ids of the MSU they belong to.
-- `MSU_paths`: a dictionary where keys are path ids and values are paths composed of MSUs instead of blocks.
+- `MSU_paths`: a dictionary where keys are path ids and values are `pypangraph.topology_utils.Walk` objects composed of MSUs instead of blocks.
- `MSU_len`: a list of the lengths of the MSUs in basepairs, i.e. the sum of consensus length of the core blocks that compose them.
We can draw a linear representation for paths in terms of the MSUs with the following code, in which each MSU is represented as a colored block of unit size. Arrows indicate inversions.
@@ -123,9 +127,9 @@ colors = defaultdict(lambda: next(color_generator))
fig, ax = plt.subplots(figsize=(8, 5))
for i, (iso, path) in enumerate(MSU_paths.items()):
- for j, node in enumerate(path.nodes):
- ax.barh(i, 1, left=j, color=colors[node.id])
- if not node.strand:
+ for j, ob in enumerate(path.oriented_blocks):
+ ax.barh(i, 1, left=j, color=colors[ob.id])
+ if not ob.strand:
ax.arrow(j + 1, i, -0.8, 0, head_width=0.2, head_length=0.2)
ax.set_yticks(range(len(MSU_paths)))
ax.set_yticklabels(list(MSU_paths.keys()))
@@ -149,15 +153,18 @@ pangraph export gfa \
And then we can export the dictionary of core-block colors with:
```python
+import csv
+
block_colors = {}
for block_id in graph.blocks.keys():
if block_id in MSU_mergers:
block_colors[block_id] = mpl.colors.to_hex(colors[MSU_mergers[block_id]])
else:
block_colors[block_id] = mpl.colors.to_hex("lightgray")
-pd.Series(block_colors, name="Colour").to_csv("block_colors.csv")
+# quote every field so Bandage matches the long numeric block ids as strings
+pd.Series(block_colors, name="Colour").to_csv("block_colors.csv", quoting=csv.QUOTE_ALL)
```
After loading the graph in Bandage and coloring the blocks we obtain the following picture:
-
\ No newline at end of file
+
diff --git a/docs/docs/pypangraph/t06-junctions-intro.md b/docs/docs/pypangraph/t06-junctions-intro.md
new file mode 100644
index 00000000..23d8956a
--- /dev/null
+++ b/docs/docs/pypangraph/t06-junctions-intro.md
@@ -0,0 +1,127 @@
+---
+sidebar_position: 8
+---
+
+# Junctions: comparing local accessory variation
+
+When comparing closely related bacterial genomes, the core genome is often largely **syntenic**: long stretches of conserved blocks appear in the same order across isolates. Between these conserved blocks, segments of **accessory** DNA (insertions, deletions, mobile elements) vary from one genome to another.
+
+The largely conserved order of core segments provides a natural **frame of reference** to meaningfully break down and compare local accessory variation across isolates. To this end we introduce the concept of a **junction**: a region of the graph delimited by two consecutive core blocks.
+
+## What is a junction?
+
+A **junction** is defined as the part of the graph "between" two consecutive core-blocks. It starts and ends with the **two flanking core blocks** that serve as anchor points, and contains all accessory diversity in-between.
+
+
+
+### Unique identifiers for junctions: core edges
+
+A graph can contain a large number of junctions, but each is uniquely defined by the two flanking core blocks. Their identity and _orientation_ determines the identity of the junction.
+
+We call a **core edge** the oriented adjacency relation of two core blocks, neglecting accessory blocks in-between. In the example below, core blocks `X` and `Y` form a core edge with orientation `[X+|Y-]`.
+
+This notation indicates that in genomes we observe first a sequence homologous to the consensus sequence of block `X` (denoted as `X+`) and then (potentially after an accessory region) a sequence homologous to the _reverse complement_ of the consensus sequence of block `Y` (denoted as `Y-`).
+
+The same junction region can occur with different strandedness on different genomes, but it should still be identified as the same junction.
+For this reason we define core-edges to be **invariant under reverse-complementation**, so `[X+|Y-]` and its reverse complement `[Y+|X-]` correspond to the same core edge.
+
+
+
+In pypangraph core edges are identified by a string of the form `X_a__Y_b`, where `X` and `Y` are large integers corresponding to the `block_id` of the two blocks, and `a` and `b` are either `f` or `r`, indicating respectively forward (`+`) or reverse (`-`) orientation.
+
+Since edges are equivalent under reverse-complementation, there are actually two possible strings that can identify an edge. To break the tie we take the lexicographically smaller. This order defines the **canonical orientation** of the junction.
+
+## The `junctions` module in pypangraph
+
+Pypangraph has a `junctions` module to facilitate junction analyses.
+
+In this part of the tutorial we will load and explore the `staph.json.gz` file, which contains a pangraph of 15 _Staphylococcus aureus_ chromosomes. You can download it from the pypangraph repository at [`packages/pypangraph/tests/data/staph.json.gz`](https://raw.githubusercontent.com/neherlab/pangraph/master/packages/pypangraph/tests/data/staph.json.gz).
+
+
+ how was `staph.json.gz` built?
+
+ The graph was built from 15 complete _S. aureus_ chromosomes downloaded from NCBI. To reproduce it, fetch the sequences and run `pangraph build`:
+
+ ```bash
+ ACCESSIONS="NZ_CP132362.1,NZ_LR822061.1,NZ_CP077852.1,NZ_CP162433.1,NZ_CP034011.1,NZ_CP092558.1,NZ_CP062358.1,NZ_CP132372.1,NZ_AP024511.1,NZ_CP181145.1,NZ_CP157420.1,NZ_CP080550.1,NZ_CP169947.1,NZ_CP022905.1,NZ_CP035791.1"
+
+ # download the sequences in FASTA format
+ curl -L "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=${ACCESSIONS}&rettype=fasta&retmode=text" > staph.fa
+
+ # build the pangraph
+ pangraph build --circular -l 100 -s 20 -b 5 staph.fa -o staph.json.gz
+ ```
+
+
+We start by loading the graph and creating a `BackboneJunctions` object. This class identifies junctions in the graph and provides methods to analyze them.
+
+```python
+import pypangraph as pp
+
+graph = pp.Pangraph.from_json("staph.json.gz")
+print(graph)
+# pangraph object with 15 paths, 664 blocks and 6817 nodes
+
+junctions = pp.junctions.BackboneJunctions(graph, L_thr=500)
+```
+
+The constructor takes only two arguments: the graph object and a length threshold `L_thr` (default 500 bp). Core blocks shorter than this threshold are too short to be reliable anchors and are treated as accessory blocks for the purpose of junction definition.
+
+### A look at edges
+
+The graph has a total of 143 valid anchor core blocks.
+
+```python
+graph.to_blockstats_df().query("len >= 500")["core"].sum()
+# 143 core blocks
+```
+
+We can get the list of edges with:
+
+```python
+print(junctions.edges())
+# ['13733442150340492168_f__17042526223432838337_f',
+# '12427448985183016017_f__13733442150340492168_f',
+# '12427448985183016017_r__14846644350907526570_f',
+# '11809679528571820295_f__14846644350907526570_r',
+# ...
+# ]
+# a total of 151 edges
+```
+
+There are 151 edges, which is only slightly more than the number of core blocks. This means that, as expected, the order of core blocks is strongly conserved and almost every core block has the same core block neighbors in every genome.
+Across the whole dataset we expect to observe only a handful of synteny changes.
+
+### Selecting a junction and a single path through it
+
+To get a feel for what a junction actually looks like, let's pick one and inspect it. Indexing into the `BackboneJunctions` object by edge id returns the `{isolate: Junction}` mapping for that edge; a second index by isolate name gives the path of a single isolate through the junction:
+
+```python
+# core edge
+edge_str = "3156970751805415521_f__4335229004353524956_f"
+# select the junction corresponding to that edge
+junction = junctions[edge_str]
+# select the path of a single isolate through that junction
+junction_path = junction["NZ_CP162433.1"]
+```
+
+Each junction path contains the left flanking core block, the center walk of accessory blocks, and the right flanking core block, accessible as attributes:
+
+```python
+print(junction_path.left)
+# {block=4335229004353524956|-}
+print(junction_path.right)
+# {block=3156970751805415521|-}
+print(junction_path.center)
+# {block=8061287138899943998|-} {block=12150994386844653378|-}
+print(len(junction_path.center))
+# 2
+```
+
+Each block occurrence is printed as `{block=|}`.
+
+Note that both flanking blocks appear on the reverse strand (`-`) and in swapped order with respect to the canonical edge ID we asked for (`3156970751805415521_f__4335229004353524956_f`). This is the reverse-complement symmetry introduced earlier: this isolate carries the junction in its reverse orientation, but the edge ID (the canonical, orientation-invariant identifier) is the same.
+
+## Next: from a single junction to summary statistics
+
+Inspecting individual junctions like the one above gives a concrete sense of what a junction is, but a graph contains _hundreds of them_, and going through one at a time quickly becomes impractical. To get an overview of the structural variation landscape, and to spot the junctions worth zooming in on, we need **summary statistics** across all junctions. This will be the focus of the [next part of the tutorial](t07-junction-stats.md).
\ No newline at end of file
diff --git a/docs/docs/pypangraph/t07-junction-stats.md b/docs/docs/pypangraph/t07-junction-stats.md
new file mode 100644
index 00000000..9748febb
--- /dev/null
+++ b/docs/docs/pypangraph/t07-junction-stats.md
@@ -0,0 +1,219 @@
+---
+sidebar_position: 9
+---
+
+# Calculating summary junction statistics
+
+Bacterial genomes can harbor hundreds of loci of accessory genome variability. Manual inspection of every locus is often impractical.
+
+[In the course of our work](https://doi.org/10.1093/molbev/msae272) we found it instructive to calculate **summary statistics** for these regions, to visualize large-scale patterns in the investigated collection. Based on these statistics, one can then pick regions of interest for more detailed inspection.
+
+Useful per-junction summaries include, for example:
+
+- the total number of unique "accessory paths" found within the junction (including the empty one). This is an indication of the **structural diversity** of the junction. We call this the number of **path categories**. In the example below, we find 4 categories across 5 paths, since one category appears in two paths.
+- the **total length of accessory genome** found in the junction. This can be calculated by summing the consensus length of all unique accessory blocks found in the junction. In the example below, the junction's total accessory length is roughly 3kb. This number gives an idea of the amount of accessory material that the junction harbors, and can help, for example, to distinguish recent changes associated with particular mobile genetic elements by typical size.
+- the number of **non-empty paths**, i.e. paths that contain at least one accessory block. In the example below, this is every path except for the last one, which only has the flanking core block. Comparing this number to the total number of genomes gives, for example, an indication of whether a junction was caused by a **recent insertion**. In this case the number of non-empty (or _occupied_) paths is expected to be very small compared to the dataset size. A recent deletion would conversely show up as a junction where the number of empty paths is very small.
+
+
+
+## Computing junction statistics with pypangraph
+
+Pypangraph provides a convenient way to quickly calculate summary statistics for all junctions. The `BackboneJunctions.stats()` method returns a `pandas.DataFrame` with one row per junction identified by the core edge:
+
+```python
+import pypangraph as pp
+
+graph = pp.Pangraph.from_json("staph.json.gz")
+junctions = pp.junctions.BackboneJunctions(graph, L_thr=500)
+
+stats = junctions.stats()
+print(stats)
+# n_isolates n_non_empty n_categories accessory_length ...
+# edge
+# 13733442150340492168_f__17042526223432838337_f 15 15 8 5773 ...
+# 12427448985183016017_f__13733442150340492168_f 15 1 2 1512 ...
+# 11809679528571820295_r__14906387308163561070_r 15 14 2 1242 ...
+```
+
+The dataframe carries nine columns (see the dropdown below for the full reference), but the three most important ones are those described above:
+
+- `n_categories`: number of distinct accessory paths observed at the junction.
+- `accessory_length`: total unique accessory content (bp) summed across all distinct accessory blocks ever seen at the junction.
+- `n_non_empty`: out of the isolates that share the edge, how many actually carry accessory content between the two flanking core blocks (the rest have the two backbone blocks sitting directly adjacent).
+
+In addition to this, the `n_isolates` column indicates in how many isolates the junction was found. Cases where this number is smaller than the total number of genomes typically indicate synteny changes in some. This is discussed further in the "_transitive junctions_" dropdown below.
+
+
+ **full column reference**
+
+ - `n_isolates`: number of isolates that have this junction. When `n_isolates` equals the total number of isolates in the graph, the junction is universal: the flanking backbone blocks appear consecutively in all genomes. Non-universal junctions are typically a sign of synteny breaks.
+ - `n_non_empty`: number of isolates where the junction is non-empty (carries at least one accessory block). The complement, `n_isolates - n_non_empty`, is the count of isolates where the two flanking backbone blocks sit directly adjacent with no accessory content in between.
+ - `n_categories`: number of distinct accessory path variants. A "category" is a unique sequence of accessory block IDs. All isolates with no accessory blocks (empty center) count as one category.
+ - `n_majority_category`: number of isolates in the most common variant. Together with `n_categories`, this tells you how diverse the junction is.
+ - `is_transitive`: `True` if `n_categories == 1`, i.e. all isolates sharing this edge have the same accessory structure (including all-empty).
+ - `is_singleton`: `True` if exactly one isolate has a different variant from all others (`n_isolates > 1` and `n_majority_category == n_isolates - 1`).
+ - `left_core_length` / `right_core_length`: consensus length (bp) of the left and right flanking backbone blocks.
+ - `accessory_length`: total unique accessory content, the sum of consensus lengths of all distinct accessory block IDs appearing in any isolate's center path for this edge. Each block is counted once even if it appears in multiple isolates.
+
+
+
+
+ **transitive junctions**
+
+ You might have noticed in the plot above that some junctions have `n_categories == 1`, i.e. we find only one unique accessory structure (sometimes empty) between the two flanking core blocks. We call these **transitive junctions**.
+
+ How can a junction be transitive? This often happens in one of two ways.
+
+ One option is the presence of fixed paralogs in the dataset. Pangraph can detect the homology of the paralogs and identify them as repeated accessory blocks. All orthologous copies of the paralog might appear in the same conserved core context, and as such in a "transitive" junction, which does not represent a locus of recent accessory variation.
+
+ A second possibility are **synteny breaks**. Changes in the order or orientation of core blocks generate edges that are present in only a subset of isolates, even if the flanking core blocks are present in every isolate.
+
+ In the scheme below, the inversion of core block `Y` generates two possible patterns: either genomes possess edges `[X+|Y+]` and `[Y+|Z+]`, or they have `[X+|Y-]` and `[Y-|Z+]`. All of these junctions, with the exception of `[Y+|Z+]`, are transitive. They are generated by the changes in synteny, and not by recent accessory genome variation between the flanking core blocks.
+
+ 
+
+
+
+## Visualizing the junction landscape
+
+A simple way to visualize the distribution of the three interesting quantities described above at a glance is a strip plot:
+
+```python
+import seaborn as sns
+import matplotlib.pyplot as plt
+
+fig, ax = plt.subplots(figsize=(7, 4.5))
+sns.stripplot(
+ data=stats,
+ x="accessory_length",
+ y="n_categories",
+ orient="y",
+ jitter=0.25,
+ alpha=0.8,
+ hue="n_non_empty",
+ palette="coolwarm",
+ log_scale=(True, False),
+ ax=ax,
+)
+ax.set_xlabel("Accessory length (bp)")
+ax.set_ylabel("Number of path categories")
+ax.invert_yaxis()
+ax.legend(title="n. non-empty", loc="upper left")
+```
+
+
+
+Reading the plot:
+
+- most variation is in junctions with a **low number of categories** (typically around 2), the loci of **limited structural variation**.
+ - amongst these, we find bands of several junctions with characteristic lengths around 1500 and 1300 bp. These are junctions with very few occupied genomes (blue dots), consistent with recent and repeated activity of mobile elements such as _Insertion Sequences_.
+- At the other end of the spectrum, in the top-right corner of the plot, we find **hotspots**. These are regions with high variability (almost every genome has a unique accessory pattern) and a vast accessory repertoire (around 100kb of unique accessory genome across 15 isolates).
+
+## Selecting and visualizing a junction
+
+With the `stats` dataframe in hand, we can easily pick a junction of interest. For example, let's select one of the 2-category junctions with length ~1300 bp:
+
+```python
+stats.query("n_categories == 2 and 1200 < accessory_length < 1400")
+
+# n_isolates n_non_empty n_categories accessory_length ...
+# edge
+# 11809679528571820295_r__14906387308163561070_r 15 14 2 1242 ...
+# 13894307413921282410_r__14205544068867539089_f 15 1 2 1324 ...
+# 10486523597117694808_f__6531151666869853507_r 15 1 2 1324 ...
+# 4535080423279022649_f__4857718550591370421_r 15 1 2 1324 ...
+# 5751814192644177414_r__8949119045531796691_f 15 2 2 1324 ...
+...
+# 15114786226276103752_r__432022604910877054_r 15 1 2 1324 ...
+# 10485686697184953244_r__1548999589339136461_f 15 1 2 1324 ...
+# 1532495113773479365_r__1534747068225797391_f 15 1 2 1324 ...
+# 1534747068225797391_f__7253571478449116197_r 15 1 2 1324 ...
+```
+
+There are several such junctions, all with a characteristic accessory length of 1324 bp. This is suggestive of the same mobile element being integrated in several locations of the genome.
+
+A linear schematic of the junction makes its structure visible at a glance. The helper `pp.plots.linear_junction_plot` draws one row per isolate and one horizontal bar per block, with width equal to the block's consensus length. Junctions are co-oriented to the canonical edge direction so the flanking core blocks line up across rows.
+
+```python
+edge = "10485686697184953244_r__1548999589339136461_f"
+
+fig, ax = plt.subplots(figsize=(10, 5))
+pp.plots.linear_junction_plot(ax, junctions, edge)
+ax.set_title(edge, fontsize=9)
+plt.tight_layout()
+plt.show()
+```
+
+
+
+
+
+ **Customizing the linear junction plot**
+
+ Once the junction decomposition is done, the code to produce this linear representation for a junction is relatively simple. Feel free to modify it and customize it to your needs.
+
+ ```python
+ from collections import defaultdict
+ import matplotlib as mpl
+ import numpy as np
+ import pypangraph as pp
+
+ # load the graph and create the junctions object
+ graph = pp.Pangraph.from_json("staph.json.gz")
+ junctions = pp.junctions.BackboneJunctions(graph, L_thr=500)
+
+ # select an edge to plot
+ edge = "10485686697184953244_r__1548999589339136461_f"
+
+ # assign a unique random color to each block
+ cmap = mpl.colormaps["rainbow"]
+ colors = defaultdict(lambda: cmap(np.random.rand()))
+
+ # consensus length for each block
+ block_len = graph.to_blockstats_df()["len"].to_dict()
+
+ # dictionary isolate -> junction
+ Js = junctions[edge]
+
+ fig, ax = plt.subplots(figsize=(10, 5))
+
+ # cycle through each isolate that has the junction
+ for row, (iso, J) in enumerate(Js.items()):
+ J = J.to_canonical() # align junction to canonical orientation
+ x = 0
+ for ob in J.oriented_blocks():
+ # consensus length and color of each block
+ length = block_len[str(ob.id)]
+ color = colors[str(ob.id)]
+ # each block is a horizontal bar
+ ax.barh(
+ row,
+ length,
+ left=x,
+ height=0.8,
+ color=color,
+ edgecolor="k",
+ linewidth=0.4,
+ )
+ x += length
+
+ isolates = list(Js.keys())
+ ax.set_yticks(range(len(isolates)))
+ ax.set_yticklabels(isolates, fontsize=8)
+ ax.set_xlabel("position along junction (bp)")
+ plt.tight_layout()
+ plt.show()
+ ```
+
+
+
+
+ **Displaying more complex junctions**
+
+ What happens if we explore the pattern of more complex junctions? Try to select junctions that have a large number of unique path categories and visualize their linear structure. Here is an example for core-edge `13256234721607664913_r__7427484406751306657_f`.
+
+ 
+
+
+
+This pattern is suggestive of an element being integrated in this specific location in a single genome. But which element? To answer this question it is useful to be able to access the **location** of these blocks on the genome, to connect it with sequence annotations, and the **sequence** itself, for specific homology search or further downstream processing. This will be the topic of the next tutorial section.
diff --git a/docs/docs/pypangraph/t08-interesting-junction.md b/docs/docs/pypangraph/t08-interesting-junction.md
new file mode 100644
index 00000000..21f63b1f
--- /dev/null
+++ b/docs/docs/pypangraph/t08-interesting-junction.md
@@ -0,0 +1,168 @@
+---
+sidebar_position: 10
+---
+
+# Investigating junctions further: positions and sequences
+
+Using summary statistics, in the [previous part](t07-junction-stats.md) of the tutorial we have singled out a junction we want to investigate further. Here we will showcase the use of two `BackboneJunctions` methods for this:
+
+- `positions()`: indicates the location of each occurrence of the junction on the genomes that carry it. Useful e.g. for cross-referencing with annotation files.
+- `sequences()`: returns the actual DNA sequence spanning the junction on every genome, as `Bio.SeqRecord` objects ready to be written to FASTA for further analysis (multiple sequence alignment, secondary pangraph construction, BLAST searches, ...).
+
+
+## The example: a candidate IS insertion
+
+As a running example we will use the same junction introduced [in the previous section](t07-junction-stats.md): the edge `10485686697184953244_r__1548999589339136461_f`. Looking up its row in the statistics dataframe:
+
+```python
+import pypangraph as pp
+
+graph = pp.Pangraph.from_json("staph.json.gz")
+junctions = pp.junctions.BackboneJunctions(graph, L_thr=500)
+stats = junctions.stats()
+
+edge = "10485686697184953244_r__1548999589339136461_f"
+print(stats.loc[edge])
+# n_isolates 15
+# n_non_empty 1
+# n_categories 2
+# accessory_length 1324
+# ...
+```
+
+We found multiple junctions with this specific accessory length and low number of non-empty paths, indicative of a putative recent _Insertion Sequence_ (IS) insertion.
+
+## Genomic coordinates of a junction
+
+The `positions()` method returns a `pandas.DataFrame` indexed by `(edge, isolate)`, with the genomic coordinates of the two flanking core blocks plus a strand flag. We slice on our edge of interest with `.loc`:
+
+```python
+all_positions = junctions.positions()
+junction_pos = all_positions.loc[edge]
+print(junction_pos)
+# left_start left_end right_start right_end strand
+# iso
+# NZ_CP132362.1 2509902 2511822 2511822 2516373 False
+# NZ_LR822061.1 1288340 1290278 1290278 1294829 False
+# NZ_CP077852.1 1661385 1665936 1665936 1667856 True
+# ...
+# NZ_CP169947.1 880245 884796 884796 886716 True
+# NZ_CP022905.1 872983 877542 878866 880786 True
+# NZ_CP035791.1 1215039 1216959 1216959 1221510 False
+```
+
+The columns are:
+
+- `left_start`, `left_end`: genomic coordinates of the left flanking core block on the genome.
+- `right_start`, `right_end`: genomic coordinates of the right flanking core block on the genome.
+- `strand`: `True` if the junction appears in canonical edge orientation on this genome, `False` if reverse-complemented.
+
+Note that the **left** / **right** labels follow each genome's own path order: the **left block** is the **first core block** of the junction encountered when walking the genome, as illustrated in the scheme below. On isolates where the junction path is inverted (`strand = False`), this means the left/right roles are swapped relative to the canonical edge direction.
+
+The coordinates on the genome should therefore always satisfy `left_start` < `left_end` < `right_start` < `right_end`, irrespective of orientation. In circular genomes there can be exceptions to this pattern, as described below.
+
+
+
+:::info circular genomes
+
+In circular genomes, a junction path might wrap around the genome origin. In such cases the left core block might be found around the end of the sequence, and the right core block at the beginning. As such, in these genomes the order of positions described will not be respected.
+
+:::
+
+The difference between the `left_end` and `right_start` coordinates gives the precise size of the accessory region. This is zero when no accessory blocks are present.
+
+```python
+junction_pos.eval("right_start - left_end")
+# iso
+# NZ_CP132362.1 0
+# NZ_LR822061.1 0
+# NZ_CP077852.1 0
+# ...
+# NZ_CP169947.1 0
+# NZ_CP022905.1 1324 <-- IS insertion
+# NZ_CP035791.1 0
+```
+
+As we saw in the linear representation of the junction in the previous tutorial section, the insertion is only present in genome `NZ_CP022905.1` between coordinates:
+
+```python
+junction_pos.loc["NZ_CP022905.1", ["left_end", "right_start"]]
+# left_end 877542
+# right_start 878866
+```
+
+Inspecting [this location in NCBI's GenBank viewer](https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP022905.1?report=graph&from=877542&to=878866#) reveals indeed the presence of an [IS256 element](https://doi.org/10.1099/00221287-133-11-3031), commonly found in _S. aureus_, and likely to have recently inserted and given rise to the junction.
+
+
+
+The location of an insertion is often not random. If inserted in a coding region, an IS element can have significant fitness effects by disrupting genes. We can easily check if this is the case for our insertion by observing whether the junction breakpoint falls in a coding region in genomes where the IS element is absent. For example, consulting the `pos` table above, we can check [genome `NZ_LR822061.1`](https://www.ncbi.nlm.nih.gov/nuccore/NZ_LR822061.1?report=graph&from=1289278&to=1291278) at coordinate `1290278`.
+
+
+
+The IS element insertion truncates the _clfA_ gene, which encodes an immunogenic surface protein involved in adhesion, a [known virulence factor](https://doi.org/10.1038/nrmicro3161).
+
+## Extracting the junction sequences
+
+For further downstream analysis you often want the nucleotide sequence of a junction. The `sequences()` method returns one `Bio.SeqRecord` per isolate, spanning **left flank + accessory center + right flank**, all co-oriented to the canonical edge direction so that no further inversion is necessary.
+
+```python
+records = junctions.sequences(edge)
+for r in records:
+ print(f"{r.id}: {r.seq[:10]}...{r.seq[-10:]} | len = {len(r.seq)} bp")
+# NZ_CP132362.1: ATTTGTAGCC...ACTCAGACAG | len = 6471 bp
+# NZ_LR822061.1: ATTTGTAGCC...ACTCAGACAG | len = 6489 bp
+# NZ_CP077852.1: ATTTGTAGCC...ACTCAGACAG | len = 6471 bp
+# ...
+# NZ_CP169947.1: ATTTGTAGCC...ACTCAGACAG | len = 6471 bp
+# NZ_CP022905.1: ATTTGTAGCC...ACTCAGACAG | len = 7803 bp
+# NZ_CP035791.1: ATTTGTAGCC...ACTCAGACAG | len = 6471 bp
+```
+
+The IS carrier (`NZ_CP022905.1`) stands out as ~1.3 kb longer than the rest. Each record has:
+
+- `id`: the isolate name.
+- `description`: the canonical edge string.
+- `seq`: the DNA from the start of the left flank to the end of the right flank.
+
+These can be conveniently exported in a fasta file using [biopython](https://biopython.org/) for further downstream processing.
+
+```python
+from Bio import SeqIO
+
+SeqIO.write(records, "junction.fa", "fasta")
+```
+
+For example a multiple sequence alignment with [MAFFT](https://mafft.cbrc.jp/alignment/software/):
+
+```bash
+mafft junction.fa > junction_aligned.fa
+```
+
+
+
+Or a second, junction-specific pangraph. This is usually useful in more complex junctions, to refine the junction structure.
+
+```bash
+pangraph build -l 100 -s 20 junction.fa -o junction.json
+pangraph export gfa junction.json -o junction.gfa
+```
+
+
+
+## More junctions to explore
+
+As a further exercise you can try to explore some more interesting junctions yourself. Here are some suggestions:
+
+- look into `13654622636097204630_f__5922465870918455593_f`.
+ - Can you locate the **prophage insertion**?
+ - Find its coordinates in the genome and visualize it on NCBI. Inspecting the annotations of the accessory region should confirm that it is likely a prophage integration. Where did it get integrated? _Note_: it wraps around the start of the sequence record.
+ - Can you find the gene immediately upstream of the prophage insertion and the first gene of the accessory region? Does this suggest a **mechanism of integration**?
+- edge `10486523597117694808_f__6531151666869853507_r` is another simple example of a junction originated by an IS element insertion.
+ - Which gene is inactivated by the insertion? Take one of the isolates without the insertion (e.g. `NZ_CP132362.1`) and check which gene is present at the location where isolate `NZ_LR822061.1` shows an insertion. You should find that the insertion likely inactivated the [_staphylocoagulase_](https://www.uniprot.org/uniprotkb/P17855/entry), a known virulence factor.
+ - This is just one example of gene inactivation; in the [junction overview plot](t07-junction-stats.md#visualizing-the-junction-landscape) we noted that several junctions originated from an IS movement. What is the general pattern? Can you automate this search and check how often a recent IS element insertion interrupted a gene? Is this more or less frequent than what is expected by chance? What could be the selective forces driving this imbalance?
+- look into `17042526223432838337_f__8287974428665837959_r`
+ - Can you spot a **duplication** in some sequences? Which genes are duplicated?
+- edge `13256234721607664913_r__7427484406751306657_f` encompasses a **highly-variable region** in these genomes.
+ - draw a linear representation to compare the variation across different genomes.
+ - focus on one genome, e.g. `NZ_CP022905.1`, and look at the annotations. Can you find mobile genetic elements? Can you find resistance determinants?
+ - the presence of the _mecA_ gene and recombinases suggests that this is a [SCCmec cassette](https://en.wikipedia.org/wiki/SCCmec) ([IWG-SCC 2009](https://doi.org/10.1128/AAC.00579-09)), a known mobile genetic element of _Staphylococcus_ bacteria that confers methicillin resistance.
diff --git a/docs/docs/pypangraph/t09-untangling-with-junctions.md b/docs/docs/pypangraph/t09-untangling-with-junctions.md
new file mode 100644
index 00000000..6b109149
--- /dev/null
+++ b/docs/docs/pypangraph/t09-untangling-with-junctions.md
@@ -0,0 +1,85 @@
+---
+sidebar_position: 11
+---
+
+# Bonus: visually untangling graph complexity using junctions
+
+When a whole pangenome graph is exported to GFA and opened in [Bandage](https://rrwick.github.io/Bandage/), the result is usually a **tangle**. As we saw in the [build tutorial](../tutorial/t01-building-pangraph.md), the same accessory or duplicated block can occur in many different genomic contexts, so a single segment ends up linked to many distant parts of the graph. These long-range links are what make the layout look like a hairball. This can be mitigated by filtering out duplicated blocks or even all accessory blocks, but at the cost of losing visualization of the accessory diversity.
+
+Junctions (introduced in [the junctions tutorial](t06-junctions-intro.md)) offer a better solution: we can keep the accessory blocks but **disentangle them by context**. The idea is to *paralog-split* each block according to the junction it sits in, so a block shared across several junctions becomes one segment per junction instead of one segment wired to all of them. The core blocks stay shared and act as anchors, and the accessory diversity is laid out as clean bubbles strung along the core-synteny backbone.
+
+This is made more concrete in the example below:
+
+
+
+Here we consider three core blocks (`X`, `Y` and `Z`) and two accessory blocks (`a` and `b`). Block `a` is found in two junction contexts (`[X|Y]` and `[Y|Z]`). As a consequence, paths need to traverse it twice, before and after core block `Y`. This generates tangles in the representation.
+
+To circumvent this, we can **make use of the junction information** to de-duplicate block `a` in two occurrences (`a1` and `a2`) that can clearly be distinguished by context. As a result of this operation, the tangle is resolved and all core junctions are clearly separated.
+
+
+## The tangled starting point
+
+We use the same `staph.json.gz` graph as the rest of the junction tutorials (15 _Staphylococcus aureus_ chromosomes). Exporting the whole graph to GFA with the CLI gives the usual tangled representation, even if we avoid exporting duplicated blocks:
+
+```bash
+pangraph export gfa --no-duplicated staph.json.gz -o staph_full.gfa
+```
+
+Opened in Bandage, this graph (494 segments, 748 links) is hard to read: accessory blocks that recur in different contexts create crossings all over the layout.
+
+
+
+## Untangling with junction context
+
+
+`pypangraph` can instead build a **junction-context resolved GFA** directly from a `BackboneJunctions` object. We load the graph, identify the junctions, and export:
+
+```python
+import pypangraph as pp
+from pypangraph.export import junction_context_gfa
+
+graph = pp.Pangraph.from_json("staph.json.gz")
+junctions = pp.junctions.BackboneJunctions(graph, L_thr=500)
+
+# paralog-split blocks by junction context (scaffold options explained below)
+gfa, prefix_map = junction_context_gfa(junctions, scaffold="consensus")
+gfa.write("staph_untangled.gfa")
+```
+
+This emits each block once *per junction context*: a block that occurs in several junctions becomes several segments, each prefixed by a junction tag (`J{n}__{block_id}`), while the core blocks -- the shared anchors -- keep their plain id and are emitted only once. Because blocks are de-duplicated by context, the export actually has **more** segments than the whole-graph GFA (1312 segments, 1742 links here) -- but each one is local to a single junction, so the long-range crossings disappear and the graph reads as a clean chain of bubbles.
+
+
+
+The second return value, `prefix_map`, maps each junction tag `J{n}` back to the core edge it came from (e.g. for writing a companion TSV). The export also records per-segment depth as a `DP:f:` tag, so Bandage can colour segments by how many isolates traverse them. The `scaffold` argument is covered in the next section.
+
+## A second source of tangle: synteny changes
+
+Paralog-splitting removes the crossings caused by accessory blocks, but there is a second -- usually minor -- source of tangle: **changes in core-genome synteny**. When the order or orientation of the core blocks themselves differs between genomes, the difference shows up as a "pinching point" in the visual representation (this is illustrated in the [marginalization tutorial](../tutorial/t04-graph-projection.md#a-look-at-the-marginalized-pangraph)).
+
+The `scaffold` argument of `junction_context_gfa` decides how this is handled, by choosing which set of core-synteny edges defines the backbone:
+
+- `"consensus"` (the default used above) keeps, for every pair of consecutive core blocks, the orientation seen in a **strict majority** of genomes. Minority rearrangements are dropped, so the synteny pinching points are smoothed out.
+- a **reference isolate name**, e.g. `scaffold="NZ_CP162433.1"`, uses that single genome's own core synteny as the backbone (useful when one assembly is your reference of record).
+- `"all"` keeps **every** junction -- the union of all synteny variants seen across the dataset. This reintroduces the alternative orderings that the consensus scaffold drops, so it is slightly more tangled, but loses no synteny variation.
+
+
+ keeping synteny changes with scaffold="all"
+
+ ```python
+ gfa_all, _ = junction_context_gfa(junctions, scaffold="all")
+ gfa_all.write("staph_all.gfa")
+ ```
+
+ For this graph `"all"` keeps 151 junctions versus the 143 of the consensus scaffold -- the extra ones are the minority synteny rearrangements.
+
+ 
+
+ Here the core blocks are coloured by [minimal synteny unit](t04-core-synteny.md#core-genome-synteny): a change of colour along the backbone marks a **synteny breakpoint** -- one of the rearrangements that the consensus scaffold smooths away.
+
+
+
+## Bonus: a much larger graph
+
+The approach scales also to **large datasets**. Below is the junction-context export (consensus scaffold) of a much larger pangraph built from **222 _E. coli_ ST131 sequences** previously analyzed in [our paper](https://doi.org/10.1093/molbev/msae272). The junction decomposition resolves the structure into an interpretable overall genomic context, whereas a standard whole-graph GFA would essentially be unreadable.
+
+
diff --git a/docs/docs/tutorial/t01-building-pangraph.md b/docs/docs/tutorial/t01-building-pangraph.md
index 97784291..9cde3691 100644
--- a/docs/docs/tutorial/t01-building-pangraph.md
+++ b/docs/docs/tutorial/t01-building-pangraph.md
@@ -149,7 +149,7 @@ pangraph export gfa \
This will create a `graph.gfa` file, which can be visualized using [Bandage](https://rrwick.github.io/Bandage/).
-The option `--no-duplications` causes the export function to avoid including duplicated blocks in the graph representation (they are instead exported as isolated blocks). In our experience this results in a less "tangled" visual representation. Below is how the Bandage visualization of this example pangraph looks like. Blocks are colored by frequency, with common blocks (appearing in many different chromosomes) in red and rare blocks (appearing in only a few chromosomes) in black.
+The option `--no-duplicated` causes the export function to avoid including duplicated blocks in the graph representation (they are instead exported as isolated blocks). In our experience this results in a less "tangled" visual representation. Below is how the Bandage visualization of this example pangraph looks like. Blocks are colored by frequency, with common blocks (appearing in many different chromosomes) in red and rare blocks (appearing in only a few chromosomes) in black.
@@ -166,4 +166,8 @@ pangraph export gfa \
The resulting graph is much simpler. The remaining crossings are due to changes in core-genome synteny. Each change in order of core blocks results in a crossing in the graph, as will be discussed in [a later tutorial section](../pypangraph/t04-core-synteny.md).
-
\ No newline at end of file
+
+
+An even better option, based on our [python library](../pypangraph/about-pypangraph.md), is discussed in [a later part of the tutorial](../pypangraph/t09-untangling-with-junctions.md): it keeps the accessory blocks but disentangles them by their core-genome context, so the graph reads clearly without sacrificing any of the accessory diversity.
+
+
\ No newline at end of file
diff --git a/docs/docs/tutorial/t04-graph-projection.md b/docs/docs/tutorial/t04-graph-projection.md
index 0092b272..0d6ad4d7 100644
--- a/docs/docs/tutorial/t04-graph-projection.md
+++ b/docs/docs/tutorial/t04-graph-projection.md
@@ -35,7 +35,9 @@ The output file can be visualized using [Bandage](https://rrwick.github.io/Banda
-Colors indicate the number of times a block occurs. Blocks that appear in red are core blocks that are found in every chromosome, while black blocks are only present in a few strains. We used the `--no-duplication` flag in the export command, which excludes duplicated blocks from the exported graph. This simplifies the visualization, which would otherwise be highly "tangled-up" by these duplications.
+Colors indicate the number of times a block occurs. Blocks that appear in red are core blocks that are found in every chromosome, while black blocks are only present in a few strains. We used the `--no-duplicated` flag in the export command, which excludes duplicated blocks from the exported graph. This simplifies the visualization, which would otherwise be highly "tangled-up" by these duplications.
+
+An even better option, based on our python library, is discussed in [a later part of the tutorial](../pypangraph/t09-untangling-with-junctions.md): rather than dropping the duplicated blocks, it disentangles them by their junction context, giving a clear graph while keeping the full accessory diversity.
## Marginalize the graph on a subset of strains
@@ -77,7 +79,7 @@ As done for the main pangraph, we can export the marginalized pangraph in `gfa`
```bash
pangraph export gfa \
- --no-duplications \
+ --no-duplicated \
--minimum-length 150 \
klebs_marginal_pangraph.json \
-o klebs_marginal_pangraph.gfa
diff --git a/packages/pypangraph/CHANGELOG.md b/packages/pypangraph/CHANGELOG.md
index 9796cbed..f6d8764c 100644
--- a/packages/pypangraph/CHANGELOG.md
+++ b/packages/pypangraph/CHANGELOG.md
@@ -1,8 +1,17 @@
+## 1.1.0
+
+Adds **backbone-junction analysis**.
+
+New `pypangraph.junctions` sub-package built around `BackboneJunctions`, which splits each path at core-block boundaries and exposes:
+
+- `stats()`: per-edge summary DataFrame (number of categories, accessory length, occupied isolates, ...).
+- `positions()`: genomic coordinates of the flanking core blocks per (edge, isolate).
+- `sequences(edge)`: co-oriented `Bio.SeqRecord`s spanning the junction, ready to write to FASTA.
+
## 1.0.1
Make pypangraph compatible with Python 3.9, by removing `match` syntax.
-
## 1.0.0
Initial release for pangraph 1.0.0.
diff --git a/packages/pypangraph/pypangraph/__init__.py b/packages/pypangraph/pypangraph/__init__.py
index 481ed21c..d0217214 100644
--- a/packages/pypangraph/pypangraph/__init__.py
+++ b/packages/pypangraph/pypangraph/__init__.py
@@ -1,6 +1,21 @@
# pangraph class
-from .class_graph import Pangraph
+from .class_graph import Pangraph, PangraphLoadError
-# sinteny analysis
-from .msu import minimal_synteny_units
+# synteny analysis
+from .minimal_synteny_units import minimal_synteny_units
from .plots import dotplot
+
+# junction analysis
+from . import junctions
+
+# export
+from . import export
+
+__all__ = [
+ "Pangraph",
+ "PangraphLoadError",
+ "minimal_synteny_units",
+ "dotplot",
+ "junctions",
+ "export",
+]
diff --git a/packages/pypangraph/pypangraph/class_alignments.py b/packages/pypangraph/pypangraph/class_alignments.py
index 9c70983f..8a061476 100644
--- a/packages/pypangraph/pypangraph/class_alignments.py
+++ b/packages/pypangraph/pypangraph/class_alignments.py
@@ -5,7 +5,6 @@
from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq
-import numpy as np
@dataclass
@@ -62,11 +61,11 @@ def apply_edits(self, cons: str) -> str:
for D in self.dels:
for idx in range(D.length):
sq[D.pos + idx] = ""
- for I in self.inss:
- if I.pos > 0:
- sq[I.pos - 1] += I.seq
- elif I.pos == 0:
- sq[0] = I.seq + sq[0]
+ for ins in self.inss:
+ if ins.pos > 0:
+ sq[ins.pos - 1] += ins.seq
+ elif ins.pos == 0:
+ sq[0] = ins.seq + sq[0]
return "".join(sq)
def aligned_seq(self, cons: str) -> str:
diff --git a/packages/pypangraph/pypangraph/class_block.py b/packages/pypangraph/pypangraph/class_block.py
index 10bf8328..fdb78f2b 100644
--- a/packages/pypangraph/pypangraph/class_block.py
+++ b/packages/pypangraph/pypangraph/class_block.py
@@ -11,13 +11,16 @@ class Block:
class for details.
"""
- def __init__(self, block_id: int, alignment: pga.Alignment):
+ def __init__(self, block_id: str, alignment: pga.Alignment):
self.id = block_id
self.alignment = alignment
@staticmethod
def from_dict(block: dict) -> "Block":
- block_id = block["id"]
+ # Block ids are u64 hashes; keep them as strings internally so they never get
+ # coerced to float when used as a pandas index/column (which silently corrupts
+ # values above 2**53).
+ block_id = str(block["id"])
alignment = pga.Alignment.from_dict(block)
return Block(block_id, alignment)
@@ -39,11 +42,11 @@ def consensus(self) -> str:
"""Returns the consensus sequence of the block"""
return self.alignment.consensus
- def to_sequences(self) -> dict[int, str]:
+ def to_sequences(self) -> dict[str, str]:
"""Returns a dictionary node_id -> sequence for the block"""
return self.alignment.generate_sequences()
- def to_alignment(self) -> dict[int, str]:
+ def to_alignment(self) -> dict[str, str]:
"""Returns a dictionary node_id -> aligned sequence for the block.
The aligned sequence does not include insertions."""
return self.alignment.generate_alignment()
@@ -63,6 +66,6 @@ class BlockCollection(IndexedCollection):
"""
def __init__(self, pan_blocks):
- ids = [block["id"] for block in pan_blocks.values()]
+ ids = [str(block["id"]) for block in pan_blocks.values()]
items = [Block.from_dict(block) for block in pan_blocks.values()]
IndexedCollection.__init__(self, ids, items)
diff --git a/packages/pypangraph/pypangraph/class_graph.py b/packages/pypangraph/pypangraph/class_graph.py
index 17b54fc6..d367097a 100644
--- a/packages/pypangraph/pypangraph/class_graph.py
+++ b/packages/pypangraph/pypangraph/class_graph.py
@@ -15,6 +15,10 @@
from .pangraph_schema import schema
+class PangraphLoadError(ValueError):
+ """Raised when a Pangraph JSON file cannot be loaded or validated."""
+
+
class Pangraph:
"""Wrapper class to load and interact with the output of the Pangraph pipeline.
The class has three main attributes:
@@ -53,22 +57,29 @@ def from_json(filename):
is_json = str(filename).endswith(".json")
is_gzjson = str(filename).endswith(".json.gz")
if not (is_json or is_gzjson):
- raise Exception(
+ raise PangraphLoadError(
f"the input file {filename} should be in .json or .json.gz format"
)
- if is_gzjson:
- with gzip.open(filename, "rt") as f:
- pan_json = json.load(f)
- else:
- with open(filename, "r") as f:
- pan_json = json.load(f)
+ try:
+ if is_gzjson:
+ with gzip.open(filename, "rt") as f:
+ pan_json = json.load(f)
+ else:
+ with open(filename, "r") as f:
+ pan_json = json.load(f)
+ except (OSError, gzip.BadGzipFile, json.JSONDecodeError) as ex:
+ raise PangraphLoadError(
+ f"failed to load pangraph from {filename}: {ex}"
+ ) from ex
try:
graph = {"pangraph": pan_json}
jsonschema.validate(instance=graph, schema=schema)
except jsonschema.exceptions.ValidationError as ex:
- print(ex)
+ raise PangraphLoadError(
+ f"invalid pangraph JSON in {filename}: {ex.message}"
+ ) from ex
pan = Pangraph(pan_json)
return pan
@@ -92,8 +103,8 @@ def to_blockstats_df(self):
- count: n. times the block occurs
- n. strains: number of strains in which the block is observed
- duplicated: whether the block is duplicated in at least one strain
- - len: average block length from pangraph.
- - core: whether a gene occurrs exactly once per strain
+ - len: block consensus length in basepairs.
+ - core: whether a block occurs exactly once per strain
"""
df = self.nodes.to_blockstats_df()
df["len"] = [len(self.blocks[bid]) for bid in df.index]
@@ -159,9 +170,8 @@ def core_genome_alignment(self, guide_strain=None):
# get order of core blocks in guide strain
if guide_strain is None:
guide_strain = self.strains()[0]
- assert guide_strain in strains, (
- f"Guide strain {guide_strain} not found in the dataset"
- )
+ if guide_strain not in strains:
+ raise ValueError(f"Guide strain {guide_strain} not found in the dataset")
guide_path = self.paths[guide_strain]
# core block ids and strandedness
@@ -177,9 +187,11 @@ def core_genome_alignment(self, guide_strain=None):
for bid, guide_strand in core_blocks:
# get block alignment
aln_dict = self.blocks[bid].to_alignment()
- assert len(aln_dict) == len(strains), (
- f"error: unexpected number of strains {bid}"
- )
+ if len(aln_dict) != len(strains):
+ raise ValueError(
+ f"core block {bid} has {len(aln_dict)} alignment rows but the graph "
+ f"has {len(strains)} strains"
+ )
# append alignment to the final alignment for each strain
aln_strains = []
@@ -193,9 +205,10 @@ def core_genome_alignment(self, guide_strain=None):
aln_strains.append(strain)
# sanity check: core-blocks are present once per strain
- assert set(strains) == set(aln_strains), (
- f"error: strain missing in block {bid}: {set(strains)} != {set(aln_strains)}"
- )
+ if set(strains) != set(aln_strains):
+ raise ValueError(
+ f"strain missing in core block {bid}: {set(strains)} != {set(aln_strains)}"
+ )
# convert to biopython alignment
records = []
diff --git a/packages/pypangraph/pypangraph/class_node.py b/packages/pypangraph/pypangraph/class_node.py
index 3043bb95..9b98a7f9 100644
--- a/packages/pypangraph/pypangraph/class_node.py
+++ b/packages/pypangraph/pypangraph/class_node.py
@@ -19,7 +19,7 @@ class Nodes:
"""Dataset with pangraph nodes information:
- id (str): node id. Index of the dataframe
- block_id (str): block id
- - path_id (str): path id
+ - path_id (int): path id
- strand (bool): strandedness of the node
- position (tuple): position of the node in the genome (start, end)
"""
@@ -28,7 +28,7 @@ def __init__(self, nodes_dict):
self.df = pd.DataFrame.from_dict(
{
node_id: {
- "block_id": node["block_id"],
+ "block_id": str(node["block_id"]),
"path_id": node["path_id"],
"strand": parse_strandedness(node["strand"]),
"start": node["position"][0],
@@ -78,11 +78,14 @@ def to_blockstats_df(self) -> pd.DataFrame:
df["core"] &= ~df["duplicated"]
return df
- def node_to_block(self, node_id: int) -> tuple[int, bool]:
- """Returns the node's block id and strandedness"""
+ def node_to_block(self, node_id: int) -> pd.Series:
+ """Returns the node's block id (str) and strandedness (bool) as a 2-element
+ pandas Series, unpackable as ``block_id, strand = node_to_block(...)``."""
return self.df.loc[str(node_id), ["block_id", "strand"]]
- def nodes_to_blocks(self, node_ids: list[int]) -> tuple[list[int], list[bool]]:
- """Returns the block id and strandedness of a list of nodes"""
+ def nodes_to_blocks(self, node_ids: list[int]) -> np.ndarray:
+ """Returns block ids and strandedness for a list of nodes as a (2, N) numpy
+ array, unpackable as ``block_ids, strands = nodes_to_blocks(...)`` (the first
+ row holds str block ids, the second the bool strands)."""
N = np.array(node_ids, dtype=str)
return self.df.loc[N, ["block_id", "strand"]].values.T
diff --git a/packages/pypangraph/pypangraph/export/__init__.py b/packages/pypangraph/pypangraph/export/__init__.py
new file mode 100644
index 00000000..ad2a2345
--- /dev/null
+++ b/packages/pypangraph/pypangraph/export/__init__.py
@@ -0,0 +1,4 @@
+from .gfa import GFA
+from .junctions_to_gfa import junction_context_gfa
+
+__all__ = ["GFA", "junction_context_gfa"]
diff --git a/packages/pypangraph/pypangraph/export/gfa.py b/packages/pypangraph/pypangraph/export/gfa.py
new file mode 100644
index 00000000..f837f595
--- /dev/null
+++ b/packages/pypangraph/pypangraph/export/gfa.py
@@ -0,0 +1,86 @@
+"""Minimal, graph-agnostic GFA1 representation.
+
+An in-memory GFA1 graph (segments + links, with optional per-segment depths)
+that knows how to serialize itself. Sequences are not stored: segment lengths are
+emitted as ``LN:i:`` tags and the sequence field is left as ``*``. The real bp
+lengths can optionally be rescaled at write time (e.g. ``log`` or ``1/l``) for
+visualization — see :meth:`GFA.write`.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Callable, Iterable
+
+
+def _orient(strand: bool) -> str:
+ """Map a strand boolean to a GFA orientation symbol (True -> '+', False -> '-')."""
+ return "+" if strand else "-"
+
+
+class GFA:
+ """An in-memory GFA1 graph: segments, links, and optional per-segment depths.
+
+ Build it directly or via a builder (e.g.
+ :func:`pypangraph.export.junction_context_gfa`), then call :meth:`write` to
+ serialize it to a GFA1 file.
+ """
+
+ def __init__(
+ self,
+ segments: dict[str, int],
+ links: Iterable[tuple[str, bool, str, bool]],
+ depths: dict[str, float] | None = None,
+ ) -> None:
+ """Store the graph.
+
+ Args:
+ segments: dict mapping segment name (str) -> length in bp (int).
+ links: iterable of (from_name, from_strand, to_name, to_strand) tuples,
+ where the strands are booleans (True -> '+', False -> '-'). A set
+ collapses duplicate links automatically.
+ depths: optional dict mapping segment name -> coverage depth, emitted
+ as a ``DP:f:`` tag (read by Bandage as node depth).
+ """
+ self.segments = segments
+ self.links = links
+ self.depths = depths or {}
+
+ def write(
+ self, filepath: str, length_transform: Callable[[int], float] | None = None
+ ) -> None:
+ """Write this graph to a minimal GFA1 file.
+
+ Emits the ``H`` header, one ``S`` line per segment (``LN:i:`` length and,
+ when present in :attr:`depths`, a ``DP:f:`` depth tag, sequence ``*``),
+ and one ``L`` line per link with a ``0M`` overlap.
+
+ Args:
+ filepath: Output path for the GFA1 file.
+ length_transform: optional ``Callable`` applied to each segment's real
+ bp length to rescale the emitted ``LN:i:`` value (useful for
+ visualization, since block lengths span orders of magnitude). The
+ result is rounded to the nearest integer and clamped to a minimum
+ of ``1``; ``None`` (default) emits the true length. The caller owns
+ the scaling: raw ``math.log(l)`` is tiny (~5-14) and raw ``1/l``
+ collapses everything to the ``1`` floor, so scale as needed, e.g.::
+
+ gfa.write("out.gfa", length_transform=lambda l: 100 * math.log(l))
+ gfa.write("out.gfa", length_transform=lambda l: 1e7 / l)
+ """
+ with open(filepath, "w") as f:
+ f.write("H\tVN:Z:1.0\n")
+ for name, length in self.segments.items():
+ ln = (
+ int(length)
+ if length_transform is None
+ else max(1, round(length_transform(length)))
+ )
+ line = f"S\t{name}\t*\tLN:i:{ln}"
+ if name in self.depths:
+ line += f"\tDP:f:{self.depths[name]}"
+ f.write(line + "\n")
+ for from_name, from_strand, to_name, to_strand in self.links:
+ f.write(
+ f"L\t{from_name}\t{_orient(from_strand)}"
+ f"\t{to_name}\t{_orient(to_strand)}\t0M\n"
+ )
diff --git a/packages/pypangraph/pypangraph/export/junctions_to_gfa.py b/packages/pypangraph/pypangraph/export/junctions_to_gfa.py
new file mode 100644
index 00000000..8042a5b8
--- /dev/null
+++ b/packages/pypangraph/pypangraph/export/junctions_to_gfa.py
@@ -0,0 +1,145 @@
+"""Junction-context GFA decomposition of a pangenome graph.
+
+Unlike the standard whole-graph GFA export, this disentangles the graph by
+**junction context**: blocks are paralog-split and de-duplicated per junction
+(the same block in two different junctions becomes two segments; shared core
+anchors stay single), and the topology is reduced to a chosen **core synteny**
+scaffold so the result is a clean, walkable graph instead of a tangle.
+
+Core flank blocks are emitted once globally (segment name = block id); accessory
+blocks get an enumerated per-junction prefix (``J{n}__{block_id}``) so blocks
+shared across junctions are emitted once per junction. Each isolate's junction is
+co-oriented to its canonical edge direction so the junction forms a clean bubble
+between its two core anchors.
+"""
+
+from __future__ import annotations
+
+from collections import Counter, defaultdict
+
+from ..minimal_synteny_units import core_paths
+from .gfa import GFA
+
+
+def _iso_core_edges(pan, L_thr: int) -> dict[str, frozenset[str]]:
+ """Map each isolate to the frozenset of canonical core-edge ids on its backbone.
+
+ ``core_paths`` purifies each genome to core blocks >= ``L_thr``; the canonical
+ ``Edge.to_str_id`` makes a reverse traversal yield identical ids. This works
+ unchanged for multi-replicon genomes: edges from several disjoint core cycles
+ just land in the same set.
+ """
+ return {
+ iso: frozenset(e.to_str_id() for e in w.edges())
+ for iso, w in core_paths(pan, L_thr).items()
+ }
+
+
+def _consensus_edge_set(iso_edges: dict[str, frozenset[str]]) -> set[str]:
+ """Core edges present in a strict majority of isolate backbones.
+
+ This is the single tunable consensus-policy point: change the predicate here
+ (e.g. weighted or threshold-based voting) without touching the rest of the
+ export.
+ """
+ n = len(iso_edges)
+ counts = Counter(e for edges in iso_edges.values() for e in edges)
+ return {e for e, c in counts.items() if c > n / 2}
+
+
+def _scaffold_edges(bj, scaffold: str) -> list[str]:
+ """Return the canonical edge ids to keep for a ``scaffold`` policy.
+
+ The result is sorted for deterministic, stable ``J{n}`` numbering. ``"all"``
+ keeps every junction (equivalently, the union of all isolate core-edge sets);
+ ``"consensus"`` keeps the per-edge majority; otherwise ``scaffold`` is an
+ isolate name and its own core-edge set is used. Edges absent from ``bj`` are
+ dropped.
+ """
+ available = set(bj.edges())
+ if scaffold == "all":
+ return sorted(available)
+
+ iso_edges = _iso_core_edges(bj.pan, bj.L_thr)
+ if scaffold == "consensus":
+ chosen = _consensus_edge_set(iso_edges)
+ elif scaffold in iso_edges:
+ chosen = iso_edges[scaffold]
+ else:
+ raise ValueError(
+ f"unknown scaffold isolate {scaffold!r}; "
+ f"expected 'consensus', 'all', or one of {sorted(iso_edges)}"
+ )
+ return sorted(chosen & available)
+
+
+def junction_context_gfa(bj, scaffold: str = "consensus") -> tuple[GFA, dict[str, str]]:
+ """Build a junction-context GFA decomposition of a graph.
+
+ Args:
+ bj: A ``BackboneJunctions`` object.
+ scaffold: Which core-synteny scaffold defines the kept junctions.
+ ``"consensus"`` (default) keeps the per-edge majority backbone; an
+ isolate name uses that genome's own core edges; ``"all"`` keeps every
+ junction (tangled escape hatch).
+
+ Returns:
+ A tuple ``(gfa, prefix_map)``:
+ - gfa: a :class:`~pypangraph.export.GFA` whose segments (name -> length in
+ bp), links and per-segment depths (coverage: the graph-wide occurrence
+ count for core anchors, the per-junction traversal count for accessory
+ copies) describe the decomposition. Serialize it with ``gfa.write(path)``.
+ - prefix_map: dict[str, str] ``"J{n}"`` -> canonical edge string id;
+ decomposition provenance the caller may persist (e.g. as a TSV).
+ """
+ bdf = bj.block_stats
+ kept = _scaffold_edges(bj, scaffold)
+
+ # dictionaries of block id -> length and count
+ len_map = bdf["len"].astype(int).to_dict()
+ count_map = bdf["count"].astype(int).to_dict()
+
+ segments = {}
+ links = set()
+ prefix_map = {}
+ acc_depth = defaultdict(int) # accessory segment name -> number of traversals
+
+ # Cycle through all kept junctions
+ for n, edge_str in enumerate(kept):
+ # junction unique prefix
+ prefix = f"J{n}"
+ prefix_map[prefix] = edge_str
+
+ # Cycle throught the value of the dictionary {iso -> junction path}
+ for junction in bj[edge_str].values():
+ # Co-orient to the canonical edge direction so all instances align.
+ jc = junction.to_canonical()
+ # Capture the two flanking core blocks
+ core_ids = {jc.left.id, jc.right.id}
+
+ # Build a list of (segment name, orientation)
+ segment_walk = []
+ for ob in jc.oriented_blocks():
+ if ob.id in core_ids:
+ # segment name for core blocks
+ name = ob.id
+ else:
+ # segment name for accessory blocks. Also increase the counter.
+ name = f"{prefix}__{ob.id}"
+ acc_depth[name] += 1
+
+ # save segment with its length, and add it to the walk
+ segments[name] = len_map[ob.id]
+ segment_walk.append((name, ob.strand))
+
+ # Link consecutive segments, carrying each block's strand into the GFA link.
+ for (n1, s1), (n2, s2) in zip(segment_walk, segment_walk[1:]):
+ links.add((n1, s1, n2, s2))
+
+ # Depth of anchor (core) blocks is given by the graph count (equal to number of isolates)
+ # while for accessory blocks is the number of occurrences in the specific junction.
+ depths = {
+ name: acc_depth[name] if name in acc_depth else count_map[name]
+ for name in segments
+ }
+ return GFA(segments, links, depths), prefix_map
diff --git a/packages/pypangraph/pypangraph/indexed_collection.py b/packages/pypangraph/pypangraph/indexed_collection.py
index 74a85990..b7f7915c 100644
--- a/packages/pypangraph/pypangraph/indexed_collection.py
+++ b/packages/pypangraph/pypangraph/indexed_collection.py
@@ -12,6 +12,10 @@ class IndexedCollection:
- queried for the list of ids with the `keys` function (like a dictionary).
The object's `len` is the lentgth of the list.
+
+ Ids are stored as strings, but lookups (`__getitem__`, `__contains__`) coerce the
+ key to str first, so a collection keyed by numeric ids (e.g. blocks) can be indexed
+ with either an int or a str.
"""
def __init__(self, ids, items):
@@ -20,8 +24,8 @@ def __init__(self, ids, items):
self.id_to_pos = {id_: n for n, id_ in enumerate(ids)}
def __contains__(self, id_):
- """Returns whether the id is in the collection"""
- return id_ in self.id_to_pos
+ """Returns whether the id is in the collection (accepts int or str)"""
+ return str(id_) in self.id_to_pos
def __iter__(self):
"""Returns an iterator over the paths"""
@@ -35,14 +39,16 @@ def __len__(self):
return len(self.list)
def __getitem__(self, id_):
- """Returns the item corresponding to the id"""
+ """Returns the item corresponding to the id (accepts int or str)"""
try:
- pos = self.id_to_pos[id_]
+ pos = self.id_to_pos[str(id_)]
return self.list[pos]
except KeyError:
+ # `from None`: the original KeyError only carries the key, which is
+ # already in this message, so suppress it for a clean traceback.
raise KeyError(
f"Id {id_} not found in collection {self.__class__.__name__}"
- )
+ ) from None
def keys(self):
"""Returns the list of ids (like a dictionary)"""
diff --git a/packages/pypangraph/pypangraph/junctions/__init__.py b/packages/pypangraph/pypangraph/junctions/__init__.py
new file mode 100644
index 00000000..1ad174e0
--- /dev/null
+++ b/packages/pypangraph/pypangraph/junctions/__init__.py
@@ -0,0 +1,3 @@
+from .backbone import BackboneJunctions
+
+__all__ = ["BackboneJunctions"]
diff --git a/packages/pypangraph/pypangraph/junctions/backbone.py b/packages/pypangraph/pypangraph/junctions/backbone.py
new file mode 100644
index 00000000..89fca86a
--- /dev/null
+++ b/packages/pypangraph/pypangraph/junctions/backbone.py
@@ -0,0 +1,144 @@
+import pandas as pd
+from Bio.SeqRecord import SeqRecord
+
+from ..topology_utils import Walk
+from .junction import JunctionNode, Junction, path_junction_split
+from .positions import junction_positions
+from .sequences import junction_sequences
+from .stats import junction_stats
+
+
+class BackboneJunctions:
+ """Backbone junction analysis for a pangenome graph.
+
+ Splits each path into junctions at backbone (core + above-threshold length)
+ block boundaries. Results are lazily computed and cached.
+
+ Args:
+ pan: A Pangraph object.
+ L_thr: Minimum block length to be considered backbone (default 500 bp).
+
+ Note:
+ Every genome must contain at least two backbone blocks, since a junction
+ is defined by its two flanking core blocks. A genome with fewer than two
+ (plausible with a high ``L_thr`` for a short or divergent genome) raises
+ ``ValueError`` when the split is computed, aborting the whole analysis.
+ To analyze the remaining genomes, lower ``L_thr`` or drop the offending
+ genome from ``pan`` beforehand.
+ """
+
+ def __init__(self, pan, L_thr: int = 500):
+ self.pan = pan
+ self.L_thr = L_thr
+ self._bdf = pan.to_blockstats_df()
+ self._junctions = None # dict[iso, list[Junction]]
+ self._edge_map = None # dict[edge_str, dict[iso, Junction]]
+
+ @property
+ def block_stats(self) -> pd.DataFrame:
+ """Block stats frame (``pan.to_blockstats_df()``) backing this analysis.
+
+ Indexed by str block id with a ``len`` and ``core`` column, among others.
+ Exposed read-only so callers (e.g. plots) need not reach into the private
+ cache.
+ """
+ return self._bdf
+
+ def _is_backbone(self, bid: str) -> bool:
+ """Determine if a block ID corresponds to a backbone block."""
+ return self._bdf.loc[bid, "core"] and self._bdf.loc[bid, "len"] >= self.L_thr
+
+ def _ensure_split(self):
+ """Ensure that paths have been split into junctions and edge maps built.
+
+ Raises:
+ ValueError: if any genome has fewer than two backbone blocks (see the
+ class docstring). The failure aborts the whole split, not just the
+ offending genome.
+ """
+ if self._junctions is not None:
+ return
+ self._junctions = {}
+ self._edge_map = {}
+ for name, path in self.pan.paths.items():
+ # create a Walk from a Path object, including accessory blocks
+ oriented_blocks = []
+ for nid in path.nodes:
+ row = self.pan.nodes.df.loc[str(nid)]
+ oriented_blocks.append(
+ JunctionNode(row["block_id"], row["strand"], nid)
+ )
+ tu_path = Walk(oriented_blocks, path.circular)
+
+ juncs = path_junction_split(tu_path, self._is_backbone)
+ self._junctions[name] = juncs
+ for j in juncs:
+ edge = j.flanking_edge()
+ if edge is None:
+ continue
+ edge_str = edge.to_str_id()
+ if edge_str not in self._edge_map:
+ self._edge_map[edge_str] = {}
+ self._edge_map[edge_str][name] = j
+
+ def __getitem__(self, edge_str: str) -> dict[str, Junction]:
+ """Return the {isolate -> Junction} mapping for a core edge.
+
+ Raises:
+ KeyError: if no junction with this edge exists in the graph.
+ """
+ self._ensure_split()
+ return self._edge_map[edge_str]
+
+ def __contains__(self, edge_str: str) -> bool:
+ """Whether a core edge with this id exists in the graph."""
+ self._ensure_split()
+ return edge_str in self._edge_map
+
+ def edges(self) -> list[str]:
+ """Return list of all edge string IDs."""
+ self._ensure_split()
+ return list(self._edge_map.keys())
+
+ def stats(self) -> pd.DataFrame:
+ """Compute per-edge junction statistics.
+
+ Returns:
+ DataFrame with edge string IDs as index and columns:
+ n_isolates, n_non_empty, n_categories, n_majority_category,
+ is_transitive, is_singleton, left_core_length, right_core_length,
+ accessory_length. Sorted by `n_isolates` descending.
+ """
+ self._ensure_split()
+ return junction_stats(self._edge_map, self._bdf)
+
+ def positions(self) -> pd.DataFrame:
+ """Find genomic positions of flanking core blocks for each junction.
+
+ Returns:
+ A DataFrame with MultiIndex (edge, iso) and columns:
+ - left_start, left_end: genomic position of the left flanking block
+ - right_start, right_end: genomic position of the right flanking block
+ - strand: True if canonical orientation, False if inverted
+ """
+ self._ensure_split()
+ return junction_positions(self._edge_map, self.pan)
+
+ def sequences(self, edge_str: str) -> list[SeqRecord]:
+ """Extract co-oriented sequences spanning a junction.
+
+ For each isolate with the given junction, returns a SeqRecord spanning
+ from the start of the left core block to the end of the right one.
+ All sequences are co-oriented: core blocks are in the same orientation
+ across isolates, with accessory sequence in between. Individual blocks
+ are reverse-complemented as needed based on their strand.
+
+ Args:
+ edge_str: The canonical edge string ID (e.g. "100_f__200_f").
+
+ Returns:
+ A list of SeqRecord objects, one per isolate. The record id is
+ the isolate name, and the description contains the edge string ID.
+ """
+ self._ensure_split()
+ return junction_sequences(self._edge_map, self.pan, edge_str)
diff --git a/packages/pypangraph/pypangraph/junctions/junction.py b/packages/pypangraph/pypangraph/junctions/junction.py
new file mode 100644
index 00000000..71626600
--- /dev/null
+++ b/packages/pypangraph/pypangraph/junctions/junction.py
@@ -0,0 +1,176 @@
+from __future__ import annotations
+
+from ..topology_utils import OrientedBlock, Walk, Edge
+
+
+class JunctionNode(OrientedBlock):
+ """OrientedBlock that also carries the unique node_id for unambiguous sequence lookup.
+
+ Inherits equality/hashing from OrientedBlock (compares block_id + strand only),
+ so it works transparently with Walk, Edge, and Junction.
+ """
+
+ def __init__(self, bid, strand: bool, node_id: int) -> None:
+ super().__init__(bid, strand)
+ self.node_id = node_id
+
+ def invert(self) -> "JunctionNode":
+ return JunctionNode(self.id, not self.strand, self.node_id)
+
+ def __repr__(self) -> str:
+ s = "+" if self.strand else "-"
+ return f"{{block={self.id}|{s}}}"
+
+
+class Junction:
+ """A junction is a segment of accessory blocks flanked by two core blocks,
+ with reverse-complement symmetry.
+
+ Attributes:
+ left: The core block on the left flank, or None for a terminal junction
+ on a linear path.
+ center: A Walk of accessory blocks between the flanks.
+ right: The core block on the right flank, or None for a terminal junction
+ on a linear path.
+ """
+
+ def __init__(
+ self, left: OrientedBlock | None, center: Walk, right: OrientedBlock | None
+ ) -> None:
+ self.left = left
+ self.center = center
+ self.right = right
+
+ def invert(self) -> "Junction":
+ """Reverse-complement the junction, swapping and inverting the flanks.
+
+ Robust to terminal junctions: a ``None`` flank stays ``None`` (it simply
+ moves to the opposite side).
+ """
+ left = self.right.invert() if self.right is not None else None
+ right = self.left.invert() if self.left is not None else None
+ return Junction(left, self.center.invert(), right)
+
+ def flanking_edge(self) -> Edge | None:
+ """Returns the Edge connecting the two flanking core blocks.
+
+ Returns None for terminal junctions on linear paths where one flank is missing.
+ """
+ if self.left is None or self.right is None:
+ return None
+ return Edge(self.left, self.right)
+
+ def is_canonical(self) -> bool:
+ """Whether the junction's flanks are in canonical edge orientation.
+
+ Delegates to the flanking edge's own canonical-direction predicate.
+
+ Raises:
+ ValueError: if the junction has a missing flank (terminal junction on
+ a linear path), for which canonical orientation is undefined.
+ """
+ edge = self.flanking_edge()
+ if edge is None:
+ raise ValueError(
+ "is_canonical() is undefined for terminal junctions with a missing flank"
+ )
+ return edge.is_canonical()
+
+ def to_canonical(self) -> "Junction":
+ """Return this junction co-oriented to its canonical edge direction.
+
+ Returns ``self`` when already canonical, otherwise the inverted junction.
+ Raises ``ValueError`` for terminal junctions with a missing flank (see
+ ``is_canonical``).
+ """
+ return self if self.is_canonical() else self.invert()
+
+ def oriented_blocks(self) -> list[OrientedBlock]:
+ """Flatten the junction into a left-to-right list of oriented blocks.
+
+ Includes the left flank (if any), then the center walk's blocks in order,
+ then the right flank (if any). Terminal junctions on linear paths may have
+ ``left`` or ``right`` equal to None; those flanks are silently omitted.
+ """
+ blocks: list[OrientedBlock] = []
+ if self.left is not None:
+ blocks.append(self.left)
+ blocks.extend(self.center.oriented_blocks)
+ if self.right is not None:
+ blocks.append(self.right)
+ return blocks
+
+ def __side_eq__(self, o: object) -> bool:
+ return self.left == o.left and self.center == o.center and self.right == o.right
+
+ def __eq__(self, o: object) -> bool:
+ if not isinstance(o, Junction):
+ return NotImplemented
+ return self.__side_eq__(o) or self.__side_eq__(o.invert())
+
+ def __side_hash__(self) -> int:
+ return hash((self.left, self.center, self.right))
+
+ def __hash__(self) -> int:
+ return self.__side_hash__() ^ self.invert().__side_hash__()
+
+ def __repr__(self) -> str:
+ return f"{self.left} <-- {self.center} --> {self.right}"
+
+
+def path_junction_split(path: Walk, is_core) -> list[Junction]:
+ """Split a path into junctions at core block boundaries.
+
+ Given a path and a boolean predicate on node ids, splits the path into a list
+ of Junctions. Each junction has flanking core blocks (for which the predicate
+ returns True) and a center path of non-core blocks.
+
+ For circular paths, the last junction wraps around the origin, connecting the
+ last core block back to the first.
+
+ For linear paths, there may be terminal junctions at the start and/or end with
+ only one flanking core block (left=None or right=None). These terminal junctions
+ have no flanking edge.
+
+ Args:
+ path: A Walk object (circular or linear).
+ is_core: A callable taking a block id and returning True if the block is core.
+
+ Returns:
+ A list of Junction objects.
+
+ Raises:
+ ValueError: if the path has fewer than two core blocks, in which case no
+ junction can be defined (and a circular path would have no flanks at all).
+ """
+ n_core = sum(1 for ob in path.oriented_blocks if is_core(ob.id))
+ if n_core < 2:
+ raise ValueError(
+ f"path has {n_core} core block(s); at least 2 are required to define a junction"
+ )
+
+ # collect here the junctions
+ junctions = []
+
+ current = []
+ left_node = None
+ for ob in path.oriented_blocks:
+ if is_core(ob.id):
+ J = Junction(left_node, Walk(current, circular=False), ob)
+ junctions.append(J)
+ left_node = ob
+ current = []
+ else:
+ current.append(ob)
+
+ if path.circular:
+ # complete periodic boundary: merge trailing non-core nodes into the first junction
+ J = junctions[0]
+ J.left = left_node
+ J.center = Walk(current + J.center.oriented_blocks, circular=False)
+ junctions[0] = J
+ elif current or left_node is not None:
+ # trailing accessory nodes after the last core block
+ junctions.append(Junction(left_node, Walk(current, circular=False), None))
+
+ return junctions
diff --git a/packages/pypangraph/pypangraph/junctions/positions.py b/packages/pypangraph/pypangraph/junctions/positions.py
new file mode 100644
index 00000000..6ea7aa7c
--- /dev/null
+++ b/packages/pypangraph/pypangraph/junctions/positions.py
@@ -0,0 +1,46 @@
+import pandas as pd
+
+
+def junction_positions(edge_map, pan) -> pd.DataFrame:
+ """Find genomic positions of the flanking core blocks for each junction.
+
+ For every (isolate, edge) the flanking blocks are looked up by their unique
+ node id, so the coordinates are unambiguous even if a block recurs in a genome.
+
+ The "left" and "right" labels follow each genome's own path order: ``left`` is
+ the flank encountered first when walking the path, ``right`` the one after the
+ accessory center. The ``strand`` column records whether that path order matches
+ the edge's canonical orientation (True) or is inverted relative to it (False).
+
+ Args:
+ edge_map: dict mapping edge string ID -> dict[isolate, Junction], as cached
+ by BackboneJunctions.
+ pan: A Pangraph object (used for node coordinate lookup).
+
+ Returns:
+ A DataFrame with MultiIndex (edge, iso) and columns:
+ - left_start, left_end: genomic position of the left flanking block
+ - right_start, right_end: genomic position of the right flanking block
+ - strand: True if the junction is in the canonical edge orientation
+ """
+ records = []
+ for edge_str, iso_junctions in edge_map.items():
+ for iso, junction in iso_junctions.items():
+ left_row = pan.nodes.df.loc[str(junction.left.node_id)]
+ right_row = pan.nodes.df.loc[str(junction.right.node_id)]
+ records.append(
+ {
+ "iso": iso,
+ "edge": edge_str,
+ "left_start": left_row["start"],
+ "left_end": left_row["end"],
+ "right_start": right_row["start"],
+ "right_end": right_row["end"],
+ "strand": junction.is_canonical(),
+ }
+ )
+
+ result = pd.DataFrame(records)
+ if result.empty:
+ return result
+ return result.set_index(["edge", "iso"])
diff --git a/packages/pypangraph/pypangraph/junctions/sequences.py b/packages/pypangraph/pypangraph/junctions/sequences.py
new file mode 100644
index 00000000..4e524d57
--- /dev/null
+++ b/packages/pypangraph/pypangraph/junctions/sequences.py
@@ -0,0 +1,50 @@
+from functools import cache
+
+from Bio.Seq import Seq
+from Bio.SeqRecord import SeqRecord
+
+
+def junction_sequences(edge_map, pan, edge_str: str) -> list[SeqRecord]:
+ """Extract co-oriented sequences spanning a junction.
+
+ For each isolate carrying the given junction, returns a SeqRecord spanning from
+ the start of the left core block to the end of the right one. All sequences are
+ co-oriented to the canonical edge direction, so the flanking core blocks line up
+ across isolates with the accessory content in between. Individual blocks are
+ reverse-complemented as needed based on their strand.
+
+ Args:
+ edge_map: dict mapping edge string ID -> dict[isolate, Junction], as cached
+ by BackboneJunctions.
+ pan: A Pangraph object (used for block sequence lookup).
+ edge_str: The canonical edge string ID (e.g. "100_f__200_f").
+
+ Returns:
+ A list of SeqRecord objects, one per isolate. The record id is the isolate
+ name and the description is the edge string ID. Empty if the edge is absent.
+ """
+ if edge_str not in edge_map:
+ return []
+
+ # to_sequences() regenerates a block's full alignment; the flanking core blocks
+ # recur in every isolate, so memoize per block id to compute each one only once.
+ @cache
+ def block_sequences(bid):
+ """Return the block's node_id -> sequence map, computing it once per block."""
+ return pan.blocks[bid].to_sequences()
+
+ records = []
+ for iso, junction in edge_map[edge_str].items():
+ oriented = junction.to_canonical()
+
+ seq_parts = []
+ for ob in oriented.oriented_blocks():
+ node_seq = block_sequences(ob.id)[str(ob.node_id)]
+ if not ob.strand:
+ node_seq = str(Seq(node_seq).reverse_complement())
+ seq_parts.append(node_seq)
+
+ record = SeqRecord(Seq("".join(seq_parts)), id=iso, description=edge_str)
+ records.append(record)
+
+ return records
diff --git a/packages/pypangraph/pypangraph/junctions/stats.py b/packages/pypangraph/pypangraph/junctions/stats.py
new file mode 100644
index 00000000..57274ce9
--- /dev/null
+++ b/packages/pypangraph/pypangraph/junctions/stats.py
@@ -0,0 +1,122 @@
+from collections import Counter
+
+import pandas as pd
+
+from ..topology_utils import Edge
+
+
+def _co_oriented_center_paths(iso_junctions):
+ """Co-orient all center paths for a given edge to canonical direction.
+
+ Junctions for the same edge may appear in opposite orientations on different
+ genomes. Before comparing center paths (e.g. for path categories), they must
+ be co-oriented so that the same accessory content is recognized regardless of
+ genomic strand.
+
+ Args:
+ iso_junctions: dict mapping isolate name -> Junction.
+
+ Returns:
+ dict mapping isolate name -> co-oriented center Walk.
+ """
+ result = {}
+ for iso, junction in iso_junctions.items():
+ result[iso] = junction.to_canonical().center
+ return result
+
+
+def _edge_stats(edge_str, iso_junctions, bdf):
+ """Compute statistics for a single edge.
+
+ Args:
+ edge_str: Canonical edge string ID.
+ iso_junctions: dict mapping isolate name -> Junction.
+ bdf: Block stats DataFrame (index=block_id, columns include 'len').
+
+ Returns:
+ dict with stat column names as keys.
+ """
+ edge = Edge.from_str_id(edge_str)
+ n_isolates = len(iso_junctions)
+ n_non_empty = sum(1 for j in iso_junctions.values() if len(j.center) > 0)
+
+ # Co-orient center paths, then group identical ones into categories. The empty
+ # path (junction with no accessory blocks) is a category in its own right, so it
+ # is counted like any other distinct center path.
+ center_paths = _co_oriented_center_paths(iso_junctions)
+ category_counts = Counter(center_paths.values())
+
+ n_categories = len(category_counts)
+ n_majority_category = max(category_counts.values())
+ is_transitive = n_categories == 1
+ is_singleton = n_isolates > 1 and n_majority_category == n_isolates - 1
+
+ # Core block lengths
+ left_core_length = bdf.loc[edge.left.id, "len"]
+ right_core_length = bdf.loc[edge.right.id, "len"]
+
+ # Unique accessory content: collect all distinct block IDs across all isolates
+ unique_block_ids = set()
+ for junction in iso_junctions.values():
+ for ob in junction.center.oriented_blocks:
+ unique_block_ids.add(ob.id)
+ accessory_length = sum(bdf.loc[bid, "len"] for bid in unique_block_ids)
+
+ return {
+ "n_isolates": n_isolates,
+ "n_non_empty": n_non_empty,
+ "n_categories": n_categories,
+ "n_majority_category": n_majority_category,
+ "is_transitive": is_transitive,
+ "is_singleton": is_singleton,
+ "left_core_length": left_core_length,
+ "right_core_length": right_core_length,
+ "accessory_length": accessory_length,
+ }
+
+
+def junction_stats(edge_map, bdf):
+ """Compute per-edge statistics for all junctions.
+
+ Args:
+ edge_map: dict mapping edge string ID -> dict[isolate, Junction].
+ Typically obtained from BackboneJunctions._edge_map.
+ bdf: Block stats DataFrame as returned by Pangraph.to_blockstats_df().
+ Must have 'len' column and block IDs as index.
+
+ Returns:
+ DataFrame with edge string IDs as index and columns:
+ - n_isolates: number of isolates with this junction
+ - n_non_empty: number of isolates whose center path has at least one
+ accessory block (`n_isolates - n_non_empty` gives the empty-junction count)
+ - n_categories: number of distinct center path variants
+ - n_majority_category: count of isolates in the most common variant
+ - is_transitive: True if only one variant exists
+ - is_singleton: True if all but one isolate share the same variant
+ - left_core_length: consensus length of left flanking core block
+ - right_core_length: consensus length of right flanking core block
+ - accessory_length: total unique accessory content (sum of distinct blocks consensus lengths)
+
+ Sorted by `n_isolates` descending.
+ """
+ records = {}
+ for edge_str, iso_junctions in edge_map.items():
+ records[edge_str] = _edge_stats(edge_str, iso_junctions, bdf)
+
+ df = pd.DataFrame.from_dict(records, orient="index")
+ df.index.name = "edge"
+ df = df.sort_values("n_isolates", ascending=False)
+
+ # Ensure integer types for count columns
+ for col in [
+ "n_isolates",
+ "n_non_empty",
+ "n_categories",
+ "n_majority_category",
+ "left_core_length",
+ "right_core_length",
+ "accessory_length",
+ ]:
+ df[col] = df[col].astype(int)
+
+ return df
diff --git a/packages/pypangraph/pypangraph/msu/msu.py b/packages/pypangraph/pypangraph/minimal_synteny_units.py
similarity index 67%
rename from packages/pypangraph/pypangraph/msu/msu.py
rename to packages/pypangraph/pypangraph/minimal_synteny_units.py
index 069a5148..2c96f60c 100644
--- a/packages/pypangraph/pypangraph/msu/msu.py
+++ b/packages/pypangraph/pypangraph/minimal_synteny_units.py
@@ -1,36 +1,38 @@
-from . import msu_utils as ut
from collections import defaultdict
+
import numpy as np
+from . import topology_utils as tu
+
def core_paths(pan, L_thr):
bdf = pan.to_blockstats_df()
- paths = ut.pangraph_to_path_dict(pan)
+ walks = tu.pangraph_to_walks(pan)
- def is_core(node_id):
- return (bdf.loc[node_id, "len"] >= L_thr) and bdf.loc[node_id, "core"]
+ def is_core(block_id):
+ return (bdf.loc[block_id, "len"] >= L_thr) and bdf.loc[block_id, "core"]
- return ut.filter_paths(paths, is_core)
+ return tu.filter_walks(walks, is_core)
def flip_msu_to_most_common_orientation(paths):
orient = defaultdict(int)
for iso, p in paths.items():
- for n in p.nodes:
- msu_id, strand = n.id, n.strand
+ for ob in p.oriented_blocks:
+ msu_id, strand = ob.id, ob.strand
orient[msu_id] += 1 if strand else -1
# flip all the ones with orient < 0
for iso, p in paths.items():
- nodes = [n.invert() if orient[n.id] < 0 else n for n in p.nodes]
- paths[iso] = ut.Path(nodes, p.circular)
+ flipped = [ob.invert() if orient[ob.id] < 0 else ob for ob in p.oriented_blocks]
+ paths[iso] = tu.Walk(flipped, p.circular)
return paths
def minimal_synteny_units(pan, L_thr: int, rotate: bool = True):
c_paths = core_paths(pan, L_thr)
- mergers = ut.find_mergers(c_paths)
+ mergers = tu.find_mergers(c_paths)
# MSU lengths
B_len = pan.to_blockstats_df()["len"].to_dict()
@@ -42,7 +44,7 @@ def minimal_synteny_units(pan, L_thr: int, rotate: bool = True):
MSU_order = sorted(MSU_len, key=MSU_len.get, reverse=True)
# simplify paths
- MSU_paths = ut.filter_paths(c_paths, lambda x: x in MSU_order)
+ MSU_paths = tu.filter_walks(c_paths, lambda x: x in MSU_order)
# rename MSUs
MSU_ids = {msu: f"MSU_{i}" for i, msu in enumerate(MSU_order)}
@@ -51,9 +53,8 @@ def minimal_synteny_units(pan, L_thr: int, rotate: bool = True):
MSU_mergers = {source: MSU_ids[sink] for source, sink in mergers.items()}
if rotate:
- assert np.all(
- [p.circular for p in MSU_paths.values()]
- ), "Only circular paths can be rotated"
+ if not np.all([p.circular for p in MSU_paths.values()]):
+ raise ValueError("Only circular paths can be rotated")
focal_block = max(MSU_len, key=MSU_len.get)
focal_strandedness = True
MSU_paths = {
diff --git a/packages/pypangraph/pypangraph/msu/__init__.py b/packages/pypangraph/pypangraph/msu/__init__.py
deleted file mode 100644
index 7a07c113..00000000
--- a/packages/pypangraph/pypangraph/msu/__init__.py
+++ /dev/null
@@ -1,11 +0,0 @@
-from .msu import minimal_synteny_units
-from .msu_utils import (
- Node,
- Edge,
- Path,
- pangraph_to_path_dict,
- find_mergers,
- filter_paths,
- path_block_count,
- path_edge_count,
-)
diff --git a/packages/pypangraph/pypangraph/msu/msu_utils.py b/packages/pypangraph/pypangraph/msu/msu_utils.py
deleted file mode 100644
index e5b89b5d..00000000
--- a/packages/pypangraph/pypangraph/msu/msu_utils.py
+++ /dev/null
@@ -1,236 +0,0 @@
-from collections import defaultdict, Counter
-
-
-class Node:
- """Combination of block id and strandedness"""
-
- def __init__(self, bid: str, strand: bool) -> None:
- self.id = bid
- self.strand = strand
-
- def invert(self) -> "Node":
- return Node(self.id, not self.strand)
-
- def __invert__(self) -> "Node":
- return self.invert()
-
- def __eq__(self, other: object) -> bool:
- return self.id == other.id and self.strand == other.strand
-
- def __hash__(self) -> int:
- return hash((self.id, self.strand))
-
- def __repr__(self) -> str:
- s = "+" if self.strand else "-"
- return f"[{self.id}|{s}]"
-
- def to_str_id(self):
- s = "f" if self.strand else "r"
- return f"{self.id}_{s}"
-
- @staticmethod
- def from_str_id(t) -> "Node":
- bid = t.split("_")[0]
- strand = True if t.split("_")[1] == "f" else False
- return Node(bid, strand)
-
-
-class Path:
- """A path is a list of nodes"""
-
- def __init__(self, nodes=[], circular=None) -> None:
- self.nodes = nodes
- self.circular = circular
-
- def add_left(self, node: Node) -> None:
- self.nodes.insert(0, node)
-
- def add_right(self, node: Node) -> None:
- self.nodes.append(node)
-
- def rotate_to(self, bid: str, strand: bool) -> "Path":
- assert self.circular, "Path is not circular"
- assert bid in [n.id for n in self.nodes], "Block not in path"
- n = Node(bid, strand)
- if n in self.nodes:
- idx = self.nodes.index(n)
- p = Path(self.nodes[idx:] + self.nodes[:idx], circular=True)
- else:
- p = self.invert()
- idx = p.nodes.index(n)
- p = Path(p.nodes[idx:] + p.nodes[:idx], circular=True)
- return p
-
- def rotate_to_node(self, node: Node) -> "Path":
- return self.rotate_to(node.id, node.strand)
-
- def invert(self) -> "Path":
- return Path([n.invert() for n in self.nodes[::-1]], circular=self.circular)
-
- def __invert__(self) -> "Path":
- return self.invert()
-
- def __eq__(self, o: object) -> bool:
- return self.nodes == o.nodes
-
- def __hash__(self) -> int:
- return hash(tuple(self.nodes))
-
- def __repr__(self) -> str:
- return "_".join([str(n) for n in self.nodes])
-
- def __len__(self) -> int:
- return len(self.nodes)
-
- def rename_bids(self, bid_dict: dict) -> "Path":
- return Path(
- [Node(bid_dict[n.id], n.strand) for n in self.nodes], circular=self.circular
- )
-
- def to_list(self):
- return [n.to_str_id() for n in self.nodes]
-
- @staticmethod
- def from_list(path_list: list[Node], circular: bool) -> "Path":
- return Path([Node.from_str_id(nid) for nid in path_list])
-
-
-class Edge:
- """Oriented link between two nodes/paths"""
-
- def __init__(self, left, right) -> None:
- self.left = left
- self.right = right
-
- def invert(self) -> "Edge":
- return Edge(self.right.invert(), self.left.invert())
-
- def __invert__(self) -> "Edge":
- return self.invert()
-
- def __side_eq__(self, o: object) -> bool:
- return self.left == o.left and self.right == o.right
-
- def __eq__(self, o: object) -> bool:
- return self.__side_eq__(o) or self.__side_eq__(o.invert())
-
- def __side_hash__(self) -> int:
- return hash((self.left, self.right))
-
- def __hash__(self) -> int:
- return self.__side_hash__() ^ self.invert().__side_hash__()
-
- def __repr__(self) -> str:
- return f"{self.left} <--> {self.right}"
-
- def __to_str_id(self) -> str:
- return "__".join([self.left.to_str_id(), self.right.to_str_id()])
-
- def to_str_id(self) -> str:
- A = self.__to_str_id()
- B = self.invert().__to_str_id()
- return A if A < B else B
-
- @staticmethod
- def from_str_id(t) -> "Edge":
- left, right = t.split("__")
- return Edge(Node.from_str_id(left), Node.from_str_id(right))
-
-
-def pangraph_to_path_dict(pan):
- """Creates a dictionary isolate -> path objects"""
- res = {}
- for name, path in pan.paths.items():
- B, S = pan.nodes.nodes_to_blocks(path.nodes)
- nodes = [Node(b, s) for b, s in zip(B, S)]
- res[name] = Path(nodes, path.circular)
- return res
-
-
-def filter_paths(paths, keep_f):
- """Given a filter function, removes nodes that fail the condition from
- the path dictionaries."""
- res = {}
- for iso, path in paths.items():
- filt_path = Path(
- [node for node in path.nodes if keep_f(node.id)],
- path.circular,
- )
- res[iso] = filt_path
- return res
-
-
-def path_categories(paths):
- """Returns a list of touples, one per non-empty path, with the following info:
- (count, path, [list of isolates])"""
- iso_list = defaultdict(list)
- n_paths = defaultdict(int)
- nodes = {}
- for iso, path in paths.items():
- if len(path.nodes) > 0:
- n_paths[path] += 1
- iso_list[path].append(iso)
- nodes[path] = path.nodes
-
- # sort by count
- path_cat = [(count, nodes[path], iso_list[path]) for path, count in n_paths.items()]
- path_cat.sort(key=lambda x: x[0], reverse=True)
- return path_cat
-
-
-def path_edge_count(paths):
- """Count internal edges of paths"""
- ct = Counter()
- for iso, p in paths.items():
- L = len(p.nodes)
- es = []
- for i in range(L - 1):
- e = Edge(p.nodes[i], p.nodes[i + 1])
- es.append(e)
- if p.circular:
- e = Edge(p.nodes[-1], p.nodes[0])
- es.append(e)
- ct.update(es)
- return dict(ct)
-
-
-def path_block_count(paths):
- """Count internal blocks of paths"""
- ct = Counter()
- for iso, p in paths.items():
- for node in p.nodes:
- ct.update([node.id])
- return dict(ct)
-
-
-def find_mergers(paths):
- """Create a dictionary source -> sinks of block-ids to be merged"""
- edge_ct = path_edge_count(paths)
- block_ct = path_block_count(paths)
-
- mergers = {}
- for e, ec in edge_ct.items():
- bl, br = e.left.id, e.right.id
- if (ec == block_ct[bl]) and (ec == block_ct[br]):
- # merge
- if bl in mergers:
- if br in mergers:
- source = mergers[br]
- sink = mergers[bl]
- for k in mergers:
- if mergers[k] == source:
- mergers[k] = sink
- else:
- mergers[br] = mergers[bl]
- elif br in mergers:
- mergers[bl] = mergers[br]
- else:
- mergers[br] = bl
- mergers[bl] = bl
-
- # add missing blocks that are not in a merger
- for bid in block_ct.keys():
- if bid not in mergers:
- mergers[bid] = bid
-
- return mergers
diff --git a/packages/pypangraph/pypangraph/plots/__init__.py b/packages/pypangraph/pypangraph/plots/__init__.py
index c17e107b..439d9c01 100644
--- a/packages/pypangraph/pypangraph/plots/__init__.py
+++ b/packages/pypangraph/pypangraph/plots/__init__.py
@@ -1 +1,4 @@
from .dotplot import dotplot
+from .linear_junction_plot import linear_junction_plot
+
+__all__ = ["dotplot", "linear_junction_plot"]
diff --git a/packages/pypangraph/pypangraph/plots/dotplot.py b/packages/pypangraph/pypangraph/plots/dotplot.py
index efb29650..1488031f 100644
--- a/packages/pypangraph/pypangraph/plots/dotplot.py
+++ b/packages/pypangraph/pypangraph/plots/dotplot.py
@@ -1,5 +1,4 @@
# %%
-import matplotlib.pyplot as plt
from collections import defaultdict
from itertools import product
from dataclasses import dataclass
@@ -54,8 +53,8 @@ def split_fraction(self, f, flip, L):
Split the segment into two segments based on a fraction of the segment length
"""
assert 0 <= f <= 1
- l = self.length(L)
- l1, l2 = f * l, (1 - f) * l
+ total_len = self.length(L)
+ l1, l2 = f * total_len, (1 - f) * total_len
if flip:
l1, l2 = l2, l1
s1 = Segment(self.start, (self.start + l1) % L, self.strand)
diff --git a/packages/pypangraph/pypangraph/plots/linear_junction_plot.py b/packages/pypangraph/pypangraph/plots/linear_junction_plot.py
new file mode 100644
index 00000000..3d52cf32
--- /dev/null
+++ b/packages/pypangraph/pypangraph/plots/linear_junction_plot.py
@@ -0,0 +1,96 @@
+"""Linear schematic plot of a junction across isolates."""
+
+from __future__ import annotations
+
+import random
+from collections import defaultdict
+
+
+def _random_color() -> tuple[float, float, float]:
+ """Return a mid-saturation RGB triple, each channel uniform in [0.3, 0.95]."""
+ return (
+ random.uniform(0.3, 0.95),
+ random.uniform(0.3, 0.95),
+ random.uniform(0.3, 0.95),
+ )
+
+
+def linear_junction_plot(
+ ax,
+ bj,
+ edge: str,
+ *,
+ isolates: list[str] | None = None,
+ color_map: dict | None = None,
+ left_flank_color="C0",
+ right_flank_color="C1",
+ highlight_inverted: bool = False,
+) -> dict:
+ """Draw a per-isolate linear schematic of one core-edge junction.
+
+ Each isolate becomes a row of horizontal bars (one per oriented block),
+ with bar width set to the block's consensus length. Junctions are
+ co-oriented via ``Junction.to_canonical()`` so the left flank lines up
+ across rows.
+
+ Args:
+ ax: A Matplotlib axes object; mutated in place.
+ bj: A ``BackboneJunctions`` instance.
+ edge: Canonical edge string ID (e.g. ``"100_f__200_f"``).
+ isolates: Optional explicit row order. Defaults to
+ ``sorted(bj[edge].keys())``. Unknown names raise ``KeyError``.
+ color_map: Optional seed mapping ``block_id -> color`` for accessory
+ blocks. Missing entries are filled in with random mid-saturated
+ colors. The input is not mutated. Flank colors always override
+ this mapping.
+ left_flank_color: Face color of the left-flank block.
+ right_flank_color: Face color of the right-flank block.
+ highlight_inverted: When ``True``, blocks with negative strand get
+ a red border; otherwise all borders are black.
+
+ Returns:
+ A plain ``dict`` mapping accessory block id to color, including any
+ newly-generated random colors. Pass it as ``color_map`` to a
+ subsequent call to keep block colors consistent across panels.
+ """
+ per_iso = bj[edge]
+
+ if isolates is None:
+ isolates = sorted(per_iso.keys())
+
+ seed = color_map or {}
+ colors = defaultdict(_random_color, seed)
+
+ block_len = bj.block_stats["len"]
+
+ for row, iso in enumerate(isolates):
+ J = per_iso[iso].to_canonical()
+ obs = J.oriented_blocks()
+ n = len(obs)
+ x = 0
+ for i, ob in enumerate(obs):
+ length = block_len.loc[ob.id]
+ if i == 0:
+ face = left_flank_color
+ elif i == n - 1:
+ face = right_flank_color
+ else:
+ face = colors[ob.id]
+ edge_c = "red" if (highlight_inverted and not ob.strand) else "black"
+ ax.barh(
+ row,
+ length,
+ left=x,
+ height=0.8,
+ color=face,
+ edgecolor=edge_c,
+ linewidth=0.4,
+ )
+ x += length
+
+ ax.set_yticks(range(len(isolates)))
+ ax.set_yticklabels(isolates, fontsize=8)
+ ax.invert_yaxis()
+ ax.set_xlabel("position along junction (bp)")
+
+ return dict(colors)
diff --git a/packages/pypangraph/pypangraph/topology_utils.py b/packages/pypangraph/pypangraph/topology_utils.py
new file mode 100644
index 00000000..c2270cc0
--- /dev/null
+++ b/packages/pypangraph/pypangraph/topology_utils.py
@@ -0,0 +1,271 @@
+from collections import Counter, defaultdict
+
+
+class OrientedBlock:
+ """Combination of block id and strandedness"""
+
+ def __init__(self, bid: str, strand: bool) -> None:
+ self.id = bid
+ self.strand = strand
+
+ def invert(self) -> "OrientedBlock":
+ return OrientedBlock(self.id, not self.strand)
+
+ def __invert__(self) -> "OrientedBlock":
+ return self.invert()
+
+ def __eq__(self, other: object) -> bool:
+ if not isinstance(other, OrientedBlock):
+ return NotImplemented
+ return self.id == other.id and self.strand == other.strand
+
+ def __hash__(self) -> int:
+ return hash((self.id, self.strand))
+
+ def __repr__(self) -> str:
+ s = "+" if self.strand else "-"
+ return f"[{self.id}|{s}]"
+
+ def to_str_id(self):
+ s = "f" if self.strand else "r"
+ return f"{self.id}_{s}"
+
+ @staticmethod
+ def from_str_id(t) -> "OrientedBlock":
+ # Block ids are kept as strings internally (they are u64 hashes; see the
+ # pypangraph data-model notes), so the id token round-trips as-is. Split
+ # from the right so that ids containing "_" (e.g. MSU-renamed "MSU_0")
+ # round-trip correctly: only the trailing "_" strand suffix is peeled.
+ bid_str, strand_str = t.rsplit("_", 1)
+ return OrientedBlock(bid_str, strand_str == "f")
+
+
+class Walk:
+ """An ordered traversal of oriented blocks (a walk through the block graph)."""
+
+ def __init__(self, oriented_blocks=None, circular=None) -> None:
+ if oriented_blocks is None:
+ oriented_blocks = []
+ self.oriented_blocks = oriented_blocks
+ self.circular = circular
+
+ def add_left(self, oriented_block: OrientedBlock) -> None:
+ self.oriented_blocks.insert(0, oriented_block)
+
+ def add_right(self, oriented_block: OrientedBlock) -> None:
+ self.oriented_blocks.append(oriented_block)
+
+ def rotate_to(self, bid: str, strand: bool) -> "Walk":
+ if not self.circular:
+ raise ValueError("Walk is not circular")
+ if bid not in [ob.id for ob in self.oriented_blocks]:
+ raise ValueError(f"Block {bid} not in walk")
+ target = OrientedBlock(bid, strand)
+ if target in self.oriented_blocks:
+ idx = self.oriented_blocks.index(target)
+ return Walk(
+ self.oriented_blocks[idx:] + self.oriented_blocks[:idx], circular=True
+ )
+ inv = self.invert()
+ idx = inv.oriented_blocks.index(target)
+ return Walk(
+ inv.oriented_blocks[idx:] + inv.oriented_blocks[:idx], circular=True
+ )
+
+ def rotate_to_oriented_block(self, oriented_block: OrientedBlock) -> "Walk":
+ return self.rotate_to(oriented_block.id, oriented_block.strand)
+
+ def invert(self) -> "Walk":
+ return Walk(
+ [ob.invert() for ob in self.oriented_blocks[::-1]], circular=self.circular
+ )
+
+ def __invert__(self) -> "Walk":
+ return self.invert()
+
+ def __eq__(self, o: object) -> bool:
+ if not isinstance(o, Walk):
+ return NotImplemented
+ return self.oriented_blocks == o.oriented_blocks
+
+ def __hash__(self) -> int:
+ return hash(tuple(self.oriented_blocks))
+
+ def __repr__(self) -> str:
+ return " ".join([str(ob) for ob in self.oriented_blocks])
+
+ def __len__(self) -> int:
+ return len(self.oriented_blocks)
+
+ def rename_bids(self, bid_dict: dict) -> "Walk":
+ return Walk(
+ [OrientedBlock(bid_dict[ob.id], ob.strand) for ob in self.oriented_blocks],
+ circular=self.circular,
+ )
+
+ def edges(self) -> list:
+ """Edges between consecutive oriented blocks (with the circular wrap)."""
+ obs = self.oriented_blocks
+ es = [Edge(a, b) for a, b in zip(obs, obs[1:])]
+ if self.circular and len(obs) > 1:
+ es.append(Edge(obs[-1], obs[0]))
+ return es
+
+
+class Edge:
+ """Oriented link between two oriented blocks."""
+
+ def __init__(self, left, right) -> None:
+ self.left = left
+ self.right = right
+
+ def invert(self) -> "Edge":
+ return Edge(self.right.invert(), self.left.invert())
+
+ def __invert__(self) -> "Edge":
+ return self.invert()
+
+ def __side_eq__(self, o: object) -> bool:
+ return self.left == o.left and self.right == o.right
+
+ def __eq__(self, o: object) -> bool:
+ if not isinstance(o, Edge):
+ return NotImplemented
+ return self.__side_eq__(o) or self.__side_eq__(o.invert())
+
+ def __side_hash__(self) -> int:
+ return hash((self.left, self.right))
+
+ def __hash__(self) -> int:
+ return self.__side_hash__() ^ self.invert().__side_hash__()
+
+ def __repr__(self) -> str:
+ return f"{self.left} <--> {self.right}"
+
+ def _natural_str_id(self) -> str:
+ """Join left and right ids in their stored order, without canonicalizing."""
+ return "__".join([self.left.to_str_id(), self.right.to_str_id()])
+
+ def is_canonical(self) -> bool:
+ """Whether (left, right) is the canonical (lex-min) orientation of this edge.
+
+ Ties (RC-palindromic edges) resolve canonical by convention.
+ """
+ return self._natural_str_id() <= self.invert()._natural_str_id()
+
+ def to_str_id(self) -> str:
+ return (
+ self._natural_str_id()
+ if self.is_canonical()
+ else self.invert()._natural_str_id()
+ )
+
+ @staticmethod
+ def from_str_id(t) -> "Edge":
+ left, right = t.split("__")
+ return Edge(OrientedBlock.from_str_id(left), OrientedBlock.from_str_id(right))
+
+
+def pangraph_to_walks(pan):
+ """Creates a dictionary isolate -> Walk objects"""
+ res = {}
+ for name, path in pan.paths.items():
+ B, S = pan.nodes.nodes_to_blocks(path.nodes)
+ oriented_blocks = [OrientedBlock(b, s) for b, s in zip(B, S)]
+ res[name] = Walk(oriented_blocks, path.circular)
+ return res
+
+
+def filter_walks(walks, keep_f):
+ """Given a filter function, removes oriented blocks that fail the condition
+ from the walk dictionary."""
+ res = {}
+ for iso, walk in walks.items():
+ filt_walk = Walk(
+ [ob for ob in walk.oriented_blocks if keep_f(ob.id)],
+ walk.circular,
+ )
+ res[iso] = filt_walk
+ return res
+
+
+def walk_categories(walks):
+ """Returns a list of tuples, one per non-empty walk, with the following info:
+ (count, oriented_blocks, [list of isolates])"""
+ iso_list = defaultdict(list)
+ n_walks = defaultdict(int)
+ walk_obs = {}
+ for iso, walk in walks.items():
+ if len(walk.oriented_blocks) > 0:
+ n_walks[walk] += 1
+ iso_list[walk].append(iso)
+ walk_obs[walk] = walk.oriented_blocks
+
+ # sort by count
+ walk_cat = [
+ (count, walk_obs[walk], iso_list[walk]) for walk, count in n_walks.items()
+ ]
+ walk_cat.sort(key=lambda x: x[0], reverse=True)
+ return walk_cat
+
+
+def walk_edge_count(walks):
+ """Count internal edges of walks"""
+ ct = Counter()
+ for iso, w in walks.items():
+ ct.update(w.edges())
+ return dict(ct)
+
+
+def walk_block_count(walks):
+ """Count occurrences of each block across walks"""
+ ct = Counter()
+ for iso, w in walks.items():
+ for ob in w.oriented_blocks:
+ ct.update([ob.id])
+ return dict(ct)
+
+
+def find_mergers(walks):
+ """Create a dictionary source -> sink of block-ids to be merged.
+
+ Two core blocks are merged when they are *always* adjacent: an edge whose
+ occurrence count equals the occurrence count of both its endpoint blocks.
+ Groups of such blocks are tracked with a union-find (disjoint-set) forest,
+ so fusing two groups is near-constant time instead of a full relabel scan.
+
+ Returns a dict mapping every block id to its group's representative block id
+ (a block in a group of one maps to itself).
+ """
+ edge_ct = walk_edge_count(walks)
+ block_ct = walk_block_count(walks)
+
+ # Union-find over block ids: every block starts in its own group,
+ # represented by itself.
+ parent = {bid: bid for bid in block_ct}
+
+ def find_root(bid):
+ """Return the representative (root) of bid's group, compressing the path."""
+ # Walk up to the root of the group.
+ root = bid
+ while parent[root] != root:
+ root = parent[root]
+ # Second pass: point every node visited directly at the root,
+ # so later lookups are fast.
+ while parent[bid] != root:
+ parent[bid], bid = root, parent[bid]
+ return root
+
+ def union(a, b):
+ """Merge the groups containing blocks a and b."""
+ parent[find_root(a)] = find_root(b)
+
+ # An edge whose count equals both endpoints' block counts means the two
+ # blocks are always adjacent -> they belong to the same merger group.
+ for e, ec in edge_ct.items():
+ bl, br = e.left.id, e.right.id
+ if (ec == block_ct[bl]) and (ec == block_ct[br]):
+ union(bl, br)
+
+ # Flatten: map every block directly to its group representative.
+ return {bid: find_root(bid) for bid in block_ct}
diff --git a/packages/pypangraph/tests/conftest.py b/packages/pypangraph/tests/conftest.py
new file mode 100644
index 00000000..006e1894
--- /dev/null
+++ b/packages/pypangraph/tests/conftest.py
@@ -0,0 +1,410 @@
+"""Shared test fixtures for pypangraph tests."""
+
+import pytest
+import pypangraph as pp
+
+
+def _make_node(node_id, block_id, path_id, strand, start, end):
+ return {
+ "id": node_id,
+ "block_id": block_id,
+ "path_id": path_id,
+ "strand": "+" if strand else "-",
+ "position": [start, end],
+ }
+
+
+def _make_block(block_id, length, node_ids):
+ return {
+ "id": block_id,
+ "consensus": "A" * length,
+ "alignments": {
+ str(nid): {"subs": [], "dels": [], "inss": []} for nid in node_ids
+ },
+ }
+
+
+def _make_path(path_id, name, node_ids, tot_len, circular=True):
+ return {
+ "id": path_id,
+ "nodes": node_ids,
+ "tot_len": tot_len,
+ "circular": circular,
+ "name": name,
+ "desc": None,
+ }
+
+
+def build_junction_pangraph_json():
+ """Build a synthetic pangraph JSON with core and accessory blocks.
+
+ Graph topology (all paths circular):
+ s1: C1+ A1+ A2+ C2+ C3+ C4+
+ s2: C1+ A1+ C2+ A3+ C3+ C4+
+ s3: C1+ A2- C3+ C2+ A3+ C4+ (C2/C3 swapped = rearrangement)
+
+ Core blocks (>=500bp, present once per strain):
+ C1=100 (1000bp), C2=200 (800bp), C3=300 (600bp), C4=400 (700bp)
+ Accessory blocks (<500bp):
+ A1=500 (200bp, s1+s2), A2=600 (150bp, s1+s3), A3=700 (300bp, s2+s3)
+ """
+ # --- nodes ---
+ # s1 (path 0): C1+(1) A1+(2) A2+(3) C2+(4) C3+(5) C4+(6)
+ # s2 (path 1): C1+(7) A1+(8) C2+(9) A3+(10) C3+(11) C4+(12)
+ # s3 (path 2): C1+(13) A2-(14) C3+(15) C2+(16) A3+(17) C4+(18)
+ nodes = {
+ # s1
+ "1": _make_node(1, 100, 0, True, 0, 1000),
+ "2": _make_node(2, 500, 0, True, 1000, 1200),
+ "3": _make_node(3, 600, 0, True, 1200, 1350),
+ "4": _make_node(4, 200, 0, True, 1350, 2150),
+ "5": _make_node(5, 300, 0, True, 2150, 2750),
+ "6": _make_node(6, 400, 0, True, 2750, 3450),
+ # s2
+ "7": _make_node(7, 100, 1, True, 0, 1000),
+ "8": _make_node(8, 500, 1, True, 1000, 1200),
+ "9": _make_node(9, 200, 1, True, 1200, 2000),
+ "10": _make_node(10, 700, 1, True, 2000, 2300),
+ "11": _make_node(11, 300, 1, True, 2300, 2900),
+ "12": _make_node(12, 400, 1, True, 2900, 3600),
+ # s3
+ "13": _make_node(13, 100, 2, True, 0, 1000),
+ "14": _make_node(14, 600, 2, False, 1000, 1150),
+ "15": _make_node(15, 300, 2, True, 1150, 1750),
+ "16": _make_node(16, 200, 2, True, 1750, 2550),
+ "17": _make_node(17, 700, 2, True, 2550, 2850),
+ "18": _make_node(18, 400, 2, True, 2850, 3550),
+ }
+
+ blocks = {
+ "100": _make_block(100, 1000, [1, 7, 13]),
+ "200": _make_block(200, 800, [4, 9, 16]),
+ "300": _make_block(300, 600, [5, 11, 15]),
+ "400": _make_block(400, 700, [6, 12, 18]),
+ "500": _make_block(500, 200, [2, 8]),
+ "600": _make_block(600, 150, [3, 14]),
+ "700": _make_block(700, 300, [10, 17]),
+ }
+
+ paths = {
+ "0": _make_path(0, "s1", [1, 2, 3, 4, 5, 6], 3450),
+ "1": _make_path(1, "s2", [7, 8, 9, 10, 11, 12], 3600),
+ "2": _make_path(2, "s3", [13, 14, 15, 16, 17, 18], 3550),
+ }
+
+ return {"paths": paths, "blocks": blocks, "nodes": nodes}
+
+
+def build_linear_pangraph_json():
+ """Build a synthetic pangraph JSON with linear (non-circular) paths.
+
+ Graph topology (all paths linear):
+ s1: c5+ C1+ A2+ C2+ C3+
+ s2: C1+ A3+ C2+ C3+ c5+
+
+ Backbone blocks (core AND >=500bp, present once per strain):
+ C1=100 (1000bp), C2=200 (800bp), C3=300 (600bp)
+ Accessory blocks (present in only one strain):
+ A2=600 (150bp, s1 only), A3=700 (300bp, s2 only)
+ Core-but-not-backbone block (lowercase 'c' to flag it):
+ c5=500 (200bp, present once in both strains)
+
+ c5 is strictly a `core` block (it occurs exactly once in every strain), but its
+ 200bp length is below the backbone threshold used by the junction/MSU analyses
+ (core AND len>=L_thr), so it cannot be a backbone block and is handled as accessory
+ content. It is the same block at a leading position in s1 and a trailing position in
+ s2, exercising terminal junctions at both ends.
+ """
+ # s1 (path 0): c5+(1) C1+(2) A2+(3) C2+(4) C3+(5)
+ # s2 (path 1): C1+(6) A3+(7) C2+(8) C3+(9) c5+(10)
+ nodes = {
+ # s1
+ "1": _make_node(1, 500, 0, True, 0, 200),
+ "2": _make_node(2, 100, 0, True, 200, 1200),
+ "3": _make_node(3, 600, 0, True, 1200, 1350),
+ "4": _make_node(4, 200, 0, True, 1350, 2150),
+ "5": _make_node(5, 300, 0, True, 2150, 2750),
+ # s2
+ "6": _make_node(6, 100, 1, True, 0, 1000),
+ "7": _make_node(7, 700, 1, True, 1000, 1300),
+ "8": _make_node(8, 200, 1, True, 1300, 2100),
+ "9": _make_node(9, 300, 1, True, 2100, 2700),
+ "10": _make_node(10, 500, 1, True, 2700, 2900),
+ }
+
+ blocks = {
+ "100": _make_block(100, 1000, [2, 6]),
+ "200": _make_block(200, 800, [4, 8]),
+ "300": _make_block(300, 600, [5, 9]),
+ "500": _make_block(500, 200, [1, 10]),
+ "600": _make_block(600, 150, [3]),
+ "700": _make_block(700, 300, [7]),
+ }
+
+ paths = {
+ "0": _make_path(0, "s1", [1, 2, 3, 4, 5], 2750, circular=False),
+ "1": _make_path(1, "s2", [6, 7, 8, 9, 10], 2900, circular=False),
+ }
+
+ return {"paths": paths, "blocks": blocks, "nodes": nodes}
+
+
+def build_tandem_accessory_pangraph_json():
+ """Pangraph where an accessory block is duplicated on one path within a junction.
+
+ Two circular strains; backbone blocks C1, C2 flank a single junction whose
+ accessory content is block A:
+ s1: C1+ A+ A+ C2+ (A appears twice -- a tandem duplication)
+ s2: C1+ A+ C2+ (A appears once)
+
+ Backbone blocks (core AND >=500bp): C1=100 (1000bp), C2=200 (800bp).
+ Accessory block (<500bp): A=500 (200bp), with two copies in s1 and one in s2,
+ so the C1-C2 junction traverses A three times in total.
+ """
+ # s1 (path 0): C1+(1) A+(2) A+(3) C2+(4)
+ # s2 (path 1): C1+(5) A+(6) C2+(7)
+ nodes = {
+ # s1
+ "1": _make_node(1, 100, 0, True, 0, 1000),
+ "2": _make_node(2, 500, 0, True, 1000, 1200),
+ "3": _make_node(3, 500, 0, True, 1200, 1400),
+ "4": _make_node(4, 200, 0, True, 1400, 2200),
+ # s2
+ "5": _make_node(5, 100, 1, True, 0, 1000),
+ "6": _make_node(6, 500, 1, True, 1000, 1200),
+ "7": _make_node(7, 200, 1, True, 1200, 2000),
+ }
+
+ blocks = {
+ "100": _make_block(100, 1000, [1, 5]),
+ "200": _make_block(200, 800, [4, 7]),
+ "500": _make_block(500, 200, [2, 3, 6]),
+ }
+
+ paths = {
+ "0": _make_path(0, "s1", [1, 2, 3, 4], 2200),
+ "1": _make_path(1, "s2", [5, 6, 7], 2000),
+ }
+
+ return {"paths": paths, "blocks": blocks, "nodes": nodes}
+
+
+def _make_block_with_edits(block_id, consensus, node_edits):
+ """Create a block with a custom consensus and per-node edits.
+
+ Args:
+ block_id: Block identifier.
+ consensus: The consensus sequence string.
+ node_edits: dict of node_id -> {"subs": [...], "dels": [...], "inss": [...]}
+ or node_id -> None for identity (no edits).
+ """
+ alignments = {}
+ for nid, edits in node_edits.items():
+ if edits is None:
+ edits = {"subs": [], "dels": [], "inss": []}
+ alignments[str(nid)] = edits
+ return {
+ "id": block_id,
+ "consensus": consensus,
+ "alignments": alignments,
+ }
+
+
+def build_sequence_pangraph_json():
+ """Minimal pangraph for testing junction sequence extraction.
+
+ Two strains (circular), sharing the same junction in opposite orientations:
+ s1: C1+(n1) A1+(n2) C2+(n3) ← forward orientation
+ s2: C2-(n4) A2-(n5) C1-(n6) ← inverted orientation
+
+ Block consensuses use distinct sequences for verification:
+ C1 (bid=10, 6bp): "AAACCC" - backbone (len >= L_thr=4)
+ C2 (bid=20, 6bp): "GGGAAA" - backbone
+ A1 (bid=30, 3bp): "TTC" - accessory (s1 only, len < L_thr=4)
+ A2 (bid=40, 3bp): "ACG" - accessory (s2 only, len < L_thr=4)
+
+ A1 and A2 are distinct strain-private accessory blocks with different consensus
+ sequences, so neither junction center block is `core`. Node 6 (C1 in s2) has a
+ substitution: pos=0 A->T, giving "TAACCC", which makes the isolates' flanking
+ sequences distinguishable too.
+ """
+ nodes = {
+ # s1
+ "1": _make_node(1, 10, 0, True, 0, 6),
+ "2": _make_node(2, 30, 0, True, 6, 9),
+ "3": _make_node(3, 20, 0, True, 9, 15),
+ # s2
+ "4": _make_node(4, 20, 1, False, 0, 6),
+ "5": _make_node(5, 40, 1, False, 6, 9),
+ "6": _make_node(6, 10, 1, False, 9, 15),
+ }
+
+ blocks = {
+ "10": _make_block_with_edits(
+ 10,
+ "AAACCC",
+ {
+ 1: None,
+ 6: {"subs": [{"pos": 0, "alt": "T"}], "dels": [], "inss": []},
+ },
+ ),
+ "20": _make_block_with_edits(
+ 20,
+ "GGGAAA",
+ {
+ 3: None,
+ 4: None,
+ },
+ ),
+ "30": _make_block_with_edits(
+ 30,
+ "TTC",
+ {
+ 2: None,
+ },
+ ),
+ "40": _make_block_with_edits(
+ 40,
+ "ACG",
+ {
+ 5: None,
+ },
+ ),
+ }
+
+ paths = {
+ "0": _make_path(0, "s1", [1, 2, 3], 15),
+ "1": _make_path(1, "s2", [4, 5, 6], 15),
+ }
+
+ return {"paths": paths, "blocks": blocks, "nodes": nodes}
+
+
+# Consensus sequences for build_inversion_pangraph_json, keyed by block-id string.
+# Distinct and non-palindromic so junction sequence extraction (incl. reverse-complement)
+# can be verified. Lengths are variable and <20bp; the seven core blocks (10-70) are >10bp.
+INVERSION_CONS = {
+ "10": "ACGTTGCAACCA", # C1, 12bp
+ "20": "TTGGAACCGGTTAC", # C2, 14bp
+ "30": "GATTACAGGCT", # C3, 11bp
+ "40": "CCAGTACGTGACATCA", # C4, 16bp
+ "50": "ACACGTGTACGTA", # C5, 13bp
+ "60": "TGTCATGCAATGCAT", # C6, 15bp
+ "70": "GGATCCGAATTCAGTCA", # C7, 17bp
+ "80": "ACGTGA", # A1, 6bp (accessory)
+ "90": "TTCAGGCA", # A2, 8bp (accessory)
+}
+
+
+def build_inversion_pangraph_json():
+ """Circular pangraph exercising reverse-complement, a single inversion, and mergers.
+
+ Three strains (all circular). C5/C6/C7 sit next to and co-oriented with C1/C2/C3, so
+ each pair always travels together and forms a merger (MSU):
+ s1 (reference): C1+ C5+ A1+ C2+ C6+ C3+ C7+ A2- C4+
+ s2 (whole-genome RC): C4- A2+ C7- C3- C6- C2- A1- C5- C1-
+ s3 (single inversion): C1+ C5+ C2+ C6+ C7- C3- C4+ (no accessory)
+
+ s2 is the exact reverse complement of s1 (order reversed, every strand flipped), so it
+ carries the same backbone edges as s1 (Edge equality is orientation-canonical) and must
+ co-orient back onto s1. s3 inverts the contiguous `C3 C7` segment: the internal C3-C7
+ adjacency is preserved (so {C3,C7} still merge) while C6-C3 and C7-C4 are broken, leaving
+ {C3,C7} as a separate, invertible MSU. A2 is on the `-` strand in s1 so the junction
+ sequence extraction exercises the center reverse-complement path. See INVERSION_CONS for
+ block consensuses (variable <20bp lengths; core blocks >10bp; tested at L_thr=10).
+
+ Core blocks (once per strain): C1=10, C2=20, C3=30, C4=40, C5=50, C6=60, C7=70.
+ Accessory blocks (s1+s2 only): A1=80, A2=90.
+ """
+ nodes = {
+ # s1: C1+ C5+ A1+ C2+ C6+ C3+ C7+ A2- C4+
+ "1": _make_node(1, 10, 0, True, 0, 12),
+ "2": _make_node(2, 50, 0, True, 12, 25),
+ "3": _make_node(3, 80, 0, True, 25, 31),
+ "4": _make_node(4, 20, 0, True, 31, 45),
+ "5": _make_node(5, 60, 0, True, 45, 60),
+ "6": _make_node(6, 30, 0, True, 60, 71),
+ "7": _make_node(7, 70, 0, True, 71, 88),
+ "8": _make_node(8, 90, 0, False, 88, 96),
+ "9": _make_node(9, 40, 0, True, 96, 112),
+ # s2: C4- A2+ C7- C3- C6- C2- A1- C5- C1-
+ "10": _make_node(10, 40, 1, False, 0, 16),
+ "11": _make_node(11, 90, 1, True, 16, 24),
+ "12": _make_node(12, 70, 1, False, 24, 41),
+ "13": _make_node(13, 30, 1, False, 41, 52),
+ "14": _make_node(14, 60, 1, False, 52, 67),
+ "15": _make_node(15, 20, 1, False, 67, 81),
+ "16": _make_node(16, 80, 1, False, 81, 87),
+ "17": _make_node(17, 50, 1, False, 87, 100),
+ "18": _make_node(18, 10, 1, False, 100, 112),
+ # s3: C1+ C5+ C2+ C6+ C7- C3- C4+
+ "19": _make_node(19, 10, 2, True, 0, 12),
+ "20": _make_node(20, 50, 2, True, 12, 25),
+ "21": _make_node(21, 20, 2, True, 25, 39),
+ "22": _make_node(22, 60, 2, True, 39, 54),
+ "23": _make_node(23, 70, 2, False, 54, 71),
+ "24": _make_node(24, 30, 2, False, 71, 82),
+ "25": _make_node(25, 40, 2, True, 82, 98),
+ }
+
+ block_nodes = {
+ "10": [1, 18, 19],
+ "20": [4, 15, 21],
+ "30": [6, 13, 24],
+ "40": [9, 10, 25],
+ "50": [2, 17, 20],
+ "60": [5, 14, 22],
+ "70": [7, 12, 23],
+ "80": [3, 16],
+ "90": [8, 11],
+ }
+ blocks = {
+ bid: _make_block_with_edits(
+ int(bid), INVERSION_CONS[bid], {nid: None for nid in nids}
+ )
+ for bid, nids in block_nodes.items()
+ }
+
+ paths = {
+ "0": _make_path(0, "s1", [1, 2, 3, 4, 5, 6, 7, 8, 9], 112),
+ "1": _make_path(1, "s2", [10, 11, 12, 13, 14, 15, 16, 17, 18], 112),
+ "2": _make_path(2, "s3", [19, 20, 21, 22, 23, 24, 25], 98),
+ }
+
+ return {"paths": paths, "blocks": blocks, "nodes": nodes}
+
+
+@pytest.fixture
+def sequence_pangraph():
+ """A minimal Pangraph for testing junction sequence extraction."""
+ return pp.Pangraph(build_sequence_pangraph_json())
+
+
+@pytest.fixture
+def junction_pangraph():
+ """A synthetic Pangraph with core/accessory blocks and non-trivial junctions."""
+ return pp.Pangraph(build_junction_pangraph_json())
+
+
+@pytest.fixture
+def tandem_accessory_pangraph():
+ """A Pangraph with an accessory block duplicated on one path inside a junction."""
+ return pp.Pangraph(build_tandem_accessory_pangraph_json())
+
+
+@pytest.fixture
+def linear_pangraph():
+ """A synthetic Pangraph with linear (non-circular) paths."""
+ return pp.Pangraph(build_linear_pangraph_json())
+
+
+@pytest.fixture
+def inversion_pangraph():
+ """A circular Pangraph with a reverse-complemented strain, a single inversion, and mergers."""
+ return pp.Pangraph(build_inversion_pangraph_json())
+
+
+@pytest.fixture
+def plasmid_pangraph():
+ """The real plasmids dataset, used by smoke tests."""
+ return pp.Pangraph.from_json("tests/data/plasmids.json")
diff --git a/packages/pypangraph/tests/data/staph.json.gz b/packages/pypangraph/tests/data/staph.json.gz
new file mode 100644
index 00000000..b21aaa48
Binary files /dev/null and b/packages/pypangraph/tests/data/staph.json.gz differ
diff --git a/packages/pypangraph/tests/test_block.py b/packages/pypangraph/tests/test_block.py
index ce1b969a..6fc1c4ff 100644
--- a/packages/pypangraph/tests/test_block.py
+++ b/packages/pypangraph/tests/test_block.py
@@ -1,5 +1,4 @@
import pytest
-import pypangraph as pp
from pypangraph.class_block import Block
from pypangraph.class_alignments import (
Alignment,
diff --git a/packages/pypangraph/tests/test_export.py b/packages/pypangraph/tests/test_export.py
new file mode 100644
index 00000000..36f4acbf
--- /dev/null
+++ b/packages/pypangraph/tests/test_export.py
@@ -0,0 +1,161 @@
+"""Tests for the GFA writer and the streamlined junction-aware GFA export."""
+
+import math
+import re
+
+import pytest
+
+from pypangraph.export import GFA, junction_context_gfa
+from pypangraph.junctions import BackboneJunctions
+
+ACCESSORY_RE = re.compile(r"^J\d+__\d+$")
+
+# The dominant backbone shared by s1 and s2 (C1-C2-C3-C4, circular); s3 carries
+# the rearranged C1-C3-C2-C4 instead. See conftest.build_junction_pangraph_json.
+DOMINANT_EDGES = {
+ "100_f__200_f",
+ "200_f__300_f",
+ "300_f__400_f",
+ "100_r__400_r",
+}
+
+
+def test_gfa_write_minimal(tmp_path):
+ """GFA.write emits H/S/L lines with LN/DP tags and correct orientations."""
+ segments = {"a": 100, "b": 50}
+ links = {("a", True, "b", False)}
+ depths = {"a": 3}
+ out = tmp_path / "tiny.gfa"
+ GFA(segments, links, depths).write(out)
+
+ lines = out.read_text().splitlines()
+ assert lines[0] == "H\tVN:Z:1.0"
+ s_lines = {ln.split("\t")[1]: ln for ln in lines if ln.startswith("S")}
+ assert "LN:i:100" in s_lines["a"] and "DP:f:3" in s_lines["a"]
+ assert "LN:i:50" in s_lines["b"] and "DP:f:" not in s_lines["b"]
+ (link,) = [ln for ln in lines if ln.startswith("L")]
+ assert link == "L\ta\t+\tb\t-\t0M"
+
+
+def _segment_lengths(gfa_path):
+ """Parse a GFA file into a dict of segment name -> emitted LN:i: value."""
+ lengths = {}
+ for ln in gfa_path.read_text().splitlines():
+ if not ln.startswith("S"):
+ continue
+ fields = ln.split("\t")
+ (tag,) = [f for f in fields if f.startswith("LN:i:")]
+ lengths[fields[1]] = int(tag.removeprefix("LN:i:"))
+ return lengths
+
+
+def test_gfa_write_length_transform(tmp_path):
+ """length_transform rescales LN:i:, rounds to nearest int, clamps to >= 1."""
+ gfa = GFA({"a": 100, "b": 1000, "c": 5}, set(), {})
+
+ # Linear rescale: round(l / 10) -> exact for 100 and 1000; 5 -> round(0.5)=0
+ # clamps up to the 1 floor.
+ out = tmp_path / "scaled.gfa"
+ gfa.write(out, length_transform=lambda length: length / 10)
+ assert _segment_lengths(out) == {"a": 10, "b": 100, "c": 1}
+
+ # log-based scaling also rounds to the nearest integer.
+ out_log = tmp_path / "log.gfa"
+ gfa.write(out_log, length_transform=lambda length: 100 * math.log(length))
+ assert _segment_lengths(out_log) == {
+ name: max(1, round(100 * math.log(length)))
+ for name, length in gfa.segments.items()
+ }
+
+ # Default path is unaffected: real lengths.
+ out_real = tmp_path / "real.gfa"
+ gfa.write(out_real)
+ assert _segment_lengths(out_real) == {"a": 100, "b": 1000, "c": 5}
+
+
+def test_consensus_gfa_structure(junction_pangraph):
+ """Consensus export: core anchors single & un-prefixed, accessory prefixed,
+ links well-formed, and depths reflect isolate coverage."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ gfa, prefix_map = junction_context_gfa(bj, scaffold="consensus")
+ segments, links, depths = gfa.segments, gfa.links, gfa.depths
+
+ core_ids = {"100", "200", "300", "400"}
+ # Every core block appears once, un-prefixed.
+ assert core_ids <= set(segments)
+ # Non-core segments are all per-junction accessory copies.
+ for name in segments:
+ if name not in core_ids:
+ assert ACCESSORY_RE.match(name), name
+
+ # Links only reference existing segments.
+ for a, _, b, _ in links:
+ assert a in segments and b in segments
+
+ # Core anchors are present in every isolate, so all share the same depth (3)
+ # regardless of which rearranged edges the consensus scaffold drops.
+ for cid in core_ids:
+ assert depths[cid] == 3
+ # A1 (id 500) sits in the C1-C2 junction of s1 and s2 -> depth 2.
+ a1 = next(n for n in segments if n.endswith("__500"))
+ assert depths[a1] == 2
+
+ # Every kept junction maps to a real edge of the graph.
+ for edge_str in prefix_map.values():
+ assert edge_str in bj
+
+
+def test_all_scaffold_is_superset_of_consensus(junction_pangraph):
+ """scaffold='all' keeps every junction, a superset of the consensus scaffold."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ gfa_c, pmap_c = junction_context_gfa(bj, scaffold="consensus")
+ gfa_a, pmap_a = junction_context_gfa(bj, scaffold="all")
+
+ edges_c = set(pmap_c.values())
+ edges_a = set(pmap_a.values())
+ assert edges_c < edges_a # strict superset (s3 rearrangement edges added)
+ assert len(gfa_a.links) >= len(gfa_c.links)
+
+
+def test_consensus_scaffold_follows_dominant_synteny(junction_pangraph):
+ """Per-edge majority keeps the s1/s2 backbone, dropping the s3 rearrangement."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ _, prefix_map = junction_context_gfa(bj, scaffold="consensus")
+ assert set(prefix_map.values()) == DOMINANT_EDGES
+
+
+def test_reference_scaffold_uses_that_genomes_edges(junction_pangraph):
+ """scaffold= keeps that genome's own core edges (s3 is rearranged)."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ _, prefix_map = junction_context_gfa(bj, scaffold="s3")
+ edges = set(prefix_map.values())
+ # s3 visits C1-C3-C2-C4, so its core edges differ from the consensus backbone.
+ assert edges != DOMINANT_EDGES
+ assert "100_f__300_f" in edges
+
+
+def test_unknown_scaffold_raises(junction_pangraph):
+ """A scaffold that is neither 'consensus'/'all' nor a known isolate raises."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ with pytest.raises(ValueError, match="unknown scaffold isolate"):
+ junction_context_gfa(bj, scaffold="not_a_genome")
+
+
+def test_accessory_depth_counts_repeats_on_same_path(tandem_accessory_pangraph):
+ """A block duplicated within a junction on the same path counts each copy.
+
+ s1 traverses accessory block 500 twice and s2 once, so the segment depth is 3
+ (total traversals), not 2 (distinct isolates carrying it)."""
+ bj = BackboneJunctions(tandem_accessory_pangraph, L_thr=500)
+ gfa, _ = junction_context_gfa(bj, scaffold="all")
+ a = next(n for n in gfa.segments if n.endswith("__500"))
+ assert gfa.depths[a] == 3
+
+
+def test_junction_context_gfa_runs_on_real_graph(plasmid_pangraph):
+ """Smoke test on a real plasmid graph: non-empty, all links valid."""
+ bj = BackboneJunctions(plasmid_pangraph, L_thr=500)
+ gfa, prefix_map = junction_context_gfa(bj, scaffold="consensus")
+ assert gfa.segments and gfa.links and prefix_map
+ for a, _, b, _ in gfa.links:
+ assert a in gfa.segments and b in gfa.segments
diff --git a/packages/pypangraph/tests/test_graph.py b/packages/pypangraph/tests/test_graph.py
index 81081e64..e9a82c09 100644
--- a/packages/pypangraph/tests/test_graph.py
+++ b/packages/pypangraph/tests/test_graph.py
@@ -20,6 +20,30 @@ def test_load_graph_gz():
assert len(pan.strains()) == 15
+def test_load_graph_invalid_extension(tmp_path):
+ fname = tmp_path / "plasmids.txt"
+ fname.write_text("{}")
+
+ with pytest.raises(pp.PangraphLoadError, match=r"\.json or \.json\.gz"):
+ pp.Pangraph.from_json(fname)
+
+
+def test_load_graph_invalid_json(tmp_path):
+ fname = tmp_path / "broken.json"
+ fname.write_text("{this is not valid json}")
+
+ with pytest.raises(pp.PangraphLoadError, match="failed to load pangraph"):
+ pp.Pangraph.from_json(fname)
+
+
+def test_load_graph_invalid_schema(tmp_path):
+ fname = tmp_path / "invalid.json"
+ fname.write_text("{}")
+
+ with pytest.raises(pp.PangraphLoadError, match="invalid pangraph JSON"):
+ pp.Pangraph.from_json(fname)
+
+
def test_paths(graph):
path = graph.paths["RCS48_p1"]
assert len(path) == 60
@@ -39,6 +63,26 @@ def test_blockstats_df(graph):
assert df["duplicated"].sum() == 10
+def test_block_ids_are_strings(graph):
+ """Block ids are stored as strings internally, and the blockstats index is str.
+
+ Block ids are u64 hashes that overflow int64; storing them as strings keeps them
+ out of any float64 coercion (which would silently corrupt values above 2**53).
+ """
+ bdf = graph.to_blockstats_df()
+ assert bdf.index.dtype == object
+ assert all(isinstance(bid, str) for bid in bdf.index)
+ assert all(isinstance(bid, str) for bid in graph.blocks.keys())
+
+
+def test_blocks_accessor_accepts_int_or_str(graph):
+ """pan.blocks[bid] resolves the same block whether bid is passed as str or int."""
+ bid = graph.blocks.keys()[0] # str
+ assert isinstance(bid, str)
+ assert graph.blocks[int(bid)] is graph.blocks[bid]
+ assert int(bid) in graph.blocks and bid in graph.blocks
+
+
def test_blockcount_df(graph):
df = graph.to_blockcount_df()
assert df.shape[0] == 137
@@ -54,7 +98,7 @@ def test_nodes_to_blocks(graph):
assert len(S) == len(nodes)
b, s = graph.nodes.node_to_block(8533989107945450583)
- assert b == 14710008249239879492
+ assert b == "14710008249239879492" # block ids are strings internally
assert s
@@ -64,6 +108,12 @@ def test_core_genome_alignment(graph):
assert core_aln.get_alignment_length() == 64989
+def test_core_genome_alignment_invalid_guide_strain(graph):
+ """An unknown guide strain is rejected with a ValueError (not a bare assert)."""
+ with pytest.raises(ValueError, match="Guide strain .* not found"):
+ graph.core_genome_alignment(guide_strain="does_not_exist")
+
+
def test_pairwise_accessory_genome_comparisons(graph):
ddf = graph.pairwise_accessory_genome_comparison()
assert ddf.shape == (225, 2)
diff --git a/packages/pypangraph/tests/test_junctions_analyses.py b/packages/pypangraph/tests/test_junctions_analyses.py
new file mode 100644
index 00000000..c6ad29aa
--- /dev/null
+++ b/packages/pypangraph/tests/test_junctions_analyses.py
@@ -0,0 +1,462 @@
+"""Tests for the BackboneJunctions analyses: positions(), sequences(), stats()."""
+
+import pandas as pd
+
+from pypangraph.junctions import BackboneJunctions
+
+
+# --- positions tests ---
+
+
+def test_junction_positions_forward_strand(junction_pangraph):
+ """Verify positions for forward-strand junctions in s1.
+
+ s1 path: C1+(0,1000) A1+(1000,1200) A2+(1200,1350) C2+(1350,2150) C3+(2150,2750) C4+(2750,3450)
+ """
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ pos = bj.positions()
+
+ # Edge C1+→C2+ in s1: both flanks on forward strand → strand=True
+ row = pos.loc[("100_f__200_f", "s1")]
+ assert row["strand"]
+ assert row["left_start"] == 0
+ assert row["left_end"] == 1000
+ assert row["right_start"] == 1350
+ assert row["right_end"] == 2150
+
+ # Edge C2+→C3+ in s1: empty junction, strand=True
+ row = pos.loc[("200_f__300_f", "s1")]
+ assert row["strand"]
+ assert row["left_start"] == 1350
+ assert row["left_end"] == 2150
+ assert row["right_start"] == 2150
+ assert row["right_end"] == 2750
+
+
+def test_junction_positions_inverted_edge(junction_pangraph):
+ """Verify positions for the universal C4→C1 edge (inverted in all strains).
+
+ The canonical edge is "100_r__400_r" = Edge(C1-, C4-). But in all strains
+ C1 and C4 are on + strand, so same_strand=False → positions are swapped.
+
+ s1: C4+(2750,3450) ... C1+(0,1000) (wraps around)
+ """
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ pos = bj.positions()
+
+ # In s1: junction is inverted, left=C4, right=C1
+ row = pos.loc[("100_r__400_r", "s1")]
+ assert not row["strand"]
+ assert row["left_start"] == 2750 # C4 start
+ assert row["left_end"] == 3450 # C4 end
+ assert row["right_start"] == 0 # C1 start
+ assert row["right_end"] == 1000 # C1 end
+
+
+def test_junction_positions_rearranged_strain(junction_pangraph):
+ """Verify positions for s3 which has a rearrangement (C2/C3 swapped).
+
+ s3 path: C1+(0,1000) A2-(1000,1150) C3+(1150,1750) C2+(1750,2550) A3+(2550,2850) C4+(2850,3550)
+ """
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ pos = bj.positions()
+
+ # Edge C1+→C3+: forward in s3
+ row = pos.loc[("100_f__300_f", "s3")]
+ assert row["strand"]
+ assert row["left_start"] == 0 # C1
+ assert row["left_end"] == 1000
+ assert row["right_start"] == 1150 # C3
+ assert row["right_end"] == 1750
+
+ # Edge "200_r__300_r" = Edge(C2-, C3-): in s3 both are on + strand → inverted
+ # Swapped: left=C3(1150,1750), right=C2(1750,2550)
+ row = pos.loc[("200_r__300_r", "s3")]
+ assert not row["strand"]
+ assert row["left_start"] == 1150 # C3
+ assert row["left_end"] == 1750
+ assert row["right_start"] == 1750 # C2
+ assert row["right_end"] == 2550
+
+ # Edge C2+→C4+: forward in s3, A3 between them
+ row = pos.loc[("200_f__400_f", "s3")]
+ assert row["strand"]
+ assert row["left_start"] == 1750 # C2
+ assert row["left_end"] == 2550
+ assert row["right_start"] == 2850 # C4
+ assert row["right_end"] == 3550
+
+
+def test_junction_positions_shape(junction_pangraph):
+ """Position DataFrame has one row per (isolate, edge) with non-NaN junction length."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ pos = bj.positions()
+
+ # one position row per (iso, edge) carrying that junction
+ n_junctions = bj.stats()["n_isolates"].sum()
+ assert len(pos) == n_junctions
+ assert list(pos.columns) == [
+ "left_start",
+ "left_end",
+ "right_start",
+ "right_end",
+ "strand",
+ ]
+
+
+def test_junction_positions_linear(linear_pangraph):
+ """Positions work correctly for linear paths.
+
+ s1: c5+(0,200) C1+(200,1200) A2+(1200,1350) C2+(1350,2150) C3+(2150,2750)
+ s2: C1+(0,1000) A3+(1000,1300) C2+(1300,2100) C3+(2100,2700) c5+(2700,2900)
+ """
+ bj = BackboneJunctions(linear_pangraph, L_thr=500)
+ pos = bj.positions()
+
+ # terminal junctions (leading c5 in s1, trailing c5 in s2) have no flanking
+ # edge, so only the two edge-bearing junctions appear
+ assert set(pos.index.get_level_values("edge")) == {"100_f__200_f", "200_f__300_f"}
+
+ # s1: C1→C2 has A2 in between
+ row = pos.loc[("100_f__200_f", "s1")]
+ assert row["left_start"] == 200 # C1
+ assert row["left_end"] == 1200
+ assert row["right_start"] == 1350 # C2
+ assert row["right_end"] == 2150
+
+ # s2: C1→C2 has A3 in between
+ row = pos.loc[("100_f__200_f", "s2")]
+ assert row["left_start"] == 0 # C1
+ assert row["left_end"] == 1000
+ assert row["right_start"] == 1300 # C2
+ assert row["right_end"] == 2100
+
+
+def test_junction_positions_smoke(plasmid_pangraph):
+ """positions() runs on the real plasmids dataset without errors."""
+ bj = BackboneJunctions(plasmid_pangraph, L_thr=500)
+ pos = bj.positions()
+
+ assert isinstance(pos, pd.DataFrame)
+ n_junctions = bj.stats()["n_isolates"].sum()
+ assert len(pos) == n_junctions
+ assert list(pos.columns) == [
+ "left_start",
+ "left_end",
+ "right_start",
+ "right_end",
+ "strand",
+ ]
+
+
+# --- sequence extraction tests ---
+
+
+def test_sequences_forward_junction(sequence_pangraph):
+ """Forward-strand junction produces correct concatenated sequence.
+
+ s1: C1+(n1) A1+(n2) C2+(n3) — forward orientation for edge 10_f__20_f
+ C1 consensus: AAACCC (no edits for n1)
+ A1 consensus: TTC (no edits for n2)
+ C2 consensus: GGGAAA (no edits for n3)
+ Expected: AAACCC + TTC + GGGAAA = AAACCCTTCGGGAAA
+ """
+ bj = BackboneJunctions(sequence_pangraph, L_thr=4)
+ records = bj.sequences("10_f__20_f")
+ by_iso = {r.id: str(r.seq) for r in records}
+ assert by_iso["s1"] == "AAACCCTTCGGGAAA"
+
+
+def test_sequences_inverted_junction(sequence_pangraph):
+ """Inverted junction is co-oriented with the forward one.
+
+ s2: C2-(n4) A2-(n5) C1-(n6) — inverted for edge 10_f__20_f
+ After inversion: C1+(n6) A2+(n5) C2+(n4)
+ C1 n6 has sub pos=0 A->T, so seq = TAACCC
+ A2 n5: ACG (no edits)
+ C2 n4: GGGAAA (no edits)
+ Expected: TAACCC + ACG + GGGAAA = TAACCCACGGGGAAA
+ """
+ bj = BackboneJunctions(sequence_pangraph, L_thr=4)
+ records = bj.sequences("10_f__20_f")
+ by_iso = {r.id: str(r.seq) for r in records}
+ assert by_iso["s2"] == "TAACCCACGGGGAAA"
+
+
+def test_sequences_co_orientation(sequence_pangraph):
+ """Both isolates' sequences start with C1 and end with C2 (co-oriented)."""
+ bj = BackboneJunctions(sequence_pangraph, L_thr=4)
+ records = bj.sequences("10_f__20_f")
+ assert len(records) == 2
+ for r in records:
+ seq = str(r.seq)
+ # All sequences should start with C1 block (6bp) and end with C2 block (6bp)
+ assert len(seq) == 15 # 6 + 3 + 6
+
+
+def test_sequences_record_metadata(sequence_pangraph):
+ """SeqRecord id is isolate name, description is edge string."""
+ bj = BackboneJunctions(sequence_pangraph, L_thr=4)
+ records = bj.sequences("10_f__20_f")
+ ids = {r.id for r in records}
+ assert ids == {"s1", "s2"}
+ for r in records:
+ assert r.description == "10_f__20_f"
+
+
+def test_sequences_reverse_complement_in_center(junction_pangraph):
+ """Accessory blocks on reverse strand are reverse-complemented.
+
+ s3 has edge 100_f__300_f: Junction(C1+(n13), [A2-(n14)], C3+(n15))
+ A2 is on reverse strand. Consensus is "A"*150, so rc = "T"*150.
+ C1 and C3 are forward: "A"*1000 and "A"*600.
+ Expected: "A"*1000 + "T"*150 + "A"*600
+ """
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ records = bj.sequences("100_f__300_f")
+ assert len(records) == 1
+ seq = str(records[0].seq)
+ assert records[0].id == "s3"
+ assert len(seq) == 1000 + 150 + 600
+ assert seq[:1000] == "A" * 1000
+ assert seq[1000:1150] == "T" * 150 # A2 reverse-complemented
+ assert seq[1150:] == "A" * 600
+
+
+def test_sequences_empty_junction(junction_pangraph):
+ """Empty junction (no accessory blocks) returns just the two core block sequences.
+
+ s1 edge 200_f__300_f: Junction(C2+, [], C3+) — no accessory blocks.
+ Expected: "A"*800 + "A"*600 = "A"*1400
+ """
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ records = bj.sequences("200_f__300_f")
+ by_iso = {r.id: str(r.seq) for r in records}
+ assert len(by_iso["s1"]) == 800 + 600
+ assert by_iso["s1"] == "A" * 1400
+
+
+def test_sequences_nonexistent_edge(junction_pangraph):
+ """Requesting sequences for a nonexistent edge returns empty list."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ assert bj.sequences("999_f__888_f") == []
+
+
+def test_sequences_smoke(plasmid_pangraph):
+ """sequences() runs on the real plasmids dataset without errors."""
+ bj = BackboneJunctions(plasmid_pangraph, L_thr=500)
+ stats_df = bj.stats()
+
+ # Test on the most frequent edge
+ top_edge = stats_df.index[0]
+ top_freq = stats_df.loc[top_edge, "n_isolates"]
+ records = bj.sequences(top_edge)
+ assert len(records) == top_freq
+ for r in records:
+ assert len(r.seq) > 0
+ assert r.description == top_edge
+
+
+# --- junction stats tests ---
+
+
+def test_junction_stats_values(junction_pangraph):
+ """Verify all per-edge statistics against expected values.
+
+ See plan for derivation of expected values.
+ """
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ sdf = bj.stats()
+
+ expected = {
+ "100_r__400_r": {
+ "n_isolates": 3,
+ "n_non_empty": 0,
+ "n_categories": 1,
+ "n_majority_category": 3,
+ "left_core_length": 1000,
+ "right_core_length": 700,
+ "accessory_length": 0,
+ },
+ "100_f__200_f": {
+ "n_isolates": 2,
+ "n_non_empty": 2,
+ "n_categories": 2,
+ "n_majority_category": 1,
+ "left_core_length": 1000,
+ "right_core_length": 800,
+ "accessory_length": 350,
+ },
+ "200_f__300_f": {
+ "n_isolates": 2,
+ "n_non_empty": 1,
+ "n_categories": 2,
+ "n_majority_category": 1,
+ "left_core_length": 800,
+ "right_core_length": 600,
+ "accessory_length": 300,
+ },
+ "300_f__400_f": {
+ "n_isolates": 2,
+ "n_non_empty": 0,
+ "n_categories": 1,
+ "n_majority_category": 2,
+ "left_core_length": 600,
+ "right_core_length": 700,
+ "accessory_length": 0,
+ },
+ "100_f__300_f": {
+ "n_isolates": 1,
+ "n_non_empty": 1,
+ "n_categories": 1,
+ "n_majority_category": 1,
+ "left_core_length": 1000,
+ "right_core_length": 600,
+ "accessory_length": 150,
+ },
+ "200_f__400_f": {
+ "n_isolates": 1,
+ "n_non_empty": 1,
+ "n_categories": 1,
+ "n_majority_category": 1,
+ "left_core_length": 800,
+ "right_core_length": 700,
+ "accessory_length": 300,
+ },
+ "200_r__300_r": {
+ "n_isolates": 1,
+ "n_non_empty": 0,
+ "n_categories": 1,
+ "n_majority_category": 1,
+ "left_core_length": 800,
+ "right_core_length": 600,
+ "accessory_length": 0,
+ },
+ }
+
+ assert set(sdf.index) == set(expected.keys())
+ for edge, vals in expected.items():
+ for col, val in vals.items():
+ assert sdf.loc[edge, col] == val, (
+ f"{edge}.{col}: {sdf.loc[edge, col]} != {val}"
+ )
+
+
+def test_junction_stats_transitive_and_singleton(junction_pangraph):
+ """Verify boolean flags for transitivity and singleton detection."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ sdf = bj.stats()
+
+ # Transitive: only one path category
+ transitive_edges = {
+ "100_r__400_r",
+ "300_f__400_f",
+ "100_f__300_f",
+ "200_f__400_f",
+ "200_r__300_r",
+ }
+ non_transitive = {"100_f__200_f", "200_f__300_f"}
+ for edge in transitive_edges:
+ assert sdf.loc[edge, "is_transitive"], f"{edge} should be transitive"
+ for edge in non_transitive:
+ assert not sdf.loc[edge, "is_transitive"], f"{edge} should not be transitive"
+
+ # Singleton: all but one isolate share the same path
+ singleton_edges = {"100_f__200_f", "200_f__300_f"}
+ non_singleton = set(sdf.index) - singleton_edges
+ for edge in singleton_edges:
+ assert sdf.loc[edge, "is_singleton"], f"{edge} should be singleton"
+ for edge in non_singleton:
+ assert not sdf.loc[edge, "is_singleton"], f"{edge} should not be singleton"
+
+
+def test_junction_stats_sorted_by_n_isolates(junction_pangraph):
+ """Stats DataFrame is sorted by `n_isolates` descending."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ sdf = bj.stats()
+ freqs = sdf["n_isolates"].values
+ assert all(freqs[i] >= freqs[i + 1] for i in range(len(freqs) - 1))
+
+
+def test_junction_stats_linear(linear_pangraph):
+ """Stats for linear paths fixture.
+
+ Two edges, both shared by s1 and s2:
+ 100_f__200_f: s1 center=[A2+], s2 center=[A3+] → 2 categories, singleton
+ 200_f__300_f: both empty → 1 category, transitive
+ """
+ bj = BackboneJunctions(linear_pangraph, L_thr=500)
+ sdf = bj.stats()
+
+ assert len(sdf) == 2
+ # terminal junctions (no flanking edge) are excluded
+ assert set(sdf.index) == {"100_f__200_f", "200_f__300_f"}
+
+ # 100_f__200_f: two different non-empty centers
+ row = sdf.loc["100_f__200_f"]
+ assert row["n_isolates"] == 2
+ assert row["n_non_empty"] == 2
+ assert row["n_categories"] == 2
+ assert row["n_majority_category"] == 1
+ assert not row["is_transitive"]
+ assert row["is_singleton"]
+ # Unique blocks: A2(bid=600, 150bp) + A3(bid=700, 300bp)
+ assert row["accessory_length"] == 150 + 300
+
+ # 200_f__300_f: both empty
+ row = sdf.loc["200_f__300_f"]
+ assert row["n_isolates"] == 2
+ assert row["n_non_empty"] == 0
+ assert row["n_categories"] == 1
+ assert row["n_majority_category"] == 2
+ assert row["is_transitive"]
+ assert not row["is_singleton"]
+ assert row["accessory_length"] == 0
+
+
+def test_backbone_stats_columns(junction_pangraph):
+ """stats() returns a DataFrame with the documented column set."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ stats_df = bj.stats()
+
+ assert isinstance(stats_df, pd.DataFrame)
+ assert list(stats_df.columns) == [
+ "n_isolates",
+ "n_non_empty",
+ "n_categories",
+ "n_majority_category",
+ "is_transitive",
+ "is_singleton",
+ "left_core_length",
+ "right_core_length",
+ "accessory_length",
+ ]
+ assert len(stats_df) == 7
+
+
+def test_junction_stats_smoke(plasmid_pangraph):
+ """junction stats on the real plasmids dataset — check invariants."""
+ bj = BackboneJunctions(plasmid_pangraph, L_thr=500)
+ sdf = bj.stats()
+
+ assert isinstance(sdf, pd.DataFrame)
+ assert len(sdf) > 0
+
+ # All isolate counts positive
+ assert (sdf["n_isolates"] > 0).all()
+ # Non-empty count is bounded by `n_isolates`
+ assert (sdf["n_non_empty"] >= 0).all()
+ assert (sdf["n_non_empty"] <= sdf["n_isolates"]).all()
+ # At least one category per edge
+ assert (sdf["n_categories"] >= 1).all()
+ # Majority can't exceed `n_isolates`
+ assert (sdf["n_majority_category"] <= sdf["n_isolates"]).all()
+ # Transitive iff exactly one category
+ assert (sdf["is_transitive"] == (sdf["n_categories"] == 1)).all()
+ # Singleton requires `n_isolates` > 1
+ assert not (sdf["is_singleton"] & (sdf["n_isolates"] <= 1)).any()
+ # Core lengths positive
+ assert (sdf["left_core_length"] > 0).all()
+ assert (sdf["right_core_length"] > 0).all()
+ # Accessory length non-negative
+ assert (sdf["accessory_length"] >= 0).all()
diff --git a/packages/pypangraph/tests/test_junctions_core.py b/packages/pypangraph/tests/test_junctions_core.py
new file mode 100644
index 00000000..968c415f
--- /dev/null
+++ b/packages/pypangraph/tests/test_junctions_core.py
@@ -0,0 +1,328 @@
+"""Tests for the junctions data structures: path splitting, Junction/JunctionNode
+primitives, canonical orientation, and the BackboneJunctions interface."""
+
+import pytest
+
+import pypangraph.topology_utils as tu
+from pypangraph.junctions import BackboneJunctions
+from pypangraph.junctions.junction import Junction, JunctionNode, path_junction_split
+
+
+# --- path_junction_split tests ---
+
+
+def test_path_junction_split(junction_pangraph):
+ """Splitting s1 at core blocks produces 4 junctions with correct flanking edges."""
+ pan = junction_pangraph
+ bdf = pan.to_blockstats_df()
+ walks = tu.pangraph_to_walks(pan)
+
+ def is_core(bid):
+ return (bdf.loc[bid, "len"] >= 500) and bdf.loc[bid, "core"]
+
+ junctions = path_junction_split(walks["s1"], is_core)
+ assert len(junctions) == 4
+
+ edges = [j.flanking_edge().to_str_id() for j in junctions]
+ assert set(edges) == {
+ "100_r__400_r",
+ "100_f__200_f",
+ "200_f__300_f",
+ "300_f__400_f",
+ }
+
+ # The junction between C1 and C2 should contain 2 accessory blocks (A1, A2)
+ j_c1_c2 = [j for j in junctions if j.flanking_edge().to_str_id() == "100_f__200_f"][
+ 0
+ ]
+ assert len(j_c1_c2.center) == 2
+
+
+def test_path_junction_split_rearranged(junction_pangraph):
+ """Splitting s3 (rearranged) produces 4 junctions with different edges than s1/s2."""
+ pan = junction_pangraph
+ bdf = pan.to_blockstats_df()
+ walks = tu.pangraph_to_walks(pan)
+
+ def is_core(bid):
+ return (bdf.loc[bid, "len"] >= 500) and bdf.loc[bid, "core"]
+
+ junctions = path_junction_split(walks["s3"], is_core)
+ assert len(junctions) == 4
+
+ edges = {j.flanking_edge().to_str_id() for j in junctions}
+ assert edges == {"100_r__400_r", "100_f__300_f", "200_r__300_r", "200_f__400_f"}
+
+
+def test_path_junction_split_requires_two_core_blocks():
+ """path_junction_split errors when a path has fewer than 2 core blocks.
+
+ A circular path with no core blocks previously crashed with IndexError; a path
+ with a single core block cannot form a junction either.
+ """
+ # no core blocks (circular) — would otherwise hit IndexError on junctions[0]
+ p0 = tu.Walk([tu.OrientedBlock(1, True), tu.OrientedBlock(2, True)], circular=True)
+ with pytest.raises(ValueError, match="at least 2"):
+ path_junction_split(p0, lambda bid: False)
+
+ # exactly one core block is still insufficient
+ p1 = tu.Walk([tu.OrientedBlock(1, True), tu.OrientedBlock(2, True)], circular=True)
+ with pytest.raises(ValueError, match="at least 2"):
+ path_junction_split(p1, lambda bid: bid == 1)
+
+
+def test_path_junction_split_linear(linear_pangraph):
+ """Splitting linear paths produces terminal junctions with None flanks.
+
+ s1: c5+ C1+ A2+ C2+ C3+
+ s2: C1+ A3+ C2+ C3+ c5+
+ """
+ pan = linear_pangraph
+ bdf = pan.to_blockstats_df()
+ walks = tu.pangraph_to_walks(pan)
+
+ def is_core(bid):
+ return (bdf.loc[bid, "len"] >= 500) and bdf.loc[bid, "core"]
+
+ # s1: c5+ C1+ A2+ C2+ C3+
+ # Expected: [None|c5+|C1+], [C1+||C2+], [C2+||C3+], no trailing (C3 is last)
+ # Wait - trailing after C3 is empty, but left_node=C3 so we get [C3+||None]
+ junctions_s1 = path_junction_split(walks["s1"], is_core)
+
+ # Leading terminal: left=None, center=[c5+], right=C1+
+ assert junctions_s1[0].left is None
+ assert len(junctions_s1[0].center) == 1
+ assert junctions_s1[0].right is not None
+
+ # Internal junctions have both flanks
+ assert junctions_s1[1].left is not None
+ assert junctions_s1[1].right is not None
+
+ # Trailing terminal: left=C3+, center=[], right=None
+ assert junctions_s1[-1].left is not None
+ assert junctions_s1[-1].right is None
+
+ # flanking_edge is None for terminal junctions
+ assert junctions_s1[0].flanking_edge() is None
+ assert junctions_s1[-1].flanking_edge() is None
+ assert junctions_s1[1].flanking_edge() is not None
+
+ # s2: C1+ A3+ C2+ C3+ c5+
+ junctions_s2 = path_junction_split(walks["s2"], is_core)
+
+ # No leading terminal (C1 is first)
+ assert junctions_s2[0].left is None
+ assert len(junctions_s2[0].center) == 0
+ assert junctions_s2[0].right is not None
+
+ # Trailing terminal: left=C3+, center=[c5+], right=None
+ assert junctions_s2[-1].left is not None
+ assert len(junctions_s2[-1].center) == 1
+ assert junctions_s2[-1].right is None
+
+
+# --- JunctionNode / Junction primitives ---
+
+
+def test_junction_node_inherits_equality():
+ """JunctionNode equality uses block_id + strand only (ignores node_id)."""
+ a = JunctionNode(100, True, 1)
+ b = JunctionNode(100, True, 2)
+ c = JunctionNode(100, False, 1)
+ assert a == b
+ assert a != c
+ assert hash(a) == hash(b)
+
+
+def test_junction_node_invert():
+ """JunctionNode.invert() flips strand but preserves node_id."""
+ n = JunctionNode(100, True, 42)
+ inv = n.invert()
+ assert inv.id == 100
+ assert inv.strand is False
+ assert inv.node_id == 42
+ assert isinstance(inv, JunctionNode)
+
+
+def test_junction_oriented_blocks():
+ """Junction.oriented_blocks() flattens left + center + right, skipping None flanks."""
+ left = tu.OrientedBlock(100, True)
+ a1 = tu.OrientedBlock(200, True)
+ a2 = tu.OrientedBlock(300, False)
+ right = tu.OrientedBlock(400, True)
+
+ # full junction: left + non-empty center + right
+ full = Junction(left, tu.Walk([a1, a2]), right)
+ assert full.oriented_blocks() == [left, a1, a2, right]
+
+ # empty center: just the two flanks
+ empty_center = Junction(left, tu.Walk([]), right)
+ assert empty_center.oriented_blocks() == [left, right]
+
+ # terminal junction with no left flank
+ no_left = Junction(None, tu.Walk([a1, a2]), right)
+ assert no_left.oriented_blocks() == [a1, a2, right]
+
+ # terminal junction with no right flank
+ no_right = Junction(left, tu.Walk([a1, a2]), None)
+ assert no_right.oriented_blocks() == [left, a1, a2]
+
+
+def test_junction_invert_terminal():
+ """invert() is robust to terminal junctions: a None flank stays None and the
+ present flank moves to the opposite side with its node_id preserved."""
+ a1 = tu.OrientedBlock(200, True)
+ a2 = tu.OrientedBlock(300, False)
+
+ # no left flank -> after inversion, no right flank
+ right = JunctionNode(400, True, 42)
+ no_left = Junction(None, tu.Walk([a1, a2]), right)
+ inv = no_left.invert()
+ assert inv.right is None
+ assert inv.left == right.invert()
+ assert inv.left.node_id == 42 # node_id preserved across inversion
+ assert inv.center == tu.Walk([a2.invert(), a1.invert()])
+ # double inversion round-trips
+ assert no_left.invert().invert() == no_left
+
+ # symmetric: no right flank -> after inversion, no left flank
+ left = JunctionNode(100, True, 7)
+ no_right = Junction(left, tu.Walk([a1, a2]), None)
+ rinv = no_right.invert()
+ assert rinv.left is None
+ assert rinv.right == left.invert()
+ assert rinv.right.node_id == 7
+
+
+# --- is_canonical / to_canonical ---
+
+
+def test_edge_is_canonical():
+ """Edge.is_canonical() is True iff (left, right) is the lex-min orientation."""
+ # natural form "100_f__200_f" is already lex-min
+ e = tu.Edge.from_str_id("100_f__200_f")
+ assert e.is_canonical()
+ assert e.to_str_id() == "100_f__200_f"
+
+ # inverting flips the orientation: now non-canonical
+ inv = e.invert()
+ assert not inv.is_canonical()
+ # but to_str_id() always returns the canonical form
+ assert inv.to_str_id() == "100_f__200_f"
+
+ # RC-palindromic edge: invert() yields the same string; tie resolves canonical
+ palindrome = tu.Edge(tu.OrientedBlock(100, True), tu.OrientedBlock(100, False))
+ assert palindrome.is_canonical()
+
+
+def test_str_id_roundtrip_underscore_ids():
+ """to_str_id / from_str_id round-trip block ids that themselves contain '_'
+ (e.g. MSU-renamed ids like 'MSU_0'); the strand suffix is peeled from the right."""
+ ob = tu.OrientedBlock("MSU_0", True)
+ assert tu.OrientedBlock.from_str_id(ob.to_str_id()) == ob
+
+ e = tu.Edge(tu.OrientedBlock("MSU_0", True), tu.OrientedBlock("MSU_12", False))
+ assert tu.Edge.from_str_id(e.to_str_id()) == e
+
+
+def test_junction_is_canonical(junction_pangraph):
+ """Junction.is_canonical() reports whether the junction's natural orientation
+ matches the lex-min canonical form of its flanking edge."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+
+ # s1 walks C1+ ... C2+: natural form "100_f__200_f" == canonical
+ assert bj["100_f__200_f"]["s1"].is_canonical()
+
+ # s3 walks C3+ then C2+ for the C2-C3 adjacency: natural form
+ # "300_f__200_f" loses lex-min to its invert "200_r__300_r"
+ assert not bj["200_r__300_r"]["s3"].is_canonical()
+
+
+def test_junction_to_canonical(junction_pangraph):
+ """to_canonical() returns self when already canonical, inverted otherwise;
+ the result always satisfies is_canonical()."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+
+ # canonical case: to_canonical() is a no-op
+ j_canon = bj["100_f__200_f"]["s1"]
+ assert j_canon.to_canonical() is j_canon
+
+ # non-canonical case: to_canonical() inverts and the result is canonical
+ j_inv = bj["200_r__300_r"]["s3"]
+ j_out = j_inv.to_canonical()
+ assert j_out is not j_inv
+ assert j_out.is_canonical()
+ assert j_out == j_inv.invert()
+
+
+def test_junction_canonical_terminal_raises():
+ """is_canonical()/to_canonical() are undefined for terminal junctions and raise,
+ so they cannot be applied accidentally to a junction with a missing flank."""
+ center = tu.Walk([tu.OrientedBlock(200, True)])
+
+ no_left = Junction(None, center, tu.OrientedBlock(400, True))
+ with pytest.raises(ValueError, match="terminal"):
+ no_left.is_canonical()
+ with pytest.raises(ValueError, match="terminal"):
+ no_left.to_canonical()
+
+ no_right = Junction(tu.OrientedBlock(100, True), center, None)
+ with pytest.raises(ValueError, match="terminal"):
+ no_right.is_canonical()
+ with pytest.raises(ValueError, match="terminal"):
+ no_right.to_canonical()
+
+
+# --- BackboneJunctions interface ---
+
+
+def test_junctions_edge_freq(junction_pangraph):
+ """Edge isolate count is correct and sorted descending."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ edge_freq = bj.stats()["n_isolates"]
+
+ # C4→C1 edge is universal (all 3 strains)
+ assert edge_freq["100_r__400_r"] == 3
+ # C1→C2, C2→C3, C3→C4 edges shared by s1 and s2
+ assert edge_freq["100_f__200_f"] == 2
+ assert edge_freq["200_f__300_f"] == 2
+ assert edge_freq["300_f__400_f"] == 2
+ # s3-private edges
+ assert edge_freq["100_f__300_f"] == 1
+ assert edge_freq["200_f__400_f"] == 1
+ assert edge_freq["200_r__300_r"] == 1
+
+ # sorted descending
+ counts = edge_freq.values
+ assert all(counts[i] >= counts[i + 1] for i in range(len(counts) - 1))
+
+
+def test_backbone_getitem(junction_pangraph):
+ """bj[edge_str] returns the {isolate -> Junction} mapping for an edge,
+ and `edge in bj` reports presence."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+
+ per_iso = bj["100_f__200_f"]
+ assert set(per_iso.keys()) == {"s1", "s2"}
+ j = per_iso["s1"]
+ assert j.flanking_edge().to_str_id() == "100_f__200_f"
+ assert len(j.center) == 2 # A1 + A2
+
+ assert "100_f__200_f" in bj
+ assert "nonexistent_edge" not in bj
+ with pytest.raises(KeyError):
+ _ = bj["nonexistent_edge"]
+
+
+def test_backbone_edges(junction_pangraph):
+ """edges() returns all distinct edge string IDs."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ assert set(bj.edges()) == {
+ "100_r__400_r",
+ "100_f__200_f",
+ "200_f__300_f",
+ "300_f__400_f",
+ "100_f__300_f",
+ "200_f__400_f",
+ "200_r__300_r",
+ }
diff --git a/packages/pypangraph/tests/test_junctions_inversion.py b/packages/pypangraph/tests/test_junctions_inversion.py
new file mode 100644
index 00000000..c727f78a
--- /dev/null
+++ b/packages/pypangraph/tests/test_junctions_inversion.py
@@ -0,0 +1,108 @@
+"""Cross-cutting junction tests on the inversion_pangraph fixture: a circular graph
+with a reverse-complemented strain and a single inversion, exercising edges,
+positions, stats, and sequences together."""
+
+from Bio.Seq import Seq
+
+import pypangraph.topology_utils as tu
+from pypangraph.junctions import BackboneJunctions
+
+
+def _edge(a, sa, b, sb):
+ """Canonical edge string id for OrientedBlock(a, sa) <-> OrientedBlock(b, sb)."""
+ return tu.Edge(tu.OrientedBlock(a, sa), tu.OrientedBlock(b, sb)).to_str_id()
+
+
+def test_inversion_edges_rc_and_private(inversion_pangraph):
+ """s1 and its whole-RC s2 share the same 7 backbone edges; s3's inversion adds 2 private ones."""
+ bj = BackboneJunctions(inversion_pangraph, L_thr=10)
+ freq = bj.stats()["n_isolates"].to_dict()
+
+ shared = {
+ _edge("10", True, "50", True), # C1-C5
+ _edge("50", True, "20", True), # C5-C2
+ _edge("20", True, "60", True), # C2-C6
+ _edge("60", True, "30", True), # C6-C3
+ _edge("30", True, "70", True), # C3-C7
+ _edge("70", True, "40", True), # C7-C4
+ _edge("40", True, "10", True), # C4-C1 (wrap)
+ }
+ private = {
+ _edge("60", True, "70", False), # C6 -> C7- (s3 inversion boundary)
+ _edge("30", False, "40", True), # C3- -> C4 (s3 inversion boundary)
+ }
+ assert set(bj.edges()) == shared | private
+
+ # 5 universal edges (freq 3), 2 broken by the inversion (freq 2), 2 private (freq 1)
+ assert sorted(freq.values(), reverse=True) == [3, 3, 3, 3, 3, 2, 2, 1, 1]
+ assert all(freq[e] == 1 for e in private)
+ assert freq[_edge("60", True, "30", True)] == 2 # C6-C3 broken by inversion
+ assert freq[_edge("70", True, "40", True)] == 2 # C7-C4 broken by inversion
+
+
+def test_inversion_positions_strand_mix(inversion_pangraph):
+ """A genome and its reverse complement disagree on canonical strand for every shared edge.
+
+ Unlike junction_pangraph (all core blocks forward, `strand` almost always True), s2 here is
+ the whole-genome RC of s1, so for every edge they share the is_canonical flag is opposite.
+ """
+ bj = BackboneJunctions(inversion_pangraph, L_thr=10)
+ pos = bj.positions()
+
+ shared = set(pos.xs("s1", level="iso").index) & set(pos.xs("s2", level="iso").index)
+ assert len(shared) == 7
+ for edge in shared:
+ assert pos.loc[(edge, "s1"), "strand"] != pos.loc[(edge, "s2"), "strand"]
+
+ # genuine mix of both orientations (not the near-constant column of junction_pangraph)
+ assert set(pos["strand"]) == {True, False}
+
+ # the inversion flips the canonical sense of the C3-C7 edge in s3 relative to the reference
+ c3c7 = _edge("30", True, "70", True)
+ assert pos.loc[(c3c7, "s1"), "strand"] != pos.loc[(c3c7, "s3"), "strand"]
+
+
+def test_inversion_stats(inversion_pangraph):
+ """Accessory-bearing edges: a variable center (C5-C2) vs a constant one (C7-C4)."""
+ bj = BackboneJunctions(inversion_pangraph, L_thr=10)
+ sdf = bj.stats()
+
+ # C5-C2: center [A1] in s1/s2, empty in s3 -> 2 categories, singleton
+ row = sdf.loc[_edge("50", True, "20", True)]
+ assert row["n_isolates"] == 3
+ assert row["n_categories"] == 2
+ assert row["n_majority_category"] == 2
+ assert not row["is_transitive"]
+ assert row["is_singleton"]
+ assert row["accessory_length"] == 6 # len(A1)
+
+ # C7-C4: center [A2] in both strains that carry it -> 1 category, transitive
+ row = sdf.loc[_edge("70", True, "40", True)]
+ assert row["n_isolates"] == 2
+ assert row["n_categories"] == 1
+ assert row["is_transitive"]
+ assert row["accessory_length"] == 8 # len(A2)
+
+
+def test_inversion_sequences_rc_equivalence(inversion_pangraph):
+ """Co-oriented junction sequences of a genome and its RC are identical; reverse-strand center
+ blocks are reverse-complemented."""
+ bj = BackboneJunctions(inversion_pangraph, L_thr=10)
+ # A1 (block 80) consensus, read back from the graph (the block has no per-node edits)
+ a1 = next(iter(inversion_pangraph.blocks["80"].to_sequences().values()))
+
+ # C5-C2 carries accessory A1 (reverse in the co-oriented frame) in s1/s2, empty in s3
+ seqs = {r.id: str(r.seq) for r in bj.sequences(_edge("50", True, "20", True))}
+ assert set(seqs) == {"s1", "s2", "s3"}
+ assert seqs["s1"] == seqs["s2"] # a genome equals its RC once co-oriented
+ assert len(seqs["s1"]) == 13 + 6 + 14 # C5 + A1 + C2
+ assert len(seqs["s3"]) == 13 + 14 # empty center
+ assert (
+ str(Seq(a1).reverse_complement()) in seqs["s1"]
+ ) # center block reverse-complemented
+ assert a1 not in seqs["s1"]
+
+ # C7-C4 (n_isolates=2, s1/s2 only) also co-orients to identical sequences
+ seqs2 = {r.id: str(r.seq) for r in bj.sequences(_edge("70", True, "40", True))}
+ assert set(seqs2) == {"s1", "s2"}
+ assert seqs2["s1"] == seqs2["s2"]
diff --git a/packages/pypangraph/tests/test_msu.py b/packages/pypangraph/tests/test_msu.py
index 1fc0f618..7dc85359 100644
--- a/packages/pypangraph/tests/test_msu.py
+++ b/packages/pypangraph/tests/test_msu.py
@@ -1,103 +1,125 @@
import pypangraph as pp
-import pypangraph.msu as msu
from collections import defaultdict
import pytest
+from pypangraph.minimal_synteny_units import (
+ minimal_synteny_units,
+ flip_msu_to_most_common_orientation,
+)
+import pypangraph.topology_utils as tu
+
class TestNode:
def test_eq(self):
- n1 = msu.Node("A", True)
- n2 = msu.Node("A", True)
- n3 = msu.Node("A", False)
- n4 = msu.Node("B", True)
+ n1 = tu.OrientedBlock("A", True)
+ n2 = tu.OrientedBlock("A", True)
+ n3 = tu.OrientedBlock("A", False)
+ n4 = tu.OrientedBlock("B", True)
assert n1 == n2
assert n1 != n3
assert n1 != n4
def test_invert(self):
- n = msu.Node("A", True)
- assert n.invert() == msu.Node("A", False)
+ n = tu.OrientedBlock("A", True)
+ assert n.invert() == tu.OrientedBlock("A", False)
class TestPath:
+ def test_constructor_default_nodes_are_not_shared(self):
+ p1 = tu.Walk()
+ p2 = tu.Walk()
+ p1.add_right(tu.OrientedBlock("A", True))
+ assert p1.oriented_blocks == [tu.OrientedBlock("A", True)]
+ assert p2.oriented_blocks == []
+
def test_eq(self):
- n1 = msu.Node("A", True)
- n2 = msu.Node("B", True)
- n3 = msu.Node("C", True)
- p1 = msu.Path([n1, n2, n3], circular=True)
- p2 = msu.Path([n1, n2, n3], circular=True)
- p3 = msu.Path([n1, ~n2, n3], circular=True)
- p4 = msu.Path([n1, n2], circular=True)
+ n1 = tu.OrientedBlock("A", True)
+ n2 = tu.OrientedBlock("B", True)
+ n3 = tu.OrientedBlock("C", True)
+ p1 = tu.Walk([n1, n2, n3], circular=True)
+ p2 = tu.Walk([n1, n2, n3], circular=True)
+ p3 = tu.Walk([n1, ~n2, n3], circular=True)
+ p4 = tu.Walk([n1, n2], circular=True)
assert p1 == p2
assert p1 != p3
assert p1 != p4
def test_invert(self):
- n1 = msu.Node("A", True)
- n2 = msu.Node("B", True)
- n3 = msu.Node("C", True)
- p = msu.Path([n1, n2, n3], circular=True)
- assert ~p == msu.Path([~n3, ~n2, ~n1], circular=True)
+ n1 = tu.OrientedBlock("A", True)
+ n2 = tu.OrientedBlock("B", True)
+ n3 = tu.OrientedBlock("C", True)
+ p = tu.Walk([n1, n2, n3], circular=True)
+ assert ~p == tu.Walk([~n3, ~n2, ~n1], circular=True)
def test_rotate_to(self):
- n1 = msu.Node("A", True)
- n2 = msu.Node("B", True)
- n3 = msu.Node("C", True)
- n4 = msu.Node("D", True)
- p = msu.Path([n1, n2, n3, n4], circular=True)
- assert p.rotate_to("B", True) == msu.Path([n2, n3, n4, n1], circular=True)
- assert p.rotate_to("D", True) == msu.Path([n4, n1, n2, n3], circular=True)
- assert p.rotate_to("B", False) == msu.Path([~n2, ~n1, ~n4, ~n3], circular=True)
- assert p.rotate_to("D", False) == msu.Path([~n4, ~n3, ~n2, ~n1], circular=True)
+ n1 = tu.OrientedBlock("A", True)
+ n2 = tu.OrientedBlock("B", True)
+ n3 = tu.OrientedBlock("C", True)
+ n4 = tu.OrientedBlock("D", True)
+ p = tu.Walk([n1, n2, n3, n4], circular=True)
+ assert p.rotate_to("B", True) == tu.Walk([n2, n3, n4, n1], circular=True)
+ assert p.rotate_to("D", True) == tu.Walk([n4, n1, n2, n3], circular=True)
+ assert p.rotate_to("B", False) == tu.Walk([~n2, ~n1, ~n4, ~n3], circular=True)
+ assert p.rotate_to("D", False) == tu.Walk([~n4, ~n3, ~n2, ~n1], circular=True)
def test_rename_bids(self):
- n1 = msu.Node("A", True)
- n2 = msu.Node("B", False)
- n3 = msu.Node("C", True)
- p = msu.Path([n1, n2, n3], circular=True)
- assert p.rename_bids({"A": "X", "B": "Y", "C": "Z"}) == msu.Path(
- [msu.Node("X", True), msu.Node("Y", False), msu.Node("Z", True)],
+ n1 = tu.OrientedBlock("A", True)
+ n2 = tu.OrientedBlock("B", False)
+ n3 = tu.OrientedBlock("C", True)
+ p = tu.Walk([n1, n2, n3], circular=True)
+ assert p.rename_bids({"A": "X", "B": "Y", "C": "Z"}) == tu.Walk(
+ [
+ tu.OrientedBlock("X", True),
+ tu.OrientedBlock("Y", False),
+ tu.OrientedBlock("Z", True),
+ ],
circular=True,
)
class TestEdge:
def test_invert(self):
- n1 = msu.Node("A", True)
- n2 = msu.Node("B", True)
- e = msu.Edge(n1, n2)
- assert e.invert() == msu.Edge(n2.invert(), n1.invert())
+ n1 = tu.OrientedBlock("A", True)
+ n2 = tu.OrientedBlock("B", True)
+ e = tu.Edge(n1, n2)
+ assert e.invert() == tu.Edge(n2.invert(), n1.invert())
def test_eq(self):
- n1 = msu.Node("A", True)
- n2 = msu.Node("B", True)
- e1 = msu.Edge(n1, n2)
- e2 = msu.Edge(n1, n2)
- e3 = msu.Edge(~n1, ~n2)
- e4 = msu.Edge(n1, ~n2)
- e5 = msu.Edge(~n2, ~n1)
+ n1 = tu.OrientedBlock("A", True)
+ n2 = tu.OrientedBlock("B", True)
+ e1 = tu.Edge(n1, n2)
+ e2 = tu.Edge(n1, n2)
+ e3 = tu.Edge(~n1, ~n2)
+ e4 = tu.Edge(n1, ~n2)
+ e5 = tu.Edge(~n2, ~n1)
assert e1 == e2
assert e1 != e3
assert e1 != e4
assert e1 == e5
+def test_new_import_paths_smoke():
+ assert callable(minimal_synteny_units)
+ assert tu.OrientedBlock("A", True).id == "A"
+ assert callable(pp.minimal_synteny_units)
+
+
@pytest.fixture
def generate_core_paths():
- A = msu.Node("A", True)
- B = msu.Node("B", True)
- C = msu.Node("C", True) # invert
- D = msu.Node("D", True) # invert
- E = msu.Node("E", True)
- F = msu.Node("F", True)
- G = msu.Node("G", True)
- H = msu.Node("H", True) # invert
- J = msu.Node("J", True)
-
- p1 = msu.Path([A, B, C, D, E, F, G, H, J], circular=True)
- p2 = msu.Path([A, B, C, D, E, F, G, H, J], circular=True)
- p3 = msu.Path([A, B, ~D, ~C, E, F, G, H, J], circular=True)
- p4 = msu.Path([A, B, ~D, ~C, E, F, G, ~H, J], circular=True)
+ A = tu.OrientedBlock("A", True)
+ B = tu.OrientedBlock("B", True)
+ C = tu.OrientedBlock("C", True) # invert
+ D = tu.OrientedBlock("D", True) # invert
+ E = tu.OrientedBlock("E", True)
+ F = tu.OrientedBlock("F", True)
+ G = tu.OrientedBlock("G", True)
+ H = tu.OrientedBlock("H", True) # invert
+ J = tu.OrientedBlock("J", True)
+
+ p1 = tu.Walk([A, B, C, D, E, F, G, H, J], circular=True)
+ p2 = tu.Walk([A, B, C, D, E, F, G, H, J], circular=True)
+ p3 = tu.Walk([A, B, ~D, ~C, E, F, G, H, J], circular=True)
+ p4 = tu.Walk([A, B, ~D, ~C, E, F, G, ~H, J], circular=True)
paths = {1: p1, 2: p2, 3: p3, 4: p4}
nodes = {n.id: n for n in [A, B, C, D, E, F, G, H, J]}
@@ -108,7 +130,7 @@ def generate_core_paths():
def test_find_mergers(generate_core_paths):
paths, nodes = generate_core_paths
- mg = msu.find_mergers(paths)
+ mg = tu.find_mergers(paths)
mg_groups = defaultdict(set)
for source, sink in mg.items():
mg_groups[sink].add(source)
@@ -126,7 +148,158 @@ def load_graph():
def test_msu(load_graph):
+ """Smoke test on the real plasmids graph with invariant checks."""
pan = load_graph
- MSU_mergers, MSU_paths, MSU_len = msu.minimal_synteny_units(
- pan, L_thr=50, rotate=True
+ MSU_mergers, MSU_paths, MSU_len = minimal_synteny_units(pan, L_thr=50, rotate=True)
+
+ # one MSU walk per strain
+ assert set(MSU_paths.keys()) == set(pan.strains())
+ # every block maps to a real MSU label
+ assert all(label in MSU_len for label in MSU_mergers.values())
+ assert all(label.startswith("MSU_") for label in MSU_mergers.values())
+ # MSU_0 is the longest unit (labels are assigned by descending length)
+ assert max(MSU_len, key=MSU_len.get) == "MSU_0"
+
+
+def test_minimal_synteny_units_circular(junction_pangraph):
+ """End-to-end MSU extraction on the circular junction fixture.
+
+ Core blocks: C1=100(1000bp), C2=200(800bp), C3=300(600bp), C4=400(700bp);
+ accessory blocks are dropped by the core filter. Core orders are
+ s1/s2 = C1 C2 C3 C4 and s3 = C1 C3 C2 C4. The only edge shared by all three
+ strains is the circular wrap C4->C1, so C1 and C4 merge into a single MSU
+ (sink=400=C4); C2 and C3 stay singletons. MSUs are labelled by descending
+ length: MSU_0={C1,C4} (1700bp), MSU_1=C2 (800bp), MSU_2=C3 (600bp). All paths
+ are rotated to MSU_0 (forward); everything is forward so the flip is a no-op,
+ and s3 keeps its C2/C3 rearrangement.
+ """
+ MSU_mergers, MSU_paths, MSU_len = minimal_synteny_units(
+ junction_pangraph, L_thr=500, rotate=True
+ )
+
+ assert MSU_len == {"MSU_0": 1700, "MSU_1": 800, "MSU_2": 600}
+ assert MSU_mergers == {
+ "100": "MSU_0",
+ "400": "MSU_0",
+ "200": "MSU_1",
+ "300": "MSU_2",
+ }
+
+ expected_s1 = tu.Walk(
+ [
+ tu.OrientedBlock("MSU_0", True),
+ tu.OrientedBlock("MSU_1", True),
+ tu.OrientedBlock("MSU_2", True),
+ ],
+ circular=True,
+ )
+ expected_s3 = tu.Walk(
+ [
+ tu.OrientedBlock("MSU_0", True),
+ tu.OrientedBlock("MSU_2", True),
+ tu.OrientedBlock("MSU_1", True),
+ ],
+ circular=True,
+ )
+ assert MSU_paths["s1"] == expected_s1
+ assert MSU_paths["s2"] == expected_s1
+ assert MSU_paths["s3"] == expected_s3
+ assert all(p.circular for p in MSU_paths.values())
+
+
+def test_minimal_synteny_units_no_rotate(linear_pangraph):
+ """MSU extraction with rotate=False on linear paths.
+
+ Core blocks C1=100(1000), C2=200(800), C3=300(600); both strains' core order is
+ C1 C2 C3. They form a single transitive chain, so all three merge into one MSU
+ (sink=100, 2400bp). Each path reduces to a single MSU node, unrotated.
+ """
+ MSU_mergers, MSU_paths, MSU_len = minimal_synteny_units(
+ linear_pangraph, L_thr=500, rotate=False
)
+
+ assert MSU_len == {"MSU_0": 2400}
+ assert MSU_mergers == {"100": "MSU_0", "200": "MSU_0", "300": "MSU_0"}
+
+ expected = tu.Walk([tu.OrientedBlock("MSU_0", True)], circular=False)
+ assert MSU_paths["s1"] == expected
+ assert MSU_paths["s2"] == expected
+ assert all(p.circular is False for p in MSU_paths.values())
+
+
+def test_minimal_synteny_units_rotate_requires_circular(linear_pangraph):
+ """rotate=True (the default) on linear paths raises a ValueError."""
+ with pytest.raises(ValueError, match="Only circular paths"):
+ minimal_synteny_units(linear_pangraph, L_thr=500)
+
+
+def test_minimal_synteny_units_rc_collapse(inversion_pangraph):
+ """A genome and its whole reverse-complement collapse onto the same canonical MSU walk.
+
+ Core blocks C1=10..C7=70; C5/C6/C7 are co-oriented neighbours of C1/C2/C3 and merge with
+ them. The universal adjacencies merge {C1,C2,C4,C5,C6} (sink 10); the s3 inversion breaks
+ C6-C3 and C7-C4, keeping {C3,C7} as a separate, invertible MSU (sink 30):
+ MSU_0 = {C1,C2,C4,C5,C6} (70bp), MSU_1 = {C3,C7} (28bp).
+ s1 and its RC s2 rotate/flip onto the identical walk (exercising rotate_to's invert branch);
+ s3's inversion survives as a flipped MSU_1 node.
+ """
+ MSU_mergers, MSU_paths, MSU_len = minimal_synteny_units(
+ inversion_pangraph, L_thr=10, rotate=True
+ )
+
+ assert MSU_len == {"MSU_0": 70, "MSU_1": 28}
+ assert MSU_mergers == {
+ "10": "MSU_0",
+ "50": "MSU_0",
+ "20": "MSU_0",
+ "60": "MSU_0",
+ "40": "MSU_0",
+ "70": "MSU_1",
+ "30": "MSU_1",
+ }
+
+ expected_ref = tu.Walk(
+ [tu.OrientedBlock("MSU_0", True), tu.OrientedBlock("MSU_1", True)],
+ circular=True,
+ )
+ expected_s3 = tu.Walk(
+ [tu.OrientedBlock("MSU_0", True), tu.OrientedBlock("MSU_1", False)],
+ circular=True,
+ )
+ assert MSU_paths["s1"] == expected_ref
+ assert MSU_paths["s2"] == expected_ref # RC collapses onto the reference
+ assert (
+ MSU_paths["s3"] == expected_s3
+ ) # the inversion survives as a flipped MSU node
+ assert all(p.circular for p in MSU_paths.values())
+
+
+def test_flip_msu_to_most_common_orientation():
+ """Blocks predominantly on the reverse strand are flipped to forward.
+
+ X is reverse in 2 of 3 walks (net negative) so every X occurrence is flipped;
+ Y is forward throughout (net positive) and left untouched. The dict passed in is
+ mutated in place and returned.
+ """
+ paths = {
+ "a": tu.Walk(
+ [tu.OrientedBlock("X", False), tu.OrientedBlock("Y", True)], circular=True
+ ),
+ "b": tu.Walk(
+ [tu.OrientedBlock("X", False), tu.OrientedBlock("Y", True)], circular=True
+ ),
+ "c": tu.Walk(
+ [tu.OrientedBlock("X", True), tu.OrientedBlock("Y", True)], circular=True
+ ),
+ }
+
+ result = flip_msu_to_most_common_orientation(paths)
+
+ assert result is paths # mutates in place and returns the same dict
+ # X net = (-1) + (-1) + (+1) = -1 < 0 -> every X occurrence flipped
+ assert paths["a"].oriented_blocks[0] == tu.OrientedBlock("X", True)
+ assert paths["b"].oriented_blocks[0] == tu.OrientedBlock("X", True)
+ assert paths["c"].oriented_blocks[0] == tu.OrientedBlock("X", False)
+ # Y net positive -> unchanged
+ assert paths["a"].oriented_blocks[1] == tu.OrientedBlock("Y", True)
+ assert paths["c"].oriented_blocks[1] == tu.OrientedBlock("Y", True)
diff --git a/packages/pypangraph/tests/test_plots.py b/packages/pypangraph/tests/test_plots.py
new file mode 100644
index 00000000..6a6e61d7
--- /dev/null
+++ b/packages/pypangraph/tests/test_plots.py
@@ -0,0 +1,17 @@
+import matplotlib
+
+matplotlib.use("Agg")
+
+import matplotlib.pyplot as plt
+
+from pypangraph.junctions import BackboneJunctions
+from pypangraph.plots import linear_junction_plot
+
+
+def test_linear_junction_plot_smoke(junction_pangraph):
+ """Verify linear_junction_plot runs end-to-end on the junction fixture."""
+ bj = BackboneJunctions(junction_pangraph, L_thr=500)
+ edge = next(iter(bj.edges()))
+ fig, ax = plt.subplots()
+ linear_junction_plot(ax, bj, edge)
+ plt.close(fig)