diff --git a/dev/design/annotate-block-compaction-open-questions.md b/dev/design/annotate-block-compaction-open-questions.md new file mode 100644 index 00000000..48cc8b00 --- /dev/null +++ b/dev/design/annotate-block-compaction-open-questions.md @@ -0,0 +1,111 @@ +# Block-compaction policy: open questions to settle before it's final + +> Scope: the **block-level compaction** layer only (`annotation/compact.rs`, +> `CoordinateConsensusStrategy`). The **node-level lift is not affected** by any of this — it is the +> lossless source of truth, and compaction is an opinionated summary built on top of it behind the +> `BlockCompactionStrategy` trait. Everything below can change (or a whole alternative strategy can be +> added) without touching the lift or the writers. Resolve these before P5.3 freezes the behaviour in +> user docs. Companion to [`annotate.md`](./annotate.md) §8/§11/§12 and the as-built notes in +> [`annotate-implementation-notes.md`](./annotate-implementation-notes.md). + +The layer already works and is regression-tested against the hard structural cases (repeated blocks, +inversions, origin-spanning nodes). What is still *provisional* is the **reconciliation policy**: how +strictly two genomes' placements must agree to be called "the same feature", and what the block-level +output should carry. These are decisions, not bugs. + +## 1. Coordinate-exact identity vs. a tolerance (the crux) + +> **Resolved (2026-07-06): ship exact-only for v1.** The five-species fragmentation study found +> coordinate fragmentation is ~1% and does not grow with divergence; a tolerance is a characterized, +> ~30-line fallback if ever needed. Details in +> [`annotate-implementation-notes.md`](./annotate-implementation-notes.md) → "Block-compaction +> open-question resolutions" → issue 1. + +**Now:** two per-genome instances cluster only when their whole crossing matches *bit-exactly* — +`(feature_type, [(block_id, cons_start, cons_end, strand_on_consensus), …])` (`ClusterKey`). One indel +near a feature boundary nudges a consensus endpoint by a few bp, which splits one biological gene into +several near-identical clusters, each scored (and possibly dropped) on its own smaller support. + +**Why it matters:** this is the single policy choice most likely to make the block view under-report +support on real, divergent genomes. The node-level table still holds everything, so nothing is *lost* +— but the compacted summary can look artificially fragmented. + +**Options:** +- **Exact only (status quo).** Predictable and simple; the risk above is real on diverse inputs. +- **Bounded coordinate tolerance.** Treat endpoints within ±δ bp as equal (cluster by rounded/binned + coordinates, or merge clusters whose termini are within δ). Needs a rule for transitivity/chaining + and a default δ. + +**Recommendation — measure before deciding.** On the klebs / *E. coli* graphs, run a throwaway script +that counts how many biological genes (grouped by GFF `ID`/name) fragment into **more than one** +block-level cluster *purely* because of boundary indels, and how many of those fragments fall below a +0.9 threshold that the whole gene would clear. If that count is small, ship exact-only for v1 and note +the limitation; if it is large, a tolerance is worth the added complexity. This is the concrete +real-data check called for by [`annotate.md`](./annotate.md) §8 ("whether to allow a small coordinate +**tolerance** … is an open refinement"). + +## 2. Per-segment support pools features that only *look* alike + +> **Resolved (2026-07-06): document as intended.** Empirically rare (≤20 placements pool >1 CDS id +> across the study's five datasets) and informational only (gating is per-crossing). Keep the key; +> document that per-segment support is a property of a placement, not a feature. + +**Now:** `n_support_segment` is tallied over `segment_support_sets`, keyed on +`(feature_type, block_id, cons_start, cons_end, strand)` with **no feature identity** +(`compact.rs`). Two genuinely different genes of the same type that happen to place a segment at +identical block-consensus coordinates are counted together. + +**Why it matters:** it can inflate the reported reproducibility of a segment. Rare in practice (needs +an exact coordinate coincidence between distinct features), but silent when it happens. + +**Decision:** either (a) document this as intended — "per-segment support is a property of a +placement, not of a feature" — or (b) add feature identity (e.g. the cluster key, or the base feature +id) to `SegmentKey` so only same-feature segments pool. (a) keeps the "shared body reports the same +support in every crossing" behaviour that motivated the current key; (b) is stricter. Pick one and +state it. + +## 3. No provenance / drill-down in the block output + +> **Resolved (2026-07-06): defer for v1.** Drill-down means going back to the node-level table; the +> nested-JSON supporters writer lands before the disagreement use case is advertised. P5.3 docs must +> state the node-table drill-down explicitly. + +**Now:** a block row reports `M`/`N` and a consensus name/attributes, but not *which* genomes support +or dissent. "Where do genomes disagree about this feature?" is one of the stated use cases, and the +CSV can't answer it. + +**Decision:** confirm the deferral is acceptable for v1, and land the nested-JSON writer (already +deferred in [`annotate.md`](./annotate.md) §8) that carries the per-genome supporters behind each +consensus row before the feature is advertised as covering the disagreement use case. Until then, note +in the docs that drill-down means going back to the node-level table. + +## 4. Genome identity is a display string, not a `PathId` + +> **Resolved (2026-07-06): keyed on `PathId` (implemented).** `LiftedAnnotation` gained a `path_id` +> field and compaction now keys its whole internal pipeline on it (the `genome` string is kept only +> for the node CSV). Removes the collision risk and the per-row `String` clones. + +**Now:** the whole compaction keys genomes on a `String` label — `path.name()` falling back to +`path.id().to_string()` — produced identically by the lift and by `genome_label` (so the `M`/`N` join +is correct today). + +**Why it matters:** it is a robustness/efficiency point rather than a policy one. Two paths with the +same name (or empty names) would silently merge into one genome, and the label is cloned per row at +genome scale. Keying the internal pipeline on `PathId` (resolving names only at output) removes the +collision risk and the allocation. Low priority; fold in only if the pipeline is touched for another +reason. + +## How to close these out + +1. Fix the float-rounding off-by-one in `min_count` — **done** (epsilon before `.ceil()`, with a + regression test). +2. Run the real-data fragmentation measurement (§1) — **done**. The five-species study + (`tmp/block-annotations/`, note `n00_block_compaction_fragmentation.typ`) decided **exact-only for + v1**; tolerance stays a ~30-line fallback in one strategy. +3. Make the calls on §2–§4 and record each in + [`annotate-implementation-notes.md`](./annotate-implementation-notes.md) — **done** (§2 document as + intended, §3 defer with node-table drill-down, §4 keyed on `PathId`, implemented). + +**All four resolved (2026-07-06).** The remaining work is P5.3: write the user docs against this +settled behaviour (exact-only clustering; per-segment support is placement-level; drill-down via the +node-level table). diff --git a/dev/design/annotate-implementation-notes.md b/dev/design/annotate-implementation-notes.md index 579c8b58..63515a56 100644 --- a/dev/design/annotate-implementation-notes.md +++ b/dev/design/annotate-implementation-notes.md @@ -14,8 +14,8 @@ | P2 | Inverse coordinate helper `consensus_coords_from_node` (+ round-trip tests) | ✅ done | | P3 | Node-level lift + `AnnotationWriter` trait (CSV impl) | ✅ done | | P5.1 | `pangraph annotate` CLI (node-level output, **exact** seqid matching) | ✅ done | -| P4 | Block-level compaction (coordinate agreement, CLI refinement level) + JSON writer | ⏳ todo (after P5.1) | -| P5.2 | Wire block-level output into the `annotate` command | ⏳ todo | +| P4 | Block-level compaction (coordinate agreement) + block CSV writer | ✅ done | +| P5.2 | Wire block-level output into the `annotate` command | ✅ done | | P5.3 | Docs page + CLI reference regeneration (both output levels) | ⏳ todo | | P6 | pypangraph consumer + visualization example | ⏳ todo | @@ -115,7 +115,9 @@ coordinate. `lift_features(grouped, graph)` runs it over the `match_features_to_ ### Per-endpoint flags use the **consensus-endpoint frame** (decision) `start_*` / `end_*` flags refer to the row's `cons_start` / `cons_end` endpoints, **not** the genome 5'/3' ends. On a reverse-strand node the feature's genome-start maps to `cons_end`; `strand_on_consensus` -(the source strand flipped on reverse nodes; `None` stays `None`) lets a consumer map back. +(the source strand flipped on reverse nodes; `None` stays `None`) lets a consumer map back. (A +companion `feature_strand` field — the *un-flipped* GFF genome strand — was added later by the +crossing-identity fix; see the P4 notes.) - `start_is_terminus`/`end_is_terminus` — the endpoint is a real feature terminus vs a fragment boundary from a node/block split. A fully-covered feature has exactly two real termini (its genome start and end); interior segment boundaries are non-termini. (GFF source-`partial` carry-through is @@ -157,8 +159,11 @@ wiring over the P1–P3 library API — no lift-logic changes. `annotate_run` do ### CLI surface (intentionally minimal) - **Graph**: positional `input` (`Option`, stdin if omitted), like `simplify`. -- **GFF(s)**: repeatable `--gff `, `required = true` (≥1); one path per occurrence; transparent - decompression by extension. +- **GFF(s)**: `--gff ...`, `required = true` (≥1), `num_args = 1..` — accepts several paths + after one flag (`--gff a.gff b.gff`, so shell globs `--gff *.gff` work) and/or the flag repeated + (`--gff a --gff b`); values accumulate into one `Vec`. Because the graph is the + positional arg, give it before the flag (or via stdin) so the variadic does not slurp it as an + extra GFF; a following flag such as `-o` terminates the list. Transparent decompression by extension. - **Output**: `-o/--output` (default `-` = stdout); CSV only; compression inferred from extension. - **No `--seqid-map`, no `--output-level`, no `--format`/`--delimiter`** — deferred to later phases (§12 / P4 / P5.2). The CSV delimiter is hard-wired to `,`. @@ -203,7 +208,247 @@ gene + its CDS); every other feature was already byte-exact (20,935 / 20,937). Re-running the (throwaway) klebs verification after this change: **20,937 / 20,937 OK, 0 skipped**. -## Public API introduced in P1–P3 +## P4 decisions & behaviours (block-level compaction) + +`annotation/compact.rs` adds the **block-level** view: collapsing the redundant node-level table into +one **consensus feature** per recurring placement. It is a configurable, **modular** layer built on +top of the node-level lift — never baked into it. + +### Modular strategy trait +`BlockCompactionStrategy::compact(&[LiftedAnnotation], &Pangraph) -> Vec` is the +single entry point. The one P4 impl is `CoordinateConsensusStrategy { min_frequency, +property_threshold }`; future name-/ortholog-based strategies implement the same trait. Provisional +config defaults live on `Default` (`min_frequency = 0.9`, `property_threshold = 0.5`) — the +user-facing defaults are owned by the CLI in **P5.2**. + +### Clustering = coordinate-consensus (decisions locked with the user) +> **Superseded by the crossing-identity fix (2026-06-16)** — the two-endpoint key and its +> mixed-orientation consequence below describe the *original* P4 design; the current behaviour is the +> whole-crossing multiset described in **"Crossing-identity fix"** further down. The bullets are kept +> for history. + +- **Cluster key** = the feature's two block-consensus **terminus** endpoints **+ `feature_type` + + `strand_on_consensus`**. A gene vs CDS, or opposite strands, at identical coordinates stay distinct + (the manifesto §8's coordinate-only key was extended with type + strand on the user's call). +- Each genome's feature instance is recovered by grouping node-level rows on + `(genome, feature_type, base id)` (`base id` = `parent_feature_id`, else the `feature_id` with its + `.seg{idx}` suffix stripped), then **reduced to exactly two terminus endpoints** via the + `*_is_terminus` flags. Instances not yielding exactly two termini (e.g. partial features) are + **excluded** from compaction — they remain in the node-level table. The two endpoints are stored in + **canonical sorted order** (`start` = lower `(block, coord)`, not necessarily the genome 5'); + `strand_on_consensus` disambiguates orientation. +- **Consequence:** a block traversed in **mixed orientations** across genomes splits the same + biological gene into a `+` and a `-` cluster (strand is part of the identity). Acceptable: the + consensus frame genuinely differs, and most blocks are traversed one way. + +### M-of-N threshold (frequency) +- `M` = number of distinct genomes sharing the exact key (a genome counts once even if duplicated). +- `N` = number of distinct genomes **traversing the cluster's block(s)** — the intersection of the + per-block path sets (just the one block for the common single-block case), via + `PangraphBlock::isolates`. *Not* the total number of genomes in the graph: a gene is not penalised + for being absent in genomes that lack the block entirely. +- A cluster is emitted when `M >= max(1, ceil(min_frequency * N))`. + +### Consensus metadata +Over the `M` supporters, `consensus_name` is the majority `name` and `consensus_attributes` are the +per-key majority attribute values, each emitted only when its support clears +`ceil(property_threshold * M)` (default 0.5). Ties break deterministically to the smaller string; +output is key-sorted. Resolved generically over **all** attribute keys, so `product` etc. fall out +for free. + +### Writer +The `AnnotationWriter` trait gained `write_block_annotations(&[BlockAnnotation])`, implemented on +`CsvAnnotationWriter` with a `BlockAnnotationCsvRow` mirroring the node-level row (`consensus_attributes` +as one JSON-string column, strand as `+`/`-`/empty, ids as numbers, headers on). **CSV only** this +phase; JSON / GFF-on-consensus writers are deferred. CLI wiring is **P5.2** (below). + +### Validation +Unit tests in `compact.rs` pin each behaviour (all-agree, below-threshold, type split, strand split, +name/attribute thresholds, multi-block endpoints, `N` = paths-traversing-block, partial-feature +exclusion). `tests/itest_annotate_compact.rs` lifts whole-core-block-node features across every +genome of `data/test_graph.json` and asserts they collapse to consensus feature(s) at `[0, L)` whose +supporters sum to all genomes, then round-trips the block CSV writer. + +### Crossing-identity fix (2026-06-16, supersedes the two-endpoint key) +Running `annotate blocks` on real *E. coli* graphs surfaced a bug on the gene `yjeM` (see the +investigation in the PR): a feature whose body sits on a core block and whose 37 bp tail sits on an +**inverted** block was split into separate `+` and `-` clusters even though the placement was +*identical* across genomes — and one genome was dropped below threshold. Root cause: the cluster key +sampled a single representative strand from `segment_idx 0` (the genome-lowest segment), and that +segment flips with the genome's global orientation across an inversion boundary. Strand was being +treated as one per-feature value when it is really **per-segment**. + +What changed: +- **Identity = the whole crossing.** The cluster key is now `(feature_type, Vec)` where + `Segment = (BlockId, cons_start, cons_end, strand_on_consensus)`, ordered **5'→3'**. Per-block + `strand_on_consensus` is invariant across genomes (the genome's `+/-` is absorbed by node + orientation), so the whole-crossing multiset is stable. This also distinguishes `X→Y→X` from `X→Y`. +- **Ordering** uses the new node-level `feature_strand` field (the original GFF genome strand, + un-flipped) to reverse the genome-ordered segments for reverse-strand features; unstranded features + (no reading direction) are canonicalized by the lexicographically smaller of the two orders. We + chose to **add `feature_strand`** rather than recompute `strand_on_consensus XOR node.strand()` at + compaction, because it keeps compaction independent of graph-consistent `node_id`s (its unit tests + build `LiftedAnnotation`s directly) and makes the node-level CSV self-document the GFF strand. +- **Output is one row per segment**, sharing a `cluster_id` and numbered 5'→3' by `segment_idx` + (`n_segments` total). A **duplicated block crossed twice** keeps both occurrences as separate rows, + distinguished by `segment_idx`. +- **`N` is structural capability** — genomes traversing every block in the crossing with the required + **multiplicity**, where multiplicity is the crossing's count of **distinct nodes** per block (not + segments), reduced to the element-wise minimum over the supporters so `N >= M` always holds. + Counting by node, not segment, matters: a whole-genome `region` feature crosses an **origin-spanning + node** that the lift splits into two coverage pieces — two segments on **one** node — so a + segment-count requirement of 2 against a genome with one such node gave `N = 0` while `M = 1` (an + `M > N` contract violation found on the E. coli data). `block_genome_counts` supplies per-genome + node counts; `reduce_instance` records per-block distinct-node use. A body-shared/tail-different + feature thus fragments into several clusters, each scored only against the genomes that could carry + it, so every variant is promoted at its true support. +- **Per-segment support** `n_support_segment / n_total_segment` is reported per row, computed over the + whole node-level table independently of clustering (`segment_support_sets`), so a body segment + shared by several crossings reports the same value in each. It is informational; gating stays + per-crossing. +- New unit tests: `test_inversion_crossing_collapses_to_one_cluster` (the regression), + `test_duplicated_block_crossed_twice`, `test_shared_body_three_tails_each_confident`; the existing + tests were updated to the per-segment shape and the block CSV header is now + `type,cluster_id,segment_idx,…`. + +### Feature-type filtering (`--only-type` / `--exclude-type`) +Motivated by the same E. coli run: whole-contig `region` records (one per genome) expanded to ~30k of +~65k block rows — correct output, but noise. Rather than special-case any type in the compactor, the +input is filtered by GFF `type` up front. Two **shared** `AnnotateCommonArgs` flags +(`packages/pangraph/src/commands/annotate/annotate_args.rs`), so both `nodes` and `blocks` inherit +them: `--only-type` (whitelist) and `--exclude-type` (blacklist), each `value_delimiter = ','` +(comma-separated and/or repeatable, accumulating into a `Vec`) and `conflicts_with` each other +(mutually exclusive). Matching is **exact and case-sensitive**. The filter — `filter_features_by_type` +in `annotation/feature.rs`, a pure `Vec -> Vec` — runs in `load_and_lift` **after** +GFF reading but **before** `match_features_to_paths`, so an excluded type never triggers a +seqid-mismatch error and never reaches the lift. Empty lists (the default) are a no-op, so existing +behaviour is unchanged. Unknown type names are silently no-ops (whitelisting an absent type yields an +empty result; excluding one keeps everything). + +## P5.2 decisions & behaviours (block output on the CLI) + +P5.2 exposes compaction on the `annotate` command. Following the `export` precedent (one command, +several output shapes as subcommands), `annotate` became a **subcommand group**: + +- `pangraph annotate nodes` — the P5.1 node-level table (unchanged behaviour). +- `pangraph annotate blocks` — runs `CoordinateConsensusStrategy::compact` then + `write_block_annotations`. + +### CLI surface (decisions locked with the user) +- **Subcommands, not a `--output-level` flag** — keeps block-only flags out of node-mode `--help`, + and matches `export`. Bare `annotate` requires a subcommand (no default mode; acceptable while + `annotate` is unreleased on the integration branch). +- **Names** = `nodes` / `blocks` (plural, consistent). +- **Shared options** (`input`, `--gff`, `-o/--output`) live in a flattened `AnnotateCommonArgs`; + both subcommands embed it via `#[clap(flatten)]` (DRY vs `export`'s per-variant duplication). +- **Block-only flags**: `--min-frequency` (default 0.9) and `--property-threshold` (default 0.5), + whose defaults are kept in sync with `CoordinateConsensusStrategy::default()`. Both are validated at + parse time by a `parse_fraction` `value_parser` rejecting anything outside `[0, 1]` (incl. `NaN`/ + `inf`) with a clap error, rather than silently emitting nothing / promoting everything downstream. + A `--strategy` selector is **deferred** until a second strategy exists — a one-option flag is noise, + and the trait already provides the modularity internally. + +### Wiring +`annotate_run` matches the `PangraphAnnotateArgs` enum and dispatches to `annotate_run_nodes` / +`annotate_run_blocks`. Both share `load_and_lift(&AnnotateCommonArgs) -> (Pangraph, Vec)` +(graph load → per-`--gff` `read_many` → `match_features_to_paths` → `lift_features`); the graph is +returned because block compaction needs it for per-cluster path totals. `tests/itest_annotate_cli.rs` +exercises both subcommands (`annotate blocks --min-frequency 0` emits the block-level header + rows) +and pins the threshold wiring: with a single annotated genome, `--min-frequency 1.0` drops the +single-support clusters that `0.0` keeps (a field swap would surface as the lenient run coming back +empty). + +## Block-compaction open-question resolutions + +Settling the provisional reconciliation policy enumerated in +[`annotate-block-compaction-open-questions.md`](./annotate-block-compaction-open-questions.md) before +P5.3 freezes the behaviour in user docs. Each subsection records the **as-built call** and its evidence. + +### Issue 1 — coordinate-exact identity vs. a tolerance → **ship exact-only for v1** (resolved 2026-07-06) + +The crux question (does exact-coordinate clustering under-report support by fragmenting one biological +gene into several near-identical clusters when a boundary indel nudges a consensus endpoint?) was +settled by the **real-data fragmentation study** the doc called for (`tmp/block-annotations/`, +gitignored Snakemake pipeline; write-up in `notes/n00_block_compaction_fragmentation.typ`). + +Method: for genes shared by ≥2 genomes (grouped by **protein identity** — CDS `parent_feature_id`, so +the same protein across genomes is one biological gene), count those that fragment into >1 exact +block-level cluster *purely* from boundary coordinate wobble (same blocks + strand, differing coords — +"coordinate fragmentation") vs. legitimately different blocks ("structural"); and of the genes that +would clear the 0.9 threshold as a single merged crossing, how many exact-only **drops** below it. The +analysis re-implements `CoordinateConsensusStrategy` in Python and was validated to reproduce the Rust +output cluster-for-cluster (10 368 == 10 368 on ecoli-15, 6 318 == 6 318 on saureus-50 at `--min-frequency 0`). + +Result across **five species spanning four phyla** (2.8–6.3 Mb; monomorphic → open pangenome) — +`coord-frag% / dropped@0.9%`: ecoli-15 0.8/1.0, saureus-50 1.1/1.4, klebsiella-50 0.8/0.9, +paeruginosa-30 1.1/1.1, mtb-30 1.0/0.3. Coordinate fragmentation is **~1% everywhere and does not grow +with divergence.** The reason is mechanistic: because genes are grouped by protein identity, divergence +surfaces as *different proteins* (different clusters, correctly) or *different blocks* (structural +fragmentation, legitimate), **not** as coordinate wobble of a fixed protein — the boundary-indel +mechanism is intrinsically rare. The Klebsiella open pangenome, named in the doc as the divergence +stress test, showed the *same* ~1%, not more. + +**Decision: exact-only for v1**, with the limitation documented (the node-level table remains the +lossless source of truth; the block summary can look marginally fragmented). A tolerance is a +characterized, cheap **fallback** if ever needed — δ ≈ 15–20 bp (single-linkage merge of clusters whose +termini agree within δ, *per structural group* so distinct crossings never chain) recovers ~65–95% of +coordinate fragments while wrongly merging ≤ 4 genuinely-distinct loci even at δ = 50. That is ~30 lines +in one strategy, exactly as the design anticipated; deferring it costs nothing because the node table +loses nothing. + +### Issue 2 — per-segment support pools look-alike placements → **document as intended** (resolved 2026-07-06) + +`n_support_segment` is tallied over `segment_support_sets`, keyed on +`(feature_type, block_id, cons_start, cons_end, strand_on_consensus)` with **no feature identity**, so +two genuinely different genes of the same type that place a segment at identical block-consensus +coordinates are counted together. The study's issue-2 side-check confirmed this is **rare** (≤ 20 +placements pool > 1 distinct CDS id across all five datasets, and those are overwhelmingly same-locus +variants, not distinct loci). The value is **informational** — gating is per-crossing, never per-segment +— so the coincidence cannot change which clusters are emitted, only an advisory column. + +**Decision: keep the key as-is** and document the semantics: *per-segment support is a property of a +placement, not of a feature.* This preserves the intended "a shared body reports the same support in +every crossing it appears in" behaviour that motivated the coordinate-only `SegmentKey`. The +rare-coincidence caveat is a doc line, not a code change. + +### Issue 3 — no per-genome provenance in the block output → **defer, document drill-down** (resolved 2026-07-06) + +A block row reports `M`/`N` and a consensus name/attributes but not *which* genomes support or dissent, +so the "where do genomes disagree about this feature?" use case can't be answered from the CSV alone. + +**Decision: acceptable for v1.** Drill-down means going back to the **node-level table** (the lossless +source of truth, joinable on `cluster_id`/block+coords). The nested-JSON writer carrying per-genome +supporters behind each consensus row (already deferred in [`annotate.md`](./annotate.md) §8) lands +**before** the feature is advertised as covering the disagreement use case. The P5.3 docs must state the +node-table drill-down explicitly rather than implying the block CSV answers disagreement. + +### Issue 4 — genome identity is a display string, not a `PathId` → **key on `PathId`** (resolved 2026-07-06) + +Compaction previously keyed every genome on a `String` label (`path.name()` → `path.id()` fallback), +so two paths sharing a name (or both unnamed) would silently merge into one genome, and the label was +cloned per row at genome scale. + +**Fix (as-built):** the node-level lift now records the stable path identity on each row, and +compaction keys its whole internal pipeline on it: +- `LiftedAnnotation` gained a `path_id: PathId` field (`lift.rs`), populated from `path.id()` alongside + the existing `genome` display name. `genome` is retained purely for the **node-level CSV** (the + output still shows the path name, resolved by the lift), so the node schema is unchanged — the CSV + writer projects an explicit `LiftedAnnotationCsvRow` that does not include `path_id`. +- In `compact.rs`, `FeatureInstance`, the supporter map (`ClusterAgg.supporters`), the feature-instance + grouping key, `block_genome_counts`, `crossing_n_total`, and `segment_support_sets` all key on + `PathId` instead of `String`. `PathId` is `Copy`, so this also drops the per-row `String` clones and + the `&String` borrows in the `N` intersection. `block_genome_counts` now takes the `PathId` straight + from `PangraphBlock::isolates` (no name lookup), and the `genome_label` helper was deleted. +- Compaction is therefore independent of the display name entirely; `M`/`N` join on identity. This + keeps the crossing-identity design's property that compaction does not depend on graph-consistent + `node_id`s — the unit tests still build `LiftedAnnotation`s directly, now passing a path index that + becomes both the `path_id` and the (irrelevant-to-logic) `genome` string. + +Behaviour is unchanged on well-formed graphs (distinct path names); the collision risk on +duplicate/empty names is removed. Full annotation unit + integration suite (incl. the klebs real-data +smoke) green. + +## Public API introduced in P1–P5.2 ```rust // annotation::feature @@ -228,14 +473,31 @@ pub fn consensus_coords_from_node_flagged( // annotation::lift (P3) pub struct LiftedAnnotation { /* feature_id, parent_feature_id, segment_idx, n_segments, genome, - block_id, node_id, strand_on_consensus, node_start/_end, cons_start/_end, start/end_is_terminus, - start/end_in_insertion, frac_covered, feature_type, name, attributes */ } + path_id, block_id, node_id, strand_on_consensus, node_start/_end, cons_start/_end, + start/end_is_terminus, start/end_in_insertion, frac_covered, feature_type, name, attributes */ } pub fn lift_feature(feature: &Feature, path: &PangraphPath, graph: &Pangraph) -> Result, Report>; pub fn lift_features(grouped: &BTreeMap>, graph: &Pangraph) -> Result, Report>; -// annotation::writer (P3) -pub trait AnnotationWriter { fn write_node_annotations(&mut self, annotations: &[LiftedAnnotation]) -> Result<(), Report>; } -pub struct CsvAnnotationWriter; // new(filepath, delimiter) ; the default CSV impl +// annotation::compact (P4) +pub struct BlockAnnotation { /* feature_type, strand_on_consensus, start_block_id, cons_start, + end_block_id, cons_end, consensus_name, consensus_attributes, n_support (M), n_total (N) */ } +pub trait BlockCompactionStrategy { + fn compact(&self, node_annotations: &[LiftedAnnotation], graph: &Pangraph) -> Result, Report>; +} +pub struct CoordinateConsensusStrategy { pub min_frequency: f64, pub property_threshold: f64 } // + Default + +// annotation::writer (P3 trait, P4 block method) +pub trait AnnotationWriter { + fn write_node_annotations(&mut self, annotations: &[LiftedAnnotation]) -> Result<(), Report>; + fn write_block_annotations(&mut self, annotations: &[BlockAnnotation]) -> Result<(), Report>; // P4 +} +pub struct CsvAnnotationWriter; // new(filepath, delimiter) ; the default CSV impl (node + block) + +// commands::annotate::annotate_args (P5.1 → P5.2) +pub enum PangraphAnnotateArgs { Nodes(PangraphAnnotateNodesArgs), Blocks(PangraphAnnotateBlocksArgs) } +pub struct AnnotateCommonArgs { input: Option, gff: Vec, output: PathBuf } +pub struct PangraphAnnotateNodesArgs { common: AnnotateCommonArgs } +pub struct PangraphAnnotateBlocksArgs { common: AnnotateCommonArgs, min_frequency: f64, property_threshold: f64 } ``` ## Real-data smoke tests & fixtures (P1.5) diff --git a/dev/design/annotate.md b/dev/design/annotate.md index 3d248642..01eb37b1 100644 --- a/dev/design/annotate.md +++ b/dev/design/annotate.md @@ -163,15 +163,25 @@ The hard part is largely solved by existing infrastructure: ## 7. Command interface -`pangraph annotate` (P5.1 as-built; the surface grows in later phases): +`pangraph annotate` (P5.2 as-built; the surface grows in later phases): +- **Mode = subcommand**: `annotate nodes` (node-level table) and `annotate blocks` (block-level + consensus features), mirroring `export`. Shared options below live in a flattened + `AnnotateCommonArgs` embedded by both (**as-built, P5.2**; P5.1 shipped only the node-level + command, with the selector deferred). - **Inputs**: a graph JSON as the **positional** argument (stdin if omitted) + one or more GFF files via a **repeatable `--gff`** flag. +- **Type filtering** (shared, both subcommands): `--only-type` (whitelist) / `--exclude-type` + (blacklist) restrict the input by GFF `type` column. Each takes a **comma-separated** list + (`--only-type gene,CDS`) and/or repeats; the two are **mutually exclusive**. Matching is **exact + and case-sensitive**, and the filter runs **before** seqid→path matching (so excluding a type also + sidesteps any seqid errors it would raise). Motivation: whole-contig `region` records (one per + genome) otherwise dominate the block-level output as giant clusters. - **Output**: `-o/--output` path (default `-` = stdout); compression inferred from the file - extension; format chosen by the `AnnotationWriter` (CSV is the only impl in P5.1; JSON arrives with - P4). -- **Mode**: P5.1 emits **node-level** only. The node-vs-block output-level selector arrives in P5.2, - once block-level output exists. + extension; format chosen by the `AnnotationWriter` (CSV is the only impl so far; JSON deferred). +- **Block tuning** (`annotate blocks` only): `--min-frequency` (default 0.9) and `--property-threshold` + (default 0.5), kept in sync with `CoordinateConsensusStrategy::default()`. A `--strategy` selector + is **deferred** to a future phase (only one strategy exists today). - **seqid → path matching**: P5.1 matches a GFF `seqid` to a path name by **exact** string equality. **ID matching is the #1 failure mode** — fail loudly (one error listing every offending seqid, not a silent drop). The `--seqid-map` escape hatch and any version-insensitive relaxation are @@ -190,7 +200,8 @@ Fields: ``` feature_id, parent_feature_id, segment_idx, genome (path name), block_id, node_id, -strand_on_consensus, +strand_on_consensus, # feature strand vs block consensus (flips with node orientation) +feature_strand, # original GFF genome strand (= strand_on_consensus XOR node orientation) node_start, node_end, # node-local coordinates cons_start, cons_end, # block-consensus coordinates start_is_terminus, end_is_terminus, # real terminus vs fragment boundary (per endpoint) @@ -207,42 +218,65 @@ them nested. ### Block-level annotation (`BlockAnnotation`) — compacted by coordinate consensus -One entry per **consensus feature**: a placement that recurs, at the *same* block-consensus -coordinates, across enough of the genomes that carry it. The node-level table stays the lossless -source of truth; this layer is an opinionated, configurable summary built **on top of** it, never -baked into the lift. +One entry per **segment** of a **consensus feature**: a crossing that recurs, at the *same* +block-consensus coordinates and per-block strands, across enough of the genomes that carry it. The +node-level table stays the lossless source of truth; this layer is an opinionated, configurable +summary built **on top of** it, never baked into the lift. -**Cluster key = the feature's block-consensus terminus endpoints.** For each genome, reduce a -feature's node-level lift to just its two real termini — the 5' start and the 3' end — as -`(block_id, cons_coord)` pairs (the endpoints whose `*_is_terminus` flag is set). The **internal** -fragment boundaries where a multi-segment feature crosses node/block boundaries are **ignored**: only -the outer start and end define identity. Two genomes carry "the same" feature when these endpoint -pairs match exactly: +**Cluster key = the feature's whole crossing, as a canonical ordered multiset.** For each genome, +reduce a feature's node-level lift to the 5'→3'-ordered list of its per-segment placements +`(block_id, cons_start, cons_end, strand_on_consensus)`. Two genomes carry "the same" feature when +their `(feature_type, crossing)` match **exactly** — every block, both coordinates, and the per-block +strand agree along the whole crossing (and, for a duplicated block crossed twice, the multiplicity +agrees too). A feature must be bounded by exactly two real termini (its 5' and 3' ends) to be +compacted; partial/truncated crossings stay only in the node-level table. ``` -key = (start_block_id, cons_start, end_block_id, cons_end) +key = (feature_type, [ (block_id, cons_start, cons_end, strand_on_consensus), … ]) # ordered 5'→3' ``` -For the common single-block feature `start_block_id == end_block_id`, so this is effectively **one -entry per (block, start, end)**. - -**M-of-N consensus.** Let N be the genomes contributing a candidate placement for a cluster (they -carry a feature whose termini fall on that block region) and M the number sharing the exact key. If -`M >= threshold`, the cluster is emitted as a single consensus annotation: a **consensus name** (the -majority feature name among the M supporters) and **only** the agreed block-consensus `cons_start` / -`cons_end` — the node-level detail is intentionally dropped. Clusters below threshold are not -promoted; they remain available in the node-level table. +**Why the whole crossing, and why per-block strand** (**as-built, the crossing-identity fix**): the +earlier design keyed only on the two *outer* terminus endpoints plus a single representative strand. +That single strand was sampled from the genome-lowest segment, which **flips with the genome's global +orientation** for a feature crossing an inversion boundary — so genomes with an *identical* placement +were split into `+` and `-` clusters (the yjeM bug). Per-block `strand_on_consensus` is invariant +across genomes (the genome's `+/-` is absorbed by node orientation), so making strand a *per-segment* +property of the whole crossing is both correct and stable. Keying on the whole crossing also +distinguishes `X→Y→X` from `X→Y`, which the two-endpoint key could not. + +**Ordering.** Node-level `segment_idx` is genome (low→high coordinate) order, already 5'→3' for a +forward feature. A reverse feature is reversed (using `feature_strand`); an **unstranded** feature +(no reading direction) is canonicalized by the lexicographically smaller of the two orders, so +homologous instances still collapse. + +**M-of-N consensus.** Let N be the genomes **structurally capable** of the crossing — those +traversing every block in it at least as many times as the crossing uses **distinct nodes** of that +block (the multiplicity-aware intersection of the per-block isolate sets, via +`PangraphBlock::isolates`; *not* the total genome count). Counting distinct nodes rather than +segments keeps an origin-spanning node — which the lift splits into two coverage pieces — from being +mistaken for a duplication. M is the genomes producing the exact crossing. If `M >= ceil(min_frequency · N)` (≥ 1) the cluster is emitted +as **one row per segment**, sharing a `cluster_id` and ordered 5'→3' by `segment_idx`. A feature whose +body is shared but whose tail differs fragments into several confident clusters (one per tail), each +scored only against the genomes that could carry it — so each variant is promoted at its true +support rather than penalised by the others. ``` -block_id, # single-block case; a multi-block feature carries both endpoints -consensus_name, # majority feature name across supporters -cons_start, cons_end, # the agreed block-consensus coordinates -n_support, # M — genomes sharing this exact start/end -n_total # N — genomes carrying a candidate placement here +feature_type, cluster_id, # cluster_id links the segments of one crossing +segment_idx, n_segments, # 5'→3' position within the feature, and total segments +block_id, cons_start, cons_end, # this segment's block and block-consensus coordinates +strand_on_consensus, # this segment's strand vs its block consensus +consensus_name, # majority feature name across supporters (>= property threshold) +consensus_attributes, # per-key majority attribute values (e.g. product) clearing threshold +n_support, n_total, # M / N — crossing-level (same on every row of the cluster) +n_support_segment, n_total_segment # M_seg / N_seg — this segment's reproducibility on its own block ``` -Default output is **long-format CSV**; an **optional nested JSON** can carry the per-genome -supporters behind each consensus row. +`n_support_segment / n_total_segment` are computed over the whole node-level table independently of +clustering, so a body segment shared by several crossings reports the same support in each. A +duplicated block crossed twice yields two rows, distinguished by `segment_idx`. + +Default output is **long-format CSV** (the only writer shipped so far); a nested-JSON writer carrying +the per-genome supporters behind each consensus row is **deferred** to a later phase. **Refinement level = the M-of-N threshold (CLI-selectable, tentative).** How strictly supporting genomes must agree before a coordinate is emitted is just the threshold `M` — an input choice, not a @@ -290,12 +324,21 @@ single documentation pass (P5.3) follow once both output levels exist. Execution **first end-to-end real-data exercise of the P3 lift** (P1.5 only matched features to paths, never lifted them). The output-level selector is **not** shipped yet — it arrives in P5.2 with block output. -- **P4** — block-level compaction into `BlockAnnotation` objects: cluster placements by their - block-consensus terminus endpoints and emit one **consensus feature** per cluster reaching the - **M-of-N threshold** (CLI-selectable; see §8), plus a second writer impl (CSV default, JSON - optional) over the same trait + tests. Designed against what the P5.1 prototype reveals. -- **P5.2** — wire the block-level output into the `annotate` command (the `block` output level over - the same args). +- **P4** ✅ — block-level compaction into `BlockAnnotation` objects via a modular + `BlockCompactionStrategy` trait (first impl: `CoordinateConsensusStrategy`): cluster placements by + their whole block-consensus crossing (+ type) and emit one row **per segment** for each cluster + reaching the **M-of-N threshold** (see §8), plus a **block CSV writer** over the same + `AnnotationWriter` trait + tests. A library layer only — the threshold/strategy CLI flags land in + P5.2; the JSON writer is deferred. *(Initially shipped a two-outer-endpoint key with a single + representative strand; **superseded by the crossing-identity fix** — see §8 — after the yjeM bug + showed that key splits identical placements across an inversion. The fix also added per-segment + output, `feature_strand` on the node-level lift, structural-capability `N`, and per-segment + support.)* +- **P5.2** ✅ — wire the block-level output into the `annotate` command. Implemented as **subcommands** + (`annotate nodes` / `annotate blocks`, mirroring `export`) rather than an `--output-level` flag, so + block-only tuning flags (`--min-frequency`, `--property-threshold`) stay out of node-mode help; + shared options live in a flattened `AnnotateCommonArgs`. A `--strategy` selector is deferred until a + second strategy exists. - **P5.3** — Docusaurus docs page + CLI reference regeneration, covering **both** output levels. Resolve the §12 punch list first so the docs describe final behaviour. - **P6 (future)** — pypangraph consumer (load the tables; feature→node→block→path joins) + a @@ -315,14 +358,20 @@ single documentation pass (P5.3) follow once both output levels exist. Execution EMBOSS `seqret`, or re-export from bakta/prokka). - **Embed in graph JSON vs separate file** — default to a separate file; revisit if a single self-contained artifact is wanted. -- **Block-level clustering policy** — current lean (§8) is **coordinate-exact identity**: cluster by - the feature's block-consensus terminus endpoints `(start_block, cons_start, end_block, cons_end)`, - using the feature **name only to label** the consensus, not to cluster. Name/product- or - ortholog-based clustering stays a possible alternative/extension; keep it configurable. -- **Block-level refinement level** — the **M-of-N threshold** on identical-coordinate support before - a consensus feature is emitted (all / fraction / any), exposed as a CLI flag; see §8. Distinct from - the clustering policy above. Threshold semantics and whether to allow a small coordinate tolerance - are TBD, informed by the P5.1 prototype. +- **Block-level clustering policy** — **as-built:** **coordinate-exact identity** keyed by the + feature's *whole crossing* — `feature_type` + the 5'→3'-ordered multiset of per-segment + `(block_id, cons_start, cons_end, strand_on_consensus)` — with the feature name used **only to + label** the consensus. (Superseded the original two-outer-endpoint key; see §8 and the P4 note in + §10.) Implemented behind the modular `BlockCompactionStrategy` trait (`CoordinateConsensusStrategy`), + so name/product- or ortholog-based clustering can be added as alternative strategies later. +- **Block-level refinement level** — **as-built:** the **M-of-N threshold** + `M >= max(1, ceil(min_frequency · N))`, with `N` = genomes **structurally capable** of the crossing + (traversing every block with the required multiplicity). Exposed on the CLI as + `annotate blocks --min-frequency` (default 0.9), with `--property-threshold` (default 0.5) gating + consensus name/attribute promotion (**as-built, P5.2**). Per-segment support + (`n_support_segment / n_total_segment`) is reported alongside but does not gate. Whether to allow a + small coordinate **tolerance** (to absorb a terminus nudged by a nearby indel), or to gate + per-segment rather than per-crossing, is still open. - **Features wholly inside an insertion** — drop, or keep node-level-only with no consensus coordinate? Lean towards keep-and-flag. - **Coordinate convention** — 0-based half-open internally; convert only at GFF I/O (GFF is 1-based diff --git a/docs/docs/reference.md b/docs/docs/reference.md index e6bb1725..da27dbd0 100644 --- a/docs/docs/reference.md +++ b/docs/docs/reference.md @@ -32,6 +32,8 @@ If you have Pangraph CLI installed, you can type `pangraph --help` to read the l * [`pangraph simplify`↴](#pangraph-simplify) * [`pangraph reconstruct`↴](#pangraph-reconstruct) * [`pangraph annotate`↴](#pangraph-annotate) +* [`pangraph annotate nodes`↴](#pangraph-annotate-nodes) +* [`pangraph annotate blocks`↴](#pangraph-annotate-blocks) * [`pangraph schema`↴](#pangraph-schema) * [`pangraph completions`↴](#pangraph-completions) * [`pangraph help-markdown`↴](#pangraph-help-markdown) @@ -335,7 +337,55 @@ Reconstruct all input fasta sequences from graph Lift genome annotations onto the pangenome graph -**Usage:** `pangraph annotate [OPTIONS] --gff [INPUT]` +**Usage:** `pangraph annotate ` + +###### **Subcommands:** + +* `nodes` — Lift annotations to per-node block-consensus coordinates (lossless, long-format CSV) +* `blocks` — Compact node-level annotations into block-level consensus features (CSV) + + + +## `pangraph annotate nodes` + +Lift annotations to per-node block-consensus coordinates (lossless, long-format CSV) + +**Usage:** `pangraph annotate nodes [OPTIONS] --gff ... [INPUT]` + +###### **Arguments:** + +* `` — Path to Pangraph JSON. + + Accepts plain or compressed file. If a compressed file is provided, it will be transparently decompressed. Supported compression formats: `gz`, `bz2`, `xz`, `zstd`. Decompressor is chosen based on file extension. + + If no input file provided, the uncompressed input is read from standard input (stdin). + +###### **Options:** + +* `--gff ` — Path(s) to GFF3 annotation file(s). + + Pass several files after a single flag (`--gff a.gff b.gff`, so shell globs like `--gff *.gff` work), and/or repeat the flag (`--gff a.gff --gff b.gff`); the values accumulate. To avoid the positional graph being slurped as an extra GFF, give it before the flag (`annotate nodes graph.json --gff *.gff`) or pipe it via stdin. + + Accepts plain or compressed files (`gz`, `bz2`, `xz`, `zstd`), chosen by file extension. At least one file is required. Annotation `seqid`s must match the graph path names exactly. +* `--only-type ` — Keep only annotations of these feature type(s) (GFF `type` column); drop all others. + + Comma-separated list (`--only-type gene,CDS`) and/or repeat the flag; values accumulate. Matching is exact and case-sensitive (`CDS`, `gene`, `region`). Mutually exclusive with `--exclude-type`. +* `--exclude-type ` — Drop annotations of these feature type(s) (GFF `type` column); keep all others. + + e.g. `--exclude-type region` removes whole-contig `region` declarations. Same comma-separated syntax as `--only-type`; mutually exclusive with it. +* `-o`, `--output ` — Path to the output annotation table (CSV). + + Will be created if it does not exist. The output is compressed if the path ends in a known compression extension (`gz`, `bz2`, `xz`, `zstd`). Use `-` to write uncompressed CSV to standard output (stdout). + + Default value: `-` + + + +## `pangraph annotate blocks` + +Compact node-level annotations into block-level consensus features (CSV) + +**Usage:** `pangraph annotate blocks [OPTIONS] --gff ... [INPUT]` ###### **Arguments:** @@ -347,14 +397,32 @@ Lift genome annotations onto the pangenome graph ###### **Options:** -* `--gff ` — Path to a GFF3 annotation file. Repeat the flag to provide multiple files. +* `--gff ` — Path(s) to GFF3 annotation file(s). + + Pass several files after a single flag (`--gff a.gff b.gff`, so shell globs like `--gff *.gff` work), and/or repeat the flag (`--gff a.gff --gff b.gff`); the values accumulate. To avoid the positional graph being slurped as an extra GFF, give it before the flag (`annotate nodes graph.json --gff *.gff`) or pipe it via stdin. Accepts plain or compressed files (`gz`, `bz2`, `xz`, `zstd`), chosen by file extension. At least one file is required. Annotation `seqid`s must match the graph path names exactly. -* `-o`, `--output ` — Path to the output node-level annotation table (CSV). +* `--only-type ` — Keep only annotations of these feature type(s) (GFF `type` column); drop all others. + + Comma-separated list (`--only-type gene,CDS`) and/or repeat the flag; values accumulate. Matching is exact and case-sensitive (`CDS`, `gene`, `region`). Mutually exclusive with `--exclude-type`. +* `--exclude-type ` — Drop annotations of these feature type(s) (GFF `type` column); keep all others. + + e.g. `--exclude-type region` removes whole-contig `region` declarations. Same comma-separated syntax as `--only-type`; mutually exclusive with it. +* `-o`, `--output ` — Path to the output annotation table (CSV). Will be created if it does not exist. The output is compressed if the path ends in a known compression extension (`gz`, `bz2`, `xz`, `zstd`). Use `-` to write uncompressed CSV to standard output (stdout). Default value: `-` +* `--min-frequency ` — Minimum frequency required to emit a block-level cluster. + + A cluster is kept when the number of supporting genomes `M >= ceil(min_frequency * N)`, where `N` is the number of paths traversing the cluster's block(s) (so a gene is not penalised for being absent in genomes that lack the block entirely). + + Default value: `0.9` +* `--property-threshold ` — Minimum supporter agreement required to promote a consensus name or attribute value. + + For each cluster, a `name`/attribute value is written only if at least this fraction of the supporting genomes agree on it; otherwise the field is left empty. + + Default value: `0.5` diff --git a/packages/pangraph/src/annotation/compact.rs b/packages/pangraph/src/annotation/compact.rs new file mode 100644 index 00000000..6a5dcd33 --- /dev/null +++ b/packages/pangraph/src/annotation/compact.rs @@ -0,0 +1,981 @@ +use crate::annotation::lift::LiftedAnnotation; +use crate::pangraph::pangraph::Pangraph; +use crate::pangraph::pangraph_block::BlockId; +use crate::pangraph::pangraph_path::PathId; +use crate::pangraph::strand::Strand; +use eyre::Report; +use serde::{Deserialize, Serialize}; +use std::collections::{BTreeMap, BTreeSet}; + +/// One segment of a feature's block-consensus crossing: a block, its consensus coordinates +/// `[cons_start, cons_end)`, and the strand relative to that block's consensus. +type Segment = (BlockId, usize, usize, Option); + +/// Identity of a block-level cluster: the feature type plus the feature's whole crossing as a +/// canonical, 5'→3'-ordered multiset of [`Segment`]s. Two per-genome instances are "the same" +/// consensus feature exactly when these match — every block, coordinate and per-block strand along +/// the crossing agrees (and, for a duplicated block crossed twice, the multiplicity agrees too). +type ClusterKey = (String, Vec); + +/// Identity of one placement on a single block, used to tally per-segment support across the whole +/// node-level table (independent of clustering): `(feature_type, block, cons_start, cons_end, strand)`. +type SegmentKey = (String, BlockId, usize, usize, Option); + +/// Per-genome metadata kept for one supporter of a cluster: its feature `name` and `attributes`. +type SupporterMeta = (Option, Vec<(String, String)>); + +/// A block-level consensus annotation: **one segment** of a feature crossing that recurs, at the +/// *same* block-consensus coordinates, across enough of the genomes that carry it. This is the +/// compacted, opinionated view built on top of the lossless node-level table (it never replaces it). +/// +/// A multi-block feature emits one row per segment, all sharing a `cluster_id` and ordered 5'→3' by +/// `segment_idx` (a duplicated block crossed twice yields two rows, distinguished by `segment_idx`). +/// `n_support`/`n_total` are **crossing-level** (identical across a cluster's rows); +/// `n_support_segment`/`n_total_segment` are **this segment's** reproducibility on its own block. +#[derive(Clone, Debug, Serialize, Deserialize, PartialEq, Eq)] +pub struct BlockAnnotation { + /// Feature type shared by the cluster, e.g. `"CDS"` or `"gene"`. + pub feature_type: String, + + /// Identifier of the consensus feature (crossing) this segment belongs to; shared by the rows of + /// one cluster, assigned `0..` in deterministic key order over the emitted clusters. + pub cluster_id: usize, + + /// 0-based position of this segment within the feature, ordered 5'→3' (or, for unstranded + /// features, by the canonical multiset order). + pub segment_idx: usize, + + /// Total number of segments in the crossing. + pub n_segments: usize, + + /// Block whose consensus this segment lies on. + pub block_id: BlockId, + + /// This segment's block-consensus coordinates `[cons_start, cons_end)`. + pub cons_start: usize, + pub cons_end: usize, + + /// Strand of this segment relative to its block consensus (`None` if unstranded). + pub strand_on_consensus: Option, + + /// Majority feature name across the supporters, when it clears the property threshold; else `None`. + pub consensus_name: Option, + + /// Per-attribute majority values that clear the property threshold (key-sorted). + pub consensus_attributes: Vec<(String, String)>, + + /// `M` — genomes sharing this exact crossing (crossing-level; same on every row of the cluster). + pub n_support: usize, + + /// `N` — genomes structurally capable of the crossing, i.e. traversing every block in it with the + /// required multiplicity (crossing-level; same on every row of the cluster). + pub n_total: usize, + + /// `M_seg` — genomes placing a feature-segment at exactly this `(block, coords, strand)`, counted + /// across the whole node-level table (so a body segment shared by several crossings reports the + /// same value in each). + pub n_support_segment: usize, + + /// `N_seg` — distinct genomes traversing this segment's block (its depth). + pub n_total_segment: usize, +} + +/// A pluggable strategy for compacting node-level annotations to the block level. +/// +/// Compaction is an opinionated, configurable summary built on top of the node-level lift; making +/// it a trait keeps the policy swappable (coordinate-consensus today; name- or ortholog-based +/// strategies could be added later) without touching the lift or the writers. +pub trait BlockCompactionStrategy { + /// Compact node-level [`LiftedAnnotation`]s into block-level [`BlockAnnotation`]s. + fn compact(&self, node_annotations: &[LiftedAnnotation], graph: &Pangraph) -> Result, Report>; +} + +/// Compact by **block-consensus coordinate agreement**. +/// +/// Each genome's feature instance is reduced to its whole crossing — the 5'→3'-ordered multiset of +/// per-segment `(block, cons_start, cons_end, strand)` placements; genomes whose `(feature_type, +/// crossing)` match exactly form a cluster. A cluster is emitted when its support `M` reaches +/// `ceil(min_frequency * N)` (at least 1), where `N` is the number of genomes structurally capable +/// of the crossing (traversing every block in it with the required multiplicity). Each emitted +/// cluster produces one row per segment. Consensus `name`/attributes are the per-value majorities +/// among the supporters that clear `property_threshold`. +pub struct CoordinateConsensusStrategy { + /// Minimum supporter fraction `M / N` (of genomes capable of the crossing) for a cluster to be emitted. + pub min_frequency: f64, + + /// Minimum supporter fraction for a metadata value (name / attribute) to be promoted to consensus. + pub property_threshold: f64, +} + +impl Default for CoordinateConsensusStrategy { + /// Provisional defaults; the user-facing defaults are owned by the CLI (P5.2). + fn default() -> Self { + Self { + min_frequency: 0.9, + property_threshold: 0.5, + } + } +} + +impl BlockCompactionStrategy for CoordinateConsensusStrategy { + fn compact(&self, node_annotations: &[LiftedAnnotation], graph: &Pangraph) -> Result, Report> { + // Aggregate per cluster key: dedup supporters by genome (a genome counts once) and reduce the + // per-block node requirement to its element-wise minimum over them. + let mut clusters: BTreeMap = BTreeMap::new(); + for inst in feature_instances(node_annotations) { + let agg = clusters.entry(inst.key).or_default(); + agg + .supporters + .entry(inst.path_id) + .or_insert((inst.name, inst.attributes)); + for (block_id, req) in inst.node_req { + agg + .node_req + .entry(block_id) + .and_modify(|cur| *cur = (*cur).min(req)) + .or_insert(req); + } + } + + // Per-block genome multiplicities (for crossing-level `N` and per-segment `N_seg`) and the + // per-segment support tally (for `M_seg`), both independent of the clustering above. + let block_depths = block_genome_counts(graph); + let segment_support = segment_support_sets(node_annotations); + + // Iterating the BTreeMap yields clusters in key order, so `cluster_id` (assigned only to emitted + // clusters) and the row order are deterministic. + let mut out = Vec::new(); + let mut cluster_id = 0; + for (key, agg) in clusters { + let (feature_type, segments) = key; + let ClusterAgg { supporters, node_req } = agg; + let m = supporters.len(); + let n = crossing_n_total(&block_depths, &node_req); + let min_support = min_count(self.min_frequency, n).max(1); + if m < min_support { + continue; // below threshold: stays only in the node-level table + } + let consensus_name = majority_name(&supporters, m, self.property_threshold); + let consensus_attributes = majority_attributes(&supporters, m, self.property_threshold); + let n_segments = segments.len(); + for (segment_idx, &(block_id, cons_start, cons_end, strand)) in segments.iter().enumerate() { + let seg_key = (feature_type.clone(), block_id, cons_start, cons_end, strand); + let n_support_segment = segment_support.get(&seg_key).map_or(0, BTreeSet::len); + let n_total_segment = block_depths.get(&block_id).map_or(0, BTreeMap::len); + out.push(BlockAnnotation { + feature_type: feature_type.clone(), + cluster_id, + segment_idx, + n_segments, + block_id, + cons_start, + cons_end, + strand_on_consensus: strand, + consensus_name: consensus_name.clone(), + consensus_attributes: consensus_attributes.clone(), + n_support: m, + n_total: n, + n_support_segment, + n_total_segment, + }); + } + cluster_id += 1; + } + Ok(out) + } +} + +/// One genome's instance of a feature, reduced to its cluster key (the whole crossing) plus the +/// metadata compaction needs (representative `name`/`attributes`, identical across a feature's +/// segments) and the **distinct nodes used per block** (for the structural-capability denominator). +struct FeatureInstance { + path_id: PathId, + key: ClusterKey, + /// Distinct nodes the crossing uses on each block — *not* segment count, so an origin-spanning + /// node split into two coverage pieces still requires only one node of its block. + node_req: BTreeMap, + name: Option, + attributes: Vec<(String, String)>, +} + +/// Aggregated supporters of one cluster, plus the per-block node requirement reduced across them. +#[derive(Default)] +struct ClusterAgg { + /// Supporters keyed by genome path id (a genome counts once), each with its representative metadata. + supporters: BTreeMap, + /// Per-block minimum distinct-node requirement over the supporters — the loosest hosting, so + /// every supporter satisfies it and `N >= M` always holds. + node_req: BTreeMap, +} + +/// Group node-level rows into per-genome feature instances and reduce each to its cluster key. +/// +/// Rows are grouped by `(path_id, feature_type, base feature id)` — keying the genome on its stable +/// `PathId`, not the display name, so paths that share a name never merge; instances that are not +/// bounded by exactly two real termini (e.g. partial / truncated features) are dropped from +/// compaction and remain in the node-level table. +fn feature_instances(node_annotations: &[LiftedAnnotation]) -> Vec { + let mut groups: BTreeMap<(PathId, String, String), Vec<&LiftedAnnotation>> = BTreeMap::new(); + for a in node_annotations { + let group_key = (a.path_id, a.feature_type.clone(), base_feature_id(a)); + groups.entry(group_key).or_default().push(a); + } + + groups + .into_iter() + .filter_map(|((path_id, feature_type, _base), rows)| reduce_instance(path_id, feature_type, &rows)) + .collect() +} + +/// Reduce one feature instance's node-level rows to its [`FeatureInstance`], or `None` when the +/// crossing is not bounded by exactly two real termini (a partial / truncated feature). +/// +/// The cluster key is the feature's whole crossing as a 5'→3'-ordered multiset of +/// `(block, cons_start, cons_end, strand_on_consensus)`. Node-level `segment_idx` is genome +/// (low→high coordinate) order, which is already 5'→3' for a forward feature; a reverse feature is +/// reversed, and an unstranded feature (no reading direction) is canonicalized by the +/// lexicographically smaller of the two orders so homologous instances still collapse. +fn reduce_instance(path_id: PathId, feature_type: String, rows: &[&LiftedAnnotation]) -> Option { + // A clean feature is bounded by exactly two real termini (its 5' and 3' ends); the internal + // boundaries where it crosses node/block splits are not termini. Anything else is partial. + let terminus_count: usize = rows + .iter() + .map(|a| usize::from(a.start_is_terminus) + usize::from(a.end_is_terminus)) + .sum(); + if terminus_count != 2 { + return None; + } + + // Segments in genome (arc) order, then oriented 5'→3'. + let mut rows_sorted: Vec<&&LiftedAnnotation> = rows.iter().collect(); + rows_sorted.sort_by_key(|a| a.segment_idx); + let mut segments: Vec = rows_sorted + .iter() + .map(|a| (a.block_id, a.cons_start, a.cons_end, a.strand_on_consensus)) + .collect(); + match rows_sorted[0].feature_strand { + Some(Strand::Reverse) => segments.reverse(), + Some(Strand::Forward) => {}, // genome order is already 5'→3' + None => { + let mut reversed = segments.clone(); + reversed.reverse(); + if reversed < segments { + segments = reversed; + } + }, + } + + // Distinct nodes used per block (not segments): an origin-spanning node split into two coverage + // pieces still needs only one node of its block, while a genuinely duplicated block needs as many. + let mut node_sets: BTreeMap> = BTreeMap::new(); + for a in rows { + node_sets.entry(a.block_id).or_default().insert(a.node_id); + } + let node_req = node_sets.into_iter().map(|(b, nodes)| (b, nodes.len())).collect(); + + // Name/attributes are identical across a feature's segments; take them from the 5'-most. + let rep = rows_sorted[0]; + Some(FeatureInstance { + path_id, + key: (feature_type, segments), + node_req, + name: rep.name.clone(), + attributes: rep.attributes.clone(), + }) +} + +/// Recover the base feature id shared across a feature's segments: its `parent_feature_id` when +/// present, else the per-row `feature_id` with any `.seg{idx}` suffix stripped. +fn base_feature_id(a: &LiftedAnnotation) -> String { + if let Some(parent) = &a.parent_feature_id { + return parent.clone(); + } + if a.n_segments > 1 { + let suffix = format!(".seg{}", a.segment_idx); + if let Some(base) = a.feature_id.strip_suffix(&suffix) { + return base.to_owned(); + } + } + a.feature_id.clone() +} + +/// Per block, the number of nodes (block instances) each genome carries of it, keyed on the genome's +/// stable `PathId`. A genome appears with a count > 1 exactly when the block is duplicated on its +/// path; `map.len()` is the block's depth (distinct genomes traversing it). Drives both crossing-level +/// `N` and per-segment `N_seg`. +fn block_genome_counts(graph: &Pangraph) -> BTreeMap> { + graph + .blocks + .iter() + .map(|(&bid, block)| { + let mut counts: BTreeMap = BTreeMap::new(); + for pid in block.isolates(graph) { + *counts.entry(pid).or_default() += 1; + } + (bid, counts) + }) + .collect() +} + +/// `N` for a crossing: genomes **structurally capable** of it — those traversing every block in the +/// crossing at least `required` times (the crossing's per-block distinct-node count, so an +/// origin-split node is not double-counted). With all requirements 1 this is the intersection of the +/// per-block genome sets. `required` is reduced to its element-wise min over the supporters, so every +/// supporter is capable and `N >= M` holds. +fn crossing_n_total( + block_depths: &BTreeMap>, + required: &BTreeMap, +) -> usize { + let empty = BTreeMap::new(); + let mut blocks = required.iter(); + let Some((first_block, &first_mult)) = blocks.next() else { + return 0; + }; + // Genomes carrying the first block enough times, then narrowed by every remaining block. + let mut candidates: BTreeSet = block_depths + .get(first_block) + .unwrap_or(&empty) + .iter() + .filter(|&(_, &count)| count >= first_mult) + .map(|(&pid, _)| pid) + .collect(); + for (block_id, &mult) in blocks { + let counts = block_depths.get(block_id).unwrap_or(&empty); + candidates.retain(|pid| counts.get(pid).is_some_and(|&count| count >= mult)); + } + candidates.len() +} + +/// Per-segment support across the whole node-level table: for each distinct placement +/// `(feature_type, block, cons_start, cons_end, strand)`, the set of genomes (by `PathId`) that place +/// a feature-segment there. Independent of clustering, so a body segment shared by several crossings +/// reports the same support in each. +fn segment_support_sets(node_annotations: &[LiftedAnnotation]) -> BTreeMap> { + let mut map: BTreeMap> = BTreeMap::new(); + for a in node_annotations { + let key = ( + a.feature_type.clone(), + a.block_id, + a.cons_start, + a.cons_end, + a.strand_on_consensus, + ); + map.entry(key).or_default().insert(a.path_id); + } + map +} + +/// Minimum supporter count to clear a fraction `f` of `n`: `ceil(f * n)`. +/// +/// The product is nudged down by a tiny epsilon before rounding up: a value that is mathematically +/// an integer can land just above it in `f64` (e.g. `0.55 * 100` is `55.000000000000007`), which a +/// naive `.ceil()` would round up, spuriously demanding one extra supporter. +fn min_count(fraction: f64, n: usize) -> usize { + const EPS: f64 = 1e-9; + (fraction * n as f64 - EPS).ceil().max(0.0) as usize +} + +/// The majority feature name across the supporters, when its support clears `threshold * M`. +/// Ties break to the lexicographically smallest name for determinism. +fn majority_name(per_genome: &BTreeMap, m: usize, threshold: f64) -> Option { + let mut counts: BTreeMap = BTreeMap::new(); + for (name, _) in per_genome.values() { + if let Some(name) = name { + *counts.entry(name.clone()).or_default() += 1; + } + } + let needed = min_count(threshold, m); + counts + .into_iter() + .max_by(|a, b| a.1.cmp(&b.1).then_with(|| b.0.cmp(&a.0))) + .filter(|(_, c)| *c >= needed) + .map(|(name, _)| name) +} + +/// Per-attribute majority values across the supporters that clear `threshold * M`, key-sorted. +/// A supporter contributes each of its `(key, value)` pairs once (duplicates within a supporter are +/// collapsed). Per key, the value with the most support wins; ties break to the smaller value. +fn majority_attributes( + per_genome: &BTreeMap, + m: usize, + threshold: f64, +) -> Vec<(String, String)> { + let mut pair_counts: BTreeMap<(String, String), usize> = BTreeMap::new(); + for (_, attrs) in per_genome.values() { + let distinct: BTreeSet<&(String, String)> = attrs.iter().collect(); + for kv in distinct { + *pair_counts.entry(kv.clone()).or_default() += 1; + } + } + + let mut by_key: BTreeMap = BTreeMap::new(); + for ((k, v), c) in pair_counts { + let best = by_key.entry(k).or_insert_with(|| (v.clone(), c)); + if c > best.1 || (c == best.1 && v < best.0) { + *best = (v, c); + } + } + + let needed = min_count(threshold, m); + by_key + .into_iter() + .filter(|(_, (_, c))| *c >= needed) + .map(|(k, (v, _))| (k, v)) + .collect() +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::pangraph::edits::Edit; + use crate::pangraph::pangraph_block::PangraphBlock; + use crate::pangraph::pangraph_node::{NodeId, PangraphNode}; + use crate::pangraph::pangraph_path::PangraphPath; + use crate::pangraph::strand::Strand::{Forward, Reverse}; + use pretty_assertions::assert_eq; + + /// Build a minimal graph: `path_names[i]` becomes `PathId(i)`; each `(block_idx, paths)` entry + /// makes `BlockId(block_idx)` traversed by those path indices (one node per listed index, so a + /// path index listed twice gives that genome two nodes of the block — a duplication). Only the + /// block→path traversal and path names matter to compaction, so positions/consensus are placeholders. + fn graph_with(path_names: &[&str], blocks: &[(usize, &[usize])]) -> Pangraph { + let mut nodes = BTreeMap::new(); + let mut blocks_map = BTreeMap::new(); + for &(bidx, path_idxs) in blocks { + let mut alignments = BTreeMap::new(); + for (occ, &pidx) in path_idxs.iter().enumerate() { + let nid = NodeId(bidx * 1000 + pidx * 10 + occ); + nodes.insert( + nid, + PangraphNode::new(Some(nid), BlockId(bidx), PathId(pidx), Forward, (0, 0)), + ); + alignments.insert(nid, Edit::empty()); + } + blocks_map.insert(BlockId(bidx), PangraphBlock::new(BlockId(bidx), "A", alignments)); + } + let paths = path_names + .iter() + .enumerate() + .map(|(i, name)| { + ( + PathId(i), + PangraphPath::new( + Some(PathId(i)), + Vec::::new(), + 0, + false, + Some((*name).to_owned()), + None, + ), + ) + }) + .collect::>(); + Pangraph { + paths, + blocks: blocks_map, + nodes, + } + } + + /// A single-segment lifted annotation (both endpoints are termini). `path` is the genome's path + /// index — it must match the `PathId` the genome has in the graph, since compaction keys on it. + fn lifted( + path: usize, + block: usize, + cons: (usize, usize), + strand: Option, + ftype: &str, + name: Option<&str>, + id: &str, + attrs: &[(&str, &str)], + ) -> LiftedAnnotation { + LiftedAnnotation { + feature_id: id.to_owned(), + parent_feature_id: Some(id.to_owned()), + segment_idx: 0, + n_segments: 1, + genome: format!("g{path}"), + path_id: PathId(path), + block_id: BlockId(block), + node_id: NodeId(0), + strand_on_consensus: strand, + // Single segment: ordering is irrelevant, so the feature strand can mirror the consensus strand. + feature_strand: strand, + node_start: cons.0, + node_end: cons.1, + cons_start: cons.0, + cons_end: cons.1, + start_is_terminus: true, + end_is_terminus: true, + start_in_insertion: false, + end_in_insertion: false, + frac_covered: 1.0, + feature_type: ftype.to_owned(), + name: name.map(str::to_owned), + attributes: attrs.iter().map(|(k, v)| ((*k).to_owned(), (*v).to_owned())).collect(), + } + } + + /// One segment of a multi-segment feature: `idx` is its node-level (genome/arc-order) index, + /// `soc` its strand on this block's consensus, and `fstrand` the feature's genome strand (the same + /// across all its segments). A reverse `fstrand` means the gene reads 5'→3' against genome order. + #[allow(clippy::too_many_arguments)] + fn seg( + path: usize, + base: &str, + idx: usize, + n: usize, + block: usize, + cons: (usize, usize), + termini: (bool, bool), + soc: Option, + fstrand: Option, + ) -> LiftedAnnotation { + LiftedAnnotation { + feature_id: if n > 1 { + format!("{base}.seg{idx}") + } else { + base.to_owned() + }, + parent_feature_id: Some(base.to_owned()), + segment_idx: idx, + n_segments: n, + genome: format!("g{path}"), + path_id: PathId(path), + block_id: BlockId(block), + // Distinct per segment so a block crossed at two segments counts as two nodes (a duplication), + // matching how the lift assigns a node per crossed block instance. + node_id: NodeId(idx), + strand_on_consensus: soc, + feature_strand: fstrand, + node_start: cons.0, + node_end: cons.1, + cons_start: cons.0, + cons_end: cons.1, + start_is_terminus: termini.0, + end_is_terminus: termini.1, + start_in_insertion: false, + end_in_insertion: false, + frac_covered: 0.5, + feature_type: "CDS".to_owned(), + name: Some("geneA".to_owned()), + attributes: vec![], + } + } + + fn strat(min_frequency: f64, property_threshold: f64) -> CoordinateConsensusStrategy { + CoordinateConsensusStrategy { + min_frequency, + property_threshold, + } + } + + #[test] + fn test_min_count_avoids_float_ceil_off_by_one() { + // Pairs where `f * n` lands just above an integer in `f64` (`0.55 * 100` is + // `55.000000000000007`), which a naive `.ceil()` would round up to one supporter too many. + assert_eq!(min_count(0.55, 100), 55); + assert_eq!(min_count(0.56, 25), 14); + assert_eq!(min_count(0.14, 50), 7); + // Genuinely fractional products still round up. + assert_eq!(min_count(0.5, 3), 2); + assert_eq!(min_count(0.9, 10), 9); + // Exact integers and the degenerate zero case are unchanged. + assert_eq!(min_count(1.0, 7), 7); + assert_eq!(min_count(0.0, 5), 0); + } + + #[test] + fn test_all_agree_single_block() { + let graph = graph_with(&["g0", "g1", "g2"], &[(1, &[0, 1, 2])]); + let anns = vec![ + lifted(0, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f0", &[]), + lifted(1, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f1", &[]), + lifted(2, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f2", &[]), + ]; + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + assert_eq!(out.len(), 1); + let b = &out[0]; + assert_eq!((b.cluster_id, b.segment_idx, b.n_segments), (0, 0, 1)); + assert_eq!((b.block_id, b.cons_start, b.cons_end), (BlockId(1), 10, 200)); + assert_eq!((b.n_support, b.n_total), (3, 3)); + assert_eq!((b.n_support_segment, b.n_total_segment), (3, 3)); + assert_eq!(b.consensus_name.as_deref(), Some("geneA")); + assert_eq!(b.strand_on_consensus, Some(Forward)); + } + + #[test] + fn test_below_threshold_dropped_but_kept_when_lenient() { + let graph = graph_with(&["g0", "g1", "g2", "g3"], &[(1, &[0, 1, 2, 3])]); + let anns = vec![ + lifted(0, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f0", &[]), + lifted(1, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f1", &[]), + ]; + // 2 of 4 -> ceil(0.9*4)=4 required -> dropped. + assert!(strat(0.9, 0.5).compact(&anns, &graph).unwrap().is_empty()); + // ceil(0.5*4)=2 required -> kept. + let out = strat(0.5, 0.5).compact(&anns, &graph).unwrap(); + assert_eq!(out.len(), 1); + assert_eq!((out[0].n_support, out[0].n_total), (2, 4)); + } + + #[test] + fn test_feature_type_splits_clusters() { + let graph = graph_with(&["g0", "g1"], &[(1, &[0, 1])]); + let anns = vec![ + lifted(0, 1, (10, 200), Some(Forward), "gene", Some("geneA"), "gene0", &[]), + lifted(0, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "cds0", &[]), + lifted(1, 1, (10, 200), Some(Forward), "gene", Some("geneA"), "gene1", &[]), + lifted(1, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "cds1", &[]), + ]; + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + assert_eq!(out.len(), 2); + let types: Vec<&str> = out.iter().map(|b| b.feature_type.as_str()).collect(); + assert_eq!(types, vec!["CDS", "gene"]); // same coords -> sorted by type + assert!(out.iter().all(|b| (b.n_support, b.n_total) == (2, 2))); + } + + #[test] + fn test_strand_splits_clusters() { + let graph = graph_with(&["g0", "g1"], &[(1, &[0, 1])]); + let anns = vec![ + lifted(0, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f0", &[]), + lifted(1, 1, (10, 200), Some(Reverse), "CDS", Some("geneA"), "f1", &[]), + ]; + let out = strat(0.0, 0.5).compact(&anns, &graph).unwrap(); + assert_eq!(out.len(), 2); + assert_eq!(out[0].strand_on_consensus, Some(Forward)); // Forward < Reverse + assert_eq!(out[1].strand_on_consensus, Some(Reverse)); + assert_eq!((out[0].n_support, out[0].n_total), (1, 2)); + } + + #[test] + fn test_name_and_attribute_majority_threshold() { + let graph = graph_with(&["g0", "g1", "g2"], &[(1, &[0, 1, 2])]); + let anns = vec![ + lifted( + 0, + 1, + (10, 200), + Some(Forward), + "CDS", + Some("geneA"), + "f0", + &[("product", "widget")], + ), + lifted( + 1, + 1, + (10, 200), + Some(Forward), + "CDS", + Some("geneA"), + "f1", + &[("product", "widget")], + ), + lifted( + 2, + 1, + (10, 200), + Some(Forward), + "CDS", + Some("geneB"), + "f2", + &[("product", "gadget")], + ), + ]; + // 2/3 majority clears ceil(0.5*3)=2. + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + assert_eq!(out.len(), 1); + assert_eq!(out[0].consensus_name.as_deref(), Some("geneA")); + assert_eq!( + out[0].consensus_attributes, + vec![("product".to_owned(), "widget".to_owned())] + ); + // 2/3 does not clear ceil(0.7*3)=3. + let out = strat(1.0, 0.7).compact(&anns, &graph).unwrap(); + assert_eq!(out[0].consensus_name, None); + assert!(out[0].consensus_attributes.is_empty()); + } + + #[test] + fn test_multi_block_feature_emits_one_row_per_segment() { + let graph = graph_with(&["g0", "g1"], &[(1, &[0, 1]), (2, &[0, 1])]); + let mut anns = Vec::new(); + for (g, base) in [(0, "f0"), (1, "f1")] { + anns.push(seg( + g, + base, + 0, + 2, + 1, + (150, 200), + (true, false), + Some(Forward), + Some(Forward), + )); + anns.push(seg( + g, + base, + 1, + 2, + 2, + (0, 50), + (false, true), + Some(Forward), + Some(Forward), + )); + } + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + // One cluster, two segments ordered 5'→3'. + assert_eq!(out.len(), 2); + assert!(out.iter().all(|b| b.cluster_id == 0 && b.n_segments == 2)); + assert_eq!( + (out[0].segment_idx, out[0].block_id, out[0].cons_start, out[0].cons_end), + (0, BlockId(1), 150, 200) + ); + assert_eq!( + (out[1].segment_idx, out[1].block_id, out[1].cons_start, out[1].cons_end), + (1, BlockId(2), 0, 50) + ); + assert!(out.iter().all(|b| (b.n_support, b.n_total) == (2, 2))); + } + + #[test] + fn test_n_total_is_paths_traversing_block_not_all_paths() { + // Block 1 is accessory: present only in g0, g1 (not g2, g3). + let graph = graph_with(&["g0", "g1", "g2", "g3"], &[(1, &[0, 1])]); + let anns = vec![ + lifted(0, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f0", &[]), + lifted(1, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f1", &[]), + ]; + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + assert_eq!(out.len(), 1); + assert_eq!((out[0].n_support, out[0].n_total), (2, 2)); // N is 2, not 4 + } + + #[test] + fn test_partial_feature_without_two_termini_is_excluded() { + let graph = graph_with(&["g0"], &[(1, &[0])]); + let mut a = lifted(0, 1, (10, 200), Some(Forward), "CDS", Some("geneA"), "f0", &[]); + a.end_is_terminus = false; // only one real terminus remains + assert!(strat(0.0, 0.5).compact(&[a], &graph).unwrap().is_empty()); + } + + /// Regression for the yjeM bug: the same gene in two genomes of opposite global orientation — + /// body on block 1 (`+` on its consensus), inverted tail on block 3 (`-` on its consensus). The + /// forward genome traverses them in arc order [B1, B3]; the reverse genome in [B3, B1] with the + /// per-block strands unchanged. Ordering 5'→3' must collapse both into one cluster (not split by a + /// strand sampled from the genome-lowest segment), giving M = N = 2. + #[test] + fn test_inversion_crossing_collapses_to_one_cluster() { + let graph = graph_with(&["gplus", "gminus"], &[(1, &[0, 1]), (3, &[0, 1])]); + let anns = vec![ + // Forward genome: gene reads 5'→3' with genome order, B1 then B3. + seg( + 0, + "fp", + 0, + 2, + 1, + (9228, 10694), + (true, false), + Some(Forward), + Some(Forward), + ), + seg( + 0, + "fp", + 1, + 2, + 3, + (705, 742), + (false, true), + Some(Reverse), + Some(Forward), + ), + // Reverse genome: arc order is B3 then B1; feature_strand Reverse flips it back to 5'→3'. + seg( + 1, + "fm", + 0, + 2, + 3, + (705, 742), + (true, false), + Some(Reverse), + Some(Reverse), + ), + seg( + 1, + "fm", + 1, + 2, + 1, + (9228, 10694), + (false, true), + Some(Forward), + Some(Reverse), + ), + ]; + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + // A single cluster with two segments, body (B1,+) before tail (B3,-). + assert_eq!(out.len(), 2); + assert!(out.iter().all(|b| b.cluster_id == 0 && b.n_segments == 2)); + assert_eq!( + ( + out[0].segment_idx, + out[0].block_id, + out[0].cons_start, + out[0].cons_end, + out[0].strand_on_consensus + ), + (0, BlockId(1), 9228, 10694, Some(Forward)) + ); + assert_eq!( + ( + out[1].segment_idx, + out[1].block_id, + out[1].cons_start, + out[1].cons_end, + out[1].strand_on_consensus + ), + (1, BlockId(3), 705, 742, Some(Reverse)) + ); + assert!(out.iter().all(|b| (b.n_support, b.n_total) == (2, 2))); + // Both genomes place each segment at the same coords, so per-segment support is also 2 of 2. + assert!(out.iter().all(|b| (b.n_support_segment, b.n_total_segment) == (2, 2))); + } + + /// A feature crossing a duplicated block twice (B5, B6, B5) keeps both occurrences as separate + /// segment rows, and `N` requires genomes to carry block 5 at least twice: g2, which has it only + /// once, is excluded from the denominator. + #[test] + fn test_duplicated_block_crossed_twice() { + let graph = graph_with(&["g0", "g1", "g2"], &[(5, &[0, 0, 1, 1, 2]), (6, &[0, 1, 2])]); + let mut anns = Vec::new(); + for (g, base) in [(0, "f0"), (1, "f1")] { + anns.push(seg( + g, + base, + 0, + 3, + 5, + (0, 50), + (true, false), + Some(Forward), + Some(Forward), + )); + anns.push(seg( + g, + base, + 1, + 3, + 6, + (0, 30), + (false, false), + Some(Forward), + Some(Forward), + )); + anns.push(seg( + g, + base, + 2, + 3, + 5, + (70, 100), + (false, true), + Some(Forward), + Some(Forward), + )); + } + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + // Three segment rows; block 5 appears twice, at its two distinct consensus coordinates. + assert_eq!(out.len(), 3); + assert!(out.iter().all(|b| b.cluster_id == 0 && b.n_segments == 3)); + let blocks: Vec<_> = out + .iter() + .map(|b| (b.segment_idx, b.block_id, b.cons_start, b.cons_end)) + .collect(); + assert_eq!( + blocks, + vec![(0, BlockId(5), 0, 50), (1, BlockId(6), 0, 30), (2, BlockId(5), 70, 100),] + ); + // N counts only genomes carrying block 5 twice (g0, g1) — not g2, which has a single copy. + assert!(out.iter().all(|b| (b.n_support, b.n_total) == (2, 2))); + } + + /// Two segments of one feature landing on the **same node** (an origin-spanning node split into + /// two coverage pieces) need only one node of the block — so `N` must still count the genome and + /// `N >= M` must hold. Regression for the `M > N` bug on whole-genome `region` features, where + /// counting required multiplicity by segment (2) instead of distinct node (1) gave `N = 0`. + #[test] + fn test_origin_split_same_node_not_double_counted() { + let graph = graph_with(&["g0"], &[(7, &[0])]); // block 7 present once on g0 + let mut s0 = seg(0, "f", 0, 2, 7, (100, 150), (true, false), Some(Forward), Some(Forward)); + let mut s1 = seg(0, "f", 1, 2, 7, (0, 40), (false, true), Some(Forward), Some(Forward)); + s0.node_id = NodeId(0); + s1.node_id = NodeId(0); // same node as s0: origin-split, not a duplication + let out = strat(1.0, 0.5).compact(&[s0, s1], &graph).unwrap(); + assert_eq!(out.len(), 2); + assert!(out.iter().all(|b| b.block_id == BlockId(7) && b.n_segments == 2)); + assert!(out.iter().all(|b| (b.n_support, b.n_total) == (1, 1))); + } + + /// Mini-yjeM: a body shared on block 1 with three different tails (blocks 2/3/4) across genomes. + /// Each tail variant is its own cluster, every one at M/N = 1.0, and the shared body reports the + /// same per-segment support (6 of 6) in all three clusters. + #[test] + fn test_shared_body_three_tails_each_confident() { + let graph = graph_with( + &["g0", "g1", "g2", "g3", "g4", "g5"], + &[(1, &[0, 1, 2, 3, 4, 5]), (2, &[0, 1, 2]), (3, &[3, 4]), (4, &[5])], + ); + let mut anns = Vec::new(); + for (g, base, tail) in [ + (0, "a0", 2), + (1, "a1", 2), + (2, "a2", 2), + (3, "a3", 3), + (4, "a4", 3), + (5, "a5", 4), + ] { + anns.push(seg( + g, + base, + 0, + 2, + 1, + (9228, 10694), + (true, false), + Some(Forward), + Some(Forward), + )); + anns.push(seg( + g, + base, + 1, + 2, + tail, + (0, 37), + (false, true), + Some(Forward), + Some(Forward), + )); + } + let out = strat(1.0, 0.5).compact(&anns, &graph).unwrap(); + // Three clusters (B1+B2, B1+B3, B1+B4), two rows each. + assert_eq!(out.len(), 6); + let cluster_ids: BTreeSet<_> = out.iter().map(|b| b.cluster_id).collect(); + assert_eq!(cluster_ids, BTreeSet::from([0, 1, 2])); + // Every cluster is fully supported among the genomes capable of it. + assert!(out.iter().all(|b| b.n_support == b.n_total)); + // The shared body segment reports the same reproducibility (6 of 6) in every cluster. + assert!( + out + .iter() + .filter(|b| b.block_id == BlockId(1)) + .all(|b| (b.n_support_segment, b.n_total_segment) == (6, 6)) + ); + } +} diff --git a/packages/pangraph/src/annotation/feature.rs b/packages/pangraph/src/annotation/feature.rs index 568cfa64..12bb6180 100644 --- a/packages/pangraph/src/annotation/feature.rs +++ b/packages/pangraph/src/annotation/feature.rs @@ -1,6 +1,7 @@ use crate::pangraph::strand::Strand; use crate::utils::interval::Interval; use serde::{Deserialize, Serialize}; +use std::collections::BTreeSet; /// A genome annotation feature, normalized from an annotation file (currently GFF) /// into a single, format-agnostic representation. @@ -47,11 +48,83 @@ pub fn interval_from_one_based_inclusive(start: usize, end: usize) -> Interval { Interval::new(start - 1, end) } +/// Keep only the features whose `feature_type` passes the type filters. +/// +/// When `only` is non-empty, retain only features whose type is in it (whitelist); independently, +/// drop any feature whose type is in `exclude` (blacklist). Empty slices are no-ops, so the default +/// (no flags) passes everything through. Matching is exact (case-sensitive). The two filters are +/// mutually exclusive at the CLI, but applying both here is order-independent and safe. +pub fn filter_features_by_type(features: Vec, only: &[String], exclude: &[String]) -> Vec { + let only: BTreeSet<&str> = only.iter().map(String::as_str).collect(); + let exclude: BTreeSet<&str> = exclude.iter().map(String::as_str).collect(); + features + .into_iter() + .filter(|f| { + (only.is_empty() || only.contains(f.feature_type.as_str())) && !exclude.contains(f.feature_type.as_str()) + }) + .collect() +} + #[cfg(test)] mod tests { use super::*; use pretty_assertions::assert_eq; + /// Build a minimal feature carrying only the `feature_type` the type-filter tests care about. + fn typed(feature_type: &str) -> Feature { + Feature { + seqid: "chr1".to_owned(), + source: None, + feature_type: feature_type.to_owned(), + interval: Interval::new(0, 1), + strand: Some(Strand::Forward), + id: None, + name: None, + attributes: vec![], + } + } + + /// Collect just the `feature_type`s after filtering, for compact assertions. + fn types_after(only: &[&str], exclude: &[&str]) -> Vec { + let features = vec![typed("gene"), typed("CDS"), typed("region")]; + let only: Vec = only.iter().map(|s| (*s).to_owned()).collect(); + let exclude: Vec = exclude.iter().map(|s| (*s).to_owned()).collect(); + filter_features_by_type(features, &only, &exclude) + .into_iter() + .map(|f| f.feature_type) + .collect() + } + + #[test] + fn filter_no_filters_passes_everything() { + assert_eq!(types_after(&[], &[]), vec!["gene", "CDS", "region"]); + } + + #[test] + fn filter_only_keeps_whitelisted_types() { + assert_eq!(types_after(&["gene"], &[]), vec!["gene"]); + assert_eq!(types_after(&["gene", "CDS"], &[]), vec!["gene", "CDS"]); + } + + #[test] + fn filter_exclude_drops_blacklisted_types() { + assert_eq!(types_after(&[], &["region"]), vec!["gene", "CDS"]); + } + + #[test] + fn filter_is_case_sensitive() { + // `cds` does not match `CDS`, so the whitelist keeps nothing. + assert!(types_after(&["cds"], &[]).is_empty()); + } + + #[test] + fn filter_unknown_type_is_a_silent_no_op() { + // Whitelisting a type that no feature has yields an empty result rather than an error. + assert!(types_after(&["mRNA"], &[]).is_empty()); + // Excluding a type that no feature has leaves everything in place. + assert_eq!(types_after(&[], &["mRNA"]), vec!["gene", "CDS", "region"]); + } + #[test] fn test_interval_from_one_based_inclusive_multibase() { assert_eq!(interval_from_one_based_inclusive(1, 3), Interval::new(0, 3)); diff --git a/packages/pangraph/src/annotation/lift.rs b/packages/pangraph/src/annotation/lift.rs index b2ad2e52..2388e4d4 100644 --- a/packages/pangraph/src/annotation/lift.rs +++ b/packages/pangraph/src/annotation/lift.rs @@ -45,6 +45,11 @@ pub struct LiftedAnnotation { /// Name of the genome (pangraph path) the feature belongs to. pub genome: String, + /// Stable identity of the genome (pangraph path). Unlike `genome` — a display name that can + /// collide when two paths share a name (or are both unnamed) — this uniquely keys the path. + /// Block-level compaction groups genomes on this rather than on the name string. + pub path_id: PathId, + /// Block whose consensus this segment is placed on. pub block_id: BlockId, @@ -55,6 +60,12 @@ pub struct LiftedAnnotation { /// node is on the reverse strand). `None` when the source feature was unstranded. pub strand_on_consensus: Option, + /// The feature's original genome strand, straight from the GFF (`None` when unstranded). Unlike + /// `strand_on_consensus`, this is *not* flipped on reverse nodes, so it is identical across a + /// feature's segments; `strand_on_consensus` is this XOR the node orientation. Block-level + /// compaction uses it to order a multi-segment feature's segments 5'→3'. + pub feature_strand: Option, + /// Consensus-oriented node-local coordinates `[node_start, node_end)`. pub node_start: usize, pub node_end: usize, @@ -311,9 +322,11 @@ pub fn lift_feature(feature: &Feature, path: &PangraphPath, graph: &Pangraph) -> segment_idx, n_segments, genome: genome.clone(), + path_id: path.id(), block_id, node_id: seg.node_id, strand_on_consensus, + feature_strand: feature.strand, node_start, node_end, cons_start, @@ -595,6 +608,10 @@ mod tests { (true, false), (false, false), ); + // `feature_strand` keeps the original genome strand (Forward), un-flipped by the reverse nodes + // and identical across both segments, even though `strand_on_consensus` is Reverse. + assert_eq!(lifted[0].feature_strand, Some(Forward)); + assert_eq!(lifted[1].feature_strand, Some(Forward)); } #[test] @@ -683,6 +700,7 @@ mod tests { let lifted = lift_feature(&feat(2, 5, None), &path, &graph).unwrap(); assert_eq!(lifted.len(), 1); assert_eq!(lifted[0].strand_on_consensus, None); + assert_eq!(lifted[0].feature_strand, None); } #[test] diff --git a/packages/pangraph/src/annotation/mod.rs b/packages/pangraph/src/annotation/mod.rs index 9a461f92..a6431b8e 100644 --- a/packages/pangraph/src/annotation/mod.rs +++ b/packages/pangraph/src/annotation/mod.rs @@ -1,3 +1,4 @@ +pub mod compact; pub mod feature; pub mod lift; pub mod matching; diff --git a/packages/pangraph/src/annotation/writer.rs b/packages/pangraph/src/annotation/writer.rs index f22cd003..f16fe579 100644 --- a/packages/pangraph/src/annotation/writer.rs +++ b/packages/pangraph/src/annotation/writer.rs @@ -1,3 +1,4 @@ +use crate::annotation::compact::BlockAnnotation; use crate::annotation::lift::LiftedAnnotation; use crate::io::file::create_file_or_stdout; use crate::pangraph::pangraph_block::BlockId; @@ -17,6 +18,9 @@ use std::path::Path; pub trait AnnotationWriter { /// Serialize a batch of node-level lifted annotations. fn write_node_annotations(&mut self, annotations: &[LiftedAnnotation]) -> Result<(), Report>; + + /// Serialize a batch of block-level (compacted) annotations. + fn write_block_annotations(&mut self, annotations: &[BlockAnnotation]) -> Result<(), Report>; } /// One CSV row per [`LiftedAnnotation`]. @@ -35,6 +39,7 @@ struct LiftedAnnotationCsvRow<'a> { block_id: BlockId, node_id: NodeId, strand_on_consensus: Option, + feature_strand: Option, node_start: usize, node_end: usize, cons_start: usize, @@ -61,6 +66,7 @@ impl<'a> LiftedAnnotationCsvRow<'a> { block_id: a.block_id, node_id: a.node_id, strand_on_consensus: a.strand_on_consensus, + feature_strand: a.feature_strand, node_start: a.node_start, node_end: a.node_end, cons_start: a.cons_start, @@ -77,6 +83,52 @@ impl<'a> LiftedAnnotationCsvRow<'a> { } } +/// One CSV row per [`BlockAnnotation`] — i.e. one segment of a compacted crossing. +/// +/// Mirrors [`LiftedAnnotationCsvRow`]: borrows from the source, renders `consensus_attributes` as a +/// single JSON-string column (a JSON array of `[key, value]` pairs), `strand_on_consensus` as +/// `+`/`-`/empty, ids as plain numbers, and `Option`s as empty cells. Rows of one crossing share a +/// `cluster_id` and are ordered 5'→3' by `segment_idx`. +#[derive(Serialize)] +struct BlockAnnotationCsvRow<'a> { + #[serde(rename = "type")] + feature_type: &'a str, + cluster_id: usize, + segment_idx: usize, + n_segments: usize, + block_id: BlockId, + cons_start: usize, + cons_end: usize, + strand_on_consensus: Option, + consensus_name: Option<&'a str>, + consensus_attributes: String, + n_support: usize, + n_total: usize, + n_support_segment: usize, + n_total_segment: usize, +} + +impl<'a> BlockAnnotationCsvRow<'a> { + fn from_block(a: &'a BlockAnnotation) -> Result { + Ok(Self { + feature_type: &a.feature_type, + cluster_id: a.cluster_id, + segment_idx: a.segment_idx, + n_segments: a.n_segments, + block_id: a.block_id, + cons_start: a.cons_start, + cons_end: a.cons_end, + strand_on_consensus: a.strand_on_consensus, + consensus_name: a.consensus_name.as_deref(), + consensus_attributes: serde_json::to_string(&a.consensus_attributes)?, + n_support: a.n_support, + n_total: a.n_total, + n_support_segment: a.n_support_segment, + n_total_segment: a.n_total_segment, + }) + } +} + /// The default [`AnnotationWriter`]: long-format CSV, one row per lifted segment. /// /// Backed by [`create_file_or_stdout`], so `-` writes to stdout and the output is transparently @@ -107,12 +159,22 @@ impl AnnotationWriter for CsvAnnotationWriter { self.writer.flush()?; Ok(()) } + + fn write_block_annotations(&mut self, annotations: &[BlockAnnotation]) -> Result<(), Report> { + for ann in annotations { + let row = BlockAnnotationCsvRow::from_block(ann)?; + self.writer.serialize(&row)?; + } + self.writer.flush()?; + Ok(()) + } } #[cfg(test)] mod tests { use super::*; use crate::io::csv::parse_csv; + use crate::pangraph::pangraph_path::PathId; use serde::Deserialize; use std::fs::read_to_string; use tempfile::tempdir; @@ -130,6 +192,7 @@ mod tests { block_id: usize, node_id: usize, strand_on_consensus: Option, + feature_strand: Option, node_start: usize, node_end: usize, cons_start: usize, @@ -152,9 +215,11 @@ mod tests { segment_idx, n_segments: 2, genome: "genomeA".to_owned(), + path_id: PathId(2), block_id: BlockId(7), node_id: NodeId(42), strand_on_consensus: Some(Strand::Reverse), + feature_strand: Some(Strand::Forward), node_start: 3, node_end: 8, cons_start: 3, @@ -196,6 +261,7 @@ mod tests { assert_eq!(r.parent_feature_id.as_deref(), Some("g1")); assert_eq!((r.block_id, r.node_id), (7, 42)); assert_eq!(r.strand_on_consensus.as_deref(), Some("-")); + assert_eq!(r.feature_strand.as_deref(), Some("+")); assert_eq!((r.cons_start, r.cons_end), (3, 9)); assert_eq!((r.start_is_terminus, r.end_is_terminus), (true, false)); assert_eq!((r.start_in_insertion, r.end_in_insertion), (false, true)); @@ -204,6 +270,75 @@ mod tests { assert_eq!(r.attributes, r#"[["ID","g1"],["Name","geneA"]]"#); } + /// Minimal owned mirror of the block-level CSV row. + #[allow(dead_code)] + #[derive(Debug, Deserialize, PartialEq)] + struct BlockRow { + #[serde(rename = "type")] + feature_type: String, + cluster_id: usize, + segment_idx: usize, + n_segments: usize, + block_id: usize, + cons_start: usize, + cons_end: usize, + strand_on_consensus: Option, + consensus_name: Option, + consensus_attributes: String, + n_support: usize, + n_total: usize, + n_support_segment: usize, + n_total_segment: usize, + } + + fn sample_block() -> BlockAnnotation { + BlockAnnotation { + feature_type: "CDS".to_owned(), + cluster_id: 4, + segment_idx: 1, + n_segments: 2, + block_id: BlockId(7), + cons_start: 3, + cons_end: 120, + strand_on_consensus: Some(Strand::Reverse), + consensus_name: Some("geneA".to_owned()), + consensus_attributes: vec![("product".to_owned(), "widget".to_owned())], + n_support: 18, + n_total: 21, + n_support_segment: 19, + n_total_segment: 25, + } + } + + #[test] + fn test_csv_writer_round_trips_block_rows() { + let dir = tempdir().unwrap(); + let path = dir.path().join("block_annotations.csv"); + + let anns = vec![sample_block()]; + { + let mut writer = CsvAnnotationWriter::new(&path, b',').unwrap(); + writer.write_block_annotations(&anns).unwrap(); + } + + let contents = read_to_string(&path).unwrap(); + assert!(contents.starts_with("type,cluster_id,segment_idx")); + assert!(contents.contains(r#"[[""product"",""widget""]]"#)); + + let rows: Vec = parse_csv(&contents).unwrap(); + assert_eq!(rows.len(), 1); + let r = &rows[0]; + assert_eq!(r.feature_type, "CDS"); + assert_eq!((r.cluster_id, r.segment_idx, r.n_segments), (4, 1, 2)); + assert_eq!(r.block_id, 7); + assert_eq!((r.cons_start, r.cons_end), (3, 120)); + assert_eq!(r.strand_on_consensus.as_deref(), Some("-")); + assert_eq!(r.consensus_name.as_deref(), Some("geneA")); + assert_eq!((r.n_support, r.n_total), (18, 21)); + assert_eq!((r.n_support_segment, r.n_total_segment), (19, 25)); + assert_eq!(r.consensus_attributes, r#"[["product","widget"]]"#); + } + #[test] fn test_csv_writer_empty_input_writes_no_header() { let dir = tempdir().unwrap(); diff --git a/packages/pangraph/src/commands/annotate/annotate_args.rs b/packages/pangraph/src/commands/annotate/annotate_args.rs index 613657f6..3a6f5b5d 100644 --- a/packages/pangraph/src/commands/annotate/annotate_args.rs +++ b/packages/pangraph/src/commands/annotate/annotate_args.rs @@ -1,18 +1,33 @@ -use clap::{Parser, ValueHint}; +use clap::{Parser, Subcommand, ValueHint}; use std::fmt::Debug; use std::path::PathBuf; /// Lift genome annotations onto the pangenome graph. /// /// Reads one or more GFF3 annotation files and places each feature on the graph node(s) it overlaps, -/// translating its coordinates into block-consensus coordinates. The result is a long-format, -/// node-level table (the lossless source of truth), written as CSV. +/// translating its coordinates into block-consensus coordinates. Choose the granularity of the +/// output with a subcommand: +/// +/// - `nodes` — the lossless, long-format node-level table (one row per feature/overlapped node). +/// - `blocks` — block-level consensus features, collapsing the redundant per-node placements that +/// most genomes share into one row per cluster. /// /// Annotation `seqid`s are matched to graph path names by exact string equality; any `seqid` that /// does not correspond to a path is a hard error (annotation seqids must match the FASTA record /// names used to build the graph). +#[derive(Subcommand, Debug)] +#[clap(verbatim_doc_comment)] +pub enum PangraphAnnotateArgs { + /// Lift annotations to per-node block-consensus coordinates (lossless, long-format CSV). + Nodes(PangraphAnnotateNodesArgs), + + /// Compact node-level annotations into block-level consensus features (CSV). + Blocks(PangraphAnnotateBlocksArgs), +} + +/// Options shared by every `annotate` subcommand: the graph, the GFF inputs, and the output path. #[derive(Parser, Debug)] -pub struct PangraphAnnotateArgs { +pub struct AnnotateCommonArgs { /// Path to Pangraph JSON. /// /// Accepts plain or compressed file. If a compressed file is provided, it will be transparently @@ -24,15 +39,38 @@ pub struct PangraphAnnotateArgs { #[clap(display_order = 1)] pub input: Option, - /// Path to a GFF3 annotation file. Repeat the flag to provide multiple files. + /// Path(s) to GFF3 annotation file(s). + /// + /// Pass several files after a single flag (`--gff a.gff b.gff`, so shell globs like `--gff + /// *.gff` work), and/or repeat the flag (`--gff a.gff --gff b.gff`); the values accumulate. To + /// avoid the positional graph being slurped as an extra GFF, give it before the flag (`annotate + /// nodes graph.json --gff *.gff`) or pipe it via stdin. /// /// Accepts plain or compressed files (`gz`, `bz2`, `xz`, `zstd`), chosen by file extension. At /// least one file is required. Annotation `seqid`s must match the graph path names exactly. - #[clap(long = "gff", required = true, value_hint = ValueHint::FilePath)] + #[clap(long = "gff", required = true, num_args = 1.., value_hint = ValueHint::FilePath)] #[clap(display_order = 2)] pub gff: Vec, - /// Path to the output node-level annotation table (CSV). + /// Keep only annotations of these feature type(s) (GFF `type` column); drop all others. + /// + /// Comma-separated list (`--only-type gene,CDS`) and/or repeat the flag; values accumulate. + /// Matching is exact and case-sensitive (`CDS`, `gene`, `region`). Mutually exclusive with + /// `--exclude-type`. + #[clap(long = "only-type", value_delimiter = ',', value_hint = ValueHint::Other)] + #[clap(conflicts_with = "exclude_type")] + #[clap(display_order = 3)] + pub only_type: Vec, + + /// Drop annotations of these feature type(s) (GFF `type` column); keep all others. + /// + /// e.g. `--exclude-type region` removes whole-contig `region` declarations. Same comma-separated + /// syntax as `--only-type`; mutually exclusive with it. + #[clap(long = "exclude-type", value_delimiter = ',', value_hint = ValueHint::Other)] + #[clap(display_order = 4)] + pub exclude_type: Vec, + + /// Path to the output annotation table (CSV). /// /// Will be created if it does not exist. The output is compressed if the path ends in a known /// compression extension (`gz`, `bz2`, `xz`, `zstd`). Use `-` to write uncompressed CSV to @@ -41,3 +79,217 @@ pub struct PangraphAnnotateArgs { #[clap(value_hint = ValueHint::AnyPath)] pub output: PathBuf, } + +/// Arguments for `annotate nodes`: produce the lossless node-level table. +#[derive(Parser, Debug)] +pub struct PangraphAnnotateNodesArgs { + #[clap(flatten)] + pub common: AnnotateCommonArgs, +} + +/// Arguments for `annotate blocks`: compact the node-level table into block-level consensus features. +/// +/// The two thresholds mirror the fields of `CoordinateConsensusStrategy`; their defaults are kept in +/// sync with that type's `Default` (0.9 / 0.5). +#[derive(Parser, Debug)] +pub struct PangraphAnnotateBlocksArgs { + #[clap(flatten)] + pub common: AnnotateCommonArgs, + + /// Minimum frequency required to emit a block-level cluster. + /// + /// A cluster is kept when the number of supporting genomes `M >= ceil(min_frequency * N)`, where + /// `N` is the number of paths traversing the cluster's block(s) (so a gene is not penalised for + /// being absent in genomes that lack the block entirely). + #[clap(long, default_value_t = 0.9, value_parser = parse_fraction)] + #[clap(value_hint = ValueHint::Other)] + pub min_frequency: f64, + + /// Minimum supporter agreement required to promote a consensus name or attribute value. + /// + /// For each cluster, a `name`/attribute value is written only if at least this fraction of the + /// supporting genomes agree on it; otherwise the field is left empty. + #[clap(long, default_value_t = 0.5, value_parser = parse_fraction)] + #[clap(value_hint = ValueHint::Other)] + pub property_threshold: f64, +} + +/// Parse a threshold given as a fraction in the closed unit interval `[0, 1]`. +/// +/// Both `annotate blocks` thresholds are fractions; a value outside `[0, 1]` is always a mistake (it +/// would silently emit nothing, or promote every value), so it is rejected at parse time rather than +/// failing quietly downstream. `NaN` and infinities fall outside the range and are rejected too. +fn parse_fraction(s: &str) -> Result { + let value: f64 = s.parse().map_err(|err| format!("`{s}` is not a valid number: {err}"))?; + if (0.0..=1.0).contains(&value) { + Ok(value) + } else { + Err(format!( + "must be a fraction between 0 and 1 (inclusive), but got `{value}`" + )) + } +} + +#[cfg(test)] +mod tests { + use super::{AnnotateCommonArgs, parse_fraction}; + use crate::commands::root_args::{PangraphArgs, PangraphCommands}; + use clap::Parser; + use std::path::PathBuf; + + /// Parse a full `pangraph annotate nodes …` invocation and return its common args (graph, GFFs, + /// output), so the `--gff` parsing behaviour is exercised through clap rather than constructed by + /// hand. + fn parse_nodes_common(argv: &[&str]) -> Result { + let args = PangraphArgs::try_parse_from(argv)?; + match args.command { + PangraphCommands::Annotate { + args: super::PangraphAnnotateArgs::Nodes(nodes), + } => Ok(nodes.common), + other => panic!("expected `annotate nodes`, got {other:?}"), + } + } + + #[test] + fn gff_accepts_multiple_values_after_one_flag() { + let common = parse_nodes_common(&["pangraph", "annotate", "nodes", "graph.json", "--gff", "a.gff", "b.gff"]) + .expect("space-separated GFFs parse"); + assert_eq!(common.input, Some(PathBuf::from("graph.json"))); + assert_eq!(common.gff, vec![PathBuf::from("a.gff"), PathBuf::from("b.gff")]); + } + + #[test] + fn gff_flag_is_still_repeatable_and_accumulates() { + let common = parse_nodes_common(&[ + "pangraph", + "annotate", + "nodes", + "graph.json", + "--gff", + "a.gff", + "--gff", + "b.gff", + ]) + .expect("repeated GFF flags parse"); + assert_eq!(common.gff, vec![PathBuf::from("a.gff"), PathBuf::from("b.gff")]); + } + + #[test] + fn gff_list_does_not_swallow_the_positional_graph_or_output() { + // The graph given before the flag stays bound to the positional input, and a flag (`-o`) + // terminates the variadic so the GFF list does not absorb the output path. + let common = parse_nodes_common(&[ + "pangraph", + "annotate", + "nodes", + "graph.json", + "--gff", + "a.gff", + "b.gff", + "-o", + "out.csv", + ]) + .expect("graph-first invocation parses"); + assert_eq!(common.input, Some(PathBuf::from("graph.json"))); + assert_eq!(common.gff, vec![PathBuf::from("a.gff"), PathBuf::from("b.gff")]); + assert_eq!(common.output, PathBuf::from("out.csv")); + } + + #[test] + fn gff_is_required() { + parse_nodes_common(&["pangraph", "annotate", "nodes", "graph.json"]).unwrap_err(); + } + + #[test] + fn only_type_splits_on_commas() { + let common = parse_nodes_common(&[ + "pangraph", + "annotate", + "nodes", + "graph.json", + "--gff", + "a.gff", + "--only-type", + "gene,CDS", + ]) + .expect("comma-separated --only-type parses"); + assert_eq!(common.only_type, vec!["gene".to_owned(), "CDS".to_owned()]); + assert!(common.exclude_type.is_empty()); + } + + #[test] + fn type_filters_accumulate_across_repeated_flags() { + let common = parse_nodes_common(&[ + "pangraph", + "annotate", + "nodes", + "graph.json", + "--gff", + "a.gff", + "--only-type", + "gene", + "--only-type", + "CDS", + ]) + .expect("repeated --only-type accumulates"); + assert_eq!(common.only_type, vec!["gene".to_owned(), "CDS".to_owned()]); + } + + #[test] + fn exclude_type_parses() { + let common = parse_nodes_common(&[ + "pangraph", + "annotate", + "nodes", + "graph.json", + "--gff", + "a.gff", + "--exclude-type", + "region", + ]) + .expect("--exclude-type parses"); + assert_eq!(common.exclude_type, vec!["region".to_owned()]); + assert!(common.only_type.is_empty()); + } + + #[test] + fn only_type_and_exclude_type_conflict() { + parse_nodes_common(&[ + "pangraph", + "annotate", + "nodes", + "graph.json", + "--gff", + "a.gff", + "--only-type", + "gene", + "--exclude-type", + "region", + ]) + .expect_err("--only-type and --exclude-type are mutually exclusive"); + } + + #[test] + fn type_filters_default_to_empty() { + let common = parse_nodes_common(&["pangraph", "annotate", "nodes", "graph.json", "--gff", "a.gff"]) + .expect("no type filters parses"); + assert!(common.only_type.is_empty()); + assert!(common.exclude_type.is_empty()); + } + + #[test] + fn parse_fraction_accepts_closed_unit_interval() { + parse_fraction("0").unwrap(); + parse_fraction("0.5").unwrap(); + parse_fraction("1").unwrap(); + } + + #[test] + fn parse_fraction_rejects_out_of_range_and_non_numeric() { + parse_fraction("-0.01").unwrap_err(); + parse_fraction("1.01").unwrap_err(); + parse_fraction("NaN").unwrap_err(); + parse_fraction("inf").unwrap_err(); + parse_fraction("abc").unwrap_err(); + } +} diff --git a/packages/pangraph/src/commands/annotate/annotate_run.rs b/packages/pangraph/src/commands/annotate/annotate_run.rs index c3b82074..f9ce81a5 100644 --- a/packages/pangraph/src/commands/annotate/annotate_run.rs +++ b/packages/pangraph/src/commands/annotate/annotate_run.rs @@ -1,35 +1,79 @@ -use crate::annotation::lift::lift_features; +use crate::annotation::compact::{BlockCompactionStrategy, CoordinateConsensusStrategy}; +use crate::annotation::feature::filter_features_by_type; +use crate::annotation::lift::{LiftedAnnotation, lift_features}; use crate::annotation::matching::match_features_to_paths; use crate::annotation::writer::{AnnotationWriter, CsvAnnotationWriter}; -use crate::commands::annotate::annotate_args::PangraphAnnotateArgs; +use crate::commands::annotate::annotate_args::{ + AnnotateCommonArgs, PangraphAnnotateArgs, PangraphAnnotateBlocksArgs, PangraphAnnotateNodesArgs, +}; use crate::io::gff::GffReader; use crate::pangraph::pangraph::Pangraph; use eyre::{Report, WrapErr}; use std::collections::BTreeMap; -/// Run the `annotate` command: lift GFF features onto graph nodes and write the node-level table. +/// Run the `annotate` command, dispatching on the chosen output granularity. /// -/// Loads the graph, reads every GFF file into the internal `Feature` model, matches each feature's -/// `seqid` to a graph path (exact match; unmatched seqids are a hard error), lifts the features to -/// block-consensus coordinates, and writes the resulting node-level annotations as CSV. +/// Both modes share the same front half (load graph, read GFFs, match seqids to paths, lift to +/// block-consensus coordinates); they differ only in what is written out. pub fn annotate_run(args: PangraphAnnotateArgs) -> Result<(), Report> { - let PangraphAnnotateArgs { input, gff, output } = args; + match args { + PangraphAnnotateArgs::Nodes(args) => annotate_run_nodes(args), + PangraphAnnotateArgs::Blocks(args) => annotate_run_blocks(args), + } +} - let graph = Pangraph::from_path(&input)?; +/// Load the graph, read every GFF, match seqids to paths, and lift to node-level annotations. +/// +/// Returns the loaded graph alongside the lifted annotations, as block-level compaction needs the +/// graph to compute per-cluster path totals. +fn load_and_lift(common: &AnnotateCommonArgs) -> Result<(Pangraph, Vec), Report> { + let graph = Pangraph::from_path(&common.input)?; let mut features = Vec::new(); - for path in &gff { + for path in &common.gff { let read = GffReader::from_path(path)? .read_many() .wrap_err_with(|| format!("When reading GFF file: {}", path.display()))?; features.extend(read); } + // Drop unwanted feature types before matching, so excluded types never trigger seqid errors. + let features = filter_features_by_type(features, &common.only_type, &common.exclude_type); + let grouped = match_features_to_paths(features, &graph, &BTreeMap::new())?; let lifted = lift_features(&grouped, &graph)?; - let mut writer = CsvAnnotationWriter::new(&output, b',')?; + Ok((graph, lifted)) +} + +/// Run `annotate nodes`: write the lossless, long-format node-level table as CSV. +fn annotate_run_nodes(args: PangraphAnnotateNodesArgs) -> Result<(), Report> { + let PangraphAnnotateNodesArgs { common } = args; + let (_graph, lifted) = load_and_lift(&common)?; + + let mut writer = CsvAnnotationWriter::new(&common.output, b',')?; writer.write_node_annotations(&lifted)?; Ok(()) } + +/// Run `annotate blocks`: compact the node-level table into block-level consensus features (CSV). +fn annotate_run_blocks(args: PangraphAnnotateBlocksArgs) -> Result<(), Report> { + let PangraphAnnotateBlocksArgs { + common, + min_frequency, + property_threshold, + } = args; + let (graph, lifted) = load_and_lift(&common)?; + + let strategy = CoordinateConsensusStrategy { + min_frequency, + property_threshold, + }; + let blocks = strategy.compact(&lifted, &graph)?; + + let mut writer = CsvAnnotationWriter::new(&common.output, b',')?; + writer.write_block_annotations(&blocks)?; + + Ok(()) +} diff --git a/packages/pangraph/src/commands/main.rs b/packages/pangraph/src/commands/main.rs index 871d761e..1c1b9a74 100644 --- a/packages/pangraph/src/commands/main.rs +++ b/packages/pangraph/src/commands/main.rs @@ -21,7 +21,7 @@ pub fn pangraph_main() -> Result<(), Report> { PangraphCommands::Export { args } => export_run(args), PangraphCommands::Simplify(args) => simplify_run(args), PangraphCommands::Reconstruct(args) => reconstruct_run(&args), - PangraphCommands::Annotate(args) => annotate_run(args), + PangraphCommands::Annotate { args } => annotate_run(args), PangraphCommands::Schema(args) => generate_schema(&args), PangraphCommands::HelpMarkdown => print_help_markdown(), PangraphCommands::Completions { shell } => generate_shell_completions(&shell), diff --git a/packages/pangraph/src/commands/root_args.rs b/packages/pangraph/src/commands/root_args.rs index 989afa0d..49fddfa5 100644 --- a/packages/pangraph/src/commands/root_args.rs +++ b/packages/pangraph/src/commands/root_args.rs @@ -76,7 +76,10 @@ pub enum PangraphCommands { Reconstruct(PangraphReconstructArgs), /// Lift genome annotations onto the pangenome graph. - Annotate(PangraphAnnotateArgs), + Annotate { + #[clap(subcommand)] + args: PangraphAnnotateArgs, + }, /// Generate JSON schema for Pangraph file format Schema(PangraphGenerateSchemaArgs), diff --git a/packages/pangraph/tests/itest_annotate_cli.rs b/packages/pangraph/tests/itest_annotate_cli.rs index 8d01f09f..575c555e 100644 --- a/packages/pangraph/tests/itest_annotate_cli.rs +++ b/packages/pangraph/tests/itest_annotate_cli.rs @@ -1,7 +1,9 @@ #[cfg(test)] mod tests { use eyre::Report; - use pangraph::commands::annotate::annotate_args::PangraphAnnotateArgs; + use pangraph::commands::annotate::annotate_args::{ + AnnotateCommonArgs, PangraphAnnotateArgs, PangraphAnnotateBlocksArgs, PangraphAnnotateNodesArgs, + }; use pangraph::commands::annotate::annotate_run::annotate_run; use pangraph::io::csv::parse_csv; use pangraph::io::file::open_file_or_stdin; @@ -13,7 +15,8 @@ mod tests { use tempfile::tempdir; const GRAPH: &str = "../../data/test_graph.json"; - const CSV_HEADER_PREFIX: &str = "feature_id,parent_feature_id,segment_idx"; + const NODE_HEADER_PREFIX: &str = "feature_id,parent_feature_id,segment_idx"; + const BLOCK_HEADER_PREFIX: &str = "type,cluster_id,segment_idx"; /// Subset of the node-level CSV columns; remaining columns are ignored on deserialization. #[derive(Deserialize)] @@ -42,6 +45,70 @@ mod tests { Ok(path) } + /// Like [`write_gff`] but prepends a whole-contig `region` feature, so the type filters have a + /// third type (`region`) to keep or drop alongside the `gene` and `CDS`. + fn write_gff_with_region(dir: &Path, seqid: &str) -> Result { + let gff = format!( + "##gff-version 3\n\ + {seqid}\ttest\tregion\t1\t500\t.\t+\t.\tID=reg1;Name=r1\n\ + {seqid}\ttest\tgene\t1001\t1200\t.\t+\t.\tID=gene1;Name=g1\n\ + {seqid}\ttest\tCDS\t101\t5000\t.\t-\t0\tID=cds1;Name=c1\n" + ); + let path = dir.join("ann_region.gff"); + write(&path, gff)?; + Ok(path) + } + + /// Build `annotate nodes` args against the test graph for the given GFFs and output path. + fn nodes_args(gff: Vec, output: PathBuf) -> PangraphAnnotateArgs { + nodes_args_with_filter(gff, output, vec![], vec![]) + } + + /// Like [`nodes_args`], with explicit `--only-type` / `--exclude-type` filter lists. + fn nodes_args_with_filter( + gff: Vec, + output: PathBuf, + only_type: Vec, + exclude_type: Vec, + ) -> PangraphAnnotateArgs { + PangraphAnnotateArgs::Nodes(PangraphAnnotateNodesArgs { + common: AnnotateCommonArgs { + input: Some(PathBuf::from(GRAPH)), + gff, + only_type, + exclude_type, + output, + }, + }) + } + + /// Build `annotate blocks` args against the test graph; `property_threshold` is fixed at 0.5. + fn blocks_args(gff: Vec, output: PathBuf, min_frequency: f64) -> PangraphAnnotateArgs { + PangraphAnnotateArgs::Blocks(PangraphAnnotateBlocksArgs { + common: AnnotateCommonArgs { + input: Some(PathBuf::from(GRAPH)), + gff, + only_type: vec![], + exclude_type: vec![], + output, + }, + min_frequency, + property_threshold: 0.5, + }) + } + + /// Count block-level data rows in a written CSV: non-empty lines minus the header (the writer + /// omits the header entirely for a zero-row result, so an empty file counts as zero rows). + fn block_row_count(path: &Path) -> Result { + let contents = read_to_string(path)?; + let non_empty = contents.lines().filter(|l| !l.is_empty()).count(); + Ok(if contents.starts_with(BLOCK_HEADER_PREFIX) { + non_empty - 1 + } else { + non_empty + }) + } + #[test] fn itest_annotate_cli_node_level_csv() -> Result<(), Report> { let name = first_path_name()?; @@ -49,14 +116,10 @@ mod tests { let gff = write_gff(dir.path(), &name)?; let out = dir.path().join("lifted.csv"); - annotate_run(PangraphAnnotateArgs { - input: Some(PathBuf::from(GRAPH)), - gff: vec![gff], - output: out.clone(), - })?; + annotate_run(nodes_args(vec![gff], out.clone()))?; let contents = read_to_string(&out)?; - assert!(contents.starts_with(CSV_HEADER_PREFIX), "header row present"); + assert!(contents.starts_with(NODE_HEADER_PREFIX), "header row present"); let rows: Vec = parse_csv(&contents)?; assert!(!rows.is_empty(), "at least one lifted segment"); @@ -70,6 +133,104 @@ mod tests { Ok(()) } + #[test] + fn itest_annotate_cli_blocks_csv() -> Result<(), Report> { + let name = first_path_name()?; + let dir = tempdir()?; + let gff = write_gff(dir.path(), &name)?; + let out = dir.path().join("blocks.csv"); + + // `min_frequency: 0.0` keeps every cluster, so the single annotated genome still produces rows. + annotate_run(blocks_args(vec![gff], out.clone(), 0.0))?; + + let contents = read_to_string(&out)?; + assert!( + contents.starts_with(BLOCK_HEADER_PREFIX), + "block-level header row present" + ); + // At least one consensus feature emitted (header line plus one or more data lines). + assert!( + contents.lines().filter(|l| !l.is_empty()).count() >= 2, + "at least one block-level annotation row" + ); + Ok(()) + } + + #[test] + fn itest_annotate_cli_blocks_min_frequency_filters() -> Result<(), Report> { + // Only one genome is annotated, so every cluster has support `M = 1` against `N` (> 1) genomes + // traversing its block(s). The threshold must therefore actually reach the strategy: a permissive + // `min_frequency` keeps the clusters, while requiring unanimity drops them. (This also pins the + // wiring: were the two threshold fields swapped, the lenient run would already come back empty.) + let name = first_path_name()?; + let dir = tempdir()?; + let gff = write_gff(dir.path(), &name)?; + + let lenient = dir.path().join("lenient.csv"); + annotate_run(blocks_args(vec![gff.clone()], lenient.clone(), 0.0))?; + let lenient_rows = block_row_count(&lenient)?; + assert!(lenient_rows >= 1, "min_frequency 0.0 keeps single-support clusters"); + + let strict = dir.path().join("strict.csv"); + annotate_run(blocks_args(vec![gff], strict.clone(), 1.0))?; + let strict_rows = block_row_count(&strict)?; + + assert!( + strict_rows < lenient_rows, + "a higher --min-frequency emits strictly fewer clusters ({strict_rows} vs {lenient_rows})" + ); + Ok(()) + } + + #[test] + fn itest_annotate_cli_only_type_keeps_whitelisted() -> Result<(), Report> { + // `--only-type gene` keeps the gene and drops the CDS and the whole-contig region. + let name = first_path_name()?; + let dir = tempdir()?; + let gff = write_gff_with_region(dir.path(), &name)?; + let out = dir.path().join("only.csv"); + + annotate_run(nodes_args_with_filter( + vec![gff], + out.clone(), + vec!["gene".to_owned()], + vec![], + ))?; + + let rows: Vec = parse_csv(&read_to_string(&out)?)?; + assert!(!rows.is_empty(), "the gene survives the whitelist"); + assert!( + rows.iter().all(|r| r.parent_feature_id == "gene1"), + "only the gene is lifted; CDS and region are filtered out" + ); + Ok(()) + } + + #[test] + fn itest_annotate_cli_exclude_type_drops_blacklisted() -> Result<(), Report> { + // `--exclude-type region` removes the whole-contig declaration but keeps gene and CDS. + let name = first_path_name()?; + let dir = tempdir()?; + let gff = write_gff_with_region(dir.path(), &name)?; + let out = dir.path().join("exclude.csv"); + + annotate_run(nodes_args_with_filter( + vec![gff], + out.clone(), + vec![], + vec!["region".to_owned()], + ))?; + + let rows: Vec = parse_csv(&read_to_string(&out)?)?; + assert!(rows.iter().any(|r| r.parent_feature_id == "gene1"), "gene kept"); + assert!(rows.iter().any(|r| r.parent_feature_id == "cds1"), "cds kept"); + assert!( + rows.iter().all(|r| r.parent_feature_id != "reg1"), + "the excluded region produces no rows" + ); + Ok(()) + } + #[test] fn itest_annotate_cli_unmatched_seqid_errors() -> Result<(), Report> { // `example.gff` has seqids `chr1`/`chr2`, which match no path in the test graph: exact-match @@ -77,12 +238,7 @@ mod tests { let dir = tempdir()?; let out = dir.path().join("lifted.csv"); - let err = annotate_run(PangraphAnnotateArgs { - input: Some(PathBuf::from(GRAPH)), - gff: vec![PathBuf::from("../../data/example.gff")], - output: out, - }) - .unwrap_err(); + let err = annotate_run(nodes_args(vec![PathBuf::from("../../data/example.gff")], out)).unwrap_err(); let msg = err.to_string(); assert!(msg.contains("chr1"), "error names chr1: {msg}"); @@ -97,17 +253,13 @@ mod tests { let gff = write_gff(dir.path(), &name)?; let out = dir.path().join("lifted.csv.gz"); - annotate_run(PangraphAnnotateArgs { - input: Some(PathBuf::from(GRAPH)), - gff: vec![gff], - output: out.clone(), - })?; + annotate_run(nodes_args(vec![gff], out.clone()))?; // The `.gz` output is real gzip: read it back through transparent decompression. let mut decompressed = String::new(); open_file_or_stdin(&Some(&out))?.read_to_string(&mut decompressed)?; assert!( - decompressed.starts_with(CSV_HEADER_PREFIX), + decompressed.starts_with(NODE_HEADER_PREFIX), "decompressed header present" ); Ok(()) diff --git a/packages/pangraph/tests/itest_annotate_compact.rs b/packages/pangraph/tests/itest_annotate_compact.rs new file mode 100644 index 00000000..9b200435 --- /dev/null +++ b/packages/pangraph/tests/itest_annotate_compact.rs @@ -0,0 +1,119 @@ +mod common; + +#[cfg(test)] +mod tests { + use eyre::Report; + use pangraph::annotation::compact::{BlockCompactionStrategy, CoordinateConsensusStrategy}; + use pangraph::annotation::feature::Feature; + use pangraph::annotation::lift::lift_features; + use pangraph::annotation::matching::match_features_to_paths; + use pangraph::annotation::writer::{AnnotationWriter, CsvAnnotationWriter}; + use pangraph::pangraph::pangraph::Pangraph; + use pangraph::pangraph::strand::Strand; + use pangraph::utils::interval::Interval; + use pretty_assertions::assert_eq; + use std::collections::BTreeMap; + use std::fs::read_to_string; + use tempfile::tempdir; + + const GRAPH: &str = "../../data/test_graph.json"; + + /// A `CDS` feature covering the full genome extent `[p0, p1)` of a node. + fn whole_node_feature(seqid: &str, p0: usize, p1: usize) -> Feature { + let id = format!("{seqid}_core"); + Feature { + seqid: seqid.to_owned(), + source: None, + feature_type: "CDS".to_owned(), + interval: Interval::new(p0, p1), + strand: Some(Strand::Forward), + id: Some(id.clone()), + name: Some("core_gene".to_owned()), + attributes: vec![("ID".to_owned(), id), ("product".to_owned(), "core_product".to_owned())], + } + } + + /// End-to-end: placing a feature over the full extent of a core block's node on every genome lifts + /// each to the same consensus interval `[0, L)`, so compaction must collapse them into consensus + /// feature(s) for that block whose supporters sum to the number of genomes carrying it. + #[test] + fn itest_compact_core_block_collapses_across_genomes() -> Result<(), Report> { + let graph = Pangraph::from_path(&Some(GRAPH))?; + let n_paths = graph.paths.len(); + + // A core block whose nodes are all non-wrapping, so a whole-node feature is a clean interval. + let bid = graph + .core_block_ids() + .find(|bid| { + graph.blocks[bid].alignment_keys().iter().all(|nid| { + let (p0, p1) = graph.nodes[nid].position(); + p0 < p1 + }) + }) + .expect("a core block with only non-wrapping nodes"); + let block_len = graph.blocks[&bid].consensus_len(); + + // One whole-node feature per genome over that block. + let mut features = Vec::new(); + for path in graph.paths.values() { + let name = path.name().as_deref().unwrap(); + let nid = path + .nodes() + .iter() + .copied() + .find(|nid| graph.nodes[nid].block_id() == bid) + .expect("core block present on every path"); + let (p0, p1) = graph.nodes[&nid].position(); + features.push(whole_node_feature(name, p0, p1)); + } + + let grouped = match_features_to_paths(features, &graph, &BTreeMap::new())?; + let lifted = lift_features(&grouped, &graph)?; + + let strategy = CoordinateConsensusStrategy { + min_frequency: 0.0, + property_threshold: 0.5, + }; + let blocks = strategy.compact(&lifted, &graph)?; + + // Every emitted annotation is coherent. + for b in &blocks { + assert!(b.n_support >= 1 && b.n_support <= b.n_total, "support within total"); + } + + // The core-block placement(s): one per traversal strand, together supported by all genomes. + let for_block: Vec<_> = blocks.iter().filter(|b| b.block_id == bid).collect(); + assert!(!for_block.is_empty(), "core block produced a consensus annotation"); + let total_support: usize = for_block.iter().map(|b| b.n_support).sum(); + assert_eq!(total_support, n_paths, "every genome supports the core-block placement"); + for b in &for_block { + assert_eq!( + (b.cons_start, b.cons_end), + (0, block_len), + "whole node maps to whole consensus" + ); + assert_eq!(b.n_total, n_paths, "N is all genomes traversing the core block"); + assert_eq!(b.feature_type, "CDS"); + assert_eq!(b.consensus_name.as_deref(), Some("core_gene")); + assert!( + b.consensus_attributes + .contains(&("product".to_owned(), "core_product".to_owned())), + "shared product is promoted to consensus" + ); + } + + // The block-level writer round-trips the real output. + let dir = tempdir()?; + let out = dir.path().join("block_annotations.csv"); + { + let mut writer = CsvAnnotationWriter::new(&out, b',')?; + writer.write_block_annotations(&blocks)?; + } + let contents = read_to_string(&out)?; + assert!( + contents.starts_with("type,cluster_id,segment_idx"), + "block CSV header present" + ); + Ok(()) + } +}