fix: verify reconstructed sequences by genome name

dev/design/merge.md

@@ -1,6 +1,6 @@
 # Graph merging — design manifesto
 
-Status: **`pangraph merge` implemented; verification and `build`-side name checks still pending (see §6)**
+Status: **`pangraph merge` implemented and verified; only user docs remain (see §6)**
 Integration branch: `feat/merge` (merged into `master` last, after all phase branches below)
 
 This document describes the design for a new `pangraph merge` command, which combines two
@@ -301,21 +301,39 @@ Reconstruction itself remains fully supported — every sequence comes back byte
 
 This has three consequences:
 
-1. **`compare_sequences` must stop comparing whole records.** It currently tests
-   `left != right` on `FastaRecord` (`reconstruct_run.rs:45-54`), and `FastaRecord` derives
-   `PartialEq` over all fields *including `index`* (`io/fasta.rs:17-24`) — despite an error message
-   that only mentions length. It should compare sequence contents and report the path name.
+1. **`compare_sequences` must stop comparing whole records.** It tested `left != right` on
+   `FastaRecord`, which derives `PartialEq` over all fields *including `index`* (`io/fasta.rs:17-24`)
+   — despite an error message that only mentioned length. It is now deleted; `verify_genome` compares
+   sequence contents and reports the genome name and the first differing position.
+
+   This was not only cosmetic. Two bugs followed from index/position pairing:
+   - `build --verify` indexed `&fastas[actual.index]`, which **panicked out of bounds** whenever the
+     input records' indices were not exactly `0..n-1` — as they are not for any programmatically
+     assembled record set, including the repo's own test helper.
+   - `reconstruct --verify` paired records positionally, so surplus records in the verification file
+     were silently ignored (exit 0), too few gave a misleading `expected length 0 but got N`, and a
+     merged graph failed spuriously with `expected length N but got N` at identical lengths.
 
 2. **Verification is keyed by path name, not by index or position.**
 
    ```rust
-   fn verify_graph_sequences(graph: &Pangraph, expected: &BTreeMap<String, Seq>) -> Result<(), Report>
+   pub fn verify_graph_sequences(graph: &Pangraph, expected: &BTreeMap<String, Seq>,
+                                 coverage: GenomeCoverage) -> Result<(), Report>
    ```
 
+   The third parameter was not in the original sketch: `build`'s intermediate clade graphs hold only
+   the genomes of their own clade, so they need `GenomeCoverage::Partial`, while the final graph of a
+   build, a merged graph and `reconstruct --verify` all require `Complete`. Making the *final* build
+   check `Complete` is a small gain — nothing previously noticed a genome going missing.
+
    `build` fills the map from the input FASTA records; `merge` fills it by reconstructing each
-   input graph *before* merging. Both are now sound because names are unique (§4.3). The map form
-   also keeps working for the intermediate graphs checked inside the build loop, which contain only
-   a subset of the paths.
+   input graph *before* merging. Both are sound because names are unique (§4.3). `reconstruct
+   --verify` does **not** build the map: it is the only one of the three that does not otherwise
+   need every genome resident, so it streams the verification file against a name → path id index
+   and reconstructs one genome at a time, sharing `verify_genome` and its error formatting.
+
+   Reconstruction itself moved out of the command module into `pangraph/reconstruct.rs`, since it is
+   a graph operation — `pangraph/pangraph.rs` was reaching up into `commands::`.
 
 3. **`pangraph reconstruct` documentation must state** that record order matches the original input
    FASTA order only for graphs produced directly by `build`, and that consumers should match records
@@ -378,16 +396,13 @@ merged into `master` last, once all phases have landed.
 | 3 | `feat/merge-cmd` | §3.4 — `relabel` / `make_disjoint_from` | landed |
 | 4 | `feat/merge-cmd` | §4.5 — the command itself, plus integration tests | landed |
 | 5 | `feat/merge-cmd` | §4.3 — duplicate genome names are an error in `build` and `merge` | landed |
-| 6 | `feat/merge-verify` | §4.4 — name-keyed verification shared with `build` | todo |
+| 6 | `feat/merge-verify` | §4.4 — name-keyed verification shared with `build` | landed |
 | 7 | `feat/merge-docs` | §9 — tutorial, `reconstruct` docs, CHANGELOG | todo |
 
 Phases 2–5 were implemented together, since a `merge` command without §3.4 panics on the first
 identifier collision and would not be testable.
 
-§4.4 is therefore the only piece left half-done: name-keyed verification exists as
-`verify_merged_sequences` inside `merge_run`, private to the merge command. `build` still verifies
-against input FASTA records by index, and `compare_sequences` still compares whole `FastaRecord`s
-(including `index`). Phase 6 unifies the two behind a single `BTreeMap<String, Seq>`-based helper.
+Phase 6 turned out to fix two latent bugs rather than merely unify style, both recorded in §4.4.
 
 ---
 

docs/docs/reference.md

@@ -370,20 +370,26 @@ Reconstruct all input fasta sequences from graph
 
 * `<INPUT_GRAPH>` — Path to a pangenome graph file in JSON format.
 
-   Accepts plain or compressed FASTA files. If a compressed fasta file is provided, it will be transparently decompressed. Supported compression formats: `gz`, `bz2`, `xz`, `zstd`. Decompressor is chosen based on file extension. If there's multiple input files, then different files can have different compression formats.
+   Accepts plain or compressed files. 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 files provided, the plain fasta input is read from standard input (stdin).
+   If no input file is provided, the plain JSON input is read from standard input (stdin).
 
 ###### **Options:**
 
 * `-o`, `--output-fasta <OUTPUT_FASTA>` — Path to output FASTA file with reconstructed sequences.
 
    If the provided file path ends with one of the supported extensions: "gz", "bz2", "xz", "zst", then the file will be written compressed. If the required directory tree does not exist, it will be created.
 
-   Use "-" to write the uncompressed data to standard output (stdout). This is the default, if the argument is not provided. See: https://en.wikipedia.org/wiki/FASTA_format
+   Use "-" to write the uncompressed data to standard output (stdout). This is the default, if the argument is not provided.
+
+   Records are written in order of path id, which reproduces the order of the original input FASTA only for graphs produced directly by `pangraph build`. A graph produced by `pangraph merge` renumbers its path ids, so consumers should match records by genome name rather than by position.
+
+   See: https://en.wikipedia.org/wiki/FASTA_format
 
   Default value: `-`
-* `-f`, `--verify <VERIFY>` — Path to the FASTA file with sequences to check the reconstructed sequences against. If this argument is provided, then the sequences are not being printed to standard output (stdout) as usual. Instead, if any differences are detected, a diff will be printed between the expected (original) sequence and reconstructed sequence.
+* `-f`, `--verify <VERIFY>` — Path to the FASTA file with sequences to check the reconstructed sequences against. If this argument is provided, then the sequences are not written out as usual: nothing is produced on success, and the first difference found is reported as an error.
+
+   Genomes are matched by name, so the order of the records is irrelevant. The file must contain exactly the genomes of the graph: a record the graph does not contain, a genome missing from the file, or a repeated name are all errors. Every path of the graph must be named.
 
    Accepts plain or compressed FASTA files. If a compressed fasta file is provided, it will be transparently decompressed. Supported compression formats: `gz`, `bz2`, `xz`, `zstd`. Decompressor is chosen based on file extension. If there's multiple input files, then different files can have different compression formats.
 

packages/pangraph/src/commands/build/build_run.rs

@@ -1,11 +1,12 @@
 use crate::align::alignment_args::check_alignment_backend_available;
 use crate::commands::build::build_args::PangraphBuildArgs;
-use crate::commands::reconstruct::reconstruct_run::{compare_sequences, reconstruct};
 use crate::io::fasta::{FastaReader, FastaRecord};
 use crate::io::json::{JsonPretty, json_write_file};
 use crate::pangraph::graph_merging::merge_graphs;
 use crate::pangraph::pangraph::Pangraph;
+use crate::pangraph::reconstruct::{GenomeCoverage, sequences_by_name, verify_graph_sequences};
 use crate::pangraph::strand::Strand::Forward;
+use crate::representation::seq::Seq;
 use crate::tree::clade::postorder;
 use crate::tree::neighbor_joining::build_tree_using_neighbor_joining;
 use crate::tree::newick::build_tree_from_newick;
@@ -15,31 +16,23 @@ use crate::{make_error, make_internal_error, make_internal_report};
 use eyre::{Report, WrapErr};
 use itertools::Itertools;
 use log::info;
-
-pub fn reconstruct_and_compare_graph_seqs(graph: &Pangraph, fastas: &[FastaRecord]) -> Result<(), Report> {
-  // Reconstruct sequences from the given graph.
-  let mut results = reconstruct(graph);
-
-  // Check that the reconstructed sequences match the original FASTA records.
-  results.try_for_each(|actual| -> Result<(), Report> {
-    let actual = actual?;
-    let expected = &fastas[actual.index];
-    compare_sequences(expected, &actual)?;
-    Ok(())
-  })?;
-
-  Ok(())
-}
-
-pub fn graph_sanity_checks(graph: &Pangraph, fastas: &[FastaRecord]) -> Result<(), Report> {
+use std::collections::BTreeMap;
+
+/// Checks that the graph is internally consistent and reconstructs the genomes it should.
+fn graph_sanity_checks(
+  graph: &Pangraph,
+  expected: &BTreeMap<String, Seq>,
+  coverage: GenomeCoverage,
+) -> Result<(), Report> {
   // check that graph internal structure (blocks, paths, nodes, edits...) is valid
   #[cfg(debug_assertions)]
   graph
     .sanity_check()
     .wrap_err("When performing sanity check on the pangraph")?;
 
-  // Reconstruct sequences from the graph and compare them with the original FASTA records.
-  reconstruct_and_compare_graph_seqs(graph, fastas)
+  // Reconstruct the genomes from the graph and compare them with the input sequences. Genomes are
+  // matched by name: path ids do not survive a merge, so they cannot pair sequences up.
+  verify_graph_sequences(graph, expected, coverage)
     .wrap_err("When comparing reconstructed sequences with original FASTA records")?;
 
   Ok(())
@@ -78,9 +71,12 @@ pub fn check_unique_sequence_names(fastas: &[FastaRecord]) -> Result<(), Report>
 pub fn build(fastas: Vec<FastaRecord>, args: &PangraphBuildArgs, verify: bool) -> Result<Pangraph, Report> {
   check_unique_sequence_names(&fastas).wrap_err("When checking the names of the input sequences")?;
 
-  // If verification is requested, we need to keep a copy of the original FASTA records
-  // to compare them with the sequences reconstructed from the graph.
-  let fasta_copy = verify.then(|| fastas.clone());
+  // If verification is requested, keep the input sequences, keyed by genome name, to compare them
+  // with the sequences reconstructed from the graph.
+  let expected = verify
+    .then(|| sequences_by_name(&fastas))
+    .transpose()
+    .wrap_err("When collecting the input sequences for verification")?;
 
   // Build singleton graphs from input sequences
   // TODO: initial graphs can potentially be constructed when initializing tree clades. This could avoid a lot of boilerplate code.
@@ -134,14 +130,12 @@ pub fn build(fastas: Vec<FastaRecord>, args: &PangraphBuildArgs, verify: bool) -
               clade.data.as_ref().unwrap().paths.len()
             );
 
-            // perform checks only in debug mode and if requested
+            // perform checks only in debug mode and if requested. An intermediate graph holds
+            // only the genomes of its own clade, hence `Partial`.
             #[cfg(debug_assertions)]
-            {
-              if verify {
-                // verify the graph if requested
-                graph_sanity_checks(clade.data.as_ref().unwrap(), fasta_copy.as_ref().unwrap())
-                  .wrap_err("When performing sanity checks on the merged graph")?;
-              }
+            if let Some(expected) = expected.as_ref() {
+              graph_sanity_checks(clade.data.as_ref().unwrap(), expected, GenomeCoverage::Partial)
+                .wrap_err("When performing sanity checks on the merged graph")?;
             }
 
             Ok(())
@@ -168,10 +162,11 @@ pub fn build(fastas: Vec<FastaRecord>, args: &PangraphBuildArgs, verify: bool) -
     .take()
     .ok_or_else(|| make_internal_report!("Root clade of the guide tree contains no graph after graph alignment"))?;
 
-  // verify the final graph if requested
-  if verify {
-    graph_sanity_checks(&graph, fasta_copy.as_ref().unwrap())
+  // verify the final graph if requested. It must hold every input genome, hence `Complete`.
+  if let Some(expected) = &expected {
+    graph_sanity_checks(&graph, expected, GenomeCoverage::Complete)
       .wrap_err("When performing sanity checks on the final pangraph")?;
+    info!("Pangraph reconstructs all {} input genomes exactly", expected.len());
   }
 
   Ok(graph)

packages/pangraph/src/commands/merge/merge_run.rs

@@ -1,15 +1,16 @@
 use crate::align::alignment_args::check_alignment_backend_available;
 use crate::commands::merge::merge_args::PangraphMergeArgs;
-use crate::commands::reconstruct::reconstruct_run::reconstruct;
 use crate::io::json::{JsonPretty, json_write_file};
 use crate::make_error;
 use crate::pangraph::graph_merging::merge_graphs;
 use crate::pangraph::pangraph::Pangraph;
 use crate::pangraph::pangraph_path::PangraphPath;
+use crate::pangraph::reconstruct::{GenomeCoverage, reconstruct_by_name, verify_graph_sequences};
 use crate::representation::seq::Seq;
 use crate::utils::collections::find_duplicates;
+use color_eyre::owo_colors::{AnsiColors, OwoColorize};
+use color_eyre::{Help, SectionExt};
 use eyre::{Report, WrapErr};
-use itertools::Itertools;
 use log::{info, warn};
 use std::collections::{BTreeMap, BTreeSet};
 use std::path::Path;
@@ -30,9 +31,8 @@ pub fn merge_run(args: &PangraphMergeArgs) -> Result<(), Report> {
   // merger. Keyed by path name: neither path ids nor record order survive a merge.
   let expected = args
     .verify
-    .then(|| expected_sequences(&left, &right))
-    .transpose()
-    .wrap_err("When reconstructing the sequences of the input graphs")?;
+    .then(|| expected_sequences(args, &left, &right))
+    .transpose()?;
 
   info!(
     "=== Graph merging start:     graph sizes {} + {}",
@@ -51,7 +51,11 @@ pub fn merge_run(args: &PangraphMergeArgs) -> Result<(), Report> {
   );
 
   if let Some(expected) = expected {
-    verify_merged_sequences(&merged, &expected).wrap_err("When verifying the sequences of the merged graph")?;
+    #[cfg(debug_assertions)]
+    merged.sanity_check().wrap_err("When checking the merged graph")?;
+
+    verify_graph_sequences(&merged, &expected, GenomeCoverage::Complete)
+      .wrap_err("When verifying the sequences of the merged graph")?;
     info!("Merged graph reconstructs all {} input genomes exactly", expected.len());
   }
 
@@ -91,17 +95,6 @@ fn merge_cmd_preliminary_checks(args: &PangraphMergeArgs, left: &Pangraph, right
     );
   }
 
-  if args.verify {
-    for (graph, filepath) in [(left, &args.left_graph), (right, &args.right_graph)] {
-      if graph.path_names().any(|name| name.is_none()) {
-        return make_error!(
-          "Graph '{}' contains genomes without a name, which cannot be verified: verification matches genomes by name. Re-run without `--verify`.",
-          filepath.display()
-        );
-      }
-    }
-  }
-
   // Circularity is a per-path property, so mixing is structurally fine. It is however most often a
   // mistake, since `build --circular` applies to all genomes of a graph at once.
   if circularity(left) != circularity(right) {
@@ -119,53 +112,24 @@ fn circularity(graph: &Pangraph) -> BTreeSet<bool> {
 }
 
 /// Reconstructs the genomes of both input graphs, keyed by genome name.
-fn expected_sequences(left: &Pangraph, right: &Pangraph) -> Result<BTreeMap<String, Seq>, Report> {
+///
+/// Cross-graph name collisions are already rejected by `merge_cmd_preliminary_checks`, so the two
+/// sets cannot overwrite each other here.
+fn expected_sequences(
+  args: &PangraphMergeArgs,
+  left: &Pangraph,
+  right: &Pangraph,
+) -> Result<BTreeMap<String, Seq>, Report> {
   let mut expected = BTreeMap::new();
-  for graph in [left, right] {
-    for record in reconstruct(graph) {
-      let record = record?;
-      expected.insert(record.seq_name, record.seq);
-    }
+  for (graph, filepath) in [(left, &args.left_graph), (right, &args.right_graph)] {
+    let genomes = reconstruct_by_name(graph)
+      .wrap_err_with(|| format!("When reconstructing the genomes of graph '{}'", filepath.display()))
+      .with_section(|| {
+        "Verification matches genomes by name. Re-run without `--verify` to skip it."
+          .color(AnsiColors::Cyan)
+          .header("Suggestion:")
+      })?;
+    expected.extend(genomes);
   }
   Ok(expected)
 }
-
-/// Checks that the merged graph reconstructs exactly the genomes of the input graphs.
-/// Genomes are matched by name: path ids are renumbered by the merger, and the order in which
-/// genomes are reconstructed is therefore not the order of either input graph.
-fn verify_merged_sequences(merged: &Pangraph, expected: &BTreeMap<String, Seq>) -> Result<(), Report> {
-  #[cfg(debug_assertions)]
-  merged.sanity_check().wrap_err("When checking the merged graph")?;
-
-  let mut remaining: BTreeSet<&String> = expected.keys().collect();
-
-  for record in reconstruct(merged) {
-    let record = record?;
-    let Some(expected_seq) = expected.get(&record.seq_name) else {
-      return make_error!(
-        "Merged graph contains genome '{}', which is not present in either input graph",
-        record.seq_name
-      );
-    };
-
-    if record.seq != *expected_seq {
-      return make_error!(
-        "Sequence mismatch for genome '{}': expected length {} but got {}",
-        record.seq_name,
-        expected_seq.len(),
-        record.seq.len()
-      );
-    }
-
-    remaining.remove(&record.seq_name);
-  }
-
-  if !remaining.is_empty() {
-    return make_error!(
-      "Merged graph is missing genomes from the input graphs: [{}]",
-      remaining.into_iter().sorted().join(", ")
-    );
-  }
-
-  Ok(())
-}