diff --git a/dev/design/merge.md b/dev/design/merge.md index 7b9e7ee2..465e865f 100644 --- a/dev/design/merge.md +++ b/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) -> Result<(), Report> + pub fn verify_graph_sequences(graph: &Pangraph, expected: &BTreeMap, + 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`-based helper. +Phase 6 turned out to fix two latent bugs rather than merely unify style, both recorded in §4.4. --- diff --git a/docs/docs/reference.md b/docs/docs/reference.md index 3451f1b9..ffa6bc09 100644 --- a/docs/docs/reference.md +++ b/docs/docs/reference.md @@ -370,9 +370,9 @@ Reconstruct all input fasta sequences from 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:** @@ -380,10 +380,16 @@ Reconstruct all input fasta sequences from graph 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 ` — 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 ` — 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. diff --git a/packages/pangraph/src/commands/build/build_run.rs b/packages/pangraph/src/commands/build/build_run.rs index 731e3bee..0f906099 100644 --- a/packages/pangraph/src/commands/build/build_run.rs +++ b/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, + 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, args: &PangraphBuildArgs, verify: bool) -> Result { 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, 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, 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) diff --git a/packages/pangraph/src/commands/merge/merge_run.rs b/packages/pangraph/src/commands/merge/merge_run.rs index 7782b040..1e4bd04e 100644 --- a/packages/pangraph/src/commands/merge/merge_run.rs +++ b/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 { } /// Reconstructs the genomes of both input graphs, keyed by genome name. -fn expected_sequences(left: &Pangraph, right: &Pangraph) -> Result, 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, 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) -> 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(()) -} diff --git a/packages/pangraph/src/commands/reconstruct/reconstruct_args.rs b/packages/pangraph/src/commands/reconstruct/reconstruct_args.rs index 50dc8261..4d5741e6 100644 --- a/packages/pangraph/src/commands/reconstruct/reconstruct_args.rs +++ b/packages/pangraph/src/commands/reconstruct/reconstruct_args.rs @@ -1,17 +1,18 @@ use clap::{Parser, ValueHint}; +use smart_default::SmartDefault; use std::fmt::Debug; use std::path::PathBuf; /// Reconstruct sequences from a pangenome graph -#[derive(Parser, Debug)] +#[derive(Parser, Debug, SmartDefault)] pub struct PangraphReconstructArgs { /// 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 + /// 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 there's multiple input files, then different files can have different compression formats. + /// 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). #[clap(value_hint = ValueHint::FilePath)] #[clap(display_order = 1)] pub input_graph: Option, @@ -21,12 +22,22 @@ pub struct PangraphReconstructArgs { /// 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. + /// + /// 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 #[clap(long, short = 'o', default_value = "-")] #[clap(value_hint = ValueHint::AnyPath)] pub output_fasta: PathBuf, - /// 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. + /// 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 diff --git a/packages/pangraph/src/commands/reconstruct/reconstruct_run.rs b/packages/pangraph/src/commands/reconstruct/reconstruct_run.rs index cf4c57d3..d61e4252 100644 --- a/packages/pangraph/src/commands/reconstruct/reconstruct_run.rs +++ b/packages/pangraph/src/commands/reconstruct/reconstruct_run.rs @@ -1,15 +1,14 @@ use crate::commands::reconstruct::reconstruct_args::PangraphReconstructArgs; use crate::io::fasta::{FastaReader, FastaRecord, FastaWriter}; use crate::io::json::json_read_file; -use crate::io::seq::reverse_complement; +use crate::make_error; use crate::pangraph::pangraph::Pangraph; -use crate::pangraph::pangraph_node::NodeId; -use crate::pangraph::pangraph_path::PangraphPath; -use crate::representation::seq::Seq; -use crate::{make_error, make_internal_report}; -use eyre::Report; +use crate::pangraph::reconstruct::{path_ids_by_name, reconstruct, reconstruct_genome, verify_genome}; +use eyre::{Report, WrapErr}; use itertools::Itertools; use log::info; +use std::collections::BTreeSet; +use std::path::Path; pub fn reconstruct_run(args: &PangraphReconstructArgs) -> Result<(), Report> { let PangraphReconstructArgs { @@ -19,21 +18,18 @@ pub fn reconstruct_run(args: &PangraphReconstructArgs) -> Result<(), Report> { } = &args; let graph: Pangraph = json_read_file(input_graph)?; - let mut results = reconstruct(&graph); if let Some(verify) = verify { info!("Verifying sequences reconstructed from pangenome graph"); - let mut reader = FastaReader::from_path(verify)?; - results.try_for_each(|actual| -> Result<(), Report> { - let actual = actual?; - let mut expected = FastaRecord::new(); - reader.read(&mut expected)?; - compare_sequences(&expected, &actual)?; - Ok(()) - })?; + let n_verified = verify_against_fasta(&graph, verify) + .wrap_err_with(|| format!("When verifying reconstructed sequences against '{}'", verify.display()))?; + info!( + "Graph reconstructs all {n_verified} genomes of '{}' exactly", + verify.display() + ); } else { let mut writer = FastaWriter::from_path(output_fasta)?; - results.try_for_each(|fasta| { + reconstruct(&graph).try_for_each(|fasta| { let fasta = fasta?; writer.write(fasta.seq_name, &fasta.desc, &fasta.seq) })?; @@ -42,86 +38,52 @@ pub fn reconstruct_run(args: &PangraphReconstructArgs) -> Result<(), Report> { Ok(()) } -pub fn compare_sequences(left: &FastaRecord, right: &FastaRecord) -> Result { - if left != right { - return make_error!( - "Sequence mismatch detected: expected length {} but got {}", - left.seq.len(), - right.seq.len() - ); - } - Ok(true) -} - -pub fn reconstruct(graph: &Pangraph) -> impl Iterator> + use<'_> { - graph - .paths - .iter() - .sorted_by_key(|(path_id, _)| **path_id) - .map(|(path_id, path)| { - let index = path_id.0; - let seq = reconstruct_path_sequence(graph, path)?; - let seq_name = path - .name() - .clone() - .unwrap_or_else(|| format!("Unknown sequence #{path_id}")); - let desc = path.desc().clone(); - Ok(FastaRecord { - seq_name, - desc, - seq, - index, - }) - }) -} - -fn reconstruct_path_sequence(graph: &Pangraph, path: &PangraphPath) -> Result { - if let Some(first_node_id) = path.nodes.first() { - let first_node_pos = graph.nodes[first_node_id].position().0; +/// Checks that the graph reconstructs exactly the genomes of the given FASTA file, and returns how +/// many were verified. +/// +/// Genomes are matched by name, so the order of the FASTA records is irrelevant: a merged graph +/// renumbers its path ids and reconstructs its genomes in an order that matches no input file. +/// The file is streamed and genomes are reconstructed one at a time, so only a single genome is +/// held in memory at once. +fn verify_against_fasta(graph: &Pangraph, verify: &Path) -> Result { + let path_ids = path_ids_by_name(graph)?; + let mut remaining: BTreeSet<&str> = path_ids.keys().copied().collect(); + + let mut reader = FastaReader::from_path(verify)?; + let mut record = FastaRecord::new(); + loop { + record.clear(); + reader.read(&mut record)?; + if record.is_empty() { + break; + } - let mut genome: Seq = path - .nodes - .iter() - .map(|node_id| reconstruct_block_sequence(graph, *node_id)) - .collect::>()?; + let Some(path_id) = path_ids.get(record.seq_name.as_str()) else { + return make_error!( + "Verification file contains genome '{}', which the graph does not contain", + record.seq_name + ); + }; - let genome_len = path.tot_len(); - if genome.len() != genome_len { + if !remaining.remove(record.seq_name.as_str()) { return make_error!( - "When reconstructing sequences, genome length mismatch: computed length {} expected {}", - genome.len(), - genome_len + "Verification file contains genome '{}' more than once. Genome names must be unique, because they identify genomes.", + record.seq_name ); } - genome.rotate_right(first_node_pos); - - Ok(genome) - } else { - Ok(Seq::new()) + let actual = reconstruct_genome(graph, *path_id) + .wrap_err_with(|| format!("When reconstructing genome '{}'", record.seq_name))?; + verify_genome(&record.seq_name, &record.seq, &actual)?; } -} - -fn reconstruct_block_sequence(graph: &Pangraph, node_id: NodeId) -> Result { - let node = graph - .nodes - .get(&node_id) - .ok_or_else(|| make_internal_report!("Node {node_id} not found in graph"))?; - - let block_id = node.block_id(); - let block = graph - .blocks - .get(&block_id) - .ok_or_else(|| make_internal_report!("Block {block_id} not found in graph"))?; - // Get edits and apply them to the consensus sequence - let edits = block.alignment(node_id); - - let mut s = edits.apply(block.consensus())?; - - // Reverse-complement if on opposite strand - if node.strand().is_reverse() { - s = reverse_complement(&s)?; + if !remaining.is_empty() { + return make_error!( + "Verification file is missing {} genome(s) that the graph contains: [{}]", + remaining.len(), + remaining.iter().join(", ") + ); } - Ok(s) + + Ok(path_ids.len()) } diff --git a/packages/pangraph/src/pangraph/mod.rs b/packages/pangraph/src/pangraph/mod.rs index 201c062e..aebbc33f 100644 --- a/packages/pangraph/src/pangraph/mod.rs +++ b/packages/pangraph/src/pangraph/mod.rs @@ -6,6 +6,7 @@ pub mod pangraph_block; pub mod pangraph_interval; pub mod pangraph_node; pub mod pangraph_path; +pub mod reconstruct; pub mod reweave; pub mod slice; pub mod split_matches; diff --git a/packages/pangraph/src/pangraph/pangraph.rs b/packages/pangraph/src/pangraph/pangraph.rs index 77e26c61..6336a629 100644 --- a/packages/pangraph/src/pangraph/pangraph.rs +++ b/packages/pangraph/src/pangraph/pangraph.rs @@ -409,11 +409,11 @@ mod tests { #![allow(non_snake_case, clippy::redundant_clone)] use super::*; - use crate::commands::reconstruct::reconstruct_run::reconstruct; use crate::o; use crate::pangraph::edits::Edit; use crate::pangraph::pangraph_node::PangraphNode; use crate::pangraph::pangraph_path::PangraphPath; + use crate::pangraph::reconstruct::reconstruct; use crate::pangraph::strand::Strand::{Forward, Reverse}; use itertools::Itertools; use maplit::btreemap; diff --git a/packages/pangraph/src/pangraph/reconstruct.rs b/packages/pangraph/src/pangraph/reconstruct.rs new file mode 100644 index 00000000..02bd8059 --- /dev/null +++ b/packages/pangraph/src/pangraph/reconstruct.rs @@ -0,0 +1,444 @@ +use crate::io::fasta::FastaRecord; +use crate::io::seq::reverse_complement; +use crate::pangraph::pangraph::Pangraph; +use crate::pangraph::pangraph_node::NodeId; +use crate::pangraph::pangraph_path::{PangraphPath, PathId}; +use crate::representation::seq::Seq; +use crate::utils::collections::find_duplicates; +use crate::utils::string::str_slice_safe; +use crate::{make_error, make_internal_report}; +use eyre::{Report, WrapErr}; +use itertools::Itertools; +use std::collections::BTreeMap; + +/// Number of genome names listed in full in an error message before the rest are elided. +const MAX_NAMES_IN_ERROR: usize = 10; + +/// Number of nucleotides shown on each side of the first difference between two genomes. +const MISMATCH_CONTEXT: usize = 10; + +/// How much of the expected genome set the graph is required to reconstruct. +/// +/// In both cases a genome that the graph contains but `expected` does not is an error: the +/// variants differ only in whether the graph is allowed to be missing expected genomes. +#[derive(Copy, Clone, Debug, Eq, PartialEq)] +pub enum GenomeCoverage { + /// The graph must reconstruct every expected genome, and nothing else. + Complete, + /// The graph may reconstruct only some of the expected genomes. Used for the intermediate graphs + /// of a build, which hold the genomes of one clade of the guide tree. + Partial, +} + +/// Reconstructs every genome of the graph as a FASTA record, ordered by path id. +/// +/// Record order and `index` reproduce the order of the original input FASTA only for graphs +/// produced directly by `build`. A merged graph renumbers its path ids, so consumers must match +/// records by genome name. +pub fn reconstruct(graph: &Pangraph) -> impl Iterator> + use<'_> { + graph + .paths + .iter() + .sorted_by_key(|(path_id, _)| **path_id) + .map(|(path_id, path)| { + let index = path_id.0; + let seq = reconstruct_path_sequence(graph, path)?; + let seq_name = path + .name() + .clone() + .unwrap_or_else(|| format!("Unknown sequence #{path_id}")); + let desc = path.desc().clone(); + Ok(FastaRecord { + seq_name, + desc, + seq, + index, + }) + }) +} + +/// Maps each genome name of the graph to the id of the path holding it. +/// +/// Errors if any path is unnamed or if two paths share a name. The genome name is the only +/// identifier that survives a merge unchanged, and is therefore the key that verification matches +/// on. Note that this inspects the paths directly rather than going through [`reconstruct`], which +/// masks unnamed paths behind a placeholder name. +pub fn path_ids_by_name(graph: &Pangraph) -> Result, Report> { + let unnamed = graph + .paths + .iter() + .filter(|(_, path)| path.name.is_none()) + .map(|(path_id, _)| path_id.to_string()) + .collect_vec(); + + if !unnamed.is_empty() { + return make_error!( + "Graph contains {} genome(s) without a name (path ids: {}). Genomes are identified by name, so every path must be named.", + unnamed.len(), + format_names(&unnamed) + ); + } + + let duplicates = find_duplicates(graph.path_names().flatten()); + if !duplicates.is_empty() { + return make_error!( + "Graph contains duplicate genome names: {}. Genome names must be unique, because they identify genomes.", + format_names(&duplicates) + ); + } + + Ok( + graph + .paths + .iter() + .filter_map(|(path_id, path)| path.name.as_deref().map(|name| (name, *path_id))) + .collect(), + ) +} + +/// Reconstructs the genome of a single path. +pub fn reconstruct_genome(graph: &Pangraph, path_id: PathId) -> Result { + let path = graph + .paths + .get(&path_id) + .ok_or_else(|| make_internal_report!("Path {path_id} not found in graph"))?; + reconstruct_path_sequence(graph, path) +} + +/// Reconstructs every genome of the graph, keyed by genome name. +/// +/// Errors if any path is unnamed or if two paths share a name. +pub fn reconstruct_by_name(graph: &Pangraph) -> Result, Report> { + path_ids_by_name(graph)? + .into_iter() + .map(|(name, path_id)| { + let seq = + reconstruct_genome(graph, path_id).wrap_err_with(|| format!("When reconstructing the genome of '{name}'"))?; + Ok((name.to_owned(), seq)) + }) + .collect() +} + +/// Collects FASTA records into genome sequences keyed by name. +/// +/// Errors on duplicate names: the name is the key that verification matches on. +pub fn sequences_by_name(fastas: &[FastaRecord]) -> Result, Report> { + let duplicates = find_duplicates(fastas.iter().map(|fasta| fasta.seq_name.as_str())); + if !duplicates.is_empty() { + return make_error!( + "Duplicate sequence names found: {}. Sequence names must be unique, because they identify genomes.", + format_names(&duplicates) + ); + } + + Ok( + fastas + .iter() + .map(|fasta| (fasta.seq_name.clone(), fasta.seq.clone())) + .collect(), + ) +} + +/// Compares one reconstructed genome against the sequence it is expected to have, reporting the +/// genome name and the first position at which the two differ. +pub fn verify_genome(name: &str, expected: &Seq, actual: &Seq) -> Result<(), Report> { + if expected == actual { + return Ok(()); + } + + if expected.len() != actual.len() { + return make_error!( + "Sequence mismatch for genome '{name}': expected length {} but got {}", + expected.len(), + actual.len() + ); + } + + let pos = expected + .iter() + .zip(actual.iter()) + .position(|(e, a)| e != a) + .ok_or_else(|| make_internal_report!("Genomes of '{name}' compare as different but share every character"))?; + + let (start, end) = (pos.saturating_sub(MISMATCH_CONTEXT), pos + MISMATCH_CONTEXT + 1); + make_error!( + "Sequence mismatch for genome '{name}' at position {pos} (length {}):\n expected: {}\n actual: {}", + expected.len(), + str_slice_safe(expected.as_str(), start, end), + str_slice_safe(actual.as_str(), start, end) + ) +} + +/// Checks that the graph reconstructs the expected genomes, matched by name. +/// +/// Genomes are never matched by position or by path id: a merge renumbers path ids, so neither the +/// order in which genomes are reconstructed nor [`FastaRecord::index`] can pair them up. Only the +/// sequence content is compared; descriptions are ignored. +/// +/// Genomes are reconstructed and dropped one at a time, so this holds no more than a single genome +/// in memory beyond `expected`. +pub fn verify_graph_sequences( + graph: &Pangraph, + expected: &BTreeMap, + coverage: GenomeCoverage, +) -> Result<(), Report> { + let path_ids = path_ids_by_name(graph)?; + + for (name, path_id) in &path_ids { + let Some(expected_seq) = expected.get(*name) else { + return make_error!("Graph contains genome '{name}', which is not among the expected genomes"); + }; + let actual = + reconstruct_genome(graph, *path_id).wrap_err_with(|| format!("When reconstructing genome '{name}'"))?; + verify_genome(name, expected_seq, &actual)?; + } + + if coverage == GenomeCoverage::Complete { + let missing = expected + .keys() + .filter(|name| !path_ids.contains_key(name.as_str())) + .cloned() + .collect_vec(); + + if !missing.is_empty() { + return make_error!( + "Graph is missing {} expected genome(s): {}", + missing.len(), + format_names(&missing) + ); + } + } + + Ok(()) +} + +/// Formats a list of genome names for an error message, eliding all but the first few. +fn format_names>(names: &[S]) -> String { + let shown = names.iter().take(MAX_NAMES_IN_ERROR).map(AsRef::as_ref).join(", "); + if names.len() > MAX_NAMES_IN_ERROR { + format!("[{shown}, ... and {} more]", names.len() - MAX_NAMES_IN_ERROR) + } else { + format!("[{shown}]") + } +} + +fn reconstruct_path_sequence(graph: &Pangraph, path: &PangraphPath) -> Result { + if let Some(first_node_id) = path.nodes.first() { + let first_node_pos = graph.nodes[first_node_id].position().0; + + let mut genome: Seq = path + .nodes + .iter() + .map(|node_id| reconstruct_block_sequence(graph, *node_id)) + .collect::>()?; + + let genome_len = path.tot_len(); + if genome.len() != genome_len { + return crate::make_error!( + "When reconstructing sequences, genome length mismatch: computed length {} expected {}", + genome.len(), + genome_len + ); + } + + genome.rotate_right(first_node_pos); + + Ok(genome) + } else { + Ok(Seq::new()) + } +} + +fn reconstruct_block_sequence(graph: &Pangraph, node_id: NodeId) -> Result { + let node = graph + .nodes + .get(&node_id) + .ok_or_else(|| make_internal_report!("Node {node_id} not found in graph"))?; + + let block_id = node.block_id(); + let block = graph + .blocks + .get(&block_id) + .ok_or_else(|| make_internal_report!("Block {block_id} not found in graph"))?; + + // Get edits and apply them to the consensus sequence + let edits = block.alignment(node_id); + + let mut s = edits.apply(block.consensus())?; + + // Reverse-complement if on opposite strand + if node.strand().is_reverse() { + s = reverse_complement(&s)?; + } + Ok(s) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::o; + use crate::pangraph::edits::Edit; + use crate::pangraph::pangraph_block::{BlockId, PangraphBlock}; + use crate::pangraph::pangraph_node::PangraphNode; + use crate::pangraph::strand::Strand::{Forward, Reverse}; + use crate::utils::error::report_to_string; + use maplit::btreemap; + use pretty_assertions::assert_eq; + use rstest::rstest; + + /// A two-genome graph: `a` is a forward single-block path, `b` a reverse one, so their + /// reconstructed sequences are hand-checkable. + fn two_genome_graph(names: [Option<&str>; 2]) -> Pangraph { + let blocks = btreemap! { + BlockId(0) => PangraphBlock::new(BlockId(0), "ACGTACGT", btreemap!{ NodeId(0) => Edit::empty() }), + BlockId(1) => PangraphBlock::new(BlockId(1), "TTTTGGGG", btreemap!{ NodeId(1) => Edit::empty() }), + }; + let nodes = btreemap! { + NodeId(0) => PangraphNode::new(Some(NodeId(0)), BlockId(0), PathId(0), Forward, (0, 8)), + NodeId(1) => PangraphNode::new(Some(NodeId(1)), BlockId(1), PathId(1), Reverse, (0, 8)), + }; + let paths = btreemap! { + PathId(0) => PangraphPath::new(Some(PathId(0)), [NodeId(0)], 8, false, names[0].map(String::from), None), + PathId(1) => PangraphPath::new(Some(PathId(1)), [NodeId(1)], 8, false, names[1].map(String::from), None), + }; + Pangraph { paths, blocks, nodes } + } + + fn graph() -> Pangraph { + two_genome_graph([Some("a"), Some("b")]) + } + + fn expected_genomes() -> BTreeMap { + btreemap! { o!("a") => Seq::from_str("ACGTACGT"), o!("b") => Seq::from_str("CCCCAAAA") } + } + + #[rstest] + fn test_reconstruct_by_name() { + assert_eq!(reconstruct_by_name(&graph()).unwrap(), expected_genomes()); + } + + #[rstest] + fn test_path_ids_by_name_rejects_unnamed_path() { + let graph = two_genome_graph([Some("a"), None]); + assert!(report_to_string(&path_ids_by_name(&graph).unwrap_err()).contains("without a name")); + assert!(report_to_string(&reconstruct_by_name(&graph).unwrap_err()).contains("without a name")); + } + + #[rstest] + fn test_path_ids_by_name_rejects_duplicate_names() { + let graph = two_genome_graph([Some("a"), Some("a")]); + assert!(report_to_string(&path_ids_by_name(&graph).unwrap_err()).contains("duplicate genome names")); + } + + #[rstest] + fn test_sequences_by_name_rejects_duplicate_names() { + let record = |name: &str| FastaRecord { + seq_name: name.to_owned(), + desc: None, + seq: Seq::from_str("ACGT"), + index: 0, + }; + let fastas = [record("a"), record("a")]; + assert!(report_to_string(&sequences_by_name(&fastas).unwrap_err()).contains("Duplicate sequence names")); + } + + #[rstest] + fn test_verify_graph_sequences_accepts_exact_match() { + verify_graph_sequences(&graph(), &expected_genomes(), GenomeCoverage::Complete).unwrap(); + } + + /// The whole point of keying by name: `index` and `desc` must not take part in the comparison. + /// The previous whole-`FastaRecord` comparison failed this, which is why verifying a merged + /// graph reported a mismatch between two sequences of identical length. + #[rstest] + fn test_verify_graph_sequences_ignores_index_and_desc() { + let fastas = [ + FastaRecord { + seq_name: o!("a"), + desc: Some(o!("some description")), + seq: Seq::from_str("ACGTACGT"), + index: 41, + }, + FastaRecord { + seq_name: o!("b"), + desc: None, + seq: Seq::from_str("CCCCAAAA"), + index: 42, + }, + ]; + let expected = sequences_by_name(&fastas).unwrap(); + verify_graph_sequences(&graph(), &expected, GenomeCoverage::Complete).unwrap(); + } + + #[rstest] + fn test_verify_graph_sequences_detects_missing_genome() { + let mut expected = expected_genomes(); + expected.insert(o!("c"), Seq::from_str("GGGG")); + + let err = verify_graph_sequences(&graph(), &expected, GenomeCoverage::Complete).unwrap_err(); + let err = report_to_string(&err); + assert!(err.contains("missing"), "unexpected error: {err}"); + assert!(err.contains('c'), "unexpected error: {err}"); + } + + #[rstest] + fn test_verify_graph_sequences_detects_extra_genome() { + let mut expected = expected_genomes(); + expected.remove("b"); + + let err = report_to_string(&verify_graph_sequences(&graph(), &expected, GenomeCoverage::Complete).unwrap_err()); + assert!( + err.contains("not among the expected genomes"), + "unexpected error: {err}" + ); + } + + /// `Partial` relaxes only one direction: the graph may hold a subset of the expected genomes, + /// but a genome the expected set does not know about is still an error. + #[rstest] + fn test_verify_graph_sequences_partial_accepts_superset() { + let mut expected = expected_genomes(); + expected.insert(o!("c"), Seq::from_str("GGGG")); + + verify_graph_sequences(&graph(), &expected, GenomeCoverage::Partial).unwrap(); + } + + #[rstest] + fn test_verify_graph_sequences_partial_still_rejects_extra_genome() { + let mut expected = expected_genomes(); + expected.remove("b"); + + let err = report_to_string(&verify_graph_sequences(&graph(), &expected, GenomeCoverage::Partial).unwrap_err()); + assert!( + err.contains("not among the expected genomes"), + "unexpected error: {err}" + ); + } + + #[rstest] + fn test_verify_genome_reports_first_difference() { + let expected = Seq::from_str("ACGTACGT"); + let actual = Seq::from_str("ACGTTCGT"); + + let err = report_to_string(&verify_genome("a", &expected, &actual).unwrap_err()); + assert!(err.contains("at position 4"), "unexpected error: {err}"); + assert!(err.contains("ACGTACGT"), "unexpected error: {err}"); + assert!(err.contains("ACGTTCGT"), "unexpected error: {err}"); + } + + #[rstest] + fn test_verify_genome_reports_length_mismatch() { + let err = report_to_string(&verify_genome("a", &Seq::from_str("ACGT"), &Seq::from_str("ACG")).unwrap_err()); + assert!(err.contains("expected length 4 but got 3"), "unexpected error: {err}"); + } + + #[rstest] + fn test_format_names_elides_long_lists() { + let names = (0..12).map(|i| format!("g{i}")).collect_vec(); + assert_eq!( + format_names(&names), + "[g0, g1, g2, g3, g4, g5, g6, g7, g8, g9, ... and 2 more]" + ); + assert_eq!(format_names(&names[..2]), "[g0, g1]"); + } +} diff --git a/packages/pangraph/tests/itest_merge.rs b/packages/pangraph/tests/itest_merge.rs index c9b50344..82de39a7 100644 --- a/packages/pangraph/tests/itest_merge.rs +++ b/packages/pangraph/tests/itest_merge.rs @@ -9,10 +9,10 @@ mod tests { use pangraph::commands::build::build_run::build; use pangraph::commands::merge::merge_args::PangraphMergeArgs; use pangraph::commands::merge::merge_run::merge_run; - use pangraph::commands::reconstruct::reconstruct_run::reconstruct; use pangraph::io::fasta::{FastaReader, FastaRecord}; use pangraph::io::json::{JsonPretty, json_write_file}; use pangraph::pangraph::pangraph::Pangraph; + use pangraph::pangraph::reconstruct::reconstruct; use pangraph::representation::seq::Seq; use pangraph::utils::error::report_to_string; use pretty_assertions::assert_eq; @@ -36,7 +36,7 @@ mod tests { circular: false, ..PangraphBuildArgs::default() }; - let graph = build(fastas, &args, false)?; + let graph = build(fastas, &args, true)?; let path = dir.path().join(name); json_write_file(&path, &graph, JsonPretty(false))?; Ok(path) @@ -79,6 +79,7 @@ mod tests { merge_run(&merge_args(left.clone(), right.clone(), output.clone()))?; let merged = read_graph(&output)?; + #[cfg(debug_assertions)] merged.sanity_check()?; // all genomes of both inputs are present, exactly once @@ -130,6 +131,7 @@ mod tests { merge_run(&merge_args(left, right, output.clone()))?; let merged = read_graph(&output)?; + #[cfg(debug_assertions)] merged.sanity_check()?; assert_eq!(merged.paths.len(), 4); @@ -155,6 +157,23 @@ mod tests { Ok(()) } + /// `build --verify` used to pair reconstructed genomes with input records by `FastaRecord::index`, + /// which is only valid when the records' indices happen to be exactly `0..n-1`. Here they are + /// `[3, 4, 5]` for a 3-record slice, which panicked with an out-of-bounds index before genomes + /// were matched by name. + /// + /// Note that the intermediate-clade half of this check only runs in debug builds. + #[rstest] + fn itest_build_verify_with_nonzero_record_indices() -> Result<(), Report> { + let (_, right_fastas) = read_and_split("../../data/ges-1.fa", 3, 6)?; + assert_eq!(right_fastas.iter().map(|f| f.index).collect_vec(), vec![3, 4, 5]); + + let graph = build(right_fastas, &PangraphBuildArgs::default(), true)?; + assert_eq!(graph.paths.len(), 3); + + Ok(()) + } + /// `build` enforces the same uniqueness invariant on its input FASTA records, so that a graph /// can never carry duplicate genome names into a later merge. #[rstest] diff --git a/packages/pangraph/tests/itest_reconstruct.rs b/packages/pangraph/tests/itest_reconstruct.rs new file mode 100644 index 00000000..4ea3b8ff --- /dev/null +++ b/packages/pangraph/tests/itest_reconstruct.rs @@ -0,0 +1,238 @@ +mod common; + +#[cfg(test)] +mod tests { + use eyre::Report; + use itertools::Itertools; + use maplit::btreemap; + use pangraph::align::alignment_args::GraphMergeParams; + use pangraph::commands::build::build_args::PangraphBuildArgs; + use pangraph::commands::build::build_run::build; + use pangraph::commands::merge::merge_args::PangraphMergeArgs; + use pangraph::commands::merge::merge_run::merge_run; + use pangraph::commands::reconstruct::reconstruct_args::PangraphReconstructArgs; + use pangraph::commands::reconstruct::reconstruct_run::reconstruct_run; + use pangraph::io::fasta::{FastaReader, FastaRecord, FastaWriter}; + use pangraph::io::json::{JsonPretty, json_write_file}; + use pangraph::pangraph::edits::Edit; + use pangraph::pangraph::pangraph::Pangraph; + use pangraph::pangraph::pangraph_block::{BlockId, PangraphBlock}; + use pangraph::pangraph::pangraph_node::{NodeId, PangraphNode}; + use pangraph::pangraph::pangraph_path::{PangraphPath, PathId}; + use pangraph::pangraph::strand::Strand::Forward; + use pangraph::representation::seq::Seq; + use pangraph::utils::error::report_to_string; + use pretty_assertions::assert_eq; + use rstest::rstest; + use std::collections::BTreeMap; + use std::path::{Path, PathBuf}; + use tempfile::{TempDir, tempdir}; + + const GES: &str = "../../data/ges-1.fa"; + + /// Builds a merged graph from two disjoint halves of a FASTA file, and returns its path together + /// with all the input records. The merger renumbers path ids, which is what makes this graph the + /// interesting case: neither record order nor `FastaRecord::index` can pair genomes up any more. + fn merged_graph(dir: &TempDir) -> Result<(PathBuf, Vec), Report> { + let mut fastas = FastaReader::from_paths(&[PathBuf::from(GES)])?.read_many()?; + fastas.truncate(6); + let right_fastas = fastas.split_off(3); + let all = fastas.iter().chain(right_fastas.iter()).cloned().collect_vec(); + + let args = PangraphBuildArgs::default(); + let left = write_graph(dir, "left.json", &build(fastas, &args, true)?)?; + let right = write_graph(dir, "right.json", &build(right_fastas, &args, true)?)?; + + let output = dir.path().join("merged.json"); + merge_run(&PangraphMergeArgs { + left_graph: left, + right_graph: right, + output_json: output.clone(), + merge_params: GraphMergeParams::default(), + verify: true, + })?; + + Ok((output, all)) + } + + fn write_graph(dir: &TempDir, name: &str, graph: &Pangraph) -> Result { + let path = dir.path().join(name); + json_write_file(&path, graph, JsonPretty(false))?; + Ok(path) + } + + fn write_fasta(dir: &TempDir, name: &str, records: &[FastaRecord]) -> Result { + let path = dir.path().join(name); + let mut writer = FastaWriter::from_path(&path)?; + for record in records { + writer.write(record.seq_name.clone(), &record.desc, &record.seq)?; + } + drop(writer); + Ok(path) + } + + fn verify_args(graph: &Path, verify: &Path) -> PangraphReconstructArgs { + PangraphReconstructArgs { + input_graph: Some(graph.to_owned()), + output_fasta: PathBuf::from("-"), + verify: Some(verify.to_owned()), + } + } + + /// The headline case: a merged graph verified against its genomes in a different order. Before + /// verification was keyed by name this failed spuriously, reporting a mismatch between two + /// sequences of identical length, because whole `FastaRecord`s were compared and `index` differs + /// once a merge has renumbered the path ids. + #[rstest] + fn itest_reconstruct_verify_merged_graph_ignores_order() -> Result<(), Report> { + let dir = tempdir()?; + let (graph, mut records) = merged_graph(&dir)?; + records.reverse(); + + let verify = write_fasta(&dir, "verify.fa", &records)?; + reconstruct_run(&verify_args(&graph, &verify))?; + + Ok(()) + } + + /// A record the graph does not contain used to be read past and silently ignored, so the command + /// exited 0 while the verification file and the graph disagreed. + #[rstest] + fn itest_reconstruct_verify_rejects_surplus_record() -> Result<(), Report> { + let dir = tempdir()?; + let (graph, mut records) = merged_graph(&dir)?; + records.push(FastaRecord { + seq_name: "not_in_the_graph".to_owned(), + desc: None, + seq: Seq::from_str("ACGTACGT"), + index: 0, + }); + + let verify = write_fasta(&dir, "verify.fa", &records)?; + let err = report_to_string(&reconstruct_run(&verify_args(&graph, &verify)).unwrap_err()); + + assert!(err.contains("not_in_the_graph"), "unexpected error: {err}"); + assert!( + err.contains("which the graph does not contain"), + "unexpected error: {err}" + ); + + Ok(()) + } + + #[rstest] + fn itest_reconstruct_verify_rejects_missing_record() -> Result<(), Report> { + let dir = tempdir()?; + let (graph, mut records) = merged_graph(&dir)?; + let dropped = records.remove(0).seq_name; + + let verify = write_fasta(&dir, "verify.fa", &records)?; + let err = report_to_string(&reconstruct_run(&verify_args(&graph, &verify)).unwrap_err()); + + assert!(err.contains(&dropped), "unexpected error: {err}"); + assert!(err.contains("missing"), "unexpected error: {err}"); + + Ok(()) + } + + #[rstest] + fn itest_reconstruct_verify_rejects_duplicate_record() -> Result<(), Report> { + let dir = tempdir()?; + let (graph, mut records) = merged_graph(&dir)?; + records.push(records[0].clone()); + + let verify = write_fasta(&dir, "verify.fa", &records)?; + let err = report_to_string(&reconstruct_run(&verify_args(&graph, &verify)).unwrap_err()); + + assert!(err.contains("more than once"), "unexpected error: {err}"); + + Ok(()) + } + + /// A real sequence difference must be reported with the genome name and the position, rather than + /// the old "expected length N but got N". + #[rstest] + fn itest_reconstruct_verify_reports_mutated_base() -> Result<(), Report> { + let dir = tempdir()?; + let (graph, mut records) = merged_graph(&dir)?; + + let mutated = records[1].seq_name.clone(); + let seq = records[1].seq.as_str().to_owned(); + let original = seq.as_bytes()[100] as char; + let replacement = if original == 'A' { 'T' } else { 'A' }; + let mut bases = seq.into_bytes(); + bases[100] = replacement as u8; + records[1].seq = Seq::from_str(std::str::from_utf8(&bases)?); + + let verify = write_fasta(&dir, "verify.fa", &records)?; + let err = report_to_string(&reconstruct_run(&verify_args(&graph, &verify)).unwrap_err()); + + assert!(err.contains(&mutated), "unexpected error: {err}"); + assert!(err.contains("at position 100"), "unexpected error: {err}"); + + Ok(()) + } + + /// Verification matches by name, so a graph with unnamed paths cannot be verified at all. + #[rstest] + fn itest_reconstruct_verify_rejects_unnamed_path() -> Result<(), Report> { + let dir = tempdir()?; + let graph = Pangraph { + blocks: btreemap! { + BlockId(0) => PangraphBlock::new(BlockId(0), "ACGTACGT", btreemap!{ NodeId(0) => Edit::empty() }), + }, + nodes: btreemap! { + NodeId(0) => PangraphNode::new(Some(NodeId(0)), BlockId(0), PathId(0), Forward, (0, 8)), + }, + paths: btreemap! { + PathId(0) => PangraphPath::new(Some(PathId(0)), [NodeId(0)], 8, false, None, None), + }, + }; + + let graph = write_graph(&dir, "unnamed.json", &graph)?; + let verify = write_fasta( + &dir, + "verify.fa", + &[FastaRecord { + seq_name: "a".to_owned(), + desc: None, + seq: Seq::from_str("ACGTACGT"), + index: 0, + }], + )?; + + let err = report_to_string(&reconstruct_run(&verify_args(&graph, &verify)).unwrap_err()); + assert!(err.contains("without a name"), "unexpected error: {err}"); + + Ok(()) + } + + /// The non-verify path is untouched: every genome is written out, and matching them by name + /// recovers the inputs exactly. Their order is deliberately not asserted. + #[rstest] + fn itest_reconstruct_writes_every_genome() -> Result<(), Report> { + let dir = tempdir()?; + let (graph, records) = merged_graph(&dir)?; + let output = dir.path().join("out.fa"); + + reconstruct_run(&PangraphReconstructArgs { + input_graph: Some(graph), + output_fasta: output.clone(), + verify: None, + })?; + + let written: BTreeMap = FastaReader::from_path(&output)? + .read_many()? + .into_iter() + .map(|record| (record.seq_name, record.seq)) + .collect(); + let expected: BTreeMap = records + .into_iter() + .map(|record| (record.seq_name, record.seq)) + .collect(); + + assert_eq!(written, expected); + + Ok(()) + } +}