>(filepath: &Option) -> Result {
let reader = open_file_or_stdin(filepath)?;
let data = read_reader_to_string(reader).wrap_err("When reading Pangraph JSON")?;
- Self::from_str(&data).wrap_err("When parsing Pangraph JSON")
+ let graph = Self::from_str(&data).wrap_err("When parsing Pangraph JSON")?;
+ graph.validate().wrap_err_with(|| match filepath {
+ Some(filepath) => format!("When validating the graph read from '{}'", filepath.as_ref().display()),
+ None => "When validating the graph read from standard input".to_owned(),
+ })?;
+ Ok(graph)
}
pub fn to_string_pretty(&self) -> Result {
@@ -65,6 +84,81 @@ impl Pangraph {
self.blocks.values().map(|block| block.consensus())
}
+ /// Returns this graph with its path ids renumbered contiguously from `offset`, preserving their
+ /// relative order.
+ ///
+ /// Blocks are untouched and node ids are preserved: only the `path_id` a node stores, and the id
+ /// of the path itself, are rewritten. Node ids are derived from [`PangraphPath::seed`] rather
+ /// than from the path id (see [`PangraphNode::with_derived_id`]), so renumbering cannot put a
+ /// node id out of step with the contents it was derived from.
+ ///
+ /// Used by `merge` to lift one graph's genomes above the other's: path ids are sequential within
+ /// each graph, so two graphs built independently always collide on them. The graph is taken by
+ /// value so that the sequences can be moved rather than copied.
+ ///
+ /// Errors if a node refers to a path the graph does not contain, which can only happen in a graph
+ /// that pangraph did not write.
+ pub fn renumber_paths(self, offset: usize) -> Result {
+ let Self {
+ paths,
+ blocks,
+ mut nodes,
+ } = self;
+
+ let path_map: BTreeMap = paths
+ .keys()
+ .enumerate()
+ .map(|(rank, &pid)| (pid, PathId(offset + rank)))
+ .collect();
+
+ // Node ids do not change, so the map is updated in place rather than rebuilt: only the path a
+ // node points at moves.
+ //
+ // `path_map` is keyed by the map keys of this graph, so looking a path's own new id up by its
+ // key cannot miss. The path id read off a node's *contents* is a different matter: it can
+ // dangle in a graph that pangraph did not write, so that lookup is fallible.
+ for (nid, node) in &mut nodes {
+ let path_id = *path_map.get(&node.path_id()).ok_or_else(|| {
+ make_report!(
+ "Node {nid} refers to path {}, which the graph does not contain",
+ node.path_id()
+ )
+ })?;
+ node.set_path_id(path_id);
+ }
+
+ let paths: BTreeMap = paths
+ .into_iter()
+ .map(|(pid, mut path)| {
+ path.id = path_map[&pid];
+ (path.id, path)
+ })
+ .collect();
+
+ Ok(Self { paths, blocks, nodes })
+ }
+
+ /// Returns true if this graph shares no block, node or path id with `other`.
+ ///
+ /// Block and node ids are derived from genome names, which `merge` requires to be distinct across
+ /// the two graphs, so this holds by construction for graphs written by pangraph 1.4 or later.
+ /// Graphs written by earlier versions derive their ids from the input order instead, and two of
+ /// those do collide, which is what `merge` uses this to detect.
+ pub fn is_id_disjoint_from(&self, other: &Self) -> bool {
+ self.blocks.keys().all(|bid| !other.blocks.contains_key(bid))
+ && self.nodes.keys().all(|nid| !other.nodes.contains_key(nid))
+ && self.paths.keys().all(|pid| !other.paths.contains_key(pid))
+ }
+
+ /// Returns the smallest path id that is above every path id of this graph.
+ ///
+ /// The offset `merge` renumbers the appended graph from, so that the two sets of path ids cannot
+ /// overlap. Not simply the number of paths: `simplify` drops paths without renumbering, so path
+ /// ids are not necessarily contiguous.
+ pub fn path_id_upper_bound(&self) -> usize {
+ self.paths.keys().map(|pid| pid.0 + 1).max().unwrap_or(0)
+ }
+
pub fn update(&mut self, u: &GraphUpdate) {
// Consistency check: node ids
let old_nodes_set_from_graph: BTreeSet = self.blocks[&u.b_old_id].alignment_keys();
@@ -131,47 +225,123 @@ impl Pangraph {
}
}
- #[cfg(any(test, debug_assertions))]
- pub fn sanity_check(&self) -> Result<(), Report> {
+ /// Checks the invariants that let the rest of pangraph resolve graph entities by id without
+ /// error handling, and that keep sequence reconstruction inside array bounds.
+ ///
+ /// Unlike [`Self::sanity_check`] this runs in release builds, because it guards against a
+ /// malformed *input*: a graph read from a file was not necessarily written by pangraph, so
+ /// nothing guarantees that its cross-references resolve or that its offsets are in range. Every
+ /// graph read through [`Self::from_path`] is validated, and that is what lets the code downstream
+ /// index the maps directly and treat a lookup that fails anyway as an internal error.
+ ///
+ /// Deliberately limited to what a bad file can break. It does *not* check that the graph is
+ /// semantically coherent (that node positions tile the genome, that edits do not overlap), since
+ /// those are symptoms of a bug in pangraph rather than of a bad input, and are covered by
+ /// [`Self::sanity_check`] in debug builds. Errors are reported as ordinary user-facing errors for
+ /// the same reason: the offending graph came from the user.
+ ///
+ /// Costs roughly a tenth of the JSON parse it follows, so it is always worth running.
+ pub fn validate(&self) -> Result<(), Report> {
+ // Each entity is stored under a map key *and* carries its own id. Serde keys the maps by the
+ // JSON object key, so the two can disagree in a graph that pangraph did not write. Everything
+ // that resolves an entity by key while reading its id off the entity itself depends on them
+ // agreeing, so check it first. One integer comparison per entity.
+ for (block_id, block) in &self.blocks {
+ if block.id() != *block_id {
+ return make_error!("Block is stored under id {block_id} but reports id {}", block.id());
+ }
+ }
+
for (node_id, node) in &self.nodes {
- if !self.blocks.contains_key(&node.block_id()) {
- return Err(eyre::eyre!("Block {} not found in graph", node.block_id()));
+ if node.id() != *node_id {
+ return make_error!("Node is stored under id {node_id} but reports id {}", node.id());
+ }
+ }
+
+ for (path_id, path) in &self.paths {
+ if path.id() != *path_id {
+ return make_error!("Path is stored under id {path_id} but reports id {}", path.id());
}
- let block = &self.blocks[&node.block_id()];
+ }
- if !self.paths.contains_key(&node.path_id()) {
- return Err(eyre::eyre!("Path {} not found in graph", node.path_id()));
+ // Which path walks each node, as claimed by the paths. Collected in one pass because
+ // `path.nodes` is a `Vec`: scanning it per node would re-read a genome's entire walk once for
+ // every node of that genome.
+ let mut walked_by: BTreeMap = BTreeMap::new();
+ for (path_id, path) in &self.paths {
+ for node_id in &path.nodes {
+ if !self.nodes.contains_key(node_id) {
+ return make_error!("Node {node_id} from path {path_id} not found in graph");
+ }
+ // A node id may repeat within the walk of one path, which happens for empty nodes, but two
+ // paths sharing a node would make the node's own `path_id` ambiguous.
+ if let Some(previous) = walked_by.insert(*node_id, *path_id) {
+ if previous != *path_id {
+ return make_error!("Node {node_id} is walked by both path {previous} and path {path_id}");
+ }
+ }
}
- let path = &self.paths[&node.path_id()];
+ }
+ for (node_id, node) in &self.nodes {
+ let Some(block) = self.blocks.get(&node.block_id()) else {
+ return make_error!("Block {} of node {node_id} not found in graph", node.block_id());
+ };
if !block.alignments().contains_key(node_id) {
- return Err(eyre::eyre!("Node {} not found in block {}", node_id, block.id()));
+ return make_error!("Node {node_id} not found in block {}", block.id());
}
- if !path.nodes.contains(node_id) {
- return Err(eyre::eyre!("Node {} not found in path {}", node_id, path.id()));
+ let Some(path) = self.paths.get(&node.path_id()) else {
+ return make_error!("Path {} of node {node_id} not found in graph", node.path_id());
+ };
+ if walked_by.get(node_id) != Some(&node.path_id()) {
+ return make_error!("Node {node_id} is not in the walk of its path {}", node.path_id());
}
- }
- for (block_id, block) in &self.blocks {
- if block.alignments().is_empty() {
- return Err(eyre::eyre!("Block {} has no nodes", block_id));
+ // Reconstruction rotates a genome by the offset of its first node, so an offset reaching past
+ // the end of that genome would rotate by more than the genome's length. Both ends are
+ // compared, and `tot_len` itself is a valid offset: the last node of a circular path ends
+ // where the genome does.
+ let (start, end) = node.position();
+ if start > path.tot_len() || end > path.tot_len() {
+ return make_error!(
+ "Node {node_id} has position ({start}, {end}), outside its path {} of total length {}",
+ path.id(),
+ path.tot_len()
+ );
}
+ }
- for node_id in block.alignments().keys() {
+ for (block_id, block) in &self.blocks {
+ let consensus_len = block.consensus().len();
+ for (node_id, edits) in block.alignments() {
if !self.nodes.contains_key(node_id) {
- return Err(eyre::eyre!("Node {} not found in graph", node_id));
+ return make_error!("Node {node_id} of block {block_id} not found in graph");
}
+ // `Edit::apply` indexes the consensus by edit position while reconstructing this node.
+ edits
+ .check_bounds(consensus_len)
+ .wrap_err_with(|| format!("When checking the alignment of node {node_id} against block {block_id}"))?;
}
}
- for (path_id, path) in &self.paths {
- for node_id in &path.nodes {
- if !self.nodes.contains_key(node_id) {
- return Err(eyre::eyre!("Node {node_id} from path {path_id} not found in graph"));
- }
+ Ok(())
+ }
+
+ #[cfg(any(test, debug_assertions))]
+ pub fn sanity_check(&self) -> Result<(), Report> {
+ // Referential integrity and the bounds that keep reconstruction in range. Shared with the
+ // release-build validation of graphs read from a file, since a graph that pangraph just built
+ // must satisfy at least as much as one it is willing to load.
+ self.validate()?;
+
+ for (block_id, block) in &self.blocks {
+ if block.alignments().is_empty() {
+ return Err(eyre::eyre!("Block {} has no nodes", block_id));
}
+ }
+ for (path_id, path) in &self.paths {
// // check that there are no duplicated node ids
// // currently disabled because this could rarely happen for empty nodes
// let mut seen = BTreeSet::new();
@@ -303,10 +473,14 @@ mod tests {
#![allow(non_snake_case, clippy::redundant_clone)]
use super::*;
- use crate::pangraph::edits::Edit;
+ use crate::o;
+ use crate::pangraph::edits::{Edit, Sub};
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 crate::utils::error::report_to_string;
+ use itertools::Itertools;
use maplit::btreemap;
use rstest::rstest;
@@ -320,14 +494,14 @@ mod tests {
// b2+ -> [b4+, b5-]
let nodes = btreemap! {
- NodeId(1) => PangraphNode::new(Some(NodeId(1)), BlockId(1), PathId(1), Forward, (0, 0)), // FIXME
- NodeId(2) => PangraphNode::new(Some(NodeId(2)), BlockId(1), PathId(3), Forward, (0, 0)), // FIXME
- NodeId(3) => PangraphNode::new(Some(NodeId(3)), BlockId(2), PathId(1), Forward, (0, 0)), // FIXME
- NodeId(4) => PangraphNode::new(Some(NodeId(4)), BlockId(2), PathId(2), Forward, (0, 0)), // FIXME
- NodeId(5) => PangraphNode::new(Some(NodeId(5)), BlockId(2), PathId(3), Reverse, (0, 0)), // FIXME
- NodeId(6) => PangraphNode::new(Some(NodeId(6)), BlockId(3), PathId(1), Forward, (0, 0)), // FIXME
- NodeId(7) => PangraphNode::new(Some(NodeId(7)), BlockId(3), PathId(2), Forward, (0, 0)), // FIXME
- NodeId(8) => PangraphNode::new(Some(NodeId(8)), BlockId(3), PathId(3), Forward, (0, 0)) // FIXME
+ NodeId(1) => PangraphNode::new(NodeId(1), BlockId(1), PathId(1), Forward, (0, 0)), // FIXME
+ NodeId(2) => PangraphNode::new(NodeId(2), BlockId(1), PathId(3), Forward, (0, 0)), // FIXME
+ NodeId(3) => PangraphNode::new(NodeId(3), BlockId(2), PathId(1), Forward, (0, 0)), // FIXME
+ NodeId(4) => PangraphNode::new(NodeId(4), BlockId(2), PathId(2), Forward, (0, 0)), // FIXME
+ NodeId(5) => PangraphNode::new(NodeId(5), BlockId(2), PathId(3), Reverse, (0, 0)), // FIXME
+ NodeId(6) => PangraphNode::new(NodeId(6), BlockId(3), PathId(1), Forward, (0, 0)), // FIXME
+ NodeId(7) => PangraphNode::new(NodeId(7), BlockId(3), PathId(2), Forward, (0, 0)), // FIXME
+ NodeId(8) => PangraphNode::new(NodeId(8), BlockId(3), PathId(3), Forward, (0, 0)) // FIXME
};
let blocks = btreemap! {
@@ -340,9 +514,9 @@ mod tests {
};
let paths = btreemap! {
- PathId(1) => PangraphPath::new(Some(PathId(1)), [NodeId(1), NodeId(3), NodeId(6)], 0, false, None, None),
- PathId(2) => PangraphPath::new(Some(PathId(2)), [NodeId(4), NodeId(7) ], 0, false, None, None),
- PathId(3) => PangraphPath::new(Some(PathId(3)), [NodeId(2), NodeId(5), NodeId(8)], 0, false, None, None),
+ PathId(1) => PangraphPath::new(PathId(1), [NodeId(1), NodeId(3), NodeId(6)], 0, false, None, None),
+ PathId(2) => PangraphPath::new(PathId(2), [NodeId(4), NodeId(7) ], 0, false, None, None),
+ PathId(3) => PangraphPath::new(PathId(3), [NodeId(2), NodeId(5), NodeId(8)], 0, false, None, None),
};
let mut G = Pangraph {
@@ -352,12 +526,12 @@ mod tests {
};
let new_nodes = btreemap! {
- NodeId(9) => PangraphNode::new(Some(NodeId(9)), BlockId(4), PathId(1), Forward, (0, 0)),
- NodeId(10) => PangraphNode::new(Some(NodeId(10)), BlockId(5), PathId(1), Reverse, (0, 0)),
- NodeId(11) => PangraphNode::new(Some(NodeId(11)), BlockId(4), PathId(2), Forward, (0, 0)),
- NodeId(12) => PangraphNode::new(Some(NodeId(12)), BlockId(5), PathId(2), Reverse, (0, 0)),
- NodeId(13) => PangraphNode::new(Some(NodeId(13)), BlockId(4), PathId(3), Reverse, (0, 0)),
- NodeId(14) => PangraphNode::new(Some(NodeId(14)), BlockId(5), PathId(3), Forward, (0, 0)),
+ NodeId(9) => PangraphNode::new(NodeId(9), BlockId(4), PathId(1), Forward, (0, 0)),
+ NodeId(10) => PangraphNode::new(NodeId(10), BlockId(5), PathId(1), Reverse, (0, 0)),
+ NodeId(11) => PangraphNode::new(NodeId(11), BlockId(4), PathId(2), Forward, (0, 0)),
+ NodeId(12) => PangraphNode::new(NodeId(12), BlockId(5), PathId(2), Reverse, (0, 0)),
+ NodeId(13) => PangraphNode::new(NodeId(13), BlockId(4), PathId(3), Reverse, (0, 0)),
+ NodeId(14) => PangraphNode::new(NodeId(14), BlockId(5), PathId(3), Forward, (0, 0)),
};
let new_blocks = btreemap! {
@@ -386,9 +560,9 @@ mod tests {
assert_eq!(G.blocks, expected_blocks);
let expected_paths = btreemap! {
- PathId(1) => PangraphPath::new(Some(PathId(1)), [NodeId(1), NodeId(9), NodeId(10), NodeId(6)], 0, false, None, None),
- PathId(2) => PangraphPath::new(Some(PathId(2)), [NodeId(11), NodeId(12), NodeId(7) ], 0, false, None, None),
- PathId(3) => PangraphPath::new(Some(PathId(3)), [NodeId(2), NodeId(14), NodeId(13), NodeId(8)], 0, false, None, None),
+ PathId(1) => PangraphPath::new(PathId(1), [NodeId(1), NodeId(9), NodeId(10), NodeId(6)], 0, false, None, None),
+ PathId(2) => PangraphPath::new(PathId(2), [NodeId(11), NodeId(12), NodeId(7) ], 0, false, None, None),
+ PathId(3) => PangraphPath::new(PathId(3), [NodeId(2), NodeId(14), NodeId(13), NodeId(8)], 0, false, None, None),
};
assert_eq!(G.paths, expected_paths);
@@ -415,14 +589,7 @@ mod tests {
.iter()
.enumerate()
.map(|(i, name)| {
- let path = PangraphPath::new(
- Some(PathId(i)),
- Vec::::new(),
- 0,
- false,
- name.map(String::from),
- None,
- );
+ let path = PangraphPath::new(PathId(i), Vec::::new(), 0, false, name.map(String::from), None);
(path.id, path)
})
.collect::>();
@@ -448,4 +615,352 @@ mod tests {
let g = pangraph_with_named_paths(names);
assert_eq!(g.newick_name(), expected);
}
+
+ /// Builds a two-genome graph the way `build` would, with block and node ids seeded from the
+ /// genome names. Two such graphs collide only on their path ids.
+ fn two_genome_graph(names: [&str; 2]) -> Pangraph {
+ let seeds = names.map(genome_seed);
+ let (b0, b1) = (BlockId(seeds[0]), BlockId(seeds[1]));
+ let (n0, n1) = (NodeId(seeds[0]), NodeId(seeds[1]));
+
+ let blocks = btreemap! {
+ b0 => PangraphBlock::new(b0, "ACGTACGT", btreemap!{ n0 => Edit::empty() }),
+ b1 => PangraphBlock::new(b1, "TTTTGGGG", btreemap!{ n1 => Edit::empty() }),
+ };
+ let nodes = btreemap! {
+ n0 => PangraphNode::new(n0, b0, PathId(0), Forward, (0, 8)),
+ n1 => PangraphNode::new(n1, b1, PathId(1), Reverse, (0, 8)),
+ };
+ let paths = btreemap! {
+ PathId(0) => PangraphPath::new(PathId(0), [n0], 8, false, Some(names[0].to_owned()), None),
+ PathId(1) => PangraphPath::new(PathId(1), [n1], 8, false, Some(names[1].to_owned()), None),
+ };
+ Pangraph { paths, blocks, nodes }
+ }
+
+ /// A graph read from JSON is keyed by the object key, while each entity also stores its own id;
+ /// serde never checks that the two agree. `renumber_paths` resolves paths by key, so a
+ /// disagreement would otherwise abort the process with a bare "no entry found for key" panic and
+ /// no indication of which file was at fault.
+ #[rstest]
+ fn test_renumber_paths_reports_node_referring_to_a_missing_path() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let nid = graph.node_ids().next().unwrap();
+ let node = &graph.nodes[&nid];
+ let detached = PangraphNode::new(nid, node.block_id(), PathId(99), node.strand(), node.position());
+ graph.nodes.insert(nid, detached);
+
+ let err = report_to_string(&graph.renumber_paths(7).unwrap_err());
+ assert!(
+ err.contains("refers to path 99, which the graph does not contain"),
+ "unexpected error: {err}"
+ );
+ }
+
+ #[rstest]
+ fn test_renumber_paths_keeps_graph_consistent() {
+ let original = two_genome_graph(["a", "b"]);
+ let graph = original.clone().renumber_paths(7).unwrap();
+
+ graph.sanity_check().unwrap();
+ assert_eq!(graph.blocks.len(), 2);
+ assert_eq!(graph.nodes.len(), 2);
+ assert_eq!(graph.paths.len(), 2);
+
+ // path ids are renumbered contiguously from the offset, in their original order
+ assert_eq!(graph.path_ids().collect_vec(), vec![PathId(7), PathId(8)]);
+ assert_eq!(
+ graph.paths.values().map(|p| p.name.clone()).collect_vec(),
+ vec![Some(o!("a")), Some(o!("b"))]
+ );
+
+ // blocks and node ids are untouched: only the path id a node stores is rewritten
+ assert_eq!(graph.blocks, original.blocks);
+ assert_eq!(graph.node_ids().collect_vec(), original.node_ids().collect_vec());
+ for path in graph.paths.values() {
+ for nid in path.nodes() {
+ assert_eq!(graph.nodes[nid].path_id(), path.id());
+ }
+ }
+ }
+
+ #[rstest]
+ fn test_renumber_paths_preserves_sequences() {
+ let original = two_genome_graph(["a", "b"]);
+ let renumbered = original.clone().renumber_paths(5).unwrap();
+
+ let seqs = |g: &Pangraph| {
+ reconstruct(g)
+ .map(|r| r.map(|r| (r.seq_name, r.seq)))
+ .collect::, Report>>()
+ .unwrap()
+ };
+
+ assert_eq!(seqs(&original), seqs(&renumbered));
+ }
+
+ /// `simplify` drops paths without renumbering the survivors, so the offset cannot be taken to be
+ /// the number of paths.
+ #[rstest]
+ fn test_path_id_upper_bound_clears_sparse_path_ids() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ assert_eq!(graph.path_id_upper_bound(), 2);
+
+ graph.paths.remove(&PathId(0));
+ assert_eq!(graph.paths.len(), 1);
+ assert_eq!(graph.path_id_upper_bound(), 2);
+ }
+
+ #[rstest]
+ fn test_path_id_upper_bound_of_an_empty_graph() {
+ let graph = Pangraph {
+ paths: btreemap! {},
+ blocks: btreemap! {},
+ nodes: btreemap! {},
+ };
+ assert_eq!(graph.path_id_upper_bound(), 0);
+ }
+
+ /// The property the whole identifier model rests on: two graphs built from differently named
+ /// genomes share no block or node id, so `merge` can join them without relabeling either.
+ #[rstest]
+ fn test_graphs_with_distinct_genome_names_are_id_disjoint() {
+ let left = two_genome_graph(["a", "b"]);
+ let right = two_genome_graph(["c", "d"])
+ .renumber_paths(left.path_id_upper_bound())
+ .unwrap();
+
+ assert!(right.is_id_disjoint_from(&left));
+ right.sanity_check().unwrap();
+
+ assert_eq!(left.path_ids().collect_vec(), vec![PathId(0), PathId(1)]);
+ assert_eq!(right.path_ids().collect_vec(), vec![PathId(2), PathId(3)]);
+
+ // joining the two does not conflict
+ let joined = crate::pangraph::graph_merging::graph_join(&left, &right);
+ joined.sanity_check().unwrap();
+ assert_eq!(joined.paths.len(), 4);
+ assert_eq!(joined.blocks.len(), 4);
+ assert_eq!(joined.nodes.len(), 4);
+ }
+
+ /// Appending to a graph that already absorbed another one. This needed a varying relabeling salt
+ /// back when block and node ids were re-derived on merge; deriving them from genome names instead
+ /// makes it hold with no bookkeeping at all.
+ #[rstest]
+ fn test_appending_to_an_already_merged_graph() {
+ let first = two_genome_graph(["a", "b"]);
+ let second = two_genome_graph(["c", "d"])
+ .renumber_paths(first.path_id_upper_bound())
+ .unwrap();
+ let joined = crate::pangraph::graph_merging::graph_join(&first, &second);
+
+ let third = two_genome_graph(["e", "f"])
+ .renumber_paths(joined.path_id_upper_bound())
+ .unwrap();
+
+ assert!(third.is_id_disjoint_from(&joined));
+ third.sanity_check().unwrap();
+ assert_eq!(third.path_ids().collect_vec(), vec![PathId(4), PathId(5)]);
+
+ // and the three of them can be joined without conflicts
+ let joined = crate::pangraph::graph_merging::graph_join(&joined, &third);
+ joined.sanity_check().unwrap();
+ assert_eq!(joined.paths.len(), 6);
+ assert_eq!(joined.blocks.len(), 6);
+ assert_eq!(joined.nodes.len(), 6);
+ }
+
+ /// Two graphs holding the same genome do collide, which is why `merge` rejects them by name
+ /// before it ever gets as far as joining them.
+ #[rstest]
+ fn test_graphs_sharing_a_genome_name_are_not_id_disjoint() {
+ let left = two_genome_graph(["a", "b"]);
+ let right = two_genome_graph(["a", "c"])
+ .renumber_paths(left.path_id_upper_bound())
+ .unwrap();
+
+ assert!(!right.is_id_disjoint_from(&left));
+ }
+
+ /// Ids no longer depend on the order the input sequences were read in, only on their names.
+ #[rstest]
+ fn test_singleton_ids_are_seeded_from_the_name_not_the_index() {
+ let singleton = |name: &str, index: usize| {
+ Pangraph::singleton(
+ FastaRecord {
+ seq_name: name.to_owned(),
+ desc: None,
+ seq: Seq::from_str("ACGTACGT"),
+ index,
+ },
+ Forward,
+ false,
+ )
+ };
+
+ // the same genome read at a different position in the input gets the same block and node ids
+ let first = singleton("a", 0);
+ let shifted = singleton("a", 7);
+ assert_eq!(first.block_ids().collect_vec(), shifted.block_ids().collect_vec());
+ assert_eq!(first.node_ids().collect_vec(), shifted.node_ids().collect_vec());
+
+ // but the path id still records where it was read, so genome order survives
+ assert_eq!(first.path_ids().collect_vec(), vec![PathId(0)]);
+ assert_eq!(shifted.path_ids().collect_vec(), vec![PathId(7)]);
+
+ // and a different genome read at the same position gets different block and node ids
+ let other = singleton("b", 0);
+ assert_ne!(first.block_ids().collect_vec(), other.block_ids().collect_vec());
+ assert_ne!(first.node_ids().collect_vec(), other.node_ids().collect_vec());
+ }
+
+ /// `validate` is the contract every graph read from a file must satisfy, and the reason the code
+ /// downstream is allowed to resolve ids by direct indexing. Each case below is a malformation
+ /// that a hand-edited or third-party JSON graph can carry, and that used to reach an indexing
+ /// panic in a release build, where `sanity_check` is compiled out.
+ #[rstest]
+ fn test_validate_accepts_a_well_formed_graph() {
+ two_genome_graph(["a", "b"]).validate().unwrap();
+ }
+
+ #[rstest]
+ fn test_validate_rejects_path_referring_to_a_missing_node() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ graph.paths.get_mut(&PathId(0)).unwrap().nodes = vec![NodeId(99)];
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains("Node 99 from path 0 not found in graph"),
+ "unexpected error: {err}"
+ );
+ }
+
+ #[rstest]
+ fn test_validate_rejects_node_referring_to_a_missing_block() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let nid = graph.paths[&PathId(0)].nodes[0];
+ let node = &graph.nodes[&nid];
+ graph.nodes.insert(
+ nid,
+ PangraphNode::new(nid, BlockId(99), node.path_id(), node.strand(), node.position()),
+ );
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains(&format!("Block 99 of node {nid} not found in graph")),
+ "unexpected error: {err}"
+ );
+ }
+
+ #[rstest]
+ fn test_validate_rejects_node_referring_to_a_missing_path() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let nid = graph.paths[&PathId(0)].nodes[0];
+ let node = &graph.nodes[&nid];
+ graph.nodes.insert(
+ nid,
+ PangraphNode::new(nid, node.block_id(), PathId(99), node.strand(), node.position()),
+ );
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains(&format!("Path 99 of node {nid} not found in graph")),
+ "unexpected error: {err}"
+ );
+ }
+
+ /// A node whose path exists but does not walk it: reconstruction would never emit this node, and
+ /// its `path_id` claims a genome it is not part of.
+ #[rstest]
+ fn test_validate_rejects_node_missing_from_the_walk_of_its_path() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let nid = graph.paths[&PathId(0)].nodes[0];
+ graph.paths.get_mut(&PathId(0)).unwrap().nodes = vec![];
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains(&format!("Node {nid} is not in the walk of its path 0")),
+ "unexpected error: {err}"
+ );
+ }
+
+ /// Node ids are unique across the graph, so a node claimed by two walks makes its own `path_id`
+ /// meaningless, and would have the node reconstructed into two different genomes.
+ #[rstest]
+ fn test_validate_rejects_node_walked_by_two_paths() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let (first, second) = (graph.paths[&PathId(0)].nodes[0], graph.paths[&PathId(1)].nodes[0]);
+ graph.paths.get_mut(&PathId(1)).unwrap().nodes = vec![first, second];
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains(&format!("Node {first} is walked by both path 0 and path 1")),
+ "unexpected error: {err}"
+ );
+ }
+
+ #[rstest]
+ fn test_validate_rejects_node_missing_from_the_alignments_of_its_block() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let nid = graph.paths[&PathId(0)].nodes[0];
+ let bid = graph.nodes[&nid].block_id();
+ graph
+ .blocks
+ .insert(bid, PangraphBlock::new(bid, "ACGTACGT", btreemap! {}));
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains(&format!("Node {nid} not found in block {bid}")),
+ "unexpected error: {err}"
+ );
+ }
+
+ /// Reconstruction rotates a genome by the offset of its first node. An offset past the end of the
+ /// genome used to reach `rotate_right`, which panics rather than reporting.
+ #[rstest]
+ fn test_validate_rejects_node_position_beyond_the_length_of_its_path() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let nid = graph.paths[&PathId(0)].nodes[0];
+ let node = &graph.nodes[&nid];
+ graph.nodes.insert(
+ nid,
+ PangraphNode::new(nid, node.block_id(), node.path_id(), node.strand(), (100, 8)),
+ );
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains(&format!(
+ "Node {nid} has position (100, 8), outside its path 0 of total length 8"
+ )),
+ "unexpected error: {err}"
+ );
+ }
+
+ /// `Edit::apply` indexes the consensus by edit position, so an out-of-range edit used to panic
+ /// while reconstructing the block, even though `apply` returns a `Result`.
+ #[rstest]
+ fn test_validate_rejects_edit_beyond_the_consensus_of_its_block() {
+ let mut graph = two_genome_graph(["a", "b"]);
+ let nid = graph.paths[&PathId(0)].nodes[0];
+ let bid = graph.nodes[&nid].block_id();
+ let edit = Edit {
+ subs: vec![Sub::new(99, 'T')],
+ dels: vec![],
+ inss: vec![],
+ };
+ graph
+ .blocks
+ .insert(bid, PangraphBlock::new(bid, "ACGTACGT", btreemap! { nid => edit }));
+
+ let err = report_to_string(&graph.validate().unwrap_err());
+ assert!(
+ err.contains("Substitution position 99 is out of bounds for sequence of length 8"),
+ "unexpected error: {err}"
+ );
+ assert!(
+ err.contains(&format!("alignment of node {nid} against block {bid}")),
+ "unexpected error: {err}"
+ );
+ }
}
diff --git a/packages/pangraph/src/pangraph/pangraph_block.rs b/packages/pangraph/src/pangraph/pangraph_block.rs
index 11a755d1..384bbf0e 100644
--- a/packages/pangraph/src/pangraph/pangraph_block.rs
+++ b/packages/pangraph/src/pangraph/pangraph_block.rs
@@ -1,5 +1,5 @@
+use crate::align::alignment_args::GraphMergeParams;
use crate::align::map_variations::{BandParameters, map_variations};
-use crate::commands::build::build_args::PangraphBuildArgs;
use crate::io::fasta::FastaRecord;
use crate::io::json::{JsonPretty, json_write_str};
use crate::io::seq::reverse_complement;
@@ -292,7 +292,7 @@ impl PangraphBlock {
/// Applies a set of edits to the block's consensus sequence and re-aligns the sequences
/// to the new consensus. Returns a new `PangraphBlock` object with the same BlockId.
- pub fn edit_consensus_and_realign(self, edits: &Edit, args: &PangraphBuildArgs) -> Result {
+ pub fn edit_consensus_and_realign(self, edits: &Edit, args: &GraphMergeParams) -> Result {
// apply the edits to the consensus
let new_consensus = edits.apply(&self.consensus)?;
debug_assert!(!new_consensus.is_empty(), "Consensus cannot be empty");
@@ -342,7 +342,7 @@ pub enum RecordNaming {
#[cfg(test)]
mod tests {
use super::*;
- use crate::commands::build::build_args::PangraphBuildArgs;
+ use crate::align::alignment_args::GraphMergeParams;
use crate::pangraph::edits::{Del, Edit, Ins, Sub};
use crate::pangraph::pangraph_node::NodeId;
use maplit::btreemap;
@@ -820,7 +820,7 @@ mod tests {
);
// Create build args with default values for testing
- let args = PangraphBuildArgs::default();
+ let args = GraphMergeParams::default();
// Apply the edits and realign
let result_block = block.edit_consensus_and_realign(&edits, &args).unwrap();
diff --git a/packages/pangraph/src/pangraph/pangraph_node.rs b/packages/pangraph/src/pangraph/pangraph_node.rs
index 225a040d..bc35593f 100644
--- a/packages/pangraph/src/pangraph/pangraph_node.rs
+++ b/packages/pangraph/src/pangraph/pangraph_node.rs
@@ -1,5 +1,5 @@
use crate::pangraph::pangraph_block::BlockId;
-use crate::pangraph::pangraph_path::PathId;
+use crate::pangraph::pangraph_path::{PangraphPath, PathId};
use crate::pangraph::strand::Strand;
use crate::utils::id::id;
use derive_more::{Display, From};
@@ -35,14 +35,13 @@ impl NodeId {
}
impl PangraphNode {
- pub fn new(
- node_id: Option,
- block_id: BlockId,
- path_id: PathId,
- strand: Strand,
- position: (usize, usize),
- ) -> Self {
- let id = node_id.unwrap_or_else(|| id((&block_id, &path_id, &strand, &position)));
+ /// Creates a node with an explicit id.
+ ///
+ /// Node ids are content-derived, and [`PangraphNode::with_derived_id`] is the one place that
+ /// derives them. This constructor takes the id it is given, so a caller that already has an id
+ /// (because the node is being rewritten rather than created) cannot accidentally seed a second,
+ /// divergent derivation scheme.
+ pub fn new(id: NodeId, block_id: BlockId, path_id: PathId, strand: Strand, position: (usize, usize)) -> Self {
Self {
id,
block_id,
@@ -52,6 +51,30 @@ impl PangraphNode {
}
}
+ /// Creates a node placing `block_id` on `path`, with an id derived from its contents.
+ ///
+ /// This is the only place a node id is derived. Genomes are discriminated by
+ /// [`PangraphPath::seed`] rather than by their path id, so the resulting id does not depend on
+ /// the order in which the input sequences were read. The path id is still what the node stores.
+ pub fn with_derived_id(block_id: BlockId, path: &PangraphPath, strand: Strand, position: (usize, usize)) -> Self {
+ Self {
+ id: id((&block_id, &path.seed(), &strand, &position)),
+ block_id,
+ path_id: path.id(),
+ strand,
+ position,
+ }
+ }
+
+ /// Moves this node onto a different path.
+ ///
+ /// The node id derives from the block, the genome *seed*, the strand and the position, none of
+ /// which this touches, so renumbering a path leaves every node id intact. That is why this is the
+ /// one field a node may have rewritten in place; everything else goes through a constructor.
+ pub(crate) fn set_path_id(&mut self, path_id: PathId) {
+ self.path_id = path_id;
+ }
+
// this is almost equivalent to checking if the node is empty
// except for an edge case: when a circular path contains only
// one node. In this case even if the node is not empty, the
diff --git a/packages/pangraph/src/pangraph/pangraph_path.rs b/packages/pangraph/src/pangraph/pangraph_path.rs
index 0c19a10a..2d35035c 100644
--- a/packages/pangraph/src/pangraph/pangraph_path.rs
+++ b/packages/pangraph/src/pangraph/pangraph_path.rs
@@ -32,23 +32,50 @@ pub struct PangraphPath {
}
impl PangraphPath {
+ /// Creates a path with an explicit id.
+ ///
+ /// Path ids are sequential rather than content-derived: they double as the ordering index of the
+ /// genomes, so `build` assigns them from the input order and `renumber_paths` shifts them during
+ /// a merge. The id is therefore always the caller's to supply, and there is no fallback that
+ /// could seed a second, divergent numbering.
pub fn new(
- path_id: Option,
+ id: PathId,
nodes: impl Into>,
tot_len: usize,
circular: bool,
name: Option,
desc: Option,
) -> Self {
- let nodes = nodes.into();
- let id = path_id.unwrap_or_else(|| id((&nodes, &tot_len, &circular, &desc, &name)));
Self {
id,
- nodes,
+ nodes: nodes.into(),
tot_len,
circular,
name,
desc,
}
}
+
+ /// Order-independent identity of the genome on this path, used as the discriminator when deriving
+ /// node ids.
+ ///
+ /// Taken from the genome name rather than from the path id, so that node ids do not depend on the
+ /// order in which the input sequences were read, and do not change when `renumber_paths` shifts
+ /// path ids during a merge. Genome names are unique within a graph, and `merge` rejects graphs
+ /// that share one, so the seed identifies a genome as well as the path id does.
+ ///
+ /// Recomputed on demand rather than stored: a stored field would have to be kept out of the JSON
+ /// and then recomputed on load anyway, and silently defaults to the same value for every path if
+ /// that is forgotten. Falls back to the path id for unnamed paths, which pangraph never writes.
+ pub fn seed(&self) -> usize {
+ self.name.as_deref().map_or(self.id.0, genome_seed)
+ }
+}
+
+/// Derives the identifier seed of a genome from its name.
+///
+/// The single definition of the seeding rule, so that [`PangraphPath::seed`] and
+/// `Pangraph::singleton`, which needs the seed before it has a path to ask, cannot drift apart.
+pub fn genome_seed(name: &str) -> usize {
+ id(name)
}
diff --git a/packages/pangraph/src/pangraph/reconstruct.rs b/packages/pangraph/src/pangraph/reconstruct.rs
new file mode 100644
index 00000000..7ab2425e
--- /dev/null
+++ b/packages/pangraph/src/pangraph/reconstruct.rs
@@ -0,0 +1,775 @@
+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_error, make_internal_report};
+use eyre::{Report, WrapErr};
+use itertools::Itertools;
+use std::collections::{BTreeMap, BTreeSet};
+
+/// 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.
+///
+/// The ordering comes for free: `paths` is a `BTreeMap` keyed by [`PathId`], so iterating it
+/// already yields ascending ids and the records stay lazy.
+///
+/// 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().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,
+ })
+ })
+}
+
+/// Checks that every path of `graphs` carries a non-empty name, and that no name occurs more than
+/// once across all `graphs` taken together.
+///
+/// The genome name is the only identifier that survives a merge unchanged: it is the key that
+/// verification matches on, and what `simplify` resolves genomes by. A graph whose genomes cannot
+/// be told apart by name is therefore rejected rather than processed. A name that is empty, or made
+/// of whitespace only, identifies a genome no better than a missing one, and is reported the same
+/// way. A name padded with whitespace is rejected for the same reason: every comparison pangraph
+/// makes on names is exact, so " a" would be a genome that cannot be addressed by the name it
+/// appears to have, and that does not collide with the "a" of another graph. Names are only tested
+/// here, never trimmed: they are stored exactly as given.
+///
+/// This is the single implementation of that invariant for graphs; [`check_sequence_names`] is its
+/// counterpart for FASTA records. Note that it inspects the paths directly rather than going
+/// through [`reconstruct`], which masks unnamed paths behind a placeholder name.
+///
+/// Path ids are only unique within one graph, so when several graphs are passed the reported ids
+/// are ambiguous on their own. Callers that pass more than one graph are expected to name the
+/// graphs in an error section.
+pub fn check_genome_names(graphs: &[&Pangraph]) -> Result<(), Report> {
+ let unnamed = graphs
+ .iter()
+ .flat_map(|graph| graph.paths.iter())
+ .filter(|(_, path)| path.name.as_deref().is_none_or(|name| name.trim().is_empty()))
+ .map(|(path_id, _)| path_id.to_string())
+ .collect_vec();
+
+ if !unnamed.is_empty() {
+ return make_error!(
+ "Found {} genome(s) with no name or an empty name (path ids: {}). Genomes are identified by name, so every path must be named.",
+ unnamed.len(),
+ format_names(&unnamed)
+ );
+ }
+
+ let padded = padded_names(graphs.iter().flat_map(|graph| graph.path_names().flatten()));
+
+ if !padded.is_empty() {
+ return make_error!(
+ "Found {} genome name(s) with leading or trailing whitespace: {}. Genome names are compared exactly, so surrounding whitespace makes a genome impossible to tell apart from, and impossible to address as, the name it appears to have.",
+ padded.len(),
+ format_names(&padded)
+ );
+ }
+
+ let duplicates = find_duplicates(graphs.iter().flat_map(|graph| graph.path_names().flatten()));
+ if !duplicates.is_empty() {
+ return make_error!(
+ "Duplicate genome names found: {}. Genome names must be unique, because they identify genomes.",
+ format_names(&duplicates)
+ );
+ }
+
+ Ok(())
+}
+
+/// 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, via [`check_genome_names`].
+pub fn path_ids_by_name(graph: &Pangraph) -> Result, Report> {
+ check_genome_names(&[graph])?;
+
+ 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)
+}
+
+/// Checks that every FASTA record carries a non-empty name, and that no two records share one.
+///
+/// The FASTA counterpart of [`check_genome_names`]: sequence names become genome names when the
+/// records are built into a graph, so `build` enforces the invariant on its input to guarantee that
+/// a graph can never carry unusable genome names into a later merge.
+///
+/// Empty names are reported before duplicates, so that a file of headers that are all empty is
+/// reported for the reason it actually has. Names padded with whitespace are rejected too, so that
+/// two records whose headers look alike cannot become two genomes.
+pub fn check_sequence_names(fastas: &[FastaRecord]) -> Result<(), Report> {
+ let empty = fastas
+ .iter()
+ .filter(|fasta| fasta.seq_name.trim().is_empty())
+ .map(|fasta| fasta.index.to_string())
+ .collect_vec();
+
+ if !empty.is_empty() {
+ return make_error!(
+ "Found {} input sequence(s) with an empty name (record indices: {}). Sequences are identified by name, so every record must have one. Note that a space between '>' and the identifier makes the identifier part of the description rather than the name.",
+ empty.len(),
+ format_names(&empty)
+ );
+ }
+
+ let padded = padded_names(fastas.iter().map(|fasta| fasta.seq_name.as_str()));
+
+ if !padded.is_empty() {
+ return make_error!(
+ "Found {} input sequence(s) whose name has leading or trailing whitespace: {}. Sequence names are compared exactly, so surrounding whitespace makes a genome impossible to tell apart from the name it appears to have. Note that whitespace between '>' and the identifier becomes part of the name, unless it is a plain space, which makes the identifier part of the description instead.",
+ padded.len(),
+ format_names(&padded)
+ );
+ }
+
+ 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(())
+}
+
+/// Collects FASTA records into genome sequences keyed by name.
+///
+/// Names are assumed to be unique already: pass the records through [`check_sequence_names`]
+/// first, or records sharing a name will silently collapse into one entry.
+pub(crate) fn sequences_by_name(fastas: &[FastaRecord]) -> BTreeMap {
+ 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()))
+ .map(String::as_str)
+ .collect_vec();
+
+ if !missing.is_empty() {
+ return make_error!(
+ "Graph is missing {} expected genome(s): {}",
+ missing.len(),
+ format_names(&missing)
+ );
+ }
+ }
+
+ Ok(())
+}
+
+/// Checks that `graph` reconstructs exactly the genomes of `sources`, matched by name.
+///
+/// Used to verify a merged graph against the graphs it was built from. Both sides of every
+/// comparison are reconstructed on demand and dropped again, so this holds two genomes at a time
+/// instead of the whole sequence content of `sources`. The saving is proportional to the total
+/// genome length and modest in practice (peak usage during a merge is dominated by the graphs
+/// themselves), but it keeps verification consistent with `reconstruct --verify`, which streams for
+/// the same reason.
+///
+/// Every genome of `sources` must appear in `graph`, and `graph` must contain nothing else. The
+/// `sources` must have disjoint genome names: a name appearing in two of them means the two graphs
+/// describe overlapping genome sets, and is reported rather than verified twice.
+pub fn verify_graph_against_graphs(graph: &Pangraph, sources: &[&Pangraph]) -> Result<(), Report> {
+ check_genome_names(sources).wrap_err("When checking the genome names of the input graphs")?;
+
+ let path_ids = path_ids_by_name(graph)?;
+ let mut verified: BTreeSet<&str> = BTreeSet::new();
+
+ for source in sources {
+ for (name, source_path_id) in path_ids_by_name(source)? {
+ let Some(path_id) = path_ids.get(name) else {
+ return make_error!("Graph is missing genome '{name}', which is present in the input graphs");
+ };
+
+ let expected = reconstruct_genome(source, source_path_id)
+ .wrap_err_with(|| format!("When reconstructing genome '{name}' from the input graphs"))?;
+ let actual =
+ reconstruct_genome(graph, *path_id).wrap_err_with(|| format!("When reconstructing genome '{name}'"))?;
+
+ verify_genome(name, &expected, &actual)?;
+ verified.insert(name);
+ }
+ }
+
+ // Every genome of `sources` has now been found in `graph` and checked, so anything left over is a
+ // genome the merged graph invented.
+ let extra = path_ids
+ .keys()
+ .filter(|name| !verified.contains(*name))
+ .copied()
+ .collect_vec();
+
+ if !extra.is_empty() {
+ return make_error!(
+ "Graph contains {} genome(s) that are not present in the input graphs: {}",
+ extra.len(),
+ format_names(&extra)
+ );
+ }
+
+ Ok(())
+}
+
+/// Collects the names that carry leading or trailing whitespace, quoted so that the padding is
+/// visible in an error message.
+///
+/// Shared by [`check_genome_names`] and [`check_sequence_names`]: both reject a padded name for the
+/// same reason, so what counts as padding is decided in one place. The two keep their own error
+/// messages, which name the container the offending records came from.
+fn padded_names>(names: impl Iterator
- ) -> Vec {
+ names
+ .filter(|name| {
+ let name = name.as_ref();
+ name.trim() != name
+ })
+ .map(|name| format!("{:?}", name.as_ref()))
+ .collect_vec()
+}
+
+/// Formats a list of genome names for an error message, eliding all but the first few.
+pub(crate) 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}]")
+ }
+}
+
+/// Reconstructs the genome of a path, by concatenating the sequences of its nodes and rotating the
+/// result so that it starts where the genome does.
+///
+/// Both lookups here are guaranteed by [`Pangraph::validate`], which every graph read from a file
+/// passes and every graph pangraph builds satisfies, so a failure is a bug rather than bad input.
+/// They are still checked, because reporting beats aborting: indexing the node map and rotating by
+/// an out-of-range offset both panic.
+fn reconstruct_path_sequence(graph: &Pangraph, path: &PangraphPath) -> Result {
+ let Some(first_node_id) = path.nodes.first() else {
+ return Ok(Seq::new());
+ };
+
+ let first_node = graph
+ .nodes
+ .get(first_node_id)
+ .ok_or_else(|| make_internal_report!("Node {first_node_id} not found in graph"))?;
+ let first_node_pos = first_node.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 make_error!(
+ "When reconstructing sequences, genome length mismatch: computed length {} expected {}",
+ genome.len(),
+ genome_len
+ );
+ }
+
+ if first_node_pos > genome.len() {
+ return make_internal_error!(
+ "When reconstructing sequences, the first node of the genome starts at position {first_node_pos}, past the end of a genome of length {}",
+ genome.len()
+ );
+ }
+ genome.rotate_right(first_node_pos);
+
+ Ok(genome)
+}
+
+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;
+ 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(NodeId(0), BlockId(0), PathId(0), Forward, (0, 8)),
+ NodeId(1) => PangraphNode::new(NodeId(1), BlockId(1), PathId(1), Reverse, (0, 8)),
+ };
+ let paths = btreemap! {
+ PathId(0) => PangraphPath::new(PathId(0), [NodeId(0)], 8, false, names[0].map(String::from), None),
+ PathId(1) => PangraphPath::new(PathId(1), [NodeId(1)], 8, false, names[1].map(String::from), None),
+ };
+ Pangraph { paths, blocks, nodes }
+ }
+
+ /// A one-genome graph, so that a pair of them stands in for the two inputs of a merge. Built
+ /// through `Pangraph::singleton`, so these are exactly the graphs `build` starts from and they
+ /// track its id derivation instead of restating it.
+ fn one_genome_graph(name: &str, consensus: &str, strand: Strand) -> Pangraph {
+ let fasta = FastaRecord {
+ seq_name: name.to_owned(),
+ desc: None,
+ seq: Seq::from_str(consensus),
+ index: 0,
+ };
+ Pangraph::singleton(fasta, strand, false)
+ }
+
+ /// The two single-genome graphs whose merger `graph()` stands for.
+ fn sources() -> (Pangraph, Pangraph) {
+ (
+ one_genome_graph("a", "ACGTACGT", Forward),
+ one_genome_graph("b", "TTTTGGGG", Reverse),
+ )
+ }
+
+ 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_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("no name or an empty name"));
+ }
+
+ /// A name that is empty, or made of whitespace only, identifies a genome no better than a missing
+ /// one: `--strains` cannot address it, and it round-trips to FASTA as a nameless record.
+ #[rstest]
+ #[case("")]
+ #[case(" ")]
+ fn test_path_ids_by_name_rejects_empty_path_name(#[case] name: &str) {
+ let graph = two_genome_graph([Some("a"), Some(name)]);
+ let err = report_to_string(&path_ids_by_name(&graph).unwrap_err());
+ assert!(err.contains("no name or an empty name"), "unexpected error: {err}");
+ }
+
+ /// Every comparison on genome names is exact, so " a" would be a genome that `--strains a`
+ /// cannot address, that does not collide with the "a" of the graph it is merged with, and that
+ /// round-trips to a FASTA header with an invisible space in it.
+ #[rstest]
+ #[case(" a")]
+ #[case("a ")]
+ #[case("\ta")]
+ fn test_path_ids_by_name_rejects_padded_path_name(#[case] name: &str) {
+ let graph = two_genome_graph([Some("b"), Some(name)]);
+ let err = report_to_string(&path_ids_by_name(&graph).unwrap_err());
+ assert!(
+ err.contains("leading or trailing whitespace"),
+ "unexpected error: {err}"
+ );
+ }
+
+ /// The pair a padded name is most likely to be mistaken for: without the check they are two
+ /// distinct genomes, since duplicates are detected on the name as given.
+ #[rstest]
+ fn test_path_ids_by_name_rejects_a_name_that_differs_only_by_padding() {
+ let graph = two_genome_graph([Some("a"), Some("a ")]);
+ let err = report_to_string(&path_ids_by_name(&graph).unwrap_err());
+ assert!(
+ err.contains("leading or trailing whitespace"),
+ "unexpected error: {err}"
+ );
+ }
+
+ #[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"));
+ }
+
+ /// A name may be unique within each graph and still collide across them. This is what makes
+ /// `merge` reject merging a graph with itself, and what stops two overlapping source graphs from
+ /// being "verified" against a merged graph that holds their genomes only once.
+ #[rstest]
+ fn test_check_genome_names_rejects_name_shared_across_graphs() {
+ let (left, _) = sources();
+ let other = one_genome_graph("a", "GGGGCCCC", Forward);
+
+ check_genome_names(&[&left]).unwrap();
+ let err = report_to_string(&check_genome_names(&[&left, &other]).unwrap_err());
+ assert!(err.contains("Duplicate genome names"), "unexpected error: {err}");
+ assert!(err.contains('a'), "unexpected error: {err}");
+ }
+
+ #[rstest]
+ fn test_check_sequence_names_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(&check_sequence_names(&fastas).unwrap_err()).contains("Duplicate sequence names"));
+ }
+
+ /// `> id`, with a space between the '>' and the identifier, parses into an empty name with the
+ /// identifier in the description, so this common header style used to yield a nameless genome.
+ /// The record index is part of the message, since the name cannot point at the offending record.
+ #[rstest]
+ #[case("")]
+ #[case(" ")]
+ fn test_check_sequence_names_rejects_empty_name(#[case] name: &str) {
+ let fastas = [FastaRecord {
+ seq_name: name.to_owned(),
+ desc: Some(o!("NC_000913.3 Escherichia coli")),
+ seq: Seq::from_str("ACGT"),
+ index: 3,
+ }];
+
+ let err = report_to_string(&check_sequence_names(&fastas).unwrap_err());
+ assert!(err.contains("empty name"), "unexpected error: {err}");
+ assert!(err.contains('3'), "unexpected error: {err}");
+ }
+
+ /// A header indented with anything other than a plain space, such as `>\tid`, keeps the whitespace
+ /// in the name, which is invisible in every place the name is later shown.
+ #[rstest]
+ #[case(" a")]
+ #[case("a ")]
+ #[case("\ta")]
+ fn test_check_sequence_names_rejects_padded_name(#[case] name: &str) {
+ let fastas = [FastaRecord {
+ seq_name: name.to_owned(),
+ desc: None,
+ seq: Seq::from_str("ACGT"),
+ index: 3,
+ }];
+
+ let err = report_to_string(&check_sequence_names(&fastas).unwrap_err());
+ assert!(
+ err.contains("leading or trailing whitespace"),
+ "unexpected error: {err}"
+ );
+ }
+
+ #[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);
+ 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_graph_against_graphs_accepts_exact_match() {
+ let (left, right) = sources();
+ verify_graph_against_graphs(&graph(), &[&left, &right]).unwrap();
+ }
+
+ #[rstest]
+ fn test_verify_graph_against_graphs_detects_missing_genome() {
+ let (left, right) = sources();
+ let extra = one_genome_graph("c", "GGGGCCCC", Forward);
+
+ let err = report_to_string(&verify_graph_against_graphs(&graph(), &[&left, &right, &extra]).unwrap_err());
+ assert!(err.contains("missing genome 'c'"), "unexpected error: {err}");
+ }
+
+ #[rstest]
+ fn test_verify_graph_against_graphs_detects_extra_genome() {
+ let (left, _) = sources();
+
+ let err = report_to_string(&verify_graph_against_graphs(&graph(), &[&left]).unwrap_err());
+ assert!(
+ err.contains("not present in the input graphs"),
+ "unexpected error: {err}"
+ );
+ assert!(err.contains('b'), "unexpected error: {err}");
+ }
+
+ /// Two sources holding the same genome name describe overlapping genome sets. Counting one
+ /// expectation per (source, name) pair used to make the totals disagree while no genome was
+ /// actually extra, reporting the nonsensical "contains 0 genome(s) that are not present ...: []".
+ #[rstest]
+ fn test_verify_graph_against_graphs_detects_name_shared_across_sources() {
+ let (left, right) = sources();
+ let duplicate = one_genome_graph("a", "ACGTACGT", Forward);
+
+ let err = report_to_string(&verify_graph_against_graphs(&graph(), &[&left, &right, &duplicate]).unwrap_err());
+ assert!(err.contains("Duplicate genome names"), "unexpected error: {err}");
+ assert!(!err.contains("0 genome(s)"), "unexpected error: {err}");
+ }
+
+ /// The name matches but the sequence does not: the mismatch is reported against the source graph
+ /// the genome came from, with the position of the first difference.
+ #[rstest]
+ fn test_verify_graph_against_graphs_detects_mutated_genome() {
+ let (_, right) = sources();
+ let mutated = one_genome_graph("a", "ACGTTCGT", Forward);
+
+ let err = report_to_string(&verify_graph_against_graphs(&graph(), &[&mutated, &right]).unwrap_err());
+ assert!(
+ err.contains("Sequence mismatch for genome 'a'"),
+ "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]");
+ }
+
+ /// `Pangraph::validate` rejects both of these before reconstruction is ever reached, so they
+ /// stand for a graph assembled in-process rather than read from a file. They are still worth
+ /// pinning: both used to be indexing panics, which abort the process instead of being reported,
+ /// and `sanity_check` is compiled out of release builds.
+ #[rstest]
+ fn test_reconstruct_reports_path_referring_to_a_missing_node() {
+ let mut graph = graph();
+ graph.paths.get_mut(&PathId(0)).unwrap().nodes = vec![NodeId(99)];
+
+ let err = report_to_string(&reconstruct_genome(&graph, PathId(0)).unwrap_err());
+ assert!(err.contains("Node 99 not found in graph"), "unexpected error: {err}");
+ }
+
+ #[rstest]
+ fn test_reconstruct_reports_first_node_starting_past_the_end_of_the_genome() {
+ let mut graph = graph();
+ let node = &graph.nodes[&NodeId(0)];
+ graph.nodes.insert(
+ NodeId(0),
+ PangraphNode::new(NodeId(0), node.block_id(), node.path_id(), node.strand(), (100, 8)),
+ );
+
+ let err = report_to_string(&reconstruct_genome(&graph, PathId(0)).unwrap_err());
+ assert!(
+ err.contains("starts at position 100, past the end of a genome of length 8"),
+ "unexpected error: {err}"
+ );
+ }
+}
diff --git a/packages/pangraph/src/pangraph/reweave.rs b/packages/pangraph/src/pangraph/reweave.rs
index 8364e183..8c83cf42 100644
--- a/packages/pangraph/src/pangraph/reweave.rs
+++ b/packages/pangraph/src/pangraph/reweave.rs
@@ -1,7 +1,7 @@
use crate::align::alignment::{Alignment, AnchorBlock, ExtractedHit};
+use crate::align::alignment_args::GraphMergeParams;
use crate::align::bam::cigar::{Side, add_flanking_indel, cigar_switch_ref_qry, invert_cigar};
use crate::align::map_variations::{BandParameters, map_variations};
-use crate::commands::build::build_args::PangraphBuildArgs;
use crate::io::seq::reverse_complement;
use crate::make_internal_error;
use crate::pangraph::edits::Edit;
@@ -37,7 +37,7 @@ impl MergePromise {
}
}
- pub fn solve_promise(&mut self, args: &PangraphBuildArgs) -> Result {
+ pub fn solve_promise(&mut self, args: &GraphMergeParams) -> Result {
// TODO: avoid re-aligning if cigar is only a single match (no indels)
// calculate the mean shift and bandwidth of the alignment due to the displacement
@@ -704,9 +704,9 @@ mod tests {
let nid2 = NodeId(2000);
let nid3 = NodeId(3000);
- let n1 = PangraphNode::new(Some(nid1), bid, PathId(100), Forward, (100, 230));
- let n2 = PangraphNode::new(Some(nid2), bid, PathId(200), Reverse, (1000, 1130));
- let n3 = PangraphNode::new(Some(nid3), bid, PathId(300), Reverse, (180, 110));
+ let n1 = PangraphNode::new(nid1, bid, PathId(100), Forward, (100, 230));
+ let n2 = PangraphNode::new(nid2, bid, PathId(200), Reverse, (1000, 1130));
+ let n3 = PangraphNode::new(nid3, bid, PathId(300), Reverse, (180, 110));
let b1 = PangraphBlock::new(
bid,
@@ -718,9 +718,9 @@ mod tests {
},
);
- let p1 = PangraphPath::new(Some(PathId(100)), [nid1], 2000, true, None, None);
- let p2 = PangraphPath::new(Some(PathId(200)), [nid2], 2000, true, None, None);
- let p3 = PangraphPath::new(Some(PathId(300)), [nid3], 200, true, None, None);
+ let p1 = PangraphPath::new(PathId(100), [nid1], 2000, true, None, None);
+ let p2 = PangraphPath::new(PathId(200), [nid2], 2000, true, None, None);
+ let p3 = PangraphPath::new(PathId(300), [nid3], 200, true, None, None);
let G = Pangraph {
paths: btreemap! {
@@ -885,20 +885,20 @@ mod tests {
fn generate_example() -> (Pangraph, Vec) {
let nodes = btreemap! {
- NodeId(1) => PangraphNode::new(Some(NodeId(1)), BlockId(10), PathId(100), Forward, (700, 885)),
- NodeId(2) => PangraphNode::new(Some(NodeId(2)), BlockId(30), PathId(100), Forward, (885, 988)),
- NodeId(3) => PangraphNode::new(Some(NodeId(3)), BlockId(30), PathId(200), Reverse, (100, 180)),
- NodeId(4) => PangraphNode::new(Some(NodeId(4)), BlockId(20), PathId(200), Reverse, (180, 555)),
- NodeId(5) => PangraphNode::new(Some(NodeId(5)), BlockId(10), PathId(200), Reverse, (555, 735)),
- NodeId(6) => PangraphNode::new(Some(NodeId(6)), BlockId(40), PathId(300), Forward, (600, 100)),
- NodeId(7) => PangraphNode::new(Some(NodeId(7)), BlockId(50), PathId(300), Forward, (100, 325)),
- NodeId(8) => PangraphNode::new(Some(NodeId(8)), BlockId(50), PathId(300), Reverse, (325, 580)),
+ NodeId(1) => PangraphNode::new(NodeId(1), BlockId(10), PathId(100), Forward, (700, 885)),
+ NodeId(2) => PangraphNode::new(NodeId(2), BlockId(30), PathId(100), Forward, (885, 988)),
+ NodeId(3) => PangraphNode::new(NodeId(3), BlockId(30), PathId(200), Reverse, (100, 180)),
+ NodeId(4) => PangraphNode::new(NodeId(4), BlockId(20), PathId(200), Reverse, (180, 555)),
+ NodeId(5) => PangraphNode::new(NodeId(5), BlockId(10), PathId(200), Reverse, (555, 735)),
+ NodeId(6) => PangraphNode::new(NodeId(6), BlockId(40), PathId(300), Forward, (600, 100)),
+ NodeId(7) => PangraphNode::new(NodeId(7), BlockId(50), PathId(300), Forward, (100, 325)),
+ NodeId(8) => PangraphNode::new(NodeId(8), BlockId(50), PathId(300), Reverse, (325, 580)),
};
let paths = btreemap! {
- PathId(100) => PangraphPath::new(Some(PathId(100)), [NodeId(1), NodeId(2)], 1000, true, None, None),
- PathId(200) => PangraphPath::new(Some(PathId(200)), [NodeId(3), NodeId(4), NodeId(5)], 1000, true, None, None),
- PathId(300) => PangraphPath::new(Some(PathId(300)), [NodeId(6), NodeId(7), NodeId(8)], 1000, true, None, None),
+ PathId(100) => PangraphPath::new(PathId(100), [NodeId(1), NodeId(2)], 1000, true, None, None),
+ PathId(200) => PangraphPath::new(PathId(200), [NodeId(3), NodeId(4), NodeId(5)], 1000, true, None, None),
+ PathId(300) => PangraphPath::new(PathId(300), [NodeId(6), NodeId(7), NodeId(8)], 1000, true, None, None),
};
#[rustfmt::skip]
diff --git a/packages/pangraph/src/pangraph/slice.rs b/packages/pangraph/src/pangraph/slice.rs
index 620f9458..21143142 100644
--- a/packages/pangraph/src/pangraph/slice.rs
+++ b/packages/pangraph/src/pangraph/slice.rs
@@ -167,16 +167,17 @@ pub fn block_slice(
old_strandedness
};
- let path_L = G.paths[&old_node.path_id()].tot_len;
+ let path = &G.paths[&old_node.path_id()];
+ let path_L = path.tot_len;
let node_coords = interval_node_coords(i, edits, block_L);
- let circular = G.paths[&old_node.path_id()].circular();
+ let circular = path.circular();
let new_pos = if circular {
new_position_circular(old_node.position(), node_coords, path_L, old_strandedness)
} else {
new_position_non_circular(old_node.position(), node_coords, old_strandedness)
};
- let new_node = PangraphNode::new(None, i.new_block_id, old_node.path_id(), new_strand, new_pos);
+ let new_node = PangraphNode::with_derived_id(i.new_block_id, path, new_strand, new_pos);
// extract edits for the slice
let new_edits = slice_edits(i, edits, block_L);
@@ -467,15 +468,15 @@ mod tests {
let (new_b, new_nodes) = block_slice(&b, &i, &G);
assert_eq!(new_b.consensus(), "TATATTTATC");
- let nn1 = PangraphNode::new(None, new_bid, PathId(1), Forward, (111, 120));
+ let nn1 = PangraphNode::with_derived_id(new_bid, &G.paths[&PathId(1)], Forward, (111, 120));
let nn1_slice = new_nodes[&NodeId(1)].as_ref().unwrap();
assert_eq!(nn1, *nn1_slice);
- let nn2 = PangraphNode::new(None, new_bid, PathId(2), Reverse, (1008, 1017));
+ let nn2 = PangraphNode::with_derived_id(new_bid, &G.paths[&PathId(2)], Reverse, (1008, 1017));
let nn2_slice = new_nodes[&NodeId(2)].as_ref().unwrap();
assert_eq!(nn2, *nn2_slice);
- let nn3 = PangraphNode::new(None, new_bid, PathId(3), Reverse, (96, 4));
+ let nn3 = PangraphNode::with_derived_id(new_bid, &G.paths[&PathId(3)], Reverse, (96, 4));
let nn3_slice = new_nodes[&NodeId(3)].as_ref().unwrap();
assert_eq!(nn3, *nn3_slice);
@@ -532,34 +533,13 @@ mod tests {
inss: vec![Ins::new(20, "T")],
};
- let n1 = PangraphNode::new(Some(NodeId(1)), bid, PathId(1), Forward, (100, 125));
- let n2 = PangraphNode::new(Some(NodeId(2)), bid, PathId(2), Reverse, (1000, 1025));
- let n3 = PangraphNode::new(Some(NodeId(3)), bid, PathId(3), Reverse, (90, 9));
-
- let p1 = PangraphPath::new(
- Some(PathId(1)),
- /*"p1"*/ [NodeId(1), NodeId(4)],
- 2000,
- true,
- None,
- None,
- );
- let p2 = PangraphPath::new(
- Some(PathId(2)),
- /*"p2"*/ [NodeId(2), NodeId(5)],
- 2000,
- true,
- None,
- None,
- );
- let p3 = PangraphPath::new(
- Some(PathId(3)),
- /*"p3"*/ [NodeId(3), NodeId(6)],
- 100,
- true,
- None,
- None,
- );
+ let n1 = PangraphNode::new(NodeId(1), bid, PathId(1), Forward, (100, 125));
+ let n2 = PangraphNode::new(NodeId(2), bid, PathId(2), Reverse, (1000, 1025));
+ let n3 = PangraphNode::new(NodeId(3), bid, PathId(3), Reverse, (90, 9));
+
+ let p1 = PangraphPath::new(PathId(1), /*"p1"*/ [NodeId(1), NodeId(4)], 2000, true, None, None);
+ let p2 = PangraphPath::new(PathId(2), /*"p2"*/ [NodeId(2), NodeId(5)], 2000, true, None, None);
+ let p3 = PangraphPath::new(PathId(3), /*"p3"*/ [NodeId(3), NodeId(6)], 100, true, None, None);
let b1 = PangraphBlock::new(
bid,
@@ -609,15 +589,15 @@ mod tests {
assert_eq!(new_b.consensus(), "TATATTTATC");
- let nn1 = PangraphNode::new(None, new_bid, PathId(1), Reverse, (111, 120));
+ let nn1 = PangraphNode::with_derived_id(new_bid, &G.paths[&PathId(1)], Reverse, (111, 120));
let nn1_slice = new_nodes[&NodeId(1)].as_ref().unwrap();
assert_eq!(nn1, *nn1_slice);
- let nn2 = PangraphNode::new(None, new_bid, PathId(2), Forward, (1008, 1017));
+ let nn2 = PangraphNode::with_derived_id(new_bid, &G.paths[&PathId(2)], Forward, (1008, 1017));
let nn2_slice = new_nodes[&NodeId(2)].as_ref().unwrap();
assert_eq!(nn2, *nn2_slice);
- let nn3 = PangraphNode::new(None, new_bid, PathId(3), Forward, (96, 4));
+ let nn3 = PangraphNode::with_derived_id(new_bid, &G.paths[&PathId(3)], Forward, (96, 4));
let nn3_slice = new_nodes[&NodeId(3)].as_ref().unwrap();
assert_eq!(nn3, *nn3_slice);
diff --git a/packages/pangraph/src/reconsensus/reconsensus.rs b/packages/pangraph/src/reconsensus/reconsensus.rs
index 8b410b45..e7d4afce 100644
--- a/packages/pangraph/src/reconsensus/reconsensus.rs
+++ b/packages/pangraph/src/reconsensus/reconsensus.rs
@@ -1,4 +1,4 @@
-use crate::commands::build::build_args::PangraphBuildArgs;
+use crate::align::alignment_args::GraphMergeParams;
use crate::make_report;
use crate::pangraph::detach_unaligned::detach_unaligned_nodes;
use crate::pangraph::edits::Edit;
@@ -32,7 +32,7 @@ struct BlockAnalysis {
pub fn reconsensus_graph(
graph: &mut Pangraph,
ids_updated_blocks: &[BlockId],
- args: &PangraphBuildArgs,
+ args: &GraphMergeParams,
) -> Result<(), Report> {
// there should be no empty nodes in the graph
debug_assert!(
@@ -406,7 +406,7 @@ mod tests {
let majority_edits = block.find_majority_edits();
assert!(majority_edits.has_indels()); // This block has indels requiring re-alignment
let block = block
- .edit_consensus_and_realign(&majority_edits, &PangraphBuildArgs::default())
+ .edit_consensus_and_realign(&majority_edits, &GraphMergeParams::default())
.unwrap();
// Check that the re-alignment produced the expected result
@@ -422,7 +422,7 @@ mod tests {
let majority_edits = block.find_majority_edits();
assert!(majority_edits.has_indels()); // This block has indels requiring re-alignment
let block = block
- .edit_consensus_and_realign(&majority_edits, &PangraphBuildArgs::default())
+ .edit_consensus_and_realign(&majority_edits, &GraphMergeParams::default())
.unwrap();
// Check that the re-alignment produced the expected result
@@ -436,18 +436,18 @@ mod tests {
let expected_block = block_1_reconsensus();
let nodes = btreemap! {
- NodeId(1) => PangraphNode::new(Some(NodeId(1)), block.id(), PathId(1), Forward, (0, 23)),
- NodeId(2) => PangraphNode::new(Some(NodeId(2)), block.id(), PathId(2), Forward, (0, 23)),
- NodeId(3) => PangraphNode::new(Some(NodeId(3)), block.id(), PathId(3), Forward, (0, 23)),
- NodeId(4) => PangraphNode::new(Some(NodeId(4)), block.id(), PathId(4), Forward, (0, 23)),
- NodeId(5) => PangraphNode::new(Some(NodeId(5)), block.id(), PathId(5), Forward, (0, 23)),
+ NodeId(1) => PangraphNode::new(NodeId(1), block.id(), PathId(1), Forward, (0, 23)),
+ NodeId(2) => PangraphNode::new(NodeId(2), block.id(), PathId(2), Forward, (0, 23)),
+ NodeId(3) => PangraphNode::new(NodeId(3), block.id(), PathId(3), Forward, (0, 23)),
+ NodeId(4) => PangraphNode::new(NodeId(4), block.id(), PathId(4), Forward, (0, 23)),
+ NodeId(5) => PangraphNode::new(NodeId(5), block.id(), PathId(5), Forward, (0, 23)),
};
let paths = btreemap! {
- PathId(1) => PangraphPath::new(Some(PathId(1)), [NodeId(1)], 23, false, None, None),
- PathId(2) => PangraphPath::new(Some(PathId(2)), [NodeId(2)], 23, false, None, None),
- PathId(3) => PangraphPath::new(Some(PathId(3)), [NodeId(3)], 23, false, None, None),
- PathId(4) => PangraphPath::new(Some(PathId(4)), [NodeId(4)], 23, false, None, None),
- PathId(5) => PangraphPath::new(Some(PathId(5)), [NodeId(5)], 23, false, None, None),
+ PathId(1) => PangraphPath::new(PathId(1), [NodeId(1)], 23, false, None, None),
+ PathId(2) => PangraphPath::new(PathId(2), [NodeId(2)], 23, false, None, None),
+ PathId(3) => PangraphPath::new(PathId(3), [NodeId(3)], 23, false, None, None),
+ PathId(4) => PangraphPath::new(PathId(4), [NodeId(4)], 23, false, None, None),
+ PathId(5) => PangraphPath::new(PathId(5), [NodeId(5)], 23, false, None, None),
};
let mut graph = Pangraph {
blocks: btreemap! {
@@ -457,7 +457,7 @@ mod tests {
paths,
};
- let result = reconsensus_graph(&mut graph, &[block_id], &PangraphBuildArgs::default());
+ let result = reconsensus_graph(&mut graph, &[block_id], &GraphMergeParams::default());
result.unwrap();
assert_eq!(graph.blocks[&block_id], expected_block);
@@ -510,20 +510,20 @@ mod tests {
// Create nodes for the block with lengths reflecting actual sequence lengths
let nodes = btreemap! {
- NodeId(1) => PangraphNode::new(Some(NodeId(1)), initial_block.id(), PathId(1), Reverse, (0, 10)), // 50 - 40 = 9 (deletes positions 0-39)
- NodeId(2) => PangraphNode::new(Some(NodeId(2)), initial_block.id(), PathId(2), Forward, (0, 35)), // 50 - 15 = 35 (deletes positions 35-49)
- NodeId(3) => PangraphNode::new(Some(NodeId(3)), initial_block.id(), PathId(3), Forward, (0, 35)), // 50 - 15 = 35 (deletes positions 35-49)
- NodeId(4) => PangraphNode::new(Some(NodeId(4)), initial_block.id(), PathId(4), Forward, (0, 35)), // 50 - 15 = 35 (deletes positions 35-49)
- NodeId(5) => PangraphNode::new(Some(NodeId(5)), initial_block.id(), PathId(5), Forward, (0, 49)), // no deletions
+ NodeId(1) => PangraphNode::new(NodeId(1), initial_block.id(), PathId(1), Reverse, (0, 10)), // 50 - 40 = 9 (deletes positions 0-39)
+ NodeId(2) => PangraphNode::new(NodeId(2), initial_block.id(), PathId(2), Forward, (0, 35)), // 50 - 15 = 35 (deletes positions 35-49)
+ NodeId(3) => PangraphNode::new(NodeId(3), initial_block.id(), PathId(3), Forward, (0, 35)), // 50 - 15 = 35 (deletes positions 35-49)
+ NodeId(4) => PangraphNode::new(NodeId(4), initial_block.id(), PathId(4), Forward, (0, 35)), // 50 - 15 = 35 (deletes positions 35-49)
+ NodeId(5) => PangraphNode::new(NodeId(5), initial_block.id(), PathId(5), Forward, (0, 49)), // no deletions
};
// Create paths
let paths = btreemap! {
- PathId(1) => PangraphPath::new(Some(PathId(1)), [NodeId(1)], 49, false, None, None),
- PathId(2) => PangraphPath::new(Some(PathId(2)), [NodeId(2)], 49, false, None, None),
- PathId(3) => PangraphPath::new(Some(PathId(3)), [NodeId(3)], 49, false, None, None),
- PathId(4) => PangraphPath::new(Some(PathId(4)), [NodeId(4)], 49, false, None, None),
- PathId(5) => PangraphPath::new(Some(PathId(5)), [NodeId(5)], 49, false, None, None),
+ PathId(1) => PangraphPath::new(PathId(1), [NodeId(1)], 49, false, None, None),
+ PathId(2) => PangraphPath::new(PathId(2), [NodeId(2)], 49, false, None, None),
+ PathId(3) => PangraphPath::new(PathId(3), [NodeId(3)], 49, false, None, None),
+ PathId(4) => PangraphPath::new(PathId(4), [NodeId(4)], 49, false, None, None),
+ PathId(5) => PangraphPath::new(PathId(5), [NodeId(5)], 49, false, None, None),
};
// Create blocks map
@@ -535,7 +535,7 @@ mod tests {
let mut graph = Pangraph { paths, blocks, nodes };
// Apply reconsensus_graph
- let result = reconsensus_graph(&mut graph, &[initial_block.id()], &PangraphBuildArgs::default());
+ let result = reconsensus_graph(&mut graph, &[initial_block.id()], &GraphMergeParams::default());
// Check that the operation succeeded
result.unwrap();
@@ -554,7 +554,7 @@ mod tests {
// check that the node was updated correctly, flipping the strandedness
let new_node1 = &graph.nodes[&NodeId(1)];
- let expected_node1 = PangraphNode::new(Some(NodeId(1)), singleton_block_exp.id(), PathId(1), Forward, (0, 10));
+ let expected_node1 = PangraphNode::new(NodeId(1), singleton_block_exp.id(), PathId(1), Forward, (0, 10));
assert_eq!(new_node1, &expected_node1);
}
}
diff --git a/packages/pangraph/src/reconsensus/remove_nodes.rs b/packages/pangraph/src/reconsensus/remove_nodes.rs
index cf9b88b3..929c5f18 100644
--- a/packages/pangraph/src/reconsensus/remove_nodes.rs
+++ b/packages/pangraph/src/reconsensus/remove_nodes.rs
@@ -87,17 +87,17 @@ mod tests {
fn create_input_graph() -> Pangraph {
let nodes = btreemap! {
- NodeId(1) => PangraphNode::new(Some(NodeId(1)), BlockId(1), PathId(0), Forward, (0, 10)),
- NodeId(2) => PangraphNode::new(Some(NodeId(2)), BlockId(1), PathId(1), Forward, (0, 10)),
- NodeId(3) => PangraphNode::new(Some(NodeId(3)), BlockId(1), PathId(2), Reverse, (0, 0)),
- NodeId(4) => PangraphNode::new(Some(NodeId(4)), BlockId(2), PathId(0), Forward, (10, 20)),
- NodeId(5) => PangraphNode::new(Some(NodeId(5)), BlockId(2), PathId(2), Forward, (0, 10)),
+ NodeId(1) => PangraphNode::new(NodeId(1), BlockId(1), PathId(0), Forward, (0, 10)),
+ NodeId(2) => PangraphNode::new(NodeId(2), BlockId(1), PathId(1), Forward, (0, 10)),
+ NodeId(3) => PangraphNode::new(NodeId(3), BlockId(1), PathId(2), Reverse, (0, 0)),
+ NodeId(4) => PangraphNode::new(NodeId(4), BlockId(2), PathId(0), Forward, (10, 20)),
+ NodeId(5) => PangraphNode::new(NodeId(5), BlockId(2), PathId(2), Forward, (0, 10)),
};
let paths = btreemap! {
- PathId(0) => PangraphPath::new(Some(PathId(1)), vec![NodeId(1), NodeId(4)], 20, false, None, None),
- PathId(1) => PangraphPath::new(Some(PathId(2)), vec![NodeId(2)], 10, false, None, None),
- PathId(2) => PangraphPath::new(Some(PathId(3)), vec![NodeId(3), NodeId(5)], 10, false, None, None),
+ PathId(0) => PangraphPath::new(PathId(1), vec![NodeId(1), NodeId(4)], 20, false, None, None),
+ PathId(1) => PangraphPath::new(PathId(2), vec![NodeId(2)], 10, false, None, None),
+ PathId(2) => PangraphPath::new(PathId(3), vec![NodeId(3), NodeId(5)], 10, false, None, None),
};
let blocks = btreemap! {
diff --git a/packages/pangraph/src/tree/neighbor_joining.rs b/packages/pangraph/src/tree/neighbor_joining.rs
index 8563d54e..76db403a 100644
--- a/packages/pangraph/src/tree/neighbor_joining.rs
+++ b/packages/pangraph/src/tree/neighbor_joining.rs
@@ -2,6 +2,7 @@
use crate::distance::mash::mash_distance::mash_distance;
use crate::distance::mash::minimizer::MinimizersParams;
+use crate::make_error;
use crate::pangraph::pangraph::Pangraph;
// use crate::tree::balance::balance;
use crate::tree::clade::Clade;
@@ -13,7 +14,17 @@ use ndarray::{Array1, Array2, Axis, s};
use ndarray_stats::QuantileExt;
/// Generate guide tree using neighbor joining method.
-pub fn build_tree_using_neighbor_joining(graphs: Vec) -> Result>>, Report> {
+pub fn build_tree_using_neighbor_joining(mut graphs: Vec) -> Result>>, Report> {
+ match graphs.len() {
+ 0 => {
+ return make_error!("When building the guide tree: expected at least one input genome, but none were provided");
+ },
+ // A single genome needs no joining: the guide tree is that genome's leaf. Handled here, before
+ // computing pairwise distances, which would be wasted work. `pop` yields exactly the `Some` we need.
+ 1 => return Ok(Lock::new(Clade::new(graphs.pop()))),
+ _ => {},
+ }
+
let mut distances = calculate_distances(&graphs);
let mut nodes = graphs
@@ -102,12 +113,56 @@ fn join_in_place(D: &mut Array2, nodes: &mut Vec>
#[cfg(test)]
mod tests {
use super::*;
+ use crate::assert_error;
+ use crate::io::fasta::FastaRecord;
+ use crate::pangraph::strand::Strand::Forward;
+ use crate::representation::seq::Seq;
use ndarray::array;
use pretty_assertions::assert_eq;
use rstest::rstest;
const INF: f64 = f64::INFINITY;
+ /// Builds a singleton graph for one named sequence, as `build` does for each input FASTA record.
+ fn singleton(name: &str, index: usize) -> Pangraph {
+ Pangraph::singleton(
+ FastaRecord {
+ seq_name: name.to_owned(),
+ desc: None,
+ seq: Seq::from_str("ACGTACGTACGTACGT"),
+ index,
+ },
+ Forward,
+ false,
+ )
+ }
+
+ #[rstest]
+ fn test_build_tree_rejects_no_graphs() {
+ assert_error!(
+ build_tree_using_neighbor_joining(vec![]),
+ "When building the guide tree: expected at least one input genome, but none were provided"
+ );
+ }
+
+ #[rstest]
+ fn test_build_tree_single_graph_is_a_leaf() {
+ let tree = build_tree_using_neighbor_joining(vec![singleton("A", 0)]).unwrap();
+ let root = tree.read();
+ assert!(root.is_leaf());
+ assert!(root.data.is_some());
+ assert_eq!(root.to_newick(), "A;");
+ }
+
+ #[rstest]
+ fn test_build_tree_two_graphs() {
+ let tree = build_tree_using_neighbor_joining(vec![singleton("A", 0), singleton("B", 1)]).unwrap();
+ let root = tree.read();
+ assert!(!root.is_leaf());
+ assert!(root.data.is_none());
+ assert_eq!(root.to_newick(), "(A,B);");
+ }
+
#[rstest]
fn test_create_Q_matrix() {
// example from wikipedia: https://en.wikipedia.org/wiki/Neighbor_joining#Example
diff --git a/packages/pangraph/src/tree/newick.rs b/packages/pangraph/src/tree/newick.rs
index 4f41f178..66c16535 100644
--- a/packages/pangraph/src/tree/newick.rs
+++ b/packages/pangraph/src/tree/newick.rs
@@ -310,6 +310,18 @@ mod tests {
assert_error!(parse_newick(input), expected);
}
+ #[rstest]
+ fn build_tree_from_newick_single_leaf() {
+ let dir = tempfile::tempdir().unwrap();
+ let path = dir.path().join("tree.nwk");
+ std::fs::write(&path, "A;").unwrap();
+
+ let tree = build_tree_from_newick(&path, vec![singleton("A", 0)]).unwrap();
+
+ assert!(tree.read().is_leaf());
+ assert_eq!(collect_leaf_names(&tree), vec!["A"]);
+ }
+
#[rstest]
fn build_tree_from_newick_attaches_graphs() {
let dir = tempfile::tempdir().unwrap();
diff --git a/packages/pangraph/src/utils/collections.rs b/packages/pangraph/src/utils/collections.rs
index add05a54..fe1f49f4 100644
--- a/packages/pangraph/src/utils/collections.rs
+++ b/packages/pangraph/src/utils/collections.rs
@@ -1,6 +1,6 @@
use crate::{make_error, make_internal_report};
use eyre::Report;
-use std::collections::HashSet;
+use std::collections::{BTreeMap, HashSet};
use std::hash::Hash;
pub fn concat_to_vec(x: &[T], y: &[T]) -> Vec {
@@ -32,3 +32,17 @@ pub fn has_duplicates>(iter: I) -> bool
let mut seen = HashSet::new();
iter.into_iter().any(|item| !seen.insert(item))
}
+
+/// Returns the values that occur more than once in the given iterator, in sorted order.
+/// Each duplicated value is reported once, no matter how many times it occurs.
+pub fn find_duplicates>(iter: I) -> Vec {
+ let mut counts: BTreeMap = BTreeMap::new();
+ for item in iter {
+ *counts.entry(item).or_insert(0) += 1;
+ }
+ counts
+ .into_iter()
+ .filter(|(_, count)| *count > 1)
+ .map(|(item, _)| item)
+ .collect()
+}
diff --git a/packages/pangraph/tests/itest_merge.rs b/packages/pangraph/tests/itest_merge.rs
new file mode 100644
index 00000000..b5d6e465
--- /dev/null
+++ b/packages/pangraph/tests/itest_merge.rs
@@ -0,0 +1,378 @@
+mod common;
+
+#[cfg(test)]
+mod tests {
+ use eyre::Report;
+ use itertools::Itertools;
+ 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::io::fasta::{FastaReader, FastaRecord};
+ use pangraph::io::json::{JsonPretty, json_write_file};
+ 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::reconstruct::reconstruct;
+ 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, BTreeSet};
+ use std::path::{Path, PathBuf};
+ use tempfile::{TempDir, tempdir};
+
+ /// Reads the first `n` records of a FASTA file.
+ fn read_records(path: &str, n: usize) -> Result, Report> {
+ let mut fastas = FastaReader::from_paths(&[PathBuf::from(path)])?.read_many()?;
+ fastas.truncate(n);
+ Ok(fastas)
+ }
+
+ /// Reads the first `n_total` records of a FASTA file and splits them in two groups, so that the
+ /// two resulting graphs are built from disjoint but homologous sets of genomes.
+ fn read_and_split(path: &str, n_left: usize, n_total: usize) -> Result<(Vec, Vec), Report> {
+ let mut left = read_records(path, n_total)?;
+ let right = left.split_off(n_left);
+ Ok((left, right))
+ }
+
+ /// Builds a graph out of the given records and writes it to a JSON file in `dir`.
+ fn build_graph_file(dir: &TempDir, name: &str, fastas: Vec) -> Result {
+ let args = PangraphBuildArgs {
+ circular: false,
+ ..PangraphBuildArgs::default()
+ };
+ let graph = build(fastas, &args, true)?;
+ let path = dir.path().join(name);
+ json_write_file(&path, &graph, JsonPretty(false))?;
+ Ok(path)
+ }
+
+ /// Reconstructs the genomes of a graph, keyed by genome name.
+ fn sequences(graph: &Pangraph) -> Result, Report> {
+ reconstruct(graph).map(|r| r.map(|r| (r.seq_name, r.seq))).collect()
+ }
+
+ fn read_graph(path: &Path) -> Result {
+ Pangraph::from_path(&Some(path))
+ }
+
+ fn merge_args(left: PathBuf, right: PathBuf, output: PathBuf) -> PangraphMergeArgs {
+ PangraphMergeArgs {
+ left_graph: left,
+ right_graph: right,
+ output_json: output,
+ merge_params: GraphMergeParams::default(),
+ verify: true,
+ }
+ }
+
+ /// Two graphs built independently from homologous genomes merge into a graph that reconstructs
+ /// every input genome exactly. This is the main end-to-end guarantee of the `merge` command.
+ #[rstest]
+ fn itest_merge_homologous_graphs() -> Result<(), Report> {
+ let dir = tempdir()?;
+ let (left_fastas, right_fastas) = read_and_split("../../data/ges-1.fa", 3, 6)?;
+ let left_names: BTreeSet = left_fastas.iter().map(|f| f.seq_name.clone()).collect();
+ let right_names: BTreeSet = right_fastas.iter().map(|f| f.seq_name.clone()).collect();
+ let expected_names = left_names.union(&right_names).cloned().collect_vec();
+
+ let left = build_graph_file(&dir, "left.json", left_fastas)?;
+ let right = build_graph_file(&dir, "right.json", right_fastas)?;
+ let output = dir.path().join("merged.json");
+
+ // `--verify` is on, so this already checks that every genome round-trips
+ merge_run(&merge_args(left.clone(), right.clone(), output.clone()))?;
+
+ let merged = read_graph(&output)?;
+
+ // all genomes of both inputs are present, exactly once
+ assert_eq!(merged.paths.len(), 6);
+ assert_eq!(
+ sequences(&merged)?.keys().cloned().sorted().collect_vec(),
+ expected_names
+ );
+
+ // and their sequences are unchanged
+ let mut expected = sequences(&read_graph(&left)?)?;
+ expected.extend(sequences(&read_graph(&right)?)?);
+ assert_eq!(sequences(&merged)?, expected);
+
+ // merging actually found homology across the two graphs: some blocks are now shared between
+ // genomes that came from different inputs. (Block counts alone say little: reweaving splits
+ // blocks as it merges them, so the total can go either way.)
+ let cross_graph_blocks = merged
+ .blocks
+ .values()
+ .filter(|block| {
+ let names = block
+ .alignment_keys()
+ .iter()
+ .filter_map(|nid| merged.paths[&merged.nodes[nid].path_id()].name.clone())
+ .collect_vec();
+ names.iter().any(|n| left_names.contains(n)) && names.iter().any(|n| right_names.contains(n))
+ })
+ .count();
+ assert!(
+ cross_graph_blocks > 0,
+ "no block is shared between the two input graphs"
+ );
+
+ Ok(())
+ }
+
+ /// Appending a single genome to an existing graph: the case the command is meant to serve.
+ #[rstest]
+ fn itest_merge_single_genome() -> Result<(), Report> {
+ let dir = tempdir()?;
+ let (left_fastas, right_fastas) = read_and_split("../../data/ges-1.fa", 3, 4)?;
+ assert_eq!(right_fastas.len(), 1);
+
+ let left = build_graph_file(&dir, "left.json", left_fastas)?;
+ let right = build_graph_file(&dir, "right.json", right_fastas)?;
+ let output = dir.path().join("merged.json");
+
+ merge_run(&merge_args(left, right, output.clone()))?;
+
+ let merged = read_graph(&output)?;
+ assert_eq!(merged.paths.len(), 4);
+
+ Ok(())
+ }
+
+ /// Appending to a graph that is itself the result of an earlier merge. Identifiers minted by the
+ /// first merge survive into its output whenever a block or node finds no homologue, so this is
+ /// the case that used to need a varying relabeling salt to keep the third graph's identifiers off
+ /// them. Deriving identifiers from genome names removes the problem at the source. The three
+ /// genomes here are mutually unrelated, so nothing aligns and every identifier survives;
+ /// homologous appends never exercised this.
+ #[rstest]
+ fn itest_merge_appends_to_an_already_merged_graph() -> Result<(), Report> {
+ let dir = tempdir()?;
+
+ let first = build_graph_file(&dir, "first.json", read_records("../../data/flu-h1.fa", 2)?)?;
+ let second = build_graph_file(&dir, "second.json", read_records("../../data/sc2.fa", 1)?)?;
+ let third = build_graph_file(&dir, "third.json", read_records("../../data/mpox.fa", 1)?)?;
+
+ let merged_once = dir.path().join("merged-once.json");
+ merge_run(&merge_args(first, second, merged_once.clone()))?;
+ assert_eq!(read_graph(&merged_once)?.paths.len(), 3);
+
+ let merged_twice = dir.path().join("merged-twice.json");
+ merge_run(&merge_args(merged_once, third, merged_twice.clone()))?;
+
+ let merged = read_graph(&merged_twice)?;
+ assert_eq!(merged.paths.len(), 4);
+
+ Ok(())
+ }
+
+ /// Maps each genome name to the ids of the blocks its path walks through, in order. Comparing the
+ /// id *sets* of two graphs would not do: `build` used to label singleton blocks `0, 1, 2, ...`
+ /// whatever the genome, so the sets matched even when every genome held a different block.
+ fn blocks_by_genome(graph: &Pangraph) -> BTreeMap> {
+ graph
+ .paths
+ .values()
+ .map(|path| {
+ let blocks = path.nodes().iter().map(|nid| graph.nodes[nid].block_id()).collect_vec();
+ (path.name().clone().unwrap_or_default(), blocks)
+ })
+ .collect()
+ }
+
+ /// Block and node identifiers are derived from genome names rather than from the order the input
+ /// sequences were read in, so the same genomes under the same guide tree produce the same graph
+ /// whichever order they arrive in. Only path ids, which deliberately record the input order, are
+ /// expected to differ.
+ #[rstest]
+ fn itest_build_ids_do_not_depend_on_input_order() -> Result<(), Report> {
+ let dir = tempdir()?;
+ let fastas = read_records("../../data/ges-1.fa", 4)?;
+ let names = fastas.iter().map(|f| f.seq_name.clone()).collect_vec();
+
+ // Pinned, so that the input order is the only thing that varies between the two builds.
+ let newick = dir.path().join("guide.nwk");
+ std::fs::write(
+ &newick,
+ format!("(({},{}),({},{}));", names[0], names[1], names[2], names[3]),
+ )?;
+
+ let build_in_order = |fastas: Vec| -> Result {
+ let args = PangraphBuildArgs {
+ circular: false,
+ guide_tree: Some(newick.clone()),
+ ..PangraphBuildArgs::default()
+ };
+ build(fastas, &args, true)
+ };
+
+ // Reversed *and* re-indexed, which is what reading the same genomes from a reordered file does.
+ let mut backwards = fastas.clone();
+ backwards.reverse();
+ for (index, record) in backwards.iter_mut().enumerate() {
+ record.index = index;
+ }
+
+ let forward = build_in_order(fastas)?;
+ let reversed = build_in_order(backwards)?;
+
+ assert_eq!(blocks_by_genome(&forward), blocks_by_genome(&reversed));
+ assert_eq!(
+ forward.node_ids().collect::>(),
+ reversed.node_ids().collect::>()
+ );
+ assert_eq!(sequences(&forward)?, sequences(&reversed)?);
+
+ // path ids still record the input order, so the genomes come back in the order they were read
+ assert_eq!(forward.path_names().flatten().collect_vec(), names);
+ assert_eq!(
+ reversed.path_names().flatten().collect_vec(),
+ names.iter().rev().collect_vec()
+ );
+
+ Ok(())
+ }
+
+ /// A graph written by pangraph 1.3 or earlier labels its first genome with block, node and path
+ /// id `0`, whatever that genome is called, so two of them collide even though their names differ.
+ /// `merge` has to say so: without the check, `graph_join` panics on the conflicting key.
+ #[rstest]
+ fn itest_merge_rejects_graphs_with_colliding_identifiers() -> Result<(), Report> {
+ let dir = tempdir()?;
+
+ // A singleton graph in the pre-1.4 identifier scheme, where ids came from the record index.
+ let legacy_graph = |name: &str, seq: &str| {
+ let (bid, nid, pid) = (BlockId(0), NodeId(0), PathId(0));
+ Pangraph {
+ blocks: BTreeMap::from([(bid, PangraphBlock::from_consensus(seq, bid, nid))]),
+ nodes: BTreeMap::from([(nid, PangraphNode::new(nid, bid, pid, Forward, (0, seq.len())))]),
+ paths: BTreeMap::from([(
+ pid,
+ PangraphPath::new(pid, [nid], seq.len(), false, Some(name.to_owned()), None),
+ )]),
+ }
+ };
+
+ let left = dir.path().join("left.json");
+ let right = dir.path().join("right.json");
+ json_write_file(&left, &legacy_graph("genome_a", "ACGTACGTAC"), JsonPretty(false))?;
+ json_write_file(&right, &legacy_graph("genome_b", "TTTTGGGGCC"), JsonPretty(false))?;
+
+ let result = merge_run(&merge_args(left, right, dir.path().join("out.json")));
+
+ let error = report_to_string(&result.unwrap_err());
+ assert!(
+ error.contains("share block or node identifiers"),
+ "unexpected error message: {error}"
+ );
+
+ Ok(())
+ }
+
+ /// Merging a graph with itself duplicates every genome name, and must be rejected.
+ #[rstest]
+ fn itest_merge_rejects_duplicate_genome_names() -> Result<(), Report> {
+ let dir = tempdir()?;
+ let fastas = read_records("../../data/ges-1.fa", 3)?;
+ let graph = build_graph_file(&dir, "graph.json", fastas)?;
+ let output = dir.path().join("merged.json");
+
+ let result = merge_run(&merge_args(graph.clone(), graph, output));
+
+ let error = report_to_string(&result.unwrap_err());
+ assert!(
+ error.contains("Duplicate genome names"),
+ "unexpected error message: {error}"
+ );
+
+ Ok(())
+ }
+
+ /// Unnamed paths used to slip past the duplicate-name guard entirely: it scanned
+ /// `path_names().flatten()`, which drops `None`, and the check that would have caught them ran
+ /// only under `--verify`. Merging such a graph with itself therefore succeeded and silently
+ /// emitted every genome twice. `verify: false` is the whole point of this test.
+ #[rstest]
+ fn itest_merge_rejects_unnamed_genomes_without_verify() -> Result<(), Report> {
+ let dir = tempdir()?;
+ let fastas = read_records("../../data/ges-1.fa", 3)?;
+
+ let mut graph = read_graph(&build_graph_file(&dir, "named.json", fastas)?)?;
+ for path in graph.paths.values_mut() {
+ path.name = None;
+ }
+ let anonymous = dir.path().join("anonymous.json");
+ json_write_file(&anonymous, &graph, JsonPretty(false))?;
+
+ let output = dir.path().join("merged.json");
+ let result = merge_run(&PangraphMergeArgs {
+ verify: false,
+ ..merge_args(anonymous.clone(), anonymous, output.clone())
+ });
+
+ let error = report_to_string(&result.unwrap_err());
+ assert!(
+ error.contains("no name or an empty name"),
+ "unexpected error message: {error}"
+ );
+ assert!(!output.exists(), "a rejected merge must not write an output graph");
+
+ 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]
+ fn itest_build_rejects_duplicate_sequence_names() -> Result<(), Report> {
+ let mut fastas = read_records("../../data/ges-1.fa", 2)?;
+ fastas[1].seq_name = fastas[0].seq_name.clone();
+
+ let result = build(fastas, &PangraphBuildArgs::default(), false);
+
+ let error = report_to_string(&result.unwrap_err());
+ assert!(
+ error.contains("Duplicate sequence names"),
+ "unexpected error message: {error}"
+ );
+
+ Ok(())
+ }
+
+ /// A FASTA header of the form `> id` leaves the record unnamed, with the identifier in the
+ /// description. Such a genome could not be addressed by name afterwards, so `build` rejects it.
+ #[rstest]
+ fn itest_build_rejects_empty_sequence_name() -> Result<(), Report> {
+ let mut fastas = read_records("../../data/ges-1.fa", 2)?;
+ fastas[1].desc = Some(fastas[1].seq_name.clone());
+ fastas[1].seq_name = String::new();
+
+ let result = build(fastas, &PangraphBuildArgs::default(), false);
+
+ let error = report_to_string(&result.unwrap_err());
+ assert!(error.contains("empty name"), "unexpected error message: {error}");
+
+ Ok(())
+ }
+}
diff --git a/packages/pangraph/tests/itest_reconstruct.rs b/packages/pangraph/tests/itest_reconstruct.rs
new file mode 100644
index 00000000..2be595f2
--- /dev/null
+++ b/packages/pangraph/tests/itest_reconstruct.rs
@@ -0,0 +1,277 @@
+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(NodeId(0), BlockId(0), PathId(0), Forward, (0, 8)),
+ },
+ paths: btreemap! {
+ PathId(0) => PangraphPath::new(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("no name or an empty 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(())
+ }
+
+ /// A graph read from a file was not necessarily written by pangraph, so nothing guarantees its
+ /// cross-references resolve. `reconstruct` used to index the node map directly, so a path naming
+ /// a node the graph does not contain aborted the process with a bare "no entry found for key"
+ /// panic, and in a release build, where `sanity_check` is compiled out, with no indication of
+ /// which file was at fault. Validation now happens where the file name is still known.
+ #[rstest]
+ fn itest_reconstruct_rejects_malformed_graph_naming_the_file() -> Result<(), Report> {
+ let dir = tempdir()?;
+
+ let mut graph = Pangraph {
+ blocks: btreemap! {
+ BlockId(0) => PangraphBlock::new(BlockId(0), "ACGTACGT", btreemap! { NodeId(0) => Edit::empty() }),
+ },
+ nodes: btreemap! {
+ NodeId(0) => PangraphNode::new(NodeId(0), BlockId(0), PathId(0), Forward, (0, 8)),
+ },
+ paths: btreemap! {
+ PathId(0) => PangraphPath::new(PathId(0), [NodeId(0)], 8, false, Some("a".to_owned()), None),
+ },
+ };
+ graph.paths.get_mut(&PathId(0)).unwrap().nodes = vec![NodeId(99)];
+ let graph = write_graph(&dir, "malformed.json", &graph)?;
+
+ let args = PangraphReconstructArgs {
+ input_graph: Some(graph),
+ output_fasta: PathBuf::from("-"),
+ verify: None,
+ };
+ let err = report_to_string(&reconstruct_run(&args).unwrap_err());
+
+ assert!(
+ err.contains("Node 99 from path 0 not found in graph"),
+ "unexpected error: {err}"
+ );
+ assert!(err.contains("malformed.json"), "unexpected error: {err}");
+
+ Ok(())
+ }
+}