diff --git a/fixtures/gfa/connected_mixed_sn.gfa b/fixtures/gfa/connected_mixed_sn.gfa new file mode 100644 index 00000000..e9d03a8d --- /dev/null +++ b/fixtures/gfa/connected_mixed_sn.gfa @@ -0,0 +1,6 @@ +H VN:Z:1.0 +S A AAAA SN:Z:chr1 +S B CCCC SN:Z:chr1 +S C TTTT SN:Z:chr2 +L A + B + * +L B + C + * diff --git a/fixtures/gfa/disjoint_graphs.gfa b/fixtures/gfa/disjoint_graphs.gfa new file mode 100644 index 00000000..f8620826 --- /dev/null +++ b/fixtures/gfa/disjoint_graphs.gfa @@ -0,0 +1,7 @@ +H VN:Z:1.0 +S 1 AAAA +S 2 CCCC +S 3 TTTT +S 4 GGGG +L 1 + 2 + * +L 3 + 4 + * diff --git a/fixtures/gfa/paths_no_ref.gfa b/fixtures/gfa/paths_no_ref.gfa new file mode 100644 index 00000000..35918a0f --- /dev/null +++ b/fixtures/gfa/paths_no_ref.gfa @@ -0,0 +1,8 @@ +H VN:Z:1.0 +S 1 AAAA +S 2 CCCC +S 3 TTTT +L 1 + 2 + * +L 2 + 3 + * +P short 1+,2+ * +P longer_path 1+,2+,3+ * diff --git a/fixtures/gfa/ref_tag_multi.gfa b/fixtures/gfa/ref_tag_multi.gfa new file mode 100644 index 00000000..32dee6b8 --- /dev/null +++ b/fixtures/gfa/ref_tag_multi.gfa @@ -0,0 +1,7 @@ +H VN:Z:1.0 +S 1 AAAA SN:Z:chr1 SO:i:0 SR:i:0 +S 2 CCCC SN:Z:chr1 SO:i:4 SR:i:0 +S 3 TTTT SN:Z:chr2 SO:i:0 SR:i:0 +S 4 GGGG SN:Z:chr2 SO:i:4 SR:i:0 +L 1 + 2 + * +L 3 + 4 + * diff --git a/fixtures/gfa/ref_tag_single.gfa b/fixtures/gfa/ref_tag_single.gfa new file mode 100644 index 00000000..3eba7162 --- /dev/null +++ b/fixtures/gfa/ref_tag_single.gfa @@ -0,0 +1,6 @@ +H VN:Z:1.0 +S 1 AAAA SN:Z:chr1 SO:i:0 SR:i:0 +S 2 CCCC SN:Z:chr1 SO:i:4 SR:i:0 +S 3 TTTT SN:Z:chr1 SO:i:8 SR:i:0 +L 1 + 2 + * +L 2 + 3 + * diff --git a/gen-python/src/python_api/repository/imports.rs b/gen-python/src/python_api/repository/imports.rs index 0d31dd73..d38c2b49 100644 --- a/gen-python/src/python_api/repository/imports.rs +++ b/gen-python/src/python_api/repository/imports.rs @@ -17,7 +17,7 @@ use gen_models::{ use pyo3::{exceptions::PyRuntimeError, prelude::*}; use super::{PyRepository, run_write}; -use crate::python_api::sequence_part::PySequencePart; +use crate::python_api::{block_group::PySequenceGraph, sequence_part::PySequencePart}; #[pymethods] impl PyRepository { @@ -80,12 +80,17 @@ impl PyRepository { filename: String, sample: Option, collection: Option, - ) -> PyResult { + ) -> PyResult> { let collection = collection.unwrap_or_else(|| self.get_default_collection()); let sample = sample.unwrap_or_else(|| Sample::DEFAULT_NAME.to_string()); run_write(&self.context, !self.in_transaction, |ctx| { import_gfa(ctx, &PathBuf::from(&filename), &collection, &sample) - .map(|_| format!("'{}' imported.", filename)) + .map(|(_, block_groups)| { + block_groups + .into_iter() + .map(|bg| self.into_py_block_group(bg)) + .collect() + }) .map_err(|e| match e { GFAImportError::OperationError(OperationError::NoChanges) => { PyRuntimeError::new_err(format!("'{}': already exists", filename)) diff --git a/gen-r/src/rust/src/lib.rs b/gen-r/src/rust/src/lib.rs index db36c0e9..72f38cd8 100644 --- a/gen-r/src/rust/src/lib.rs +++ b/gen-r/src/rust/src/lib.rs @@ -1073,7 +1073,7 @@ impl Repository { filename: String, sample: String, collection: Nullable, - ) -> std::result::Result { + ) -> std::result::Result { let collection_name = resolve_collection_name( self.context.operations().conn(), nullable_string_to_option(collection), @@ -1086,9 +1086,13 @@ impl Repository { &collection_name, &sample, ) { - Ok(_) => { + Ok((_, block_groups)) => { end_transactions(&self.context).map_err(Error::Other)?; - Ok("GFA imported.".to_string()) + let values = block_groups + .into_iter() + .map(|bg| r!(self.into_sequence_graph(bg))) + .collect::>(); + Ok(List::from_values(values)) } Err(r#gen::imports::gfa::GFAImportError::OperationError(OperationError::NoChanges)) => { rollback_transactions(&self.context); diff --git a/gen-r/tests/testthat/test-basic.R b/gen-r/tests/testthat/test-basic.R index 83f70663..58e9f561 100644 --- a/gen-r/tests/testthat/test-basic.R +++ b/gen-r/tests/testthat/test-basic.R @@ -76,10 +76,10 @@ test_that("GFA and GenBank import/export work", { gfa_out <- tempfile(fileext = ".gfa") gb_out <- tempfile(fileext = ".gb") - expect_match( - repo$import_gfa(fixture_path("simple.gfa"), sample = "sample-a"), - "imported", ignore.case = TRUE - ) + result <- repo$import_gfa(fixture_path("simple.gfa"), sample = "sample-a") + expect_true(is.list(result)) + expect_gt(length(result), 0) + expect_true(all(sapply(result, inherits, "SequenceGraph"))) expect_binding_result(try(repo$export_gfa(gfa_out, sample = "sample-a"), silent = TRUE)) expect_binding_result(try( diff --git a/src/commands/import/gfa.rs b/src/commands/import/gfa.rs index a0107548..1a5f4322 100644 --- a/src/commands/import/gfa.rs +++ b/src/commands/import/gfa.rs @@ -57,10 +57,15 @@ pub fn execute(cli_context: &CliContext, cmd: Command) -> Result<()> { )?; } match import_gfa(context, &PathBuf::from(cmd.path.clone()), name, sample_name) { - Ok(_) => { - println!("GFA imported."); + Ok((_, block_groups)) => { conn.execute("END TRANSACTION;", []).unwrap(); operation_conn.execute("END TRANSACTION;", []).unwrap(); + let names: Vec<&str> = block_groups.iter().map(|bg| bg.name.as_str()).collect(); + println!( + "Imported {} sequence graph(s): {}", + names.len(), + names.join(", ") + ); Ok(()) } Err(GFAImportError::OperationError(OperationError::NoChanges)) => { diff --git a/src/diffs/gfa.rs b/src/diffs/gfa.rs index 9fcdbd63..0e7cbd90 100644 --- a/src/diffs/gfa.rs +++ b/src/diffs/gfa.rs @@ -911,11 +911,15 @@ mod tests { Sample::DEFAULT_NAME, ); - let new_block_group = Collection::get_block_groups(conn, "test collection 3") - .pop() - .unwrap(); - let all_sequences = - BlockGroup::get_all_sequences(conn, &new_block_group.id, false).unwrap(); + let all_sequences: HashSet = + Collection::get_block_groups(conn, "test collection 3") + .iter() + .flat_map(|bg| { + BlockGroup::get_all_sequences(conn, &bg.id, false) + .unwrap() + .into_iter() + }) + .collect(); assert_eq!( all_sequences, @@ -1081,11 +1085,15 @@ mod tests { Sample::DEFAULT_NAME, ); - let new_block_group = Collection::get_block_groups(conn, "test collection 3") - .pop() - .unwrap(); - let all_sequences = - BlockGroup::get_all_sequences(conn, &new_block_group.id, false).unwrap(); + let all_sequences: HashSet = + Collection::get_block_groups(conn, "test collection 3") + .iter() + .flat_map(|bg| { + BlockGroup::get_all_sequences(conn, &bg.id, false) + .unwrap() + .into_iter() + }) + .collect(); assert_eq!( all_sequences, diff --git a/src/exports/gfa.rs b/src/exports/gfa.rs index acbaae75..51f3960e 100644 --- a/src/exports/gfa.rs +++ b/src/exports/gfa.rs @@ -609,7 +609,8 @@ mod tests { let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "m123", None); let all_sequences = BlockGroup::get_all_sequences(conn, &block_group_id, false).unwrap(); let temp_dir = tempdir().expect("Couldn't get handle to temp directory"); @@ -652,7 +653,12 @@ mod tests { let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = BlockGroup::get_id( + &collection_name, + Sample::DEFAULT_NAME, + "BBa_J23119#0#BBa_J23119", + None, + ); let all_sequences = BlockGroup::get_all_sequences(conn, &block_group_id, false).unwrap(); let temp_dir = tempdir().expect("Couldn't get handle to temp directory"); @@ -695,7 +701,8 @@ mod tests { let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "123", None); let all_sequences = BlockGroup::get_all_sequences(conn, &block_group_id, false).unwrap(); let temp_dir = tempdir().expect("Couldn't get handle to temp directory"); diff --git a/src/imports/gfa.rs b/src/imports/gfa.rs index cf1bd0e3..07d87358 100644 --- a/src/imports/gfa.rs +++ b/src/imports/gfa.rs @@ -1,10 +1,14 @@ -use std::{collections::HashMap, path::Path as FilePath}; +use std::{ + collections::{HashMap, HashSet}, + fs::File, + io::{BufRead, BufReader}, + path::Path as FilePath, +}; use gen_core::{ - HashId, NO_CHROMOSOME_INDEX, PATH_END_NODE_ID, PATH_START_NODE_ID, Strand, is_end_node, - is_start_node, + HashId, NO_CHROMOSOME_INDEX, PATH_END_NODE_ID, PATH_START_NODE_ID, Strand, calculate_hash, + is_end_node, is_start_node, }; -use gen_graph::{GraphEdge, GraphNode}; use gen_models::{ block_group::{BlockGroup, NewBlockGroup}, block_group_edge::{BlockGroupEdge, BlockGroupEdgeData}, @@ -26,7 +30,10 @@ use gen_models::{ }; use indexmap::IndexSet; use itertools::Itertools; -use petgraph::{algo::kosaraju_scc, prelude::UnGraphMap, visit::Dfs}; +use petgraph::{ + prelude::UnGraphMap, + visit::{Bfs, Dfs}, +}; use thiserror::Error; use crate::{ @@ -58,7 +65,7 @@ pub fn import_gfa( gfa_path: &FilePath, collection_name: &str, sample_name: &str, -) -> Result { +) -> Result<(Operation, Vec), GFAImportError> { let conn = context.graph().conn(); let progress_bar = get_handler(); let mut session = start_operation(conn); @@ -79,21 +86,124 @@ pub fn import_gfa( return Err(GFAImportError::SampleError(e)); } } - let block_group = BlockGroup::create( - conn, - NewBlockGroup { - collection_name, - sample_name, - name: "", - ..Default::default() - }, - )?; + let bar = progress_bar.add(get_time_elapsed_bar()); bar.set_message("Parsing GFA"); let gfa: Gfa = Gfa::parse_gfa_file(gfa_path.to_str().unwrap()); + let sn_tags = read_sn_tags(gfa_path); + bar.finish(); + + // Build segment connectivity graph from GFA links + let mut seg_graph: UnGraphMap<&str, ()> = UnGraphMap::new(); + for seg in &gfa.segments { + seg_graph.add_node(seg.id.as_str()); + } + for link in &gfa.links { + seg_graph.add_edge(link.from.as_str(), link.to.as_str(), ()); + } + + // Find connected components via BFS + let mut visited: HashSet<&str> = HashSet::new(); + let mut components: Vec> = vec![]; + for seg in &gfa.segments { + if visited.contains(seg.id.as_str()) { + continue; + } + let mut component = vec![]; + let mut bfs = Bfs::new(&seg_graph, seg.id.as_str()); + while let Some(node) = bfs.next(&seg_graph) { + visited.insert(node); + component.push(node); + } + components.push(component); + } + + // Assign one name per connected component: SN:Z: tags > longest path/walk > filename stem + let filename_stem = gfa_path + .file_stem() + .and_then(|s| s.to_str()) + .unwrap_or("") + .to_string(); + let raw_names: Vec = components + .iter() + .map(|segs| { + let seg_set: HashSet<&str> = segs.iter().copied().collect(); + if !sn_tags.is_empty() { + // majority SN tag value among segments; tie-break alphabetically first + let mut counts: HashMap<&str, usize> = HashMap::new(); + for seg_id in segs { + if let Some(tag) = sn_tags.get(*seg_id) { + *counts.entry(tag.as_str()).or_default() += 1; + } + } + counts + .into_iter() + .max_by(|(a, ca), (b, cb)| ca.cmp(cb).then(b.cmp(a))) + .map(|(tag, _)| tag.to_string()) + .unwrap_or_default() + } else if !gfa.paths.is_empty() || !gfa.walk.is_empty() { + gfa.paths + .iter() + .filter(|p| p.segments.iter().any(|s| seg_set.contains(s.as_str()))) + .map(|p| (p.name.as_str(), p.segments.len())) + .chain( + gfa.walk + .iter() + .filter(|w| w.segments.iter().any(|s| seg_set.contains(s.as_str()))) + .map(|w| (w.sample_id.as_str(), w.segments.len())), + ) + .max_by_key(|&(_, len)| len) + .map(|(name, _)| name.to_string()) + .unwrap_or_default() + } else { + filename_stem.clone() + } + }) + .collect(); + + // Disambiguate duplicate base names with _N suffix + let mut name_counts: HashMap = HashMap::new(); + for n in &raw_names { + *name_counts.entry(n.clone()).or_default() += 1; + } + let mut name_seen: HashMap = HashMap::new(); + let component_names: Vec = raw_names + .into_iter() + .map(|n| { + if name_counts[&n] > 1 { + let idx = name_seen.entry(n.clone()).or_default(); + *idx += 1; + format!("{}_{}", n, idx) + } else { + n + } + }) + .collect(); + + let bg_name_by_segment: HashMap = components + .iter() + .zip(component_names.iter()) + .flat_map(|(segs, name)| segs.iter().map(move |seg| (seg.to_string(), name.clone()))) + .collect(); + + // Create one block group per unique name + let unique_bg_names: IndexSet = bg_name_by_segment.values().cloned().collect(); + let mut block_groups_by_name: HashMap = HashMap::new(); + for bg_name in &unique_bg_names { + let bg = BlockGroup::create( + conn, + NewBlockGroup { + collection_name, + sample_name, + name: bg_name, + ..Default::default() + }, + )?; + block_groups_by_name.insert(bg_name.clone(), bg); + } + let mut sequences_by_segment_id: HashMap<&String, Sequence> = HashMap::new(); let mut node_ids_by_segment_id: HashMap<&String, HashId> = HashMap::new(); - bar.finish(); let bar = progress_bar.add(get_progress_bar(gfa.segments.len() as u64)); bar.set_message("Parsing Segments"); @@ -104,14 +214,29 @@ pub fn import_gfa( .sequence(input_sequence) .save(conn)?; sequences_by_segment_id.insert(&segment.id, sequence.clone()); - // TODO: Node hash is always new, it's sorted by insert time via being a v7 uuid but maybe want to - // define the hash itself for idempotency? - let node_id = Node::create(conn, &sequence.hash, &HashId::uuid7())?; + let node_hash = HashId(calculate_hash(&format!( + "{collection_name}.{seg_id}:{seq_hash}", + seg_id = segment.id, + seq_hash = sequence.hash + ))); + let node_id = Node::create(conn, &sequence.hash, &node_hash)?; node_ids_by_segment_id.insert(&segment.id, node_id); bar.inc(1); } bar.finish(); + // Map node_id → block_group_id for routing edges to the correct block group + let mut bg_id_by_node_id: HashMap = HashMap::new(); + for (seg_id, node_id) in &node_ids_by_segment_id { + if let Some(bg_id) = bg_name_by_segment + .get(*seg_id) + .and_then(|name| block_groups_by_name.get(name)) + .map(|bg| bg.id) + { + bg_id_by_node_id.insert(*node_id, bg_id); + } + } + let mut edges = IndexSet::new(); let bar = progress_bar.add(get_progress_bar(gfa.links.len() as u64)); let mut source_refs_in_links = IndexSet::new(); @@ -233,7 +358,8 @@ pub fn import_gfa( let edge_ids = Edge::bulk_create(conn, &edges.into_iter().collect::>()); let saved_edges = Edge::query_by_ids(conn, &edge_ids); - let mut edge_ids_by_data = HashMap::new(); + let mut edge_ids_by_data: HashMap = HashMap::new(); + let mut edge_data_by_id: HashMap = HashMap::new(); for edge in saved_edges { let key = edge_data_from_fields( edge.source_node_id, @@ -243,6 +369,7 @@ pub fn import_gfa( edge.target_strand, ); edge_ids_by_data.insert(key, edge.id); + edge_data_by_id.insert(edge.id, key); } let mut created_blockgroup_edges: IndexSet = IndexSet::new(); @@ -281,19 +408,26 @@ pub fn import_gfa( path_edge_ids.push(edge_id); created_blockgroup_edges.extend(path_edge_ids.iter()); + let first_seg = &input_path.segments[0]; + let first_node_id = *node_ids_by_segment_id.get(first_seg).unwrap(); + let bg_id = bg_id_by_node_id + .get(&first_node_id) + .copied() + .unwrap_or_else(|| block_groups_by_name.values().next().unwrap().id); + BlockGroupEdge::bulk_create( conn, &path_edge_ids .iter() .map(|id| BlockGroupEdgeData { - block_group_id: block_group.id, + block_group_id: bg_id, edge_id: *id, chromosome_index: NO_CHROMOSOME_INDEX, phased: 0, }) .collect::>(), ); - Path::create(conn, path_name, &block_group.id, &path_edge_ids)?; + Path::create(conn, path_name, &bg_id, &path_edge_ids)?; } for input_walk in &gfa.walk { @@ -330,129 +464,138 @@ pub fn import_gfa( path_edge_ids.push(edge_id); created_blockgroup_edges.extend(path_edge_ids.iter()); + let first_seg = &input_walk.segments[0]; + let first_node_id = *node_ids_by_segment_id.get(first_seg).unwrap(); + let bg_id = bg_id_by_node_id + .get(&first_node_id) + .copied() + .unwrap_or_else(|| block_groups_by_name.values().next().unwrap().id); + BlockGroupEdge::bulk_create( conn, &path_edge_ids .iter() .map(|id| BlockGroupEdgeData { - block_group_id: block_group.id, + block_group_id: bg_id, edge_id: *id, chromosome_index: NO_CHROMOSOME_INDEX, phased: 0, }) .collect::>(), ); - Path::create(conn, path_name, &block_group.id, &path_edge_ids)?; + Path::create(conn, path_name, &bg_id, &path_edge_ids)?; } // make any block group edges not in paths or walks - BlockGroupEdge::bulk_create( - conn, - &edge_ids - .iter() - .filter_map(|id| { - if !created_blockgroup_edges.contains(id) { - Some(BlockGroupEdgeData { - block_group_id: block_group.id, - edge_id: *id, - chromosome_index: NO_CHROMOSOME_INDEX, - phased: 0, - }) - } else { - None - } + let leftover_bge: Vec = edge_ids + .iter() + .filter_map(|id| { + if created_blockgroup_edges.contains(id) { + return None; + } + let edge_data = edge_data_by_id.get(id)?; + let node_key = if is_start_node(edge_data.source_node_id) { + edge_data.target_node_id + } else { + edge_data.source_node_id + }; + let bg_id = bg_id_by_node_id + .get(&node_key) + .copied() + .unwrap_or_else(|| block_groups_by_name.values().next().unwrap().id); + Some(BlockGroupEdgeData { + block_group_id: bg_id, + edge_id: *id, + chromosome_index: NO_CHROMOSOME_INDEX, + phased: 0, }) - .collect::>(), - ); + }) + .collect(); + BlockGroupEdge::bulk_create(conn, &leftover_bge); - // check the graph for cycles and make start/end nodes if so + // check each block group graph for cycles and wire up start/end nodes let bar = progress_bar.add(get_progress_bar(None)); bar.set_message("Breaking cycles"); let message_bar = progress_bar.add(get_message_bar()); - let graph = BlockGroup::get_graph(conn, &block_group.id)?; - let mut undirected_graph: UnGraphMap = UnGraphMap::new(); - for node in graph.nodes() { - undirected_graph.add_node(node); - } - for (src, dst, weights) in graph.all_edges() { - undirected_graph.add_edge(src, dst, weights[0]); - } - let connected_components = kosaraju_scc(&undirected_graph); - let mut new_edges = vec![]; - for subgraph in connected_components.iter() { - if subgraph.len() >= 3 { - let mut has_start = false; - let mut has_end = false; - for node in subgraph.iter() { - if !has_start && is_start_node(node.node_id) { - has_start = true; - } else if !has_end && is_end_node(node.node_id) { - has_end = true; - }; - if has_start && has_end { - break; - } + let mut all_new_cycle_edges: Vec<(EdgeData, HashId)> = vec![]; + for bg in block_groups_by_name.values() { + let graph = BlockGroup::get_graph(conn, &bg.id)?; + if graph.node_count() < 3 { + continue; + } + let has_start = graph + .all_edges() + .any(|(src, _, _)| is_start_node(src.node_id)); + let has_end = graph + .all_edges() + .any(|(_, dst, _)| is_end_node(dst.node_id)); + if !has_start && !has_end { + // Cycle with no explicit entry/exit — wire start/end via DFS ordering + let first = graph.nodes().min_by_key(|n| n.node_id).unwrap(); + let mut order = vec![]; + let mut dfs = Dfs::new(&graph, first); + while let Some(nx) = dfs.next(&graph) { + order.push(nx); } - // For graphs with just one enter/exit point, we log a message - if !has_start && !has_end { - // from the subgraph, we want to find a deterministic sort of ordered elements. - // Kosaraju returns nodes in arbitrary order. We use DFS and then rotate the vector - // so the first node_id starts the list for consistency. If a node in the DFS is in - // a known start node for a path, we use that one. - let mut order = vec![]; - let mut dfs = Dfs::new(&graph, subgraph[0]); - while let Some(nx) = dfs.next(&graph) { - order.push(nx); - } - let min_index = order.iter().enumerate().min_set_by_key(|(_, k)| k.node_id)[0].0; - order.rotate_left(min_index); - bar.inc(1); - new_edges.push(edge_data_from_fields( + let min_index = order.iter().enumerate().min_set_by_key(|(_, k)| k.node_id)[0].0; + order.rotate_left(min_index); + bar.inc(1); + let last_node = *order.last().unwrap(); + all_new_cycle_edges.push(( + edge_data_from_fields( PATH_START_NODE_ID, 0, Strand::Forward, order[0].node_id, Strand::Forward, - )); - let last_node = order.last().unwrap(); - new_edges.push(edge_data_from_fields( + ), + bg.id, + )); + all_new_cycle_edges.push(( + edge_data_from_fields( last_node.node_id, last_node.sequence_end, Strand::Forward, PATH_END_NODE_ID, Strand::Forward, - )); - new_edges.push(edge_data_from_fields( + ), + bg.id, + )); + all_new_cycle_edges.push(( + edge_data_from_fields( PATH_END_NODE_ID, 0, Strand::Forward, PATH_START_NODE_ID, Strand::Forward, - )); - } else if has_start && has_end { - // there's a cycle, but has a start/end already. At some point we should track this - // so we know ahead of time where the cycles are - } else { - message_bar.set_message("Path encountered with cycle after start/end node, no cycle breaking will apply."); - } + ), + bg.id, + )); + } else if !has_start || !has_end { + message_bar.set_message( + "Path encountered with cycle after start/end node, no cycle breaking will apply.", + ); } } message_bar.finish(); - let new_edge_ids = Edge::bulk_create(conn, &new_edges.into_iter().collect::>()); - BlockGroupEdge::bulk_create( - conn, - &new_edge_ids - .iter() - .map(|id| BlockGroupEdgeData { - block_group_id: block_group.id, - edge_id: *id, - chromosome_index: NO_CHROMOSOME_INDEX, - phased: 0, - }) - .collect::>(), - ); + let new_edge_data: Vec = all_new_cycle_edges.iter().map(|(e, _)| *e).collect(); + let new_edge_ids = Edge::bulk_create(conn, &new_edge_data); + let cycle_bge: Vec = new_edge_ids + .iter() + .zip(all_new_cycle_edges.iter()) + .map(|(id, (_, bg_id))| BlockGroupEdgeData { + block_group_id: *bg_id, + edge_id: *id, + chromosome_index: NO_CHROMOSOME_INDEX, + phased: 0, + }) + .collect(); + BlockGroupEdge::bulk_create(conn, &cycle_bge); bar.finish(); + let mut block_groups: Vec = block_groups_by_name.into_values().collect(); + block_groups.sort_by(|a, b| a.name.cmp(&b.name)); + let op = end_operation( context, &mut session, @@ -468,7 +611,32 @@ pub fn import_gfa( ) .map_err(GFAImportError::OperationError); gen_bar.finish(); - op + op.map(|op| (op, block_groups)) +} + +fn read_sn_tags(gfa_path: &FilePath) -> HashMap { + let Ok(file) = File::open(gfa_path) else { + return HashMap::new(); + }; + let mut sn_tags = HashMap::new(); + for line in BufReader::new(file).lines().map_while(Result::ok) { + if !line.starts_with('S') { + continue; + } + let mut fields = line.splitn(10, '\t'); + fields.next(); // 'S' + let Some(seg_id) = fields.next() else { + continue; + }; + fields.next(); // sequence + for opt in fields { + if let Some(sn_value) = opt.strip_prefix("SN:Z:") { + sn_tags.insert(seg_id.to_string(), sn_value.to_string()); + break; + } + } + } + sn_tags } fn edge_data_from_fields( @@ -509,7 +677,8 @@ mod tests { track_database(conn, context.operations().conn()).unwrap(); let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "m123", None); let path = Path::query( conn, "select * from paths where block_group_id = ?1 AND name = ?2", @@ -540,6 +709,39 @@ mod tests { ); } + #[test] + fn test_double_import_gfa_is_idempotent() { + let mut gfa_path = PathBuf::from(env!("CARGO_MANIFEST_DIR")); + gfa_path.push("fixtures/simple.gfa"); + let collection_name = "test".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + + track_database(conn, context.operations().conn()).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let node_count_after_first = + Node::query(conn, "select * from nodes", rusqlite::params!()).len(); + + let second_result = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + assert!( + matches!( + second_result, + Err(GFAImportError::OperationError( + gen_models::errors::OperationError::NoChanges + )) + ), + "expected NoChanges on duplicate import, got {second_result:?}" + ); + + let node_count_after_second = + Node::query(conn, "select * from nodes", rusqlite::params!()).len(); + assert_eq!( + node_count_after_first, node_count_after_second, + "duplicate import must not create new nodes" + ); + } + #[test] fn test_import_no_path_gfa() { let mut gfa_path = PathBuf::from(env!("CARGO_MANIFEST_DIR")); @@ -551,7 +753,8 @@ mod tests { track_database(conn, context.operations().conn()).unwrap(); let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "no_path", None); let all_sequences = BlockGroup::get_all_sequences(conn, &block_group_id, false).unwrap(); assert_eq!( all_sequences, @@ -573,7 +776,8 @@ mod tests { track_database(conn, context.operations().conn()).unwrap(); let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "291344", None); let path = Path::query( conn, "select * from paths where block_group_id = ?1 AND name = ?2", @@ -599,7 +803,8 @@ mod tests { track_database(conn, context.operations().conn()).unwrap(); let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "123", None); let path = Path::query( conn, "select * from paths where block_group_id = ?1 AND name = ?2", @@ -629,7 +834,12 @@ mod tests { let paths = Path::query_for_collection(conn, &collection_name); assert_eq!(paths.len(), 20); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = BlockGroup::get_id( + &collection_name, + Sample::DEFAULT_NAME, + "BBa_J23119#0#BBa_J23119", + None, + ); let path = Path::query( conn, "select * from paths where block_group_id = ?1 AND name = ?2", @@ -732,7 +942,8 @@ mod tests { track_database(conn, op_conn).unwrap(); let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "123", None); let path = Path::query( conn, "select * from paths where block_group_id = ?1 AND name = ?2", @@ -762,7 +973,12 @@ mod tests { track_database(conn, op_conn).unwrap(); let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = BlockGroup::get_id( + &collection_name, + Sample::DEFAULT_NAME, + "cycle_no_path", + None, + ); let all_sequences = BlockGroup::get_all_sequences(conn, &block_group_id, false).unwrap(); assert_eq!( @@ -771,6 +987,41 @@ mod tests { ); } + #[test] + fn test_import_disjoint_graphs() { + // Two separate linear chains with no cross-links; each connected component becomes its own block group. + let gfa_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/gfa/disjoint_graphs.gfa"); + let collection_name = "disjoint".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + + track_database(conn, op_conn).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let bg1_id = BlockGroup::get_id( + &collection_name, + Sample::DEFAULT_NAME, + "disjoint_graphs_1", + None, + ); + let bg2_id = BlockGroup::get_id( + &collection_name, + Sample::DEFAULT_NAME, + "disjoint_graphs_2", + None, + ); + assert_eq!( + BlockGroup::get_all_sequences(conn, &bg1_id, false).unwrap(), + HashSet::from_iter(vec!["AAAACCCC".to_string()]) + ); + assert_eq!( + BlockGroup::get_all_sequences(conn, &bg2_id, false).unwrap(), + HashSet::from_iter(vec!["TTTTGGGG".to_string()]) + ); + } + #[test] fn test_breaks_cycle_using_path_node() { // here the fixture has a path indicting the cycle starts in the middle of where it would @@ -785,7 +1036,8 @@ mod tests { track_database(conn, op_conn).unwrap(); let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); - let block_group_id = BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "", None); + let block_group_id = + BlockGroup::get_id(&collection_name, Sample::DEFAULT_NAME, "m123", None); let all_sequences = BlockGroup::get_all_sequences(conn, &block_group_id, false).unwrap(); assert_eq!( @@ -793,4 +1045,148 @@ mod tests { HashSet::from_iter(vec!["TTTGGGACTCTAAAACCC".to_string()]) ); } + + // --- Chris' proposed naming cases --- + + #[test] + fn test_import_gfa_ref_tag_names_block_group() { + // SN:Z: tag present on all segments; block group should be named by that value. + let gfa_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/gfa/ref_tag_single.gfa"); + let collection_name = "ref_tag_single".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + + track_database(conn, op_conn).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let block_groups = BlockGroup::query( + conn, + "select * from block_groups where collection_name = ?1", + params![collection_name], + ); + assert_eq!(block_groups.len(), 1); + assert_eq!(block_groups[0].name, "chr1"); + } + + #[test] + fn test_import_gfa_multi_ref_tags_create_multiple_block_groups() { + // Different SN:Z: values on disjoint subgraphs; one block group per unique SN value. + let gfa_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/gfa/ref_tag_multi.gfa"); + let collection_name = "ref_tag_multi".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + + track_database(conn, op_conn).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let block_groups = BlockGroup::query( + conn, + "select * from block_groups where collection_name = ?1", + params![collection_name], + ); + let names: HashSet = block_groups.iter().map(|bg| bg.name.clone()).collect(); + assert_eq!( + names, + HashSet::from_iter(vec!["chr1".to_string(), "chr2".to_string()]) + ); + } + + #[test] + fn test_import_gfa_longest_path_names_block_group() { + // No SN:Z: tags; block group named by the path with the most segments. + let gfa_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/gfa/paths_no_ref.gfa"); + let collection_name = "paths_no_ref".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + + track_database(conn, op_conn).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let block_groups = BlockGroup::query( + conn, + "select * from block_groups where collection_name = ?1", + params![collection_name], + ); + assert_eq!(block_groups.len(), 1); + assert_eq!(block_groups[0].name, "longer_path"); + } + + #[test] + fn test_import_gfa_single_graph_no_ref_uses_filename() { + // No SN:Z: tags, no paths, single connected graph; block group named by GFA filename stem. + let gfa_path = PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/no_path.gfa"); + let collection_name = "no_path_filename".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + + track_database(conn, op_conn).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let block_groups = BlockGroup::query( + conn, + "select * from block_groups where collection_name = ?1", + params![collection_name], + ); + assert_eq!(block_groups.len(), 1); + assert_eq!(block_groups[0].name, "no_path"); + } + + #[test] + fn test_import_gfa_multiple_graphs_no_ref_uses_filename() { + // No SN:Z: tags, no paths, multiple disjoint subgraphs; one block group per component, names disambiguated with _N. + let gfa_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/gfa/disjoint_graphs.gfa"); + let collection_name = "disjoint_filename".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + + track_database(conn, op_conn).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let block_groups = BlockGroup::query( + conn, + "select * from block_groups where collection_name = ?1", + params![collection_name], + ); + let names: HashSet = block_groups.iter().map(|bg| bg.name.clone()).collect(); + assert_eq!(block_groups.len(), 2); + assert_eq!( + names, + HashSet::from_iter(vec![ + "disjoint_graphs_1".to_string(), + "disjoint_graphs_2".to_string() + ]) + ); + } + + #[test] + fn test_import_gfa_connected_mixed_sn_tags() { + // All segments connected (A→B→C) but with mixed SN:Z: tags (chr1, chr1, chr2). + // Connectivity takes priority: one block group named by majority SN tag (chr1, 2 vs 1). + let gfa_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/gfa/connected_mixed_sn.gfa"); + let collection_name = "mixed_sn".to_string(); + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + + track_database(conn, op_conn).unwrap(); + let _ = import_gfa(&context, &gfa_path, &collection_name, Sample::DEFAULT_NAME); + + let block_groups = BlockGroup::query( + conn, + "select * from block_groups where collection_name = ?1", + params![collection_name], + ); + assert_eq!(block_groups.len(), 1); + assert_eq!(block_groups[0].name, "chr1"); + } }