diff --git a/.github/workflows/cli.yml b/.github/workflows/cli.yml index 7e2cc2b5..8a2796c5 100644 --- a/.github/workflows/cli.yml +++ b/.github/workflows/cli.yml @@ -387,7 +387,7 @@ jobs: - name: "Check that the git diff is clean" run: | - git -c color.ui=always diff --exit-code 'docs/docs/reference.md' || (echo "Looks like command-line interface has changed, but the autogenerated CLI reference documentation at 'docs/docs/reference.md' is not up-to-date. Please build the fresh version of pangraph, then run 'cd docs && ./generate-reference-docs docs/docs/reference.md', then verify and commit changes to the file docs/docs/reference.md." >&2; exit 1) + git -c color.ui=always diff --exit-code 'docs/docs/reference.md' || (echo "Looks like the command-line interface has changed, but the autogenerated CLI reference at 'docs/docs/reference.md' is not up-to-date. Please build pangraph, then from the repository root run 'cd docs && ./generate-reference-docs docs/reference.md' (the executable path is relative to docs/, e.g. ../target/release/pangraph), then verify and commit changes to the file docs/docs/reference.md." >&2; exit 1) publish-to-github-releases: diff --git a/data/klebs_annotations/NC_017540.gff.gz b/data/klebs_annotations/NC_017540.gff.gz index 14217c6b..5b9ed919 100644 Binary files a/data/klebs_annotations/NC_017540.gff.gz and b/data/klebs_annotations/NC_017540.gff.gz differ diff --git a/data/klebs_annotations/NZ_CP013711.gff.gz b/data/klebs_annotations/NZ_CP013711.gff.gz index 6cbd214b..54c280c5 100644 Binary files a/data/klebs_annotations/NZ_CP013711.gff.gz and b/data/klebs_annotations/NZ_CP013711.gff.gz differ diff --git a/dev/design/annotate-implementation-notes.md b/dev/design/annotate-implementation-notes.md index a1987952..579c8b58 100644 --- a/dev/design/annotate-implementation-notes.md +++ b/dev/design/annotate-implementation-notes.md @@ -13,10 +13,16 @@ | P1.5 | Real-data smoke tests (klebs graph + NCBI GFF annotations) | ✅ done | | P2 | Inverse coordinate helper `consensus_coords_from_node` (+ round-trip tests) | ✅ done | | P3 | Node-level lift + `AnnotationWriter` trait (CSV impl) | ✅ done | -| P4 | Block-level compaction (coordinate agreement) + JSON writer | ⏳ todo | -| P5 | `pangraph annotate` CLI command + `--seqid-map` + docs/CLI reference | ⏳ todo | +| 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 | +| P5.3 | Docs page + CLI reference regeneration (both output levels) | ⏳ todo | | P6 | pypangraph consumer + visualization example | ⏳ todo | +> **Ordering note (2026-06-12):** P4 (block compaction) is **deferred until after** a working +> node-level CLI (P5.1), so the still-undecided compaction policy can be designed against real +> node-level output. See the roadmap in [`annotate.md`](./annotate.md) §10. + ## Module layout - `packages/pangraph/src/annotation/` — domain logic @@ -102,10 +108,9 @@ coordinate. `lift_features(grouped, graph)` runs it over the `match_features_to_ - `node_start`/`node_end` are **consensus-oriented** node-local coordinates (already strand-flipped), i.e. exactly the input to `consensus_coords_from_node`; `cons_start`/`cons_end` are block-consensus coordinates. All half-open, 0-based. -- A feature interval is assumed **non-wrapping** (`f_s < f_e`); origin-spanning features are separate - GFF records. **Nodes**, however, may wrap the circular origin and that is handled - (`node_coverage_pieces` splits a wrapping node into its `[p0, tot_len)` and `[0, p1)` pieces, and a - whole-circle node `p0 == p1`). +- Both **nodes** and **features** may wrap the circular origin. Nodes are split by + `node_coverage_pieces` (`[p0, tot_len)` + `[0, p1)`, or a whole-circle node `p0 == p1`); features + by `feature_pieces` (see the origin-spanning section below). ### 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 @@ -142,6 +147,62 @@ origin-wrapping node in both strands) and asserting the reassembled segment base substring. Unit tests in `lift.rs` pin exact coordinates for each edge case (strand, deletion, insertion, in-insertion, multi-segment termini, circular wrap, unstranded). +## P5.1 decisions & behaviours (the `annotate` command) + +The `annotate` command (`packages/pangraph/src/commands/annotate/`, mirroring `simplify`) is pure +wiring over the P1–P3 library API — no lift-logic changes. `annotate_run` does: `Pangraph::from_path` +→ for each `--gff` file `GffReader::from_path(..).read_many()` → +`match_features_to_paths(.., &BTreeMap::new())` → `lift_features` → +`CsvAnnotationWriter::new(out, b',').write_node_annotations(..)`. + +### 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. +- **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 `,`. + +### seqid matching = exact (decision) +Matching passes an **empty** `seqid_map`, so `match_features_to_paths` matches `seqid == path.name` +exactly and aggregates **all** unmatched seqids into one hard error (existing P1 behaviour). On real +NCBI data (versioned `NZ_*.1` seqids vs bare `NZ_*` paths) this **errors by design** — the concrete +signal motivating the §12 relaxation decision. (The P1 matcher error text was reworded to drop the +"provide an explicit seqid-to-path mapping" suggestion — a capability the CLI does not expose yet — +and now simply states that seqids must match the path names.) + +### Tests +`packages/pangraph/tests/itest_annotate_cli.rs` drives `annotate_run` end-to-end against +`data/test_graph.json` with a GFF written at test time using a **real** path name (exercising the +`GffReader` path that `itest_annotate_lift.rs` skips): asserts the CSV header + that both features +lift + the genome column; a mismatch case (`data/example.gff`, seqids `chr1`/`chr2`) asserts the loud +error; and a `.csv.gz` output is read back through transparent decompression. + +## Origin-spanning features (circular paths) + +Real NCBI GFF encodes a feature crossing the replicon origin as a **single record with +`end > sequence length`** (e.g. `5278484..5279188` on a 5,278,493-bp chromosome → wraps to +`[0, 695)`), per +. +The P5.1 verification on real klebs data found exactly two such features (the chromosome's origin +gene + its CDS); every other feature was already byte-exact (20,935 / 20,937). + +**Behaviour (`lift.rs`):** +- A feature with `f_e > tot_len` on a **circular** path is decomposed by `feature_pieces` into its + arc pieces `[f_s, tot_len)` + `[0, f_e − tot_len)` (each tagged with an `arc_off`), intersected with + the node pieces, and emitted as a **single feature** (shared `parent_feature_id`). +- `merge_wrapped_segments` folds consecutive same-node, node-contiguous raw segments back together, so + a wrap landing on **one** origin-wrapping (or whole-circle) node yields **one** segment — not an + artificial head/tail split. A genuine multi-block wrap stays multiple segments, ordered 5'→3'. +- Terminus flags and `frac_covered` are computed from **arc position** along the feature + (`arc_start == 0` / `arc_end == len`), identical to the old genome-coordinate test for non-wrapping + features. +- **Skipped with a `warn!`** (returns no segments, so the user is told): a wrap on a **non-circular** + path, a feature **longer than the genome**, or a **start beyond the genome**. Warnings go through the + `log` crate, surfaced by `--verbosity`. + +Re-running the (throwaway) klebs verification after this change: **20,937 / 20,937 OK, 0 skipped**. + ## Public API introduced in P1–P3 ```rust @@ -187,15 +248,23 @@ data: accession (the FASTA record id), e.g. `NZ_CP013711`. - **`data/klebs_annotations/{NZ_CP013711,NC_017540}.gff.gz`** — RefSeq GFF annotations for two of those genomes (NCBI sviewer `report=gff3`). Only 2 kept to bound fixture size. -- Tests: parse each GFF (thousands of features, CDS/strand/name present) and **match both against - the graph**, building the seqid map from the files themselves. +- Tests: parse each GFF (thousands of features, CDS/strand/name present) and **match + lift both + against the graph** by **exact** seqid equality (no seqid map), asserting the lifted segments stay + in-bounds (P5.1 real-data lift smoke). -**Key real-world finding — seqid version mismatch.** Annotation seqids are the **versioned** -accession (`NZ_CP013711.1`); graph path names are the **bare** accession (`NZ_CP013711`). The smoke -test bridges this with a `version → bare` seqid map. This strongly motivates **version-insensitive -matching** (or a documented `--seqid-map`) in P5, since it is the default situation for NCBI data. +**Key real-world finding — seqid version mismatch.** As downloaded, annotation seqids are the +**versioned** accession (`NZ_CP013711.1`) while graph path names are the **bare** accession +(`NZ_CP013711`). The committed fixtures have since been **normalized** (the `.N` version stripped from +the seqid column only, leaving attribute values intact), so the smoke test matches by exact equality +without a map. The underlying NCBI reality still holds for *user* data, so it remains the motivation +for the §12 version-insensitive-matching decision — it is the default situation for NCBI downloads. ## Carried-forward items for later phases / user docs + +> Several of these are now consolidated into the **post-prototype, pre-docs punch list** in +> [`annotate.md`](./annotate.md) §12 (notably seqid↔path matching and the sequence-version-drift +> guard) — settle them there before P5.3. + - Document the `--seqid-map` file format (P5) and the "seqids must match FASTA record names" rule. - **Version-insensitive seqid matching (P5):** strongly consider auto-stripping the `.N` version so `NZ_CP013711.1` matches a `NZ_CP013711` path without an explicit map (see finding above). diff --git a/dev/design/annotate.md b/dev/design/annotate.md index 920d988c..3d248642 100644 --- a/dev/design/annotate.md +++ b/dev/design/annotate.md @@ -114,8 +114,12 @@ partial in the source GFF are carried through as non-termini too.) - **Deletions** — consensus positions absent from the node are skipped by the inverse map. - **Substitutions** — do not move coordinates (identity differs, position does not); irrelevant to the lift. - **Reverse strand** — offset flip + feature-strand flip (§4 step 3). -- **Circular wrap** — feature and/or node crossing the origin; use `PangraphPath.tot_len` + - `circular` and the modular pattern from `new_position_circular`. +- **Circular wrap** — both nodes and features may cross the origin. NCBI encodes an origin-spanning + feature as a single record with `end > tot_len` (it wraps into `[f_s, tot_len) ∪ [0, f_e−tot_len)`, + [NCBI ref](https://www.ncbi.nlm.nih.gov/datasets/docs/v2/reference-docs/file-formats/annotation-files/about-ncbi-gff3/#origin-spanning-features)). + On a **circular** path it is decomposed into arc pieces and lifted as a **single feature** — landing + on one origin-wrapping node it stays **one** segment (no head/tail split). A wrap on a non-circular + path, a feature longer than the genome, or a start beyond the genome is **skipped with a `warn!`**. - **Multi-block features** — split into segments with `parent_feature_id` + `segment_idx`; interior boundaries are fragment ends. - **Partial / truncated features** — recorded via the per-endpoint terminus flags and `frac_covered`. @@ -159,15 +163,19 @@ The hard part is largely solved by existing infrastructure: ## 7. Command interface -`pangraph annotate`: +`pangraph annotate` (P5.1 as-built; the surface grows in later phases): -- **Inputs**: a graph JSON + one or more GFF files. -- **Mode**: select node-level vs block-level output. -- **Output**: path (default `-` = stdout); compression inferred from the file extension; format - chosen by the `AnnotationWriter` (CSV default, JSON optional). -- **`--seqid-map`** escape hatch: map a GFF `seqid` to a PanGraph path name. - **ID matching is the #1 failure mode** — fail loudly (error, not silent drop) on any annotation - whose seqid does not match a path. +- **Inputs**: a graph JSON as the **positional** argument (stdin if omitted) + one or more GFF files + via a **repeatable `--gff`** flag. +- **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. +- **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 + **deferred** to the §12 post-prototype decision. Wire it by mirroring `simplify` (args/run), then register in `root_args.rs` + `main.rs`. @@ -197,24 +205,53 @@ one column; the **JSON writer** keeps it nested. Same `LiftedAnnotation` objects `AnnotationWriter` impl — the CSV writer renders one row per annotation, the JSON writer renders them nested. -### Block-level annotation (`BlockAnnotation`) — compacted +### Block-level annotation (`BlockAnnotation`) — compacted by coordinate consensus -One entry per (block, clustered feature). Carries a consensus-coordinate summary plus the -**coordinate-agreement** information across the genomes that carry the feature: +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. + +**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: + +``` +key = (start_block_id, cons_start, end_block_id, cons_end) +``` + +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. ``` -block_id, feature (cluster key: name/product), -cons_start_median, cons_start_min, cons_start_max, -cons_end_median, cons_end_min, cons_end_max, -n_support, # genomes carrying this feature on this block -n_total, # genomes traversing this block -coords_agree # all supporting genomes agree on start/end (easy case) vs disagree +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 ``` -Default output is **long-format CSV**; an **optional nested JSON** carries the full per-genome -coordinate distribution. Clustering "the same annotation" across genomes (by name/product + -consensus overlap) is **opinionated** — keep it a clearly documented, configurable layer **on top -of** the node-level lift, never baked into it. +Default output is **long-format CSV**; an **optional nested JSON** can carry the per-genome +supporters behind each consensus row. + +**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 +fixed policy. Natural settings, chosen by a CLI flag: **all** (`M == N`, only unanimous placements), +**a fraction** (`M >= ceil(f · N)` for a user fraction `f`), or **any** (`M == 1`, every distinct +placement). Exact-coordinate equality is the baseline identity test; whether to allow a small +coordinate **tolerance** (to absorb a terminus nudged by a nearby indel) is an open refinement. The +exact flag surface and fraction semantics are still **to be decided** — the P5.1 node-level prototype +on real data is meant to inform this choice (see §10). ## 9. Validation strategy @@ -239,10 +276,28 @@ Land incrementally, one PR/commit per phase: - **P3** ✅ — per-feature node-level lift (overlap, segmentation, strand, per-endpoint terminus flags) producing `LiftedAnnotation` objects + the `AnnotationWriter` trait with a CSV impl + integration test on `data/test_graph.json`. -- **P4** — block-level compaction (clustering + coordinate-agreement) into `BlockAnnotation` objects - + a second writer impl (CSV default, JSON optional) over the same trait + tests. -- **P5** — `annotate` command wiring (args/run/register/dispatch), Docusaurus docs page, CLI - reference regeneration. +The original "P4 then P5" ordering is **deliberately inverted** below: we land a working node-level +CLI **first** (P5.1) so we can run `pangraph annotate` on real data and use the actual node-level +output to inform the still-undecided block-compaction policy (P4). Block CLI wiring (P5.2) and a +single documentation pass (P5.3) follow once both output levels exist. Execution order is therefore +**P5.1 → P4 → P5.2 → P5.3 → P6**. + +- **P5.1** ✅ — `annotate` command wiring for the **node-level** lift (args/run/register/dispatch, + mirroring `simplify`): graph as the positional input, GFF(s) via a repeatable `--gff` flag, + `-o/--output` with transparent compression, CSV output. seqid→path matching is **exact** with a + hard error on any mismatch; `--seqid-map` and version-insensitive relaxation are **deferred** to + the §12 decision. A working prototype: the lossless node-level table on real graphs. Also the + **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). +- **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 visualization example. Out of scope for this Rust manifesto beyond this forward pointer. @@ -260,9 +315,53 @@ Land incrementally, one PR/commit per phase: 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** — by feature name/product, by consensus overlap, or by an - external ortholog grouping; make it configurable. +- **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. - **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 closed) to avoid off-by-one creep. + +## 12. Punch list — decide after the P5.1 prototype, before docs (P5.3) + +Improvements and design decisions to settle once the node-level prototype has been **run and tested +on real data**, but **before** the user-facing documentation is written — so the docs describe final +behaviour, not a moving target. These make the §11 open questions concrete; record the chosen +behaviour in the implementation notes as each is resolved. + +- **seqid ↔ path identity matching** — *how do we connect an annotation file's `seqid` to a graph + path name?* This is the **#1 user failure mode**: P1.5 found NCBI annotation seqids carry the + **versioned** accession (`NZ_CP013711.1`) while graph paths use the **bare** accession + (`NZ_CP013711`). **P5.1 ships the exact-match baseline** (hard error on any mismatch, no + `--seqid-map`); this item decides how — and whether — to relax it. Options, not mutually exclusive: + - **Exact match only** — require `seqid == path.name`, hard error otherwise. Most predictable and + explicit; breaks out-of-the-box on the common NCBI version-suffix case. + - **Relaxed / version-insensitive** — normalize before comparing, e.g. strip the trailing `.N` + accession version (possibly other canonicalizations) so `NZ_CP013711.1` matches `NZ_CP013711`. + Convenient default for NCBI data; must define behaviour when two seqids normalize to the **same** + path (ambiguity → error) and accept that it can mask a genuine wrong-file mismatch. + - **Explicit `--seqid-map`** — user-supplied mapping; always available as the escape hatch. + + *Decisions:* is relaxation **opt-out (default on)** or **opt-in (a flag)**? In what **order** are + the strategies tried (explicit map → relaxed → exact)? Keep the existing "aggregate all unmatched + seqids into one hard error" behaviour regardless. *Lean:* explicit map first, then + version-insensitive matching on by default, falling back to a loud aggregated error — but confirm + against what the prototype shows on real files. + +- **Sequence-version drift guard** — same accession but **different length** between the annotation's + source record and the graph's genome silently mis-lifts coordinates past the divergence point + (P1.5: `NZ_CP011582` 43433 bp vs graph 45279 bp). Decide whether `annotate` runs a length/identity + sanity check, and whether a mismatch is a **warning** or a **hard error**. + +- **Unstranded / partial features in the output** — confirm and document how `strand = None` (GFF + `.`/`?`) and any source-`partial` features render in the node-level table. + +- **Output ergonomics** — once real output has been eyeballed: decide the default column set/order, + whether the `attributes` JSON-string flattening is configurable, and whether a column selector is + worth adding. Only act on this if the default proves unwieldy in practice. diff --git a/docs/docs/reference.md b/docs/docs/reference.md index 24133e6d..e6bb1725 100644 --- a/docs/docs/reference.md +++ b/docs/docs/reference.md @@ -31,6 +31,7 @@ If you have Pangraph CLI installed, you can type `pangraph --help` to read the l * [`pangraph export core-genome`↴](#pangraph-export-core-genome) * [`pangraph simplify`↴](#pangraph-simplify) * [`pangraph reconstruct`↴](#pangraph-reconstruct) +* [`pangraph annotate`↴](#pangraph-annotate) * [`pangraph schema`↴](#pangraph-schema) * [`pangraph completions`↴](#pangraph-completions) * [`pangraph help-markdown`↴](#pangraph-help-markdown) @@ -58,6 +59,7 @@ Questions, ideas, bug reports: https://github.com/neherlab/pangraph/issues * `export` — Export a pangraph to a chosen file format(s) * `simplify` — Generates a simplified graph that only contains a subset of the input genomes * `reconstruct` — Reconstruct all input fasta sequences from graph +* `annotate` — Lift genome annotations onto the pangenome graph * `schema` — Generate JSON schema for Pangraph file format * `completions` — Generate shell completions * `help-markdown` — Print command-line reference documentation in Markdown format @@ -329,6 +331,33 @@ Reconstruct all input fasta sequences from graph +## `pangraph annotate` + +Lift genome annotations onto the pangenome graph + +**Usage:** `pangraph annotate [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 to a GFF3 annotation file. Repeat the flag to provide multiple files. + + 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). + + 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 schema` Generate JSON schema for Pangraph file format diff --git a/packages/pangraph/src/annotation/lift.rs b/packages/pangraph/src/annotation/lift.rs index eaab0b20..b2ad2e52 100644 --- a/packages/pangraph/src/annotation/lift.rs +++ b/packages/pangraph/src/annotation/lift.rs @@ -7,6 +7,7 @@ use crate::pangraph::pangraph_path::{PangraphPath, PathId}; use crate::pangraph::slice::consensus_coords_from_node_flagged; use crate::pangraph::strand::Strand; use eyre::Report; +use log::warn; use serde::{Deserialize, Serialize}; use std::collections::BTreeMap; @@ -86,15 +87,15 @@ pub struct LiftedAnnotation { pub attributes: Vec<(String, String)>, } -/// A feature's overlap with one node, in genome-oriented coordinates, before strand handling. +/// A feature's overlap with one node, before strand handling. struct RawSegment { node_id: NodeId, - /// Genome coordinates of the overlap `[genome_start, genome_end)`. - genome_start: usize, - genome_end: usize, /// Genome-oriented node-local coordinates of the overlap `[node_a, node_b)`. node_a: usize, node_b: usize, + /// Offset of this overlap along the feature (5'→3' from `f_s`), `[arc_start, arc_end)`. + arc_start: usize, + arc_end: usize, } /// Genome coverage piece(s) of a node on its path, each as `(genome_start, genome_end, @@ -121,13 +122,59 @@ fn node_coverage_pieces(p0: usize, p1: usize, tot_len: usize, circular: bool) -> pieces } -/// Build the per-segment `feature_id` from the parent id (or a coordinate-based fallback when -/// the source had no `ID`), suffixing multi-segment features so each row is unique. -fn make_feature_id(parent: Option<&str>, feature: &Feature, segment_idx: usize, n_segments: usize) -> String { - let base = parent.map_or_else( - || format!("{}:{}-{}", feature.seqid, feature.interval.start, feature.interval.end), - ToOwned::to_owned, - ); +/// Genome piece(s) of a feature `[f_s, f_e)` on its path, each as `(genome_start, genome_end, +/// arc_off)` where `arc_off` is the offset along the feature (5'→3' from `f_s`) at `genome_start`. +/// +/// A feature normally has a single non-wrapping piece. An **origin-spanning** feature on a circular +/// path is encoded by NCBI as one record with `f_e > tot_len`; it wraps into `[f_s, tot_len)` then +/// `[0, f_e - tot_len)`. Callers must validate circularity/length first (see `lift_feature`). +fn feature_pieces(f_s: usize, f_e: usize, tot_len: usize) -> Vec<(usize, usize, usize)> { + if f_e <= tot_len { + return vec![(f_s, f_e, 0)]; + } + vec![(f_s, tot_len, 0), (0, f_e - tot_len, tot_len - f_s)] +} + +/// Merge consecutive (arc-ordered) raw segments on the same node that are contiguous in node-local +/// coordinates. An origin wrap staying within a single (origin-wrapping or whole-circle) node must +/// remain one segment, not an artificial head/tail split at the origin. +fn merge_wrapped_segments(segments: Vec) -> Vec { + let mut merged: Vec = Vec::with_capacity(segments.len()); + for seg in segments { + let extend = merged.last().is_some_and(|last| { + let same_node = last.node_id == seg.node_id; + let abuts = last.node_b == seg.node_a; // prev node-local end meets next node-local start + same_node && abuts + }); + if extend { + let last = merged.last_mut().unwrap(); + // Node-contiguity (the `extend` check) implies arc-contiguity given how `feature_pieces` + // and `node_coverage_pieces` decompose the wrap; pin that invariant in debug builds. + debug_assert_eq!( + last.arc_end, seg.arc_start, + "merged wrapped segments must be arc-contiguous" + ); + last.node_b = seg.node_b; + last.arc_end = seg.arc_end; + } else { + merged.push(seg); + } + } + merged +} + +/// A human-readable label for a feature in log messages: its `ID`, else `seqid:start-end`. +fn feature_label(feature: &Feature) -> String { + feature + .id + .clone() + .unwrap_or_else(|| format!("{}:{}-{}", feature.seqid, feature.interval.start, feature.interval.end)) +} + +/// Build the per-segment `feature_id` (the [`feature_label`], suffixed for multi-segment features so +/// each row is unique). +fn make_feature_id(feature: &Feature, segment_idx: usize, n_segments: usize) -> String { + let base = feature_label(feature); if n_segments > 1 { format!("{base}.seg{segment_idx}") } else { @@ -137,10 +184,13 @@ fn make_feature_id(parent: Option<&str>, feature: &Feature, segment_idx: usize, /// Lift a single annotation [`Feature`] onto the nodes of the path it belongs to. /// -/// Returns one [`LiftedAnnotation`] per node the feature overlaps, in genome (low→high -/// coordinate) order. The feature interval is assumed to be a non-wrapping `[start, end)` -/// (origin-spanning features are encoded as separate GFF records); nodes, however, may wrap the -/// circular origin, which is handled here. +/// Returns one [`LiftedAnnotation`] per node the feature overlaps, ordered 5'→3' along the feature. +/// An **origin-spanning** feature on a *circular* path — encoded by NCBI as a single record whose +/// `end` exceeds the sequence length, so `[f_s, f_e)` wraps through the origin — is supported and +/// stays a single feature: if its wrapped span lands on one origin-wrapping node it yields one +/// segment, not an artificial head/tail split. Features that cannot be placed (a wrap on a +/// non-circular path, a feature longer than the genome, or a start beyond the genome) are skipped +/// with a `warn!` so the user is aware. pub fn lift_feature(feature: &Feature, path: &PangraphPath, graph: &Pangraph) -> Result, Report> { let tot_len = path.tot_len(); let circular = path.circular(); @@ -152,7 +202,37 @@ pub fn lift_feature(feature: &Feature, path: &PangraphPath, graph: &Pangraph) -> return Ok(vec![]); // empty feature: nothing to lift } - // 1+2. Collect the feature's overlap with each node, then order by genome coordinate. + // Validate placement; skip-and-warn on features we cannot put on this path. + if f_s >= tot_len { + warn!( + "Skipping feature {} on path {genome}: start {f_s} is beyond genome length {tot_len}", + feature_label(feature) + ); + return Ok(vec![]); + } + let len = f_e - f_s; + let wraps = f_e > tot_len; // NCBI encodes an origin-spanning feature as end > sequence length + if wraps { + if !circular { + warn!( + "Skipping origin-spanning feature {} on non-circular path {genome}: [{f_s}, {f_e}) exceeds genome length {tot_len}", + feature_label(feature) + ); + return Ok(vec![]); + } + if len > tot_len { + warn!( + "Skipping feature {} on path {genome}: length {len} exceeds genome length {tot_len}", + feature_label(feature) + ); + return Ok(vec![]); + } + } + + // 1+2. Collect the feature's overlap with each node — decomposing both the feature and the nodes + // into genome pieces so an origin wrap is handled on each side — ordered 5'→3' along the feature; + // then merge an origin wrap that stays within a single node back into one segment. + let f_pieces = feature_pieces(f_s, f_e, tot_len); let mut raw_segments: Vec = Vec::new(); for &node_id in path.nodes() { let node = graph @@ -160,25 +240,30 @@ pub fn lift_feature(feature: &Feature, path: &PangraphPath, graph: &Pangraph) -> .get(&node_id) .ok_or_else(|| make_internal_report!("When lifting feature: node {node_id} not found in graph"))?; let (p0, p1) = node.position(); - for (g_start, g_end, node_off) in node_coverage_pieces(p0, p1, tot_len, circular) { - let ov_s = f_s.max(g_start); - let ov_e = f_e.min(g_end); - if ov_s < ov_e { - raw_segments.push(RawSegment { - node_id, - genome_start: ov_s, - genome_end: ov_e, - node_a: node_off + (ov_s - g_start), - node_b: node_off + (ov_e - g_start), - }); + let node_pieces = node_coverage_pieces(p0, p1, tot_len, circular); + for &(fp_s, fp_e, arc_off) in &f_pieces { + for &(g_start, g_end, node_off) in &node_pieces { + let ov_s = fp_s.max(g_start); + let ov_e = fp_e.min(g_end); + if ov_s < ov_e { + raw_segments.push(RawSegment { + node_id, + node_a: node_off + (ov_s - g_start), + node_b: node_off + (ov_e - g_start), + arc_start: arc_off + (ov_s - fp_s), + arc_end: arc_off + (ov_e - fp_s), + }); + } } } } - raw_segments.sort_by_key(|s| s.genome_start); + raw_segments.sort_by_key(|s| s.arc_start); + if wraps { + raw_segments = merge_wrapped_segments(raw_segments); + } let n_segments = raw_segments.len(); - let feature_len = (f_e - f_s) as f64; - let parent = feature.id.as_deref(); + let feature_len = len as f64; let mut out = Vec::with_capacity(n_segments); for (segment_idx, seg) in raw_segments.into_iter().enumerate() { @@ -206,19 +291,19 @@ pub fn lift_feature(feature: &Feature, path: &PangraphPath, graph: &Pangraph) -> let ((cons_start, start_in_insertion), (cons_end, end_in_insertion)) = consensus_coords_from_node_flagged((node_start, node_end), edits, block_l); - // Terminus flags, in the consensus-endpoint frame. The genome-low feature end (f_s) lives in - // the segment whose genome overlap starts at f_s; the genome-high end (f_e) in the segment - // whose overlap ends at f_e. Reverse-strand nodes swap which consensus endpoint each maps to. - let holds_feature_start = seg.genome_start == f_s; - let holds_feature_end = seg.genome_end == f_e; + // Terminus flags, in the consensus-endpoint frame, by arc position along the feature: the + // feature's 5' end is in the segment at arc 0, its 3' end in the segment ending at arc `len`. + // Reverse-strand nodes swap which consensus endpoint each maps to. + let holds_feature_start = seg.arc_start == 0; + let holds_feature_end = seg.arc_end == len; let (start_is_terminus, end_is_terminus) = if reverse { (holds_feature_end, holds_feature_start) } else { (holds_feature_start, holds_feature_end) }; - let frac_covered = (seg.genome_end - seg.genome_start) as f64 / feature_len; - let feature_id = make_feature_id(parent, feature, segment_idx, n_segments); + let frac_covered = (seg.arc_end - seg.arc_start) as f64 / feature_len; + let feature_id = make_feature_id(feature, segment_idx, n_segments); out.push(LiftedAnnotation { feature_id, @@ -454,25 +539,7 @@ mod tests { fn test_lift_multi_node_segments_termini_and_frac() { // A feature spanning two forward nodes splits into two segments; only the outer endpoints are // termini and frac_covered sums to 1. - let (graph, path) = build_graph( - "p", - 20, - false, - vec![ - NodeSpec { - consensus: "ACGTACGTAC", - edits: Edit::empty(), - strand: Forward, - position: (0, 10), - }, - NodeSpec { - consensus: "TGCATGCATG", - edits: Edit::empty(), - strand: Forward, - position: (10, 20), - }, - ], - ); + let (graph, path) = two_node_forward_graph(false); let lifted = lift_feature(&feat(5, 15, Some(Forward)), &path, &graph).unwrap(); assert_eq!(lifted.len(), 2); @@ -575,6 +642,30 @@ mod tests { ) } + /// Path of length 20 split at the midpoint into two forward nodes: node 0 `(0,10)` and node 1 + /// `(10,20)`. `circular` toggles origin wrapping. + fn two_node_forward_graph(circular: bool) -> (Pangraph, PangraphPath) { + build_graph( + "p", + 20, + circular, + vec![ + NodeSpec { + consensus: "ACGTACGTAC", + edits: Edit::empty(), + strand: Forward, + position: (0, 10), + }, + NodeSpec { + consensus: "TGCATGCATG", + edits: Edit::empty(), + strand: Forward, + position: (10, 20), + }, + ], + ) + } + #[test] fn test_lift_unstranded_feature_stays_unstranded() { // An unstranded feature stays unstranded on the consensus, even on a reverse node. @@ -631,4 +722,94 @@ mod tests { let lifted = lift_features(&grouped, &graph).unwrap(); assert_eq!(lifted.len(), 2); } + + #[test] + fn test_lift_origin_spanning_single_wrapping_node_one_segment() { + // Origin-spanning feature [18,23) (end > tot_len) on a circular path lands entirely on the + // origin-wrapping node 0 -> a single merged segment, not a head/tail split at the origin. + let (graph, path) = wrapping_graph(); + let lifted = lift_feature(&feat(18, 23, Some(Forward)), &path, &graph).unwrap(); + assert_eq!(lifted.len(), 1); + let a = &lifted[0]; + assert_eq!(a.node_id, NodeId(0)); + assert_eq!((a.segment_idx, a.n_segments), (0, 1)); + assert_eq!(a.feature_id, "g1"); + assert_seg(a, (3, 8), (3, 8), Some(Forward), (true, true), (false, false)); + assert!((a.frac_covered - 1.0).abs() < 1e-9); + } + + #[test] + fn test_lift_origin_spanning_reverse_single_node_one_segment() { + // Same wrap on a reverse-strand origin-wrapping node: one segment, strand flipped, node coords + // measured from the other end (genome-oriented [3,8) -> reverse flip on len_node 10 -> [2,7)). + let (graph, path) = build_graph( + "p", + 20, + true, + vec![ + NodeSpec { + consensus: "ACGTACGTAC", + edits: Edit::empty(), + strand: Reverse, + position: (15, 5), + }, + NodeSpec { + consensus: "TGCATGCATG", + edits: Edit::empty(), + strand: Forward, + position: (5, 15), + }, + ], + ); + let lifted = lift_feature(&feat(18, 23, Some(Forward)), &path, &graph).unwrap(); + assert_eq!(lifted.len(), 1); + let a = &lifted[0]; + assert_eq!((a.segment_idx, a.n_segments), (0, 1)); + assert_seg(a, (2, 7), (2, 7), Some(Reverse), (true, true), (false, false)); + assert!((a.frac_covered - 1.0).abs() < 1e-9); + } + + #[test] + fn test_lift_origin_spanning_across_two_nodes_two_segments() { + // Origin at a node boundary: feature [18,23) wraps across node 1 (high side) then node 0 (low + // side) -> two segments in arc (5'->3') order; only the outer endpoints are termini. + let (graph, path) = two_node_forward_graph(true); + let lifted = lift_feature(&feat(18, 23, Some(Forward)), &path, &graph).unwrap(); + assert_eq!(lifted.len(), 2); + + let s0 = &lifted[0]; + assert_eq!(s0.node_id, NodeId(1)); + assert_eq!((s0.segment_idx, s0.n_segments), (0, 2)); + assert_eq!(s0.feature_id, "g1.seg0"); + assert_seg(s0, (8, 10), (8, 10), Some(Forward), (true, false), (false, false)); + assert!((s0.frac_covered - 0.4).abs() < 1e-9); + + let s1 = &lifted[1]; + assert_eq!(s1.node_id, NodeId(0)); + assert_eq!((s1.segment_idx, s1.n_segments), (1, 2)); + assert_eq!(s1.feature_id, "g1.seg1"); + assert_seg(s1, (0, 3), (0, 3), Some(Forward), (false, true), (false, false)); + assert!((s1.frac_covered - 0.6).abs() < 1e-9); + } + + #[test] + fn test_lift_wrapping_feature_on_noncircular_path_is_skipped() { + let (graph, path) = two_node_forward_graph(false); + let lifted = lift_feature(&feat(18, 23, Some(Forward)), &path, &graph).unwrap(); + assert!(lifted.is_empty(), "wrap on a non-circular path is skipped"); + } + + #[test] + fn test_lift_feature_longer_than_genome_is_skipped() { + let (graph, path) = wrapping_graph(); // tot_len 20, circular + let lifted = lift_feature(&feat(5, 30, Some(Forward)), &path, &graph).unwrap(); + assert!(lifted.is_empty(), "feature longer than the genome is skipped"); + } + + #[test] + fn test_lift_feature_start_beyond_genome_is_skipped() { + let (graph, path) = wrapping_graph(); // tot_len 20 + let lifted = lift_feature(&feat(25, 28, Some(Forward)), &path, &graph).unwrap(); + assert!(lifted.is_empty(), "start beyond the genome is skipped"); + } } diff --git a/packages/pangraph/src/annotation/matching.rs b/packages/pangraph/src/annotation/matching.rs index ed0ec08b..15f83f7e 100644 --- a/packages/pangraph/src/annotation/matching.rs +++ b/packages/pangraph/src/annotation/matching.rs @@ -41,8 +41,8 @@ pub fn match_features_to_paths( let names = unmatched.iter().cloned().collect::>().join(", "); return make_error!( "Could not match {} annotation sequence id(s) to any pangraph path: {names}. \ - Check that annotation seqids correspond to the FASTA record names used to build the \ - graph, or provide an explicit seqid-to-path mapping.", + Annotation seqids must exactly match the FASTA record names (the pangraph path names) \ + used to build the graph.", unmatched.len() ); } diff --git a/packages/pangraph/src/commands/annotate/annotate_args.rs b/packages/pangraph/src/commands/annotate/annotate_args.rs new file mode 100644 index 00000000..613657f6 --- /dev/null +++ b/packages/pangraph/src/commands/annotate/annotate_args.rs @@ -0,0 +1,43 @@ +use clap::{Parser, 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. +/// +/// 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(Parser, Debug)] +pub struct PangraphAnnotateArgs { + /// 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). + #[clap(value_hint = ValueHint::FilePath)] + #[clap(display_order = 1)] + pub input: Option, + + /// Path to a GFF3 annotation file. Repeat the flag to provide multiple files. + /// + /// 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(display_order = 2)] + pub gff: Vec, + + /// Path to the output node-level 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). + #[clap(long, short = 'o', default_value = "-")] + #[clap(value_hint = ValueHint::AnyPath)] + pub output: PathBuf, +} diff --git a/packages/pangraph/src/commands/annotate/annotate_run.rs b/packages/pangraph/src/commands/annotate/annotate_run.rs new file mode 100644 index 00000000..c3b82074 --- /dev/null +++ b/packages/pangraph/src/commands/annotate/annotate_run.rs @@ -0,0 +1,35 @@ +use crate::annotation::lift::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::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. +/// +/// 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. +pub fn annotate_run(args: PangraphAnnotateArgs) -> Result<(), Report> { + let PangraphAnnotateArgs { input, gff, output } = args; + + let graph = Pangraph::from_path(&input)?; + + let mut features = Vec::new(); + for path in &gff { + let read = GffReader::from_path(path)? + .read_many() + .wrap_err_with(|| format!("When reading GFF file: {}", path.display()))?; + features.extend(read); + } + + let grouped = match_features_to_paths(features, &graph, &BTreeMap::new())?; + let lifted = lift_features(&grouped, &graph)?; + + let mut writer = CsvAnnotationWriter::new(&output, b',')?; + writer.write_node_annotations(&lifted)?; + + Ok(()) +} diff --git a/packages/pangraph/src/commands/annotate/mod.rs b/packages/pangraph/src/commands/annotate/mod.rs new file mode 100644 index 00000000..2fab28ee --- /dev/null +++ b/packages/pangraph/src/commands/annotate/mod.rs @@ -0,0 +1,2 @@ +pub mod annotate_args; +pub mod annotate_run; diff --git a/packages/pangraph/src/commands/main.rs b/packages/pangraph/src/commands/main.rs index 4d94f27c..871d761e 100644 --- a/packages/pangraph/src/commands/main.rs +++ b/packages/pangraph/src/commands/main.rs @@ -1,3 +1,4 @@ +use crate::commands::annotate::annotate_run::annotate_run; use crate::commands::build::build_run::build_run; use crate::commands::export::export_run::export_run; use crate::commands::md_help::print_help_markdown::print_help_markdown; @@ -20,6 +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::Schema(args) => generate_schema(&args), PangraphCommands::HelpMarkdown => print_help_markdown(), PangraphCommands::Completions { shell } => generate_shell_completions(&shell), diff --git a/packages/pangraph/src/commands/mod.rs b/packages/pangraph/src/commands/mod.rs index 0c11d129..20b72ba8 100644 --- a/packages/pangraph/src/commands/mod.rs +++ b/packages/pangraph/src/commands/mod.rs @@ -1,3 +1,4 @@ +pub mod annotate; pub mod build; pub mod export; pub mod main; diff --git a/packages/pangraph/src/commands/root_args.rs b/packages/pangraph/src/commands/root_args.rs index 5020a3cf..989afa0d 100644 --- a/packages/pangraph/src/commands/root_args.rs +++ b/packages/pangraph/src/commands/root_args.rs @@ -1,5 +1,6 @@ #![allow(unused_qualifications)] +use crate::commands::annotate::annotate_args::PangraphAnnotateArgs; use crate::commands::build::build_args::PangraphBuildArgs; use crate::commands::export::export_args::PangraphExportArgs; use crate::commands::reconstruct::reconstruct_args::PangraphReconstructArgs; @@ -74,6 +75,9 @@ pub enum PangraphCommands { /// Reconstruct all input fasta sequences from graph Reconstruct(PangraphReconstructArgs), + /// Lift genome annotations onto the pangenome graph. + Annotate(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 new file mode 100644 index 00000000..8d01f09f --- /dev/null +++ b/packages/pangraph/tests/itest_annotate_cli.rs @@ -0,0 +1,115 @@ +#[cfg(test)] +mod tests { + use eyre::Report; + use pangraph::commands::annotate::annotate_args::PangraphAnnotateArgs; + use pangraph::commands::annotate::annotate_run::annotate_run; + use pangraph::io::csv::parse_csv; + use pangraph::io::file::open_file_or_stdin; + use pangraph::pangraph::pangraph::Pangraph; + use serde::Deserialize; + use std::fs::{read_to_string, write}; + use std::io::Read; + use std::path::{Path, PathBuf}; + use tempfile::tempdir; + + const GRAPH: &str = "../../data/test_graph.json"; + const CSV_HEADER_PREFIX: &str = "feature_id,parent_feature_id,segment_idx"; + + /// Subset of the node-level CSV columns; remaining columns are ignored on deserialization. + #[derive(Deserialize)] + struct Row { + genome: String, + parent_feature_id: String, + } + + /// Name of the first path in the test graph (a valid `seqid` to match against). + fn first_path_name() -> Result { + let graph = Pangraph::from_path(&Some(GRAPH))?; + let path = graph.paths.values().next().expect("at least one path"); + Ok(path.name().as_deref().unwrap().to_owned()) + } + + /// Write a small GFF3 file (1-based inclusive coords) with two features on `seqid`: a short + /// forward gene (single node) and a long reverse CDS (spans block boundaries). + fn write_gff(dir: &Path, seqid: &str) -> Result { + let gff = format!( + "##gff-version 3\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.gff"); + write(&path, gff)?; + Ok(path) + } + + #[test] + fn itest_annotate_cli_node_level_csv() -> Result<(), Report> { + let name = first_path_name()?; + let dir = tempdir()?; + 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(), + })?; + + let contents = read_to_string(&out)?; + assert!(contents.starts_with(CSV_HEADER_PREFIX), "header row present"); + + let rows: Vec = parse_csv(&contents)?; + assert!(!rows.is_empty(), "at least one lifted segment"); + // Both source features were lifted (parent_feature_id carries the GFF ID). + assert!(rows.iter().any(|r| r.parent_feature_id == "gene1"), "gene1 lifted"); + assert!(rows.iter().any(|r| r.parent_feature_id == "cds1"), "cds1 lifted"); + // Every row belongs to the genome we annotated. + for r in &rows { + assert_eq!(r.genome, name, "genome column is the path name"); + } + 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 + // policy must fail loudly and name every offending seqid. + 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 msg = err.to_string(); + assert!(msg.contains("chr1"), "error names chr1: {msg}"); + assert!(msg.contains("chr2"), "error names chr2: {msg}"); + Ok(()) + } + + #[test] + fn itest_annotate_cli_gz_output_roundtrips() -> Result<(), Report> { + let name = first_path_name()?; + let dir = tempdir()?; + 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(), + })?; + + // 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 header present" + ); + Ok(()) + } +} diff --git a/packages/pangraph/tests/itest_klebs_annotations.rs b/packages/pangraph/tests/itest_klebs_annotations.rs index 0f1d5875..2605b867 100644 --- a/packages/pangraph/tests/itest_klebs_annotations.rs +++ b/packages/pangraph/tests/itest_klebs_annotations.rs @@ -3,17 +3,19 @@ mod common; #[cfg(test)] mod tests { use eyre::Report; + use pangraph::annotation::lift::lift_features; use pangraph::annotation::matching::match_features_to_paths; use pangraph::io::gff::GffReader; use pangraph::pangraph::pangraph::Pangraph; use pangraph::pangraph::strand::Strand; use rstest::rstest; - use std::collections::BTreeMap; + use std::collections::{BTreeMap, BTreeSet}; // Real Klebsiella annotations (downloaded from NCBI) for two of the genomes in - // `data/klebs_graph.json.gz`. The path name is the bare accession (the FASTA record id), - // while the annotation files carry the versioned accession as their seqid. - // path name (bare) gff path + // `data/klebs_graph.json.gz`. The GFF seqid column has been normalized to the bare accession (the + // FASTA record id / pangraph path name) by stripping the `.N` version, so annotations match the + // graph paths by exact string equality — no seqid map needed. + // bare accession gff path const KLEBS: &[(&str, &str)] = &[ ("NZ_CP013711", "../../data/klebs_annotations/NZ_CP013711.gff.gz"), ("NC_017540", "../../data/klebs_annotations/NC_017540.gff.gz"), @@ -23,7 +25,7 @@ mod tests { #[case(0)] #[case(1)] fn smoke_parse_real_klebs_gff(#[case] i: usize) -> Result<(), Report> { - let (_, gff) = KLEBS[i]; + let (name, gff) = KLEBS[i]; let features = GffReader::from_path(gff)?.read_many()?; assert!( features.len() > 1000, @@ -33,30 +35,45 @@ mod tests { assert!(features.iter().any(|f| f.feature_type == "CDS")); assert!(features.iter().any(|f| f.strand == Some(Strand::Reverse))); assert!(features.iter().any(|f| f.name.is_some())); + // Seqids are the bare accession (version stripped), matching the graph path name exactly. + assert!( + features.iter().all(|f| f.seqid == name), + "all seqids should be the bare accession {name}" + ); Ok(()) } #[test] - fn smoke_match_klebs_annotations_to_graph() -> Result<(), Report> { + fn smoke_lift_klebs_annotations_to_graph() -> Result<(), Report> { let graph = Pangraph::from_path(&Some("../../data/klebs_graph.json.gz"))?; let mut all_features = Vec::new(); - let mut seqid_map = BTreeMap::new(); - for (name, gff) in KLEBS { - let features = GffReader::from_path(gff)?.read_many()?; - // The annotation seqid is the versioned accession (e.g. `NZ_CP013711.1`); map it onto - // the bare-accession path name used when the graph was built. - let seqid = features.first().expect("annotation has features").seqid.clone(); - seqid_map.insert(seqid, (*name).to_owned()); - all_features.extend(features); + for (_, gff) in KLEBS { + all_features.extend(GffReader::from_path(gff)?.read_many()?); } - let grouped = match_features_to_paths(all_features, &graph, &seqid_map)?; + // Exact seqid -> path matching (empty seqid map): both annotated genomes group cleanly. + let grouped = match_features_to_paths(all_features, &graph, &BTreeMap::new())?; assert_eq!(grouped.len(), 2); for (name, _) in KLEBS { let pid = graph.path_id_by_name(name)?; assert!(grouped.get(&pid).is_some_and(|features| !features.is_empty())); } + + // Lifting thousands of real features onto the real graph succeeds and stays in-bounds. + let lifted = lift_features(&grouped, &graph)?; + assert!( + lifted.len() > 1000, + "expected many lifted segments, got {}", + lifted.len() + ); + let names: BTreeSet<&str> = KLEBS.iter().map(|(n, _)| *n).collect(); + for ann in &lifted { + assert!(names.contains(ann.genome.as_str()), "unexpected genome {}", ann.genome); + assert!(ann.cons_start <= ann.cons_end, "consensus coords ordered"); + let block = graph.blocks.get(&ann.block_id).expect("block exists"); + assert!(ann.cons_end <= block.consensus_len(), "consensus end within block"); + } Ok(()) } }