diff --git a/Cargo.lock b/Cargo.lock index cbe113e1..214e6481 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -775,6 +775,18 @@ version = "0.4.32" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "cc14f565cf027a105f7a44ccf9e5b424348421a1d8952a8fc9d499d313107789" +[[package]] +name = "console" +version = "0.15.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "054ccb5b10f9f2cbf51eb355ca1d05c2d279ce1804688d0db74b4733a5aeafd8" +dependencies = [ + "encode_unicode", + "libc", + "once_cell", + "windows-sys 0.59.0", +] + [[package]] name = "console" version = "0.16.3" @@ -1669,6 +1681,7 @@ dependencies = [ "gen-diff", "gen-graph", "gen-models", + "gen-parsers", "gen-tui", "getrandom 0.3.4", "html-escape", @@ -1805,6 +1818,17 @@ dependencies = [ "url", ] +[[package]] +name = "gen-parsers" +version = "0.1.31" +dependencies = [ + "gen-core", + "gen-graph", + "petgraph 0.6.5", + "similar-asserts", + "thiserror 2.0.18", +] + [[package]] name = "gen-sugiyama" version = "0.1.31" @@ -2322,7 +2346,7 @@ version = "0.18.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "25470f23803092da7d239834776d653104d551bc4d7eacaf31e6837854b8e9eb" dependencies = [ - "console", + "console 0.16.3", "portable-atomic", "unicode-width", "unit-prefix", @@ -2354,7 +2378,7 @@ version = "1.47.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7b4a6248eb93a4401ed2f37dfe8ea592d3cf05b7cf4f8efa867b6895af7e094e" dependencies = [ - "console", + "console 0.16.3", "once_cell", "similar", "tempfile", @@ -5135,6 +5159,20 @@ name = "similar" version = "2.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "bbbb5d9659141646ae647b42fe094daf6c6192d1620870b449d9557f748b2daa" +dependencies = [ + "bstr", + "unicode-segmentation", +] + +[[package]] +name = "similar-asserts" +version = "1.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5b441962c817e33508847a22bd82f03a30cff43642dc2fae8b050566121eb9a" +dependencies = [ + "console 0.15.11", + "similar", +] [[package]] name = "simple_asn1" @@ -6308,6 +6346,15 @@ dependencies = [ "windows-targets 0.52.6", ] +[[package]] +name = "windows-sys" +version = "0.59.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b" +dependencies = [ + "windows-targets 0.52.6", +] + [[package]] name = "windows-sys" version = "0.60.2" diff --git a/Cargo.toml b/Cargo.toml index 4b29c137..5c0a902c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,8 +21,8 @@ name = "gen" path = "src/main.rs" [workspace] -members = [".", "gen-core", "gen-models", "gen-graph", "gen-diff", "gen-tui", "gen-sugiyama", "gen-annotations", "gen-capnp-schemas"] -default-members = [".", "gen-core", "gen-models", "gen-graph", "gen-diff", "gen-capnp-schemas", "gen-annotations"] +members = [".", "gen-core", "gen-models", "gen-graph", "gen-diff", "gen-tui", "gen-sugiyama", "gen-annotations", "gen-capnp-schemas", "gen-parsers"] +default-members = [".", "gen-core", "gen-models", "gen-graph", "gen-diff", "gen-capnp-schemas", "gen-annotations", "gen-parsers"] exclude = ["gen-python", "gen-r/src/rust"] [features] @@ -37,6 +37,7 @@ gen-core = { path = "gen-core", version = "0.1.31" } gen-diff = { path = "gen-diff", version = "0.1.31" } gen-models = { path = "gen-models", version = "0.1.31" } gen-graph = { path = "gen-graph", version = "0.1.31" } +gen-parsers = { path = "gen-parsers", version = "0.1.31" } gen-tui = { path = "gen-tui", version = "0.1.31" } gen-annotations = { path = "gen-annotations", version = "0.1.31" } gen-capnp-schemas = { path = "gen-capnp-schemas", version = "0.1.31" } diff --git a/fixtures/alignments/plain_names.aln b/fixtures/alignments/plain_names.aln new file mode 100644 index 00000000..5ddf0fbb --- /dev/null +++ b/fixtures/alignments/plain_names.aln @@ -0,0 +1,5 @@ +CLUSTAL W (2.1) multiple sequence alignment + +Base ACGT +Child A-GT + * ** diff --git a/fixtures/alignments/simple.aln b/fixtures/alignments/simple.aln new file mode 100644 index 00000000..f4970928 --- /dev/null +++ b/fixtures/alignments/simple.aln @@ -0,0 +1,6 @@ +CLUSTAL W (2.1) multiple sequence alignment + +SeqA.chr1 ACGT +SeqB.chr1 A-GT +SeqC.chr1 ACG- + * ** diff --git a/gen-models/src/block_group.rs b/gen-models/src/block_group.rs index 0d8944aa..16d6d806 100644 --- a/gen-models/src/block_group.rs +++ b/gen-models/src/block_group.rs @@ -525,7 +525,7 @@ impl BlockGroup { ) } - fn copy_contents_from( + pub fn copy_contents_from( &self, conn: &GraphConnection, source_block_group: &BlockGroup, diff --git a/gen-parsers/Cargo.toml b/gen-parsers/Cargo.toml new file mode 100644 index 00000000..e98ba252 --- /dev/null +++ b/gen-parsers/Cargo.toml @@ -0,0 +1,17 @@ +[package] +name = "gen-parsers" +description = "Bioinformatic file format parsers for Gen." +version = "0.1.31" +edition = "2024" +repository = "https://github.com/genhub-bio/gen" +homepage = "https://genhub.bio" +license = "Apache-2.0" + +[dependencies] +gen-core = { path = "../gen-core", version = "0.1.31" } +gen-graph = { path = "../gen-graph", version = "0.1.31" } +thiserror = "2.0.18" + +[dev-dependencies] +petgraph = "0.6.5" +similar-asserts = "1.7.0" diff --git a/gen-parsers/src/blast.rs b/gen-parsers/src/blast.rs new file mode 100644 index 00000000..8d895155 --- /dev/null +++ b/gen-parsers/src/blast.rs @@ -0,0 +1,332 @@ +use std::{collections::HashMap, io::BufRead, marker::PhantomData}; + +use crate::{ParseError, ParsedAlignment, build_graph, validate_lengths}; + +pub struct BlastParser { + alignments: std::vec::IntoIter>, + _reader: PhantomData, +} + +impl BlastParser +where + R: BufRead, +{ + pub fn new(reader: R) -> Self { + let alignments = parse_blast_reader(reader); + Self { + alignments: alignments.into_iter(), + _reader: PhantomData, + } + } +} + +impl Iterator for BlastParser { + type Item = Result; + + fn next(&mut self) -> Option { + self.alignments.next() + } +} + +fn parse_blast_reader(reader: R) -> Vec> +where + R: BufRead, +{ + let lines = match reader.lines().collect::, _>>() { + Ok(lines) => lines, + Err(error) => return vec![Err(ParseError::Read(error.to_string()))], + }; + + let mut query_name = String::from("Query"); + let mut subject_name = String::from("Sbjct"); + let mut query_parts = Vec::new(); + let mut subject_parts = Vec::new(); + let mut results = Vec::new(); + let mut saw_score = false; + + for line in lines { + let trimmed = line.trim(); + if let Some(name) = trimmed.strip_prefix("Query=") { + query_name = name.trim().to_string(); + continue; + } + if let Some(name) = trimmed.strip_prefix('>') { + flush_blast_hsp( + &mut results, + &query_name, + &subject_name, + &mut query_parts, + &mut subject_parts, + ); + subject_name = name + .split_whitespace() + .next() + .unwrap_or("Sbjct") + .to_string(); + saw_score = false; + continue; + } + if trimmed.starts_with("Score =") { + flush_blast_hsp( + &mut results, + &query_name, + &subject_name, + &mut query_parts, + &mut subject_parts, + ); + saw_score = true; + continue; + } + if let Some(fragment) = blast_alignment_fragment(trimmed, "Query") { + query_parts.push(fragment.to_ascii_uppercase()); + continue; + } + if let Some(fragment) = blast_alignment_fragment(trimmed, "Sbjct") { + subject_parts.push(fragment.to_ascii_uppercase()); + } + } + + flush_blast_hsp( + &mut results, + &query_name, + &subject_name, + &mut query_parts, + &mut subject_parts, + ); + + if results.is_empty() && saw_score { + results.push(Err(ParseError::MissingBlastSubject)); + } else if results.is_empty() { + results.push(Err(ParseError::NoBlastAlignments)); + } + + results +} + +fn flush_blast_hsp( + results: &mut Vec>, + query_name: &str, + subject_name: &str, + query_parts: &mut Vec, + subject_parts: &mut Vec, +) { + if query_parts.is_empty() && subject_parts.is_empty() { + return; + } + if subject_parts.is_empty() { + results.push(Err(ParseError::MissingBlastSubject)); + query_parts.clear(); + return; + } + + let query_sequence = query_parts.join(""); + let subject_sequence = subject_parts.join(""); + query_parts.clear(); + subject_parts.clear(); + + let sequence_order = vec![query_name.to_string(), subject_name.to_string()]; + let mut aligned_sequences = HashMap::new(); + aligned_sequences.insert(query_name.to_string(), query_sequence); + aligned_sequences.insert(subject_name.to_string(), subject_sequence); + + let result = validate_lengths(&sequence_order, &aligned_sequences).map(|()| ParsedAlignment { + base_name: query_name.to_string(), + sequence_order: sequence_order.clone(), + graph: build_graph(query_name, &sequence_order, &aligned_sequences), + aligned_sequences, + }); + results.push(result); +} + +fn blast_alignment_fragment<'a>(line: &'a str, label: &str) -> Option<&'a str> { + let mut fields = line.split_whitespace(); + if fields.next()? != label { + return None; + } + let _start = fields.next()?; + fields.next() +} + +#[cfg(test)] +mod tests { + use gen_core::{PATH_END_NODE_ID, PATH_START_NODE_ID}; + use gen_graph::GraphNode; + use petgraph::visit::{EdgeRef as _, IntoEdgeReferences as _}; + use similar_asserts::assert_eq; + + use crate::{BlastParser, ParsedAlignment}; + + fn parse_all(input: &str) -> Vec { + BlastParser::new(input.as_bytes()) + .collect::, _>>() + .expect("should parse BLAST pairwise output") + } + + fn non_terminal_nodes(alignment: &ParsedAlignment) -> Vec { + let mut nodes = alignment + .graph + .nodes() + .filter(|node| node.node_id != PATH_START_NODE_ID && node.node_id != PATH_END_NODE_ID) + .collect::>(); + nodes.sort(); + nodes + } + + fn has_edge_for_sequence( + alignment: &ParsedAlignment, + source: GraphNode, + target: GraphNode, + index: i64, + ) -> bool { + alignment.graph.edge_references().any(|edge| { + edge.source() == source + && edge.target() == target + && edge + .weight() + .iter() + .any(|metadata| metadata.chromosome_index == index) + }) + } + + #[test] + fn parses_default_pairwise_hsp_into_base_relative_graph() { + let input = "\ +BLASTN 2.16.0+ + +Query= SeqA + +>SeqB +Length=28 + + Score = 52.8 bits (28), Expect = 2e-09 + Identities = 26/28 (93%), Gaps = 2/28 (7%) + Strand=Plus/Plus + +Query 1 MAGAASAVAAALAAA--AAGAAATAAAG 26 + ||||||||||||||| ||||||||||| +Sbjct 1 MAGAASAVAAALAAAGAAAGAAATAAAG 28 +"; + + let alignments = parse_all(input); + assert_eq!(alignments.len(), 1); + let alignment = &alignments[0]; + + assert_eq!(alignment.base_name, "SeqA"); + assert_eq!(alignment.sequence_order, vec!["SeqA", "SeqB"]); + assert_eq!( + alignment + .aligned_sequences + .get("SeqA") + .expect("should contain query"), + "MAGAASAVAAALAAA--AAGAAATAAAG" + ); + assert_eq!( + alignment.ungapped_sequence("SeqA"), + Some("MAGAASAVAAALAAAAAGAAATAAAG".to_string()) + ); + + let mut lengths = non_terminal_nodes(alignment) + .iter() + .map(GraphNode::length) + .collect::>(); + lengths.sort(); + assert_eq!(lengths, vec![2, 11, 15]); + } + + #[test] + fn parses_multiple_hsps_as_multiple_alignments() { + let input = "\ +BLASTN 2.16.0+ + +Query= QueryOne + +>SubjectOne +Length=10 + + Score = 10.0 bits (5), Expect = 0.01 + +Query 1 ACGT 4 + |||| +Sbjct 5 ACGT 8 + + Score = 8.0 bits (4), Expect = 0.02 + +Query 7 GG-T 9 + || | +Sbjct 1 GGAT 4 +"; + + let alignments = parse_all(input); + + assert_eq!(alignments.len(), 2); + assert_eq!(alignments[0].base_name, "QueryOne"); + assert_eq!(alignments[0].sequence_order, vec!["QueryOne", "SubjectOne"]); + assert_eq!( + alignments[1] + .aligned_sequences + .get("QueryOne") + .expect("should contain query"), + "GG-T" + ); + } + + #[test] + fn deletion_paths_skip_query_variant_nodes() { + let input = "\ +BLASTN 2.16.0+ + +Query= Base + +>Sample + + Score = 8.0 bits (4), Expect = 0.02 + +Query 1 ACGT 4 + | || +Sbjct 1 A-GT 3 +"; + + let alignments = parse_all(input); + let alignment = &alignments[0]; + let mut nodes = non_terminal_nodes(alignment); + assert_eq!(nodes.len(), 3); + nodes.sort_by_key(|node| node.sequence_start); + let source = nodes[0]; + let deleted_query_node = nodes[1]; + let target = nodes[2]; + + assert!( + has_edge_for_sequence(alignment, source, deleted_query_node, 0), + "query path should include the deleted query node" + ); + assert!( + has_edge_for_sequence(alignment, source, target, 1), + "subject path should skip the deleted query node" + ); + } + + #[test] + fn rejects_hsp_missing_subject_sequence() { + let input = "\ +BLASTN 2.16.0+ + +Query= SeqA + +>SeqB + + Score = 10.0 bits (5), Expect = 0.01 + +Query 1 ACGT 4 +"; + + let mut parser = BlastParser::new(input.as_bytes()); + let error = parser + .next() + .expect("should emit parse result") + .expect_err("should reject incomplete HSP"); + assert_eq!( + error.to_string(), + "BLAST HSP is missing a Sbjct sequence row" + ); + } +} diff --git a/gen-parsers/src/clustal.rs b/gen-parsers/src/clustal.rs new file mode 100644 index 00000000..f54152a7 --- /dev/null +++ b/gen-parsers/src/clustal.rs @@ -0,0 +1,340 @@ +use std::{collections::HashMap, io::BufRead}; + +use crate::{ParseError, ParsedAlignment, build_graph, validate_lengths}; + +pub struct ClustalwParser { + reader: Option, +} + +impl ClustalwParser { + pub const fn new(reader: R) -> Self { + Self { + reader: Some(reader), + } + } +} + +impl Iterator for ClustalwParser +where + R: BufRead, +{ + type Item = Result; + + fn next(&mut self) -> Option { + self.reader.take().map(parse_reader) + } +} + +fn parse_reader(reader: R) -> Result +where + R: BufRead, +{ + let mut lines = reader.lines(); + let Some(header) = lines.next() else { + return Err(ParseError::MissingHeader); + }; + let header = header.map_err(|error| ParseError::Read(error.to_string()))?; + if !is_clustal_header(&header) { + return Err(ParseError::MissingHeader); + } + + let mut sequence_order = Vec::new(); + let mut aligned_sequences = HashMap::new(); + let mut block = Vec::new(); + + for (line_index, line) in lines.enumerate() { + let line_number = line_index + 2; + let line = line.map_err(|error| ParseError::Read(error.to_string()))?; + if line.trim().is_empty() { + flush_block(&mut block, &mut aligned_sequences)?; + continue; + } + if is_conservation_line(&line) { + continue; + } + + let mut fields = line.split_whitespace(); + let name = fields + .next() + .ok_or(ParseError::MissingSequenceFragment { line_number })?; + let fragment = fields + .next() + .ok_or(ParseError::MissingSequenceFragment { line_number })?; + + let name = name.to_string(); + if !aligned_sequences.contains_key(&name) { + sequence_order.push(name.clone()); + } + + let count = fields.next().and_then(|field| field.parse::().ok()); + block.push(BlockRow { + name: name.clone(), + fragment: fragment.to_ascii_uppercase(), + count, + }); + } + flush_block(&mut block, &mut aligned_sequences)?; + + if sequence_order.is_empty() { + return Err(ParseError::NoSequences); + } + + validate_lengths(&sequence_order, &aligned_sequences)?; + let base_name = sequence_order[0].clone(); + let graph = build_graph(&base_name, &sequence_order, &aligned_sequences); + + Ok(ParsedAlignment { + base_name, + sequence_order, + aligned_sequences, + graph, + }) +} + +fn flush_block( + block: &mut Vec, + aligned_sequences: &mut HashMap, +) -> Result<(), ParseError> { + if block.is_empty() { + return Ok(()); + } + + let expected = block + .iter() + .map(|row| row.fragment.len()) + .max() + .expect("should have block rows"); + let has_counts = block.iter().any(|row| row.count.is_some()); + for row in block.drain(..) { + let actual = row.fragment.len(); + if actual != expected && !has_counts { + return Err(ParseError::MismatchedAlignedLength { + name: row.name, + actual, + expected, + }); + } + + let sequence = aligned_sequences.entry(row.name).or_default(); + sequence.push_str(&row.fragment); + for _ in actual..expected { + sequence.push('-'); + } + } + + Ok(()) +} + +fn is_clustal_header(header: &str) -> bool { + let header = header.trim_start(); + header.starts_with("CLUSTAL") +} + +fn is_conservation_line(line: &str) -> bool { + line.starts_with(char::is_whitespace) + && line + .trim() + .chars() + .all(|character| matches!(character, '*' | ':' | '.' | ' ')) +} + +#[derive(Debug)] +struct BlockRow { + name: String, + fragment: String, + count: Option, +} + +#[cfg(test)] +mod tests { + use gen_core::{PATH_END_NODE_ID, PATH_START_NODE_ID}; + use gen_graph::GraphNode; + use petgraph::visit::{EdgeRef as _, IntoEdgeReferences as _}; + use similar_asserts::assert_eq; + + use crate::{ClustalwParser, ParsedAlignment}; + + fn parse_single(input: &str) -> ParsedAlignment { + let mut parser = ClustalwParser::new(input.as_bytes()); + let alignment = parser + .next() + .expect("should emit one alignment") + .expect("should parse CLUSTAL alignment"); + assert!( + parser.next().is_none(), + "CLUSTAL parser should emit one alignment per file" + ); + alignment + } + + fn non_terminal_nodes(alignment: &ParsedAlignment) -> Vec { + let mut nodes = alignment + .graph + .nodes() + .filter(|node| node.node_id != PATH_START_NODE_ID && node.node_id != PATH_END_NODE_ID) + .collect::>(); + nodes.sort(); + nodes + } + + fn has_edge_for_sequence( + alignment: &ParsedAlignment, + source: GraphNode, + target: GraphNode, + index: i64, + ) -> bool { + alignment.graph.edge_references().any(|edge| { + edge.source() == source + && edge.target() == target + && edge + .weight() + .iter() + .any(|metadata| metadata.chromosome_index == index) + }) + } + + #[test] + fn parses_clustal_header_without_w() { + let input = "\ +CLUSTAL 2.1 multiple sequence alignment + +SeqA ACGT +SeqB ACGT + **** +"; + + let alignment = parse_single(input); + + assert_eq!(alignment.base_name, "SeqA"); + assert_eq!(alignment.sequence_order, vec!["SeqA", "SeqB"]); + assert_eq!( + alignment.ungapped_sequence("SeqB"), + Some("ACGT".to_string()) + ); + } + + #[test] + fn parses_pairwise_alignment_into_base_relative_graph() { + let input = "\ +CLUSTAL W (2.1) multiple sequence alignment + +SeqA MAGAASAVAAALAAA--AAGAAATAAAG +SeqB MAGAASAVAAALAAAGAAAGAAATAAAG + ************* ************ +"; + + let alignment = parse_single(input); + + assert_eq!(alignment.base_name, "SeqA"); + assert_eq!(alignment.sequence_order, vec!["SeqA", "SeqB"]); + assert_eq!( + alignment + .aligned_sequences + .get("SeqA") + .expect("should contain SeqA"), + "MAGAASAVAAALAAA--AAGAAATAAAG" + ); + assert_eq!( + alignment.ungapped_sequence("SeqA"), + Some("MAGAASAVAAALAAAAAGAAATAAAG".to_string()) + ); + + let nodes = non_terminal_nodes(&alignment); + let mut lengths = nodes.iter().map(GraphNode::length).collect::>(); + lengths.sort(); + assert_eq!(lengths, vec![2, 11, 15]); + } + + #[test] + fn parses_multi_block_alignment_and_keeps_sequence_order() { + let input = "\ +CLUSTAL W (1.83) multiple sequence alignment + +Seq1 ATGCCTAGCTAGCTAGCATCGATCGATCGATCGATCGTACGATCGATCGATCGATCGATC 60 +Seq2 ATGCCTAGCTAGCTAGC---TCGATCGATCGATCGATCGTACGATCGATCGATCGATC---- 53 +Seq3 ATGCCT---------GC---TCGATCGATCGATC----GTACGATCGATCGATCGATCGATC 50 + ****** ** ****** *** ********************** + +Seq1 TACGATCGATCGTACGTA 78 +Seq2 ----ATCGATCGTACGTA 69 +Seq3 TACGATCGATCGTACGTA 68 + ********** +"; + + let alignment = parse_single(input); + + assert_eq!(alignment.sequence_order, vec!["Seq1", "Seq2", "Seq3"]); + assert_eq!(alignment.base_name, "Seq1"); + assert_eq!( + alignment + .ungapped_sequence("Seq1") + .expect("should contain Seq1") + .len(), + 78 + ); + assert_eq!( + alignment + .ungapped_sequence("Seq2") + .expect("should contain Seq2") + .len(), + 69 + ); + assert_eq!( + alignment + .ungapped_sequence("Seq3") + .expect("should contain Seq3") + .len(), + 64 + ); + } + + #[test] + fn deletion_paths_skip_base_variant_nodes() { + let input = "\ +CLUSTALW multiple sequence alignment + +Base ACGT +Sample A-GT + * ** +"; + + let alignment = parse_single(input); + let nodes = non_terminal_nodes(&alignment); + assert_eq!(nodes.len(), 3); + + let mut base_nodes = nodes; + base_nodes.sort_by_key(|node| node.sequence_start); + let source = base_nodes[0]; + let deleted_base_node = base_nodes[1]; + let target = base_nodes[2]; + + assert!( + has_edge_for_sequence(&alignment, source, deleted_base_node, 0), + "base path should include the deleted base node" + ); + assert!( + has_edge_for_sequence(&alignment, source, target, 1), + "sample path should skip the deleted base node" + ); + } + + #[test] + fn rejects_mismatched_aligned_lengths() { + let input = "\ +CLUSTAL W multiple sequence alignment + +SeqA ACGT +SeqB ACG +"; + + let mut parser = ClustalwParser::new(input.as_bytes()); + let error = parser + .next() + .expect("should emit parse result") + .expect_err("should reject malformed alignment"); + assert_eq!( + error.to_string(), + "aligned sequence SeqB has length 3 but expected 4" + ); + } +} diff --git a/gen-parsers/src/lib.rs b/gen-parsers/src/lib.rs new file mode 100644 index 00000000..94cf31a9 --- /dev/null +++ b/gen-parsers/src/lib.rs @@ -0,0 +1,335 @@ +use std::collections::HashMap; + +use gen_core::{HashId, PATH_END_NODE_ID, PATH_START_NODE_ID, Strand}; +use gen_graph::{GenGraph, GraphEdge, GraphNode}; +use thiserror::Error; + +pub mod blast; +pub mod clustal; +pub mod maf; +pub mod paf; +pub mod psl; + +pub use blast::BlastParser; +pub use clustal::ClustalwParser; +pub use maf::MafParser; +pub use paf::PafParser; +pub use psl::PslParser; + +#[derive(Debug, Error, PartialEq, Eq)] +pub enum ParseError { + #[error("failed to read alignment input: {0}")] + Read(String), + #[error("CLUSTAL alignment must start with CLUSTAL W or CLUSTALW")] + MissingHeader, + #[error("CLUSTAL alignment does not contain sequence rows")] + NoSequences, + #[error("alignment row {line_number} is missing a sequence fragment")] + MissingSequenceFragment { line_number: usize }, + #[error("BLAST HSP is missing a Sbjct sequence row")] + MissingBlastSubject, + #[error("BLAST output does not contain pairwise HSP alignments")] + NoBlastAlignments, + #[error("PAF line {line_number} has {actual} fields but expected at least 12")] + PafTooFewFields { line_number: usize, actual: usize }, + #[error("PAF line {line_number} is missing required cg:Z CIGAR tag")] + MissingPafCigar { line_number: usize }, + #[error("PAF line {line_number} has invalid integer in field {field}")] + InvalidPafInteger { line_number: usize, field: usize }, + #[error("PAF line {line_number} has invalid strand {strand}")] + InvalidPafStrand { line_number: usize, strand: String }, + #[error("PAF line {line_number} has invalid CIGAR string {cigar}")] + InvalidCigar { line_number: usize, cigar: String }, + #[error("PSL line {line_number} has {actual} fields but expected at least 21")] + PslTooFewFields { line_number: usize, actual: usize }, + #[error("PSL line {line_number} has invalid integer in field {field}")] + InvalidPslInteger { line_number: usize, field: usize }, + #[error("PSL line {line_number} has invalid strand {strand}")] + InvalidPslStrand { line_number: usize, strand: String }, + #[error( + "PSL line {line_number} blockCount is {block_count} but blockSizes, qStarts, and tStarts contain {block_sizes}, {query_starts}, and {target_starts} entries" + )] + MismatchedPslBlocks { + line_number: usize, + block_count: usize, + block_sizes: usize, + query_starts: usize, + target_starts: usize, + }, + #[error("MAF block ending at line {line_number} does not contain sequence rows")] + MafBlockWithoutSequences { line_number: usize }, + #[error("MAF line {line_number} has {actual} fields but expected at least 7")] + MafTooFewFields { line_number: usize, actual: usize }, + #[error("MAF line {line_number} has invalid integer in field {field}")] + InvalidMafInteger { line_number: usize, field: usize }, + #[error("MAF line {line_number} has invalid strand {strand}")] + InvalidMafStrand { line_number: usize, strand: String }, + #[error( + "MAF line {line_number} sequence {name} declares size {declared} but has {actual} non-gap bases" + )] + MafSizeMismatch { + line_number: usize, + name: String, + declared: i64, + actual: usize, + }, + #[error("aligned sequence {name} has length {actual} but expected {expected}")] + MismatchedAlignedLength { + name: String, + actual: usize, + expected: usize, + }, +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum CigarOp { + Match(i64), + Insertion(i64), + Deletion(i64), + ReferenceSkip(i64), + SoftClip(i64), + HardClip(i64), + Equal(i64), + Difference(i64), +} + +#[derive(Debug)] +pub struct ParsedAlignment { + pub base_name: String, + pub sequence_order: Vec, + pub aligned_sequences: HashMap, + pub graph: GenGraph, +} + +impl ParsedAlignment { + pub fn ungapped_sequence(&self, name: &str) -> Option { + self.aligned_sequences + .get(name) + .map(|sequence| sequence.chars().filter(|base| *base != '-').collect()) + } +} + +#[derive(Debug)] +pub struct ParsedMapping { + pub query_name: String, + pub query_len: i64, + pub query_start: i64, + pub query_end: i64, + pub strand: Strand, + pub target_name: String, + pub target_len: i64, + pub target_start: i64, + pub target_end: i64, + pub matching_bases: i64, + pub block_len: i64, + pub mapping_quality: i64, + pub cigar: Vec, + pub tags: HashMap, + pub graph: GenGraph, +} + +pub(crate) fn validate_lengths( + sequence_order: &[String], + aligned_sequences: &HashMap, +) -> Result<(), ParseError> { + let expected = aligned_sequences + .get(&sequence_order[0]) + .expect("should have base sequence") + .len(); + + for name in sequence_order.iter().skip(1) { + let actual = aligned_sequences + .get(name) + .expect("should have sequence") + .len(); + if actual != expected { + return Err(ParseError::MismatchedAlignedLength { + name: name.clone(), + actual, + expected, + }); + } + } + + Ok(()) +} + +pub(crate) fn build_graph( + base_name: &str, + sequence_order: &[String], + aligned_sequences: &HashMap, +) -> GenGraph { + let base_aligned = aligned_sequences + .get(base_name) + .expect("should have base sequence"); + let base_sequence = ungapped(base_aligned); + let base_node_id = HashId::convert_str(&format!("clustalw:base:{base_name}:{base_sequence}")); + let runs = alignment_runs(base_aligned, sequence_order, aligned_sequences); + + let mut graph = GenGraph::new(); + let mut previous_by_sequence = vec![terminal_node(PATH_START_NODE_ID); sequence_order.len()]; + let end_node = terminal_node(PATH_END_NODE_ID); + let mut base_offset = 0_i64; + + for run in runs { + let base_text = ungapped(&base_aligned[run.start..run.end]); + let base_node = (!base_text.is_empty()).then(|| { + let node = GraphNode { + node_id: base_node_id, + sequence_start: base_offset, + sequence_end: base_offset + base_text.len() as i64, + }; + base_offset += base_text.len() as i64; + node + }); + + for (sequence_index, name) in sequence_order.iter().enumerate() { + let aligned = aligned_sequences.get(name).expect("should have sequence"); + let text = ungapped(&aligned[run.start..run.end]); + let next = if run.kind == RunKind::Invariant || sequence_index == 0 { + base_node + } else if text.is_empty() { + None + } else if text == base_text { + base_node + } else { + Some(GraphNode { + node_id: HashId::convert_str(&format!( + "clustalw:variant:{name}:{}:{}:{text}", + run.start, run.end + )), + sequence_start: 0, + sequence_end: text.len() as i64, + }) + }; + + if let Some(next_node) = next { + add_path_edge( + &mut graph, + previous_by_sequence[sequence_index], + next_node, + sequence_index as i64, + ); + previous_by_sequence[sequence_index] = next_node; + } + } + } + + for (sequence_index, previous) in previous_by_sequence.iter().enumerate() { + add_path_edge(&mut graph, *previous, end_node, sequence_index as i64); + } + + graph +} + +fn alignment_runs( + base_aligned: &str, + sequence_order: &[String], + aligned_sequences: &HashMap, +) -> Vec { + let mut runs = Vec::new(); + let mut start = 0; + let mut current_kind = None; + + for column in 0..base_aligned.len() { + let kind = column_kind(column, base_aligned, sequence_order, aligned_sequences); + if current_kind.is_some_and(|current| current != kind) { + runs.push(AlignmentRun { + start, + end: column, + kind: current_kind.expect("should have run kind"), + }); + start = column; + } + current_kind = Some(kind); + } + + if let Some(kind) = current_kind { + runs.push(AlignmentRun { + start, + end: base_aligned.len(), + kind, + }); + } + + runs +} + +fn column_kind( + column: usize, + base_aligned: &str, + sequence_order: &[String], + aligned_sequences: &HashMap, +) -> RunKind { + let base = base_aligned.as_bytes()[column]; + if base == b'-' { + return RunKind::Variant; + } + + let invariant = sequence_order.iter().all(|name| { + aligned_sequences + .get(name) + .expect("should have sequence") + .as_bytes()[column] + == base + }); + + if invariant { + RunKind::Invariant + } else { + RunKind::Variant + } +} + +pub(crate) fn add_path_edge( + graph: &mut GenGraph, + source: GraphNode, + target: GraphNode, + sequence_index: i64, +) { + let edge = GraphEdge { + edge_id: HashId::convert_str(&format!( + "clustalw:edge:{sequence_index}:{}:{}:{}:{}", + source.node_id, source.sequence_start, target.node_id, target.sequence_start + )), + source_strand: Strand::Forward, + target_strand: Strand::Forward, + chromosome_index: sequence_index, + phased: 0, + created_on: 0, + }; + + if let Some(edges) = graph.edge_weight_mut(source, target) { + edges.push(edge); + } else { + graph.add_edge(source, target, vec![edge]); + } +} + +pub(crate) fn terminal_node(node_id: HashId) -> GraphNode { + GraphNode { + node_id, + sequence_start: 0, + sequence_end: 0, + } +} + +pub(crate) fn ungapped(sequence: &str) -> String { + sequence + .chars() + .filter(|character| *character != '-') + .collect() +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +struct AlignmentRun { + start: usize, + end: usize, + kind: RunKind, +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum RunKind { + Invariant, + Variant, +} diff --git a/gen-parsers/src/maf.rs b/gen-parsers/src/maf.rs new file mode 100644 index 00000000..7818a314 --- /dev/null +++ b/gen-parsers/src/maf.rs @@ -0,0 +1,354 @@ +use std::{collections::HashMap, io::BufRead, marker::PhantomData}; + +use gen_core::Strand; + +use crate::{ParseError, ParsedAlignment, build_graph, validate_lengths}; + +pub struct MafParser { + alignments: std::vec::IntoIter>, + _reader: PhantomData, +} + +impl MafParser +where + R: BufRead, +{ + pub fn new(reader: R) -> Self { + let alignments = parse_maf_reader(reader); + Self { + alignments: alignments.into_iter(), + _reader: PhantomData, + } + } +} + +impl Iterator for MafParser { + type Item = Result; + + fn next(&mut self) -> Option { + self.alignments.next() + } +} + +fn parse_maf_reader(reader: R) -> Vec> +where + R: BufRead, +{ + let mut results = Vec::new(); + let mut rows = Vec::new(); + let mut in_alignment = false; + + for (line_index, line) in reader.lines().enumerate() { + let line_number = line_index + 1; + let line = match line { + Ok(line) => line, + Err(error) => { + results.push(Err(ParseError::Read(error.to_string()))); + continue; + } + }; + let trimmed = line.trim(); + + if trimmed.is_empty() { + flush_maf_block(&mut results, &mut rows, line_number, &mut in_alignment); + continue; + } + if should_skip_maf_line(trimmed) { + continue; + } + if trimmed.starts_with("a ") || trimmed == "a" { + flush_maf_block(&mut results, &mut rows, line_number, &mut in_alignment); + in_alignment = true; + continue; + } + if trimmed.starts_with("s ") { + match parse_maf_sequence_row(trimmed, line_number) { + Ok(row) => rows.push(row), + Err(error) => results.push(Err(error)), + } + } + } + + flush_maf_block(&mut results, &mut rows, 0, &mut in_alignment); + results +} + +fn should_skip_maf_line(line: &str) -> bool { + line.starts_with('#') + || line.starts_with("track ") + || line.starts_with("i ") + || line.starts_with("e ") + || line.starts_with("q ") +} + +fn flush_maf_block( + results: &mut Vec>, + rows: &mut Vec, + line_number: usize, + in_alignment: &mut bool, +) { + if !*in_alignment { + return; + } + if rows.is_empty() { + results.push(Err(ParseError::MafBlockWithoutSequences { line_number })); + *in_alignment = false; + return; + } + + let result = build_maf_alignment(rows); + results.push(result); + rows.clear(); + *in_alignment = false; +} + +fn parse_maf_sequence_row(line: &str, line_number: usize) -> Result { + let fields = line.split_whitespace().collect::>(); + if fields.len() < 7 { + return Err(ParseError::MafTooFewFields { + line_number, + actual: fields.len(), + }); + } + + let name = fields[1].to_string(); + let _start = parse_maf_i64(fields[2], line_number, 3)?; + let size = parse_maf_i64(fields[3], line_number, 4)?; + let _strand = match fields[4] { + "+" => Strand::Forward, + "-" => Strand::Reverse, + strand => { + return Err(ParseError::InvalidMafStrand { + line_number, + strand: strand.to_string(), + }); + } + }; + let _source_size = parse_maf_i64(fields[5], line_number, 6)?; + let text = fields[6].to_ascii_uppercase(); + let actual = text.chars().filter(|character| *character != '-').count(); + if actual != size as usize { + return Err(ParseError::MafSizeMismatch { + line_number, + name, + declared: size, + actual, + }); + } + + Ok(MafSequenceRow { name, text }) +} + +fn parse_maf_i64(field: &str, line_number: usize, field_number: usize) -> Result { + field.parse().map_err(|_| ParseError::InvalidMafInteger { + line_number, + field: field_number, + }) +} + +fn build_maf_alignment(rows: &[MafSequenceRow]) -> Result { + let base = rows.first().expect("should have MAF sequence row"); + let mut sequence_order = Vec::new(); + let mut aligned_sequences = HashMap::new(); + + for row in rows { + sequence_order.push(row.name.clone()); + let text = if row.name == base.name { + row.text.clone() + } else { + expand_maf_dots(&base.text, &row.text) + }; + aligned_sequences.insert(row.name.clone(), text); + } + + validate_lengths(&sequence_order, &aligned_sequences)?; + let graph = build_graph(&base.name, &sequence_order, &aligned_sequences); + Ok(ParsedAlignment { + base_name: base.name.clone(), + sequence_order, + graph, + aligned_sequences, + }) +} + +fn expand_maf_dots(base: &str, text: &str) -> String { + base.chars() + .zip(text.chars()) + .map(|(base_character, character)| { + if character == '.' { + base_character + } else { + character + } + }) + .collect() +} + +#[derive(Clone, Debug)] +struct MafSequenceRow { + name: String, + text: String, +} + +#[cfg(test)] +mod tests { + use gen_core::{PATH_END_NODE_ID, PATH_START_NODE_ID}; + use gen_graph::GraphNode; + use petgraph::visit::{EdgeRef as _, IntoEdgeReferences as _}; + use similar_asserts::assert_eq; + + use crate::{MafParser, ParsedAlignment}; + + fn parse_all(input: &str) -> Vec { + MafParser::new(input.as_bytes()) + .collect::, _>>() + .expect("should parse MAF alignments") + } + + fn non_terminal_nodes(alignment: &ParsedAlignment) -> Vec { + let mut nodes = alignment + .graph + .nodes() + .filter(|node| node.node_id != PATH_START_NODE_ID && node.node_id != PATH_END_NODE_ID) + .collect::>(); + nodes.sort(); + nodes + } + + fn has_edge_for_sequence( + alignment: &ParsedAlignment, + source: GraphNode, + target: GraphNode, + index: i64, + ) -> bool { + alignment.graph.edge_references().any(|edge| { + edge.source() == source + && edge.target() == target + && edge + .weight() + .iter() + .any(|metadata| metadata.chromosome_index == index) + }) + } + + #[test] + fn parses_maf_block_into_base_relative_graph() { + let input = "\ +##maf version=1 scoring=tba.v8 +# generated example + +a score=23262.0 +s hg16.chr7 27707221 13 + 158545518 gcagctgaaaaca +s panTro1.chr6 28869787 13 + 161576975 gcagctgaaaaca +s mm4.chr6 53310102 13 + 151104725 ACAGCTGAAAATA +"; + + let alignments = parse_all(input); + assert_eq!(alignments.len(), 1); + let alignment = &alignments[0]; + + assert_eq!(alignment.base_name, "hg16.chr7"); + assert_eq!( + alignment.sequence_order, + vec!["hg16.chr7", "panTro1.chr6", "mm4.chr6"] + ); + assert_eq!( + alignment + .aligned_sequences + .get("hg16.chr7") + .expect("should contain reference"), + "GCAGCTGAAAACA" + ); + + let mut lengths = non_terminal_nodes(alignment) + .iter() + .map(GraphNode::length) + .collect::>(); + lengths.sort(); + assert_eq!(lengths, vec![1, 1, 1, 1, 1, 10]); + } + + #[test] + fn parses_multiple_blocks_and_expands_dot_notation() { + let input = "\ +track name=sample +##maf version=1 + +a score=1 +s ref.chr1 0 4 + 100 ACGT +s alt.chr1 5 4 + 100 .... + +a score=2 +s ref.chr1 10 3 + 100 GG-T +s alt.chr1 20 4 + 100 ..AT +"; + + let alignments = parse_all(input); + + assert_eq!(alignments.len(), 2); + assert_eq!( + alignments[0] + .aligned_sequences + .get("alt.chr1") + .expect("should contain first alternate"), + "ACGT" + ); + assert_eq!( + alignments[1] + .aligned_sequences + .get("alt.chr1") + .expect("should contain second alternate"), + "GGAT" + ); + } + + #[test] + fn deletion_paths_skip_base_variant_nodes() { + let input = "\ +##maf version=1 + +a score=3 +s ref 0 4 + 10 ACGT +s alt 0 3 + 10 A-GT +"; + + let alignments = parse_all(input); + let alignment = &alignments[0]; + let mut nodes = non_terminal_nodes(alignment); + assert_eq!(nodes.len(), 3); + nodes.sort_by_key(|node| node.sequence_start); + let source = nodes[0]; + let deleted_base_node = nodes[1]; + let target = nodes[2]; + + assert!( + has_edge_for_sequence(alignment, source, deleted_base_node, 0), + "base path should include deleted base interval" + ); + assert!( + has_edge_for_sequence(alignment, source, target, 1), + "alternate path should skip deleted base interval" + ); + } + + #[test] + fn rejects_maf_size_mismatch() { + let input = "\ +##maf version=1 + +a score=4 +s ref 0 4 + 10 ACGT +s alt 0 4 + 10 A-GT +"; + + let mut parser = MafParser::new(input.as_bytes()); + let error = parser + .next() + .expect("should emit parse result") + .expect_err("should reject malformed MAF size"); + assert_eq!( + error.to_string(), + "MAF line 5 sequence alt declares size 4 but has 3 non-gap bases" + ); + } +} diff --git a/gen-parsers/src/paf.rs b/gen-parsers/src/paf.rs new file mode 100644 index 00000000..88de7337 --- /dev/null +++ b/gen-parsers/src/paf.rs @@ -0,0 +1,402 @@ +use std::{collections::HashMap, io::BufRead, marker::PhantomData}; + +use gen_core::{HashId, PATH_END_NODE_ID, PATH_START_NODE_ID, Strand}; +use gen_graph::{GenGraph, GraphNode}; + +use crate::{CigarOp, ParseError, ParsedMapping, add_path_edge, terminal_node}; + +pub struct PafParser { + mappings: std::vec::IntoIter>, + _reader: PhantomData, +} + +impl PafParser +where + R: BufRead, +{ + pub fn new(reader: R) -> Self { + let mappings = parse_paf_reader(reader); + Self { + mappings: mappings.into_iter(), + _reader: PhantomData, + } + } +} + +impl Iterator for PafParser { + type Item = Result; + + fn next(&mut self) -> Option { + self.mappings.next() + } +} + +fn parse_paf_reader(reader: R) -> Vec> +where + R: BufRead, +{ + reader + .lines() + .enumerate() + .filter_map(|(line_index, line)| { + let line_number = line_index + 1; + let line = match line { + Ok(line) => line, + Err(error) => return Some(Err(ParseError::Read(error.to_string()))), + }; + if line.trim().is_empty() { + return None; + } + Some(parse_paf_line(&line, line_number)) + }) + .collect() +} + +fn parse_paf_line(line: &str, line_number: usize) -> Result { + let fields = line.split('\t').collect::>(); + if fields.len() < 12 { + return Err(ParseError::PafTooFewFields { + line_number, + actual: fields.len(), + }); + } + + let query_name = fields[0].to_string(); + let query_len = parse_paf_i64(fields[1], line_number, 2)?; + let query_start = parse_paf_i64(fields[2], line_number, 3)?; + let query_end = parse_paf_i64(fields[3], line_number, 4)?; + let strand = match fields[4] { + "+" => Strand::Forward, + "-" => Strand::Reverse, + strand => { + return Err(ParseError::InvalidPafStrand { + line_number, + strand: strand.to_string(), + }); + } + }; + let target_name = fields[5].to_string(); + let target_len = parse_paf_i64(fields[6], line_number, 7)?; + let target_start = parse_paf_i64(fields[7], line_number, 8)?; + let target_end = parse_paf_i64(fields[8], line_number, 9)?; + let matching_bases = parse_paf_i64(fields[9], line_number, 10)?; + let block_len = parse_paf_i64(fields[10], line_number, 11)?; + let mapping_quality = parse_paf_i64(fields[11], line_number, 12)?; + + let mut tags = HashMap::new(); + let mut cigar = None; + for field in fields.iter().skip(12) { + let mut parts = field.splitn(3, ':'); + let Some(key) = parts.next() else { + continue; + }; + let Some(tag_type) = parts.next() else { + continue; + }; + let Some(value) = parts.next() else { + continue; + }; + tags.insert(key.to_string(), format!("{tag_type}:{value}")); + if key == "cg" && tag_type == "Z" { + cigar = Some(parse_cigar(value, line_number)?); + } + } + + let Some(cigar) = cigar else { + return Err(ParseError::MissingPafCigar { line_number }); + }; + + let mut mapping = ParsedMapping { + query_name, + query_len, + query_start, + query_end, + strand, + target_name, + target_len, + target_start, + target_end, + matching_bases, + block_len, + mapping_quality, + cigar, + tags, + graph: GenGraph::new(), + }; + mapping.graph = build_paf_graph(&mapping); + Ok(mapping) +} + +fn parse_paf_i64(field: &str, line_number: usize, field_number: usize) -> Result { + field.parse().map_err(|_| ParseError::InvalidPafInteger { + line_number, + field: field_number, + }) +} + +fn parse_cigar(cigar: &str, line_number: usize) -> Result, ParseError> { + let mut ops = Vec::new(); + let mut len = 0_i64; + let mut has_digits = false; + + for character in cigar.chars() { + if let Some(digit) = character.to_digit(10) { + has_digits = true; + len = len * 10 + i64::from(digit); + continue; + } + + if !has_digits || len == 0 { + return Err(ParseError::InvalidCigar { + line_number, + cigar: cigar.to_string(), + }); + } + + let op = match character { + 'M' => CigarOp::Match(len), + 'I' => CigarOp::Insertion(len), + 'D' => CigarOp::Deletion(len), + 'N' => CigarOp::ReferenceSkip(len), + 'S' => CigarOp::SoftClip(len), + 'H' => CigarOp::HardClip(len), + '=' => CigarOp::Equal(len), + 'X' => CigarOp::Difference(len), + _ => { + return Err(ParseError::InvalidCigar { + line_number, + cigar: cigar.to_string(), + }); + } + }; + ops.push(op); + len = 0; + has_digits = false; + } + + if has_digits || ops.is_empty() { + return Err(ParseError::InvalidCigar { + line_number, + cigar: cigar.to_string(), + }); + } + + Ok(ops) +} + +fn build_paf_graph(mapping: &ParsedMapping) -> GenGraph { + let mut graph = GenGraph::new(); + let mut target_previous = terminal_node(PATH_START_NODE_ID); + let mut query_previous = terminal_node(PATH_START_NODE_ID); + let end_node = terminal_node(PATH_END_NODE_ID); + let mut target_offset = mapping.target_start; + let mut query_offset = mapping.query_start; + let target_node_id = HashId::convert_str(&format!( + "paf:target:{}:{}", + mapping.target_name, mapping.target_len + )); + let query_node_id = HashId::convert_str(&format!( + "paf:query:{}:{}", + mapping.query_name, mapping.query_len + )); + + for op in &mapping.cigar { + match *op { + CigarOp::Match(len) | CigarOp::Equal(len) | CigarOp::Difference(len) => { + let node = GraphNode { + node_id: target_node_id, + sequence_start: target_offset, + sequence_end: target_offset + len, + }; + add_path_edge(&mut graph, target_previous, node, 0); + add_path_edge(&mut graph, query_previous, node, 1); + target_previous = node; + query_previous = node; + target_offset += len; + query_offset += len; + } + CigarOp::Insertion(len) => { + let node = GraphNode { + node_id: query_node_id, + sequence_start: query_offset, + sequence_end: query_offset + len, + }; + add_path_edge(&mut graph, query_previous, node, 1); + query_previous = node; + query_offset += len; + } + CigarOp::Deletion(len) | CigarOp::ReferenceSkip(len) => { + let node = GraphNode { + node_id: target_node_id, + sequence_start: target_offset, + sequence_end: target_offset + len, + }; + add_path_edge(&mut graph, target_previous, node, 0); + target_previous = node; + target_offset += len; + } + CigarOp::SoftClip(len) => { + query_offset += len; + } + CigarOp::HardClip(_) => {} + } + } + + add_path_edge(&mut graph, target_previous, end_node, 0); + add_path_edge(&mut graph, query_previous, end_node, 1); + graph +} + +#[cfg(test)] +mod tests { + use gen_core::{PATH_END_NODE_ID, PATH_START_NODE_ID, Strand}; + use gen_graph::GraphNode; + use petgraph::visit::{EdgeRef as _, IntoEdgeReferences as _}; + use similar_asserts::assert_eq; + + use crate::{CigarOp, PafParser, ParsedMapping}; + + fn parse_all(input: &str) -> Vec { + PafParser::new(input.as_bytes()) + .collect::, _>>() + .expect("should parse PAF mappings") + } + + fn non_terminal_nodes(mapping: &ParsedMapping) -> Vec { + let mut nodes = mapping + .graph + .nodes() + .filter(|node| node.node_id != PATH_START_NODE_ID && node.node_id != PATH_END_NODE_ID) + .collect::>(); + nodes.sort(); + nodes + } + + fn has_edge_for_sequence( + mapping: &ParsedMapping, + source: GraphNode, + target: GraphNode, + index: i64, + ) -> bool { + mapping.graph.edge_references().any(|edge| { + edge.source() == source + && edge.target() == target + && edge + .weight() + .iter() + .any(|metadata| metadata.chromosome_index == index) + }) + } + + #[test] + fn parses_paf_with_cigar_into_target_relative_graph() { + let input = + "query1\t100\t10\t24\t+\ttarget1\t200\t50\t63\t12\t14\t60\tcg:Z:4M2I3M1D5M\ttp:A:P\n"; + + let mappings = parse_all(input); + assert_eq!(mappings.len(), 1); + let mapping = &mappings[0]; + + assert_eq!(mapping.query_name, "query1"); + assert_eq!(mapping.query_len, 100); + assert_eq!(mapping.query_start, 10); + assert_eq!(mapping.query_end, 24); + assert_eq!(mapping.strand, Strand::Forward); + assert_eq!(mapping.target_name, "target1"); + assert_eq!(mapping.target_len, 200); + assert_eq!(mapping.target_start, 50); + assert_eq!(mapping.target_end, 63); + assert_eq!(mapping.matching_bases, 12); + assert_eq!(mapping.block_len, 14); + assert_eq!(mapping.mapping_quality, 60); + assert_eq!( + mapping.cigar, + vec![ + CigarOp::Match(4), + CigarOp::Insertion(2), + CigarOp::Match(3), + CigarOp::Deletion(1), + CigarOp::Match(5), + ] + ); + assert_eq!( + mapping.tags.get("tp").expect("should retain optional tag"), + "A:P" + ); + + let mut lengths = non_terminal_nodes(mapping) + .iter() + .map(GraphNode::length) + .collect::>(); + lengths.sort(); + assert_eq!(lengths, vec![1, 2, 3, 4, 5]); + } + + #[test] + fn query_path_uses_insertions_and_skips_target_deletions() { + let input = "query1\t100\t10\t24\t+\ttarget1\t200\t50\t63\t12\t14\t60\tcg:Z:4M2I3M1D5M\n"; + let mappings = parse_all(input); + let mapping = &mappings[0]; + let mut target_nodes = non_terminal_nodes(mapping) + .into_iter() + .filter(|node| node.sequence_start >= 50) + .collect::>(); + target_nodes.sort_by_key(|node| node.sequence_start); + + let target_match_left = target_nodes[0]; + let target_match_middle = target_nodes[1]; + let target_deletion = target_nodes[2]; + let target_match_right = target_nodes[3]; + let query_insertion = non_terminal_nodes(mapping) + .into_iter() + .find(|node| node.sequence_start == 14 && node.sequence_end == 16) + .expect("should contain query insertion node"); + + assert!( + has_edge_for_sequence(mapping, target_match_left, query_insertion, 1), + "query path should enter inserted query interval" + ); + assert!( + has_edge_for_sequence(mapping, query_insertion, target_match_middle, 1), + "query path should leave inserted query interval" + ); + assert!( + has_edge_for_sequence(mapping, target_match_middle, target_deletion, 0), + "target path should include deleted target interval" + ); + assert!( + has_edge_for_sequence(mapping, target_match_middle, target_match_right, 1), + "query path should skip deleted target interval" + ); + } + + #[test] + fn parses_multiple_paf_lines_and_reverse_strand() { + let input = "\ +query1\t100\t0\t4\t+\ttarget1\t200\t0\t4\t4\t4\t60\tcg:Z:4M +query2\t80\t5\t10\t-\ttarget2\t90\t30\t35\t5\t5\t255\tcg:Z:5= +"; + + let mappings = parse_all(input); + + assert_eq!(mappings.len(), 2); + assert_eq!(mappings[0].strand, Strand::Forward); + assert_eq!(mappings[1].strand, Strand::Reverse); + assert_eq!(mappings[1].cigar, vec![CigarOp::Equal(5)]); + } + + #[test] + fn rejects_paf_without_cigar_tag() { + let input = "query1\t100\t10\t24\t+\ttarget1\t200\t50\t63\t12\t14\t60\n"; + + let mut parser = PafParser::new(input.as_bytes()); + let error = parser + .next() + .expect("should emit parse result") + .expect_err("should reject PAF without CIGAR"); + assert_eq!( + error.to_string(), + "PAF line 1 is missing required cg:Z CIGAR tag" + ); + } +} diff --git a/gen-parsers/src/psl.rs b/gen-parsers/src/psl.rs new file mode 100644 index 00000000..e9cd13bd --- /dev/null +++ b/gen-parsers/src/psl.rs @@ -0,0 +1,432 @@ +use std::{collections::HashMap, io::BufRead, marker::PhantomData}; + +use gen_core::{HashId, PATH_END_NODE_ID, PATH_START_NODE_ID, Strand}; +use gen_graph::{GenGraph, GraphNode}; + +use crate::{ParseError, ParsedMapping, add_path_edge, terminal_node}; + +pub struct PslParser { + mappings: std::vec::IntoIter>, + _reader: PhantomData, +} + +impl PslParser +where + R: BufRead, +{ + pub fn new(reader: R) -> Self { + let mappings = parse_psl_reader(reader); + Self { + mappings: mappings.into_iter(), + _reader: PhantomData, + } + } +} + +impl Iterator for PslParser { + type Item = Result; + + fn next(&mut self) -> Option { + self.mappings.next() + } +} + +fn parse_psl_reader(reader: R) -> Vec> +where + R: BufRead, +{ + reader + .lines() + .enumerate() + .filter_map(|(line_index, line)| { + let line_number = line_index + 1; + let line = match line { + Ok(line) => line, + Err(error) => return Some(Err(ParseError::Read(error.to_string()))), + }; + if should_skip_psl_line(&line) { + return None; + } + Some(parse_psl_line(&line, line_number)) + }) + .collect() +} + +fn should_skip_psl_line(line: &str) -> bool { + let trimmed = line.trim(); + trimmed.is_empty() + || trimmed.starts_with("psLayout") + || trimmed.starts_with("match ") + || trimmed.starts_with("-----") + || trimmed.starts_with("browser ") + || trimmed.starts_with("track ") +} + +fn parse_psl_line(line: &str, line_number: usize) -> Result { + let fields = line.split_whitespace().collect::>(); + if fields.len() < 21 { + return Err(ParseError::PslTooFewFields { + line_number, + actual: fields.len(), + }); + } + + let matches = parse_psl_i64(fields[0], line_number, 1)?; + let mismatches = parse_psl_i64(fields[1], line_number, 2)?; + let rep_matches = parse_psl_i64(fields[2], line_number, 3)?; + let n_count = parse_psl_i64(fields[3], line_number, 4)?; + let query_insert_count = parse_psl_i64(fields[4], line_number, 5)?; + let query_insert_bases = parse_psl_i64(fields[5], line_number, 6)?; + let target_insert_count = parse_psl_i64(fields[6], line_number, 7)?; + let target_insert_bases = parse_psl_i64(fields[7], line_number, 8)?; + let strand = parse_psl_strand(fields[8], line_number)?; + let query_name = fields[9].to_string(); + let query_len = parse_psl_i64(fields[10], line_number, 11)?; + let query_start = parse_psl_i64(fields[11], line_number, 12)?; + let query_end = parse_psl_i64(fields[12], line_number, 13)?; + let target_name = fields[13].to_string(); + let target_len = parse_psl_i64(fields[14], line_number, 15)?; + let target_start = parse_psl_i64(fields[15], line_number, 16)?; + let target_end = parse_psl_i64(fields[16], line_number, 17)?; + let block_count = parse_psl_usize(fields[17], line_number, 18)?; + let block_sizes = parse_psl_i64_list(fields[18], line_number, 19)?; + let query_starts = parse_psl_i64_list(fields[19], line_number, 20)?; + let target_starts = parse_psl_i64_list(fields[20], line_number, 21)?; + + if block_sizes.len() != block_count + || query_starts.len() != block_count + || target_starts.len() != block_count + { + return Err(ParseError::MismatchedPslBlocks { + line_number, + block_count, + block_sizes: block_sizes.len(), + query_starts: query_starts.len(), + target_starts: target_starts.len(), + }); + } + + let blocks = block_sizes + .iter() + .zip(query_starts.iter()) + .zip(target_starts.iter()) + .map(|((size, query_start), target_start)| PslBlock { + size: *size, + query_start: normalize_psl_query_start(*query_start, *size, query_len, strand), + target_start: *target_start, + }) + .collect::>(); + let mut tags = HashMap::new(); + tags.insert("psl:misMatches".to_string(), mismatches.to_string()); + tags.insert("psl:repMatches".to_string(), rep_matches.to_string()); + tags.insert("psl:nCount".to_string(), n_count.to_string()); + tags.insert("psl:qNumInsert".to_string(), query_insert_count.to_string()); + tags.insert( + "psl:qBaseInsert".to_string(), + query_insert_bases.to_string(), + ); + tags.insert( + "psl:tNumInsert".to_string(), + target_insert_count.to_string(), + ); + tags.insert( + "psl:tBaseInsert".to_string(), + target_insert_bases.to_string(), + ); + + let mut mapping = ParsedMapping { + query_name, + query_len, + query_start, + query_end, + strand, + target_name, + target_len, + target_start, + target_end, + matching_bases: matches, + block_len: block_sizes.iter().sum(), + mapping_quality: 0, + cigar: Vec::new(), + tags, + graph: GenGraph::new(), + }; + mapping.graph = build_psl_graph(&mapping, &blocks); + Ok(mapping) +} + +fn parse_psl_i64(field: &str, line_number: usize, field_number: usize) -> Result { + field.parse().map_err(|_| ParseError::InvalidPslInteger { + line_number, + field: field_number, + }) +} + +fn parse_psl_usize( + field: &str, + line_number: usize, + field_number: usize, +) -> Result { + field.parse().map_err(|_| ParseError::InvalidPslInteger { + line_number, + field: field_number, + }) +} + +fn parse_psl_i64_list( + field: &str, + line_number: usize, + field_number: usize, +) -> Result, ParseError> { + field + .trim_end_matches(',') + .split(',') + .filter(|value| !value.is_empty()) + .map(|value| parse_psl_i64(value, line_number, field_number)) + .collect() +} + +fn parse_psl_strand(field: &str, line_number: usize) -> Result { + match field.chars().next() { + Some('+') => Ok(Strand::Forward), + Some('-') => Ok(Strand::Reverse), + _ => Err(ParseError::InvalidPslStrand { + line_number, + strand: field.to_string(), + }), + } +} + +fn normalize_psl_query_start(query_start: i64, size: i64, query_len: i64, strand: Strand) -> i64 { + if strand == Strand::Reverse { + query_len - query_start - size + } else { + query_start + } +} + +fn build_psl_graph(mapping: &ParsedMapping, blocks: &[PslBlock]) -> GenGraph { + let mut graph = GenGraph::new(); + let mut target_previous = terminal_node(PATH_START_NODE_ID); + let mut query_previous = terminal_node(PATH_START_NODE_ID); + let end_node = terminal_node(PATH_END_NODE_ID); + let target_node_id = HashId::convert_str(&format!( + "psl:target:{}:{}", + mapping.target_name, mapping.target_len + )); + let query_node_id = HashId::convert_str(&format!( + "psl:query:{}:{}", + mapping.query_name, mapping.query_len + )); + + for (index, block) in blocks.iter().enumerate() { + if let Some(previous_block) = index + .checked_sub(1) + .and_then(|previous| blocks.get(previous)) + { + let target_gap_start = previous_block.target_start + previous_block.size; + let target_gap_end = block.target_start; + let (query_gap_start, query_gap_end) = + psl_query_gap(previous_block, block, mapping.strand); + + if target_gap_end > target_gap_start { + let node = GraphNode { + node_id: target_node_id, + sequence_start: target_gap_start, + sequence_end: target_gap_end, + }; + add_path_edge(&mut graph, target_previous, node, 0); + target_previous = node; + } + + if query_gap_end > query_gap_start { + let node = GraphNode { + node_id: query_node_id, + sequence_start: query_gap_start, + sequence_end: query_gap_end, + }; + add_path_edge(&mut graph, query_previous, node, 1); + query_previous = node; + } + } + + let node = GraphNode { + node_id: target_node_id, + sequence_start: block.target_start, + sequence_end: block.target_start + block.size, + }; + add_path_edge(&mut graph, target_previous, node, 0); + add_path_edge(&mut graph, query_previous, node, 1); + target_previous = node; + query_previous = node; + } + + add_path_edge(&mut graph, target_previous, end_node, 0); + add_path_edge(&mut graph, query_previous, end_node, 1); + graph +} + +fn psl_query_gap(previous_block: &PslBlock, block: &PslBlock, strand: Strand) -> (i64, i64) { + if strand == Strand::Reverse { + (block.query_start + block.size, previous_block.query_start) + } else { + ( + previous_block.query_start + previous_block.size, + block.query_start, + ) + } +} + +#[derive(Clone, Copy, Debug)] +struct PslBlock { + size: i64, + query_start: i64, + target_start: i64, +} + +#[cfg(test)] +mod tests { + use gen_core::{PATH_END_NODE_ID, PATH_START_NODE_ID, Strand}; + use gen_graph::GraphNode; + use petgraph::visit::{EdgeRef as _, IntoEdgeReferences as _}; + use similar_asserts::assert_eq; + + use crate::{ParsedMapping, PslParser}; + + fn parse_all(input: &str) -> Vec { + PslParser::new(input.as_bytes()) + .collect::, _>>() + .expect("should parse PSL mappings") + } + + fn non_terminal_nodes(mapping: &ParsedMapping) -> Vec { + let mut nodes = mapping + .graph + .nodes() + .filter(|node| node.node_id != PATH_START_NODE_ID && node.node_id != PATH_END_NODE_ID) + .collect::>(); + nodes.sort(); + nodes + } + + fn has_edge_for_sequence( + mapping: &ParsedMapping, + source: GraphNode, + target: GraphNode, + index: i64, + ) -> bool { + mapping.graph.edge_references().any(|edge| { + edge.source() == source + && edge.target() == target + && edge + .weight() + .iter() + .any(|metadata| metadata.chromosome_index == index) + }) + } + + #[test] + fn parses_psl_blocks_into_target_relative_graph() { + let input = "\ +psLayout version 3 + +match mis- rep. N's Q gap Q gap T gap T gap strand Q Q Q Q T T T T block blockSizes qStarts tStarts +----- ---- ---- --- ----- ----- ----- ----- ------ -------- ---- ----- --- -------- ---- ----- --- ----- ---------- -------- -------- +59 9 0 0 1 823 1 96 + query1 1200 10 33 target1 5000 100 128 2 10,13, 10,20, 100,115, +"; + + let mappings = parse_all(input); + assert_eq!(mappings.len(), 1); + let mapping = &mappings[0]; + + assert_eq!(mapping.query_name, "query1"); + assert_eq!(mapping.query_len, 1200); + assert_eq!(mapping.query_start, 10); + assert_eq!(mapping.query_end, 33); + assert_eq!(mapping.strand, Strand::Forward); + assert_eq!(mapping.target_name, "target1"); + assert_eq!(mapping.target_len, 5000); + assert_eq!(mapping.target_start, 100); + assert_eq!(mapping.target_end, 128); + assert_eq!(mapping.matching_bases, 59); + assert_eq!(mapping.block_len, 23); + assert_eq!(mapping.mapping_quality, 0); + assert_eq!(mapping.tags.get("psl:misMatches"), Some(&"9".to_string())); + + let mut lengths = non_terminal_nodes(mapping) + .iter() + .map(GraphNode::length) + .collect::>(); + lengths.sort(); + assert_eq!(lengths, vec![5, 10, 13]); + } + + #[test] + fn query_path_uses_query_gap_and_skips_target_gap() { + let input = + "23 0 0 0 1 2 1 5 + query1 100 10 35 target1 200 50 78 2 10,13, 10,22, 50,65,\n"; + + let mappings = parse_all(input); + let mapping = &mappings[0]; + let mut target_nodes = non_terminal_nodes(mapping) + .into_iter() + .filter(|node| node.sequence_start >= 50) + .collect::>(); + target_nodes.sort_by_key(|node| node.sequence_start); + + let target_left = target_nodes[0]; + let target_gap = target_nodes[1]; + let target_right = target_nodes[2]; + let query_gap = non_terminal_nodes(mapping) + .into_iter() + .find(|node| node.sequence_start == 20 && node.sequence_end == 22) + .expect("should contain query insertion interval"); + + assert!( + has_edge_for_sequence(mapping, target_left, target_gap, 0), + "target path should include target-only gap interval" + ); + assert!( + has_edge_for_sequence(mapping, target_left, query_gap, 1), + "query path should include query-only gap interval" + ); + assert!( + has_edge_for_sequence(mapping, query_gap, target_right, 1), + "query path should rejoin the next aligned target block" + ); + } + + #[test] + fn normalizes_negative_query_block_starts() { + let input = + "38 0 0 0 1 14 1 17 - query1 61 4 56 target1 1000 100 155 2 20,18, 5,39, 100,137,\n"; + + let mappings = parse_all(input); + let mapping = &mappings[0]; + + assert_eq!(mapping.strand, Strand::Reverse); + assert_eq!(mapping.query_start, 4); + assert_eq!(mapping.query_end, 56); + assert!( + non_terminal_nodes(mapping) + .iter() + .any(|node| node.sequence_start == 22 && node.sequence_end == 36), + "minus-strand query gap should use normalized forward query coordinates" + ); + } + + #[test] + fn rejects_psl_with_mismatched_block_lists() { + let input = "23 0 0 0 0 0 0 0 + query1 100 10 33 target1 200 50 73 2 10,13, 10, 50,65,\n"; + + let mut parser = PslParser::new(input.as_bytes()); + let error = parser + .next() + .expect("should emit parse result") + .expect_err("should reject malformed PSL block lists"); + assert_eq!( + error.to_string(), + "PSL line 1 blockCount is 2 but blockSizes, qStarts, and tStarts contain 2, 1, and 2 entries" + ); + } +} diff --git a/src/commands/import/alignment.rs b/src/commands/import/alignment.rs new file mode 100644 index 00000000..5d1c26c5 --- /dev/null +++ b/src/commands/import/alignment.rs @@ -0,0 +1,58 @@ +use anyhow::Result; +use clap::Args; +use gen_models::errors::OperationError; + +use crate::{ + commands::{cli_context::CliContext, get_default_collection}, + imports::alignment::{AlignmentImportError, import_alignment_aln}, +}; + +/// Import an alignment file +#[derive(Debug, Args, Clone)] +pub struct Command { + /// Alignment file path + #[clap(index = 1)] + pub path: String, + /// The name of the collection to store the alignment under + #[arg(short, long)] + name: Option, + /// Override the sample name for all imported alignment rows + #[arg(short, long)] + sample: Option, +} + +pub fn execute(cli_context: &CliContext, cmd: Command) -> Result<()> { + println!("Alignment import called"); + + let context = cli_context.context; + let operation_conn = context.operations().conn(); + let conn = context.graph().conn(); + + conn.execute("BEGIN TRANSACTION", []).unwrap(); + operation_conn.execute("BEGIN TRANSACTION", []).unwrap(); + + let name = &cmd + .name + .clone() + .unwrap_or_else(|| get_default_collection(operation_conn)); + + match import_alignment_aln(context, &cmd.path.clone(), name, cmd.sample.as_deref()) { + Ok(_) => { + println!("Alignment imported."); + conn.execute("END TRANSACTION;", []).unwrap(); + operation_conn.execute("END TRANSACTION;", []).unwrap(); + Ok(()) + } + Err(AlignmentImportError::OperationError(OperationError::NoChanges)) => { + conn.execute("ROLLBACK TRANSACTION;", []).unwrap(); + operation_conn.execute("ROLLBACK TRANSACTION;", []).unwrap(); + println!("Alignment contents already exist."); + Ok(()) + } + Err(e) => { + conn.execute("ROLLBACK TRANSACTION;", []).unwrap(); + operation_conn.execute("ROLLBACK TRANSACTION;", []).unwrap(); + Err(e.into()) + } + } +} diff --git a/src/commands/import/mod.rs b/src/commands/import/mod.rs index aa47de9d..6703807e 100644 --- a/src/commands/import/mod.rs +++ b/src/commands/import/mod.rs @@ -3,6 +3,7 @@ use clap::{Args, Subcommand}; use crate::commands::cli_context::CliContext; +mod alignment; mod fasta; mod genbank; mod gfa; @@ -17,6 +18,8 @@ pub struct Command { /// Import commands #[derive(Clone, Debug, Subcommand)] pub enum Commands { + /// Import alignment + Alignment(alignment::Command), /// Import fasta Fasta(fasta::Command), /// Import genbank @@ -29,6 +32,7 @@ pub enum Commands { pub fn execute(ctx: &CliContext, command: Command) -> anyhow::Result<()> { match command.command { + Commands::Alignment(cmd) => crate::commands::import::alignment::execute(ctx, cmd), Commands::Fasta(cmd) => crate::commands::import::fasta::execute(ctx, cmd), Commands::Genbank(cmd) => crate::commands::import::genbank::execute(ctx, cmd), Commands::Gfa(cmd) => crate::commands::import::gfa::execute(ctx, cmd), diff --git a/src/imports.rs b/src/imports.rs index ebdf137a..e5115a89 100644 --- a/src/imports.rs +++ b/src/imports.rs @@ -1,3 +1,4 @@ +pub mod alignment; pub mod fasta; pub mod genbank; pub mod gfa; diff --git a/src/imports/alignment.rs b/src/imports/alignment.rs new file mode 100644 index 00000000..8dea8360 --- /dev/null +++ b/src/imports/alignment.rs @@ -0,0 +1,749 @@ +use std::{ + collections::HashSet, + io::{BufRead, BufReader}, + path::PathBuf, +}; + +use flate2::read::MultiGzDecoder; +use gen_core::{HashId, PATH_END_NODE_ID, PATH_START_NODE_ID, Strand}; +use gen_graph::GraphNode; +use gen_models::{ + assets::AssetUri, + block_group::{BlockGroup, NewBlockGroup}, + block_group_edge::{BlockGroupEdge, BlockGroupEdgeData}, + collection::Collection, + db::DbContext, + edge::Edge, + errors::{ + BlockGroupError, CollectionError, EdgeError, FileAdditionError, NodeError, OperationError, + PathError, SampleError, SequenceError, + }, + file_types::FileTypes, + node::Node, + operations::{Operation, OperationFile, OperationInfo}, + path::Path, + sample::{NewSample, Sample}, + sample_lineage::SampleLineage, + sequence::Sequence, + session_operations::{end_operation, start_operation}, + traits::Query as _, +}; +use gen_parsers::{ClustalwParser, ParseError, ParsedAlignment}; +use noodles::bgzf; +use petgraph::visit::{EdgeRef as _, IntoEdgeReferences as _}; +use thiserror::Error; + +use crate::progress_bar::{add_saving_operation_bar, get_handler, get_progress_bar}; + +#[derive(Debug, Error)] +pub enum AlignmentImportError { + #[error("Operation Error: {0}")] + OperationError(#[from] OperationError), + #[error("Asset Error: {0}")] + FileAdditionError(#[from] FileAdditionError), + #[error("IO Error: {0}")] + IOError(#[from] std::io::Error), + #[error("Collection creation error: {0}")] + CollectionError(#[from] CollectionError), + #[error("Sample creation error: {0}")] + SampleError(#[from] SampleError), + #[error("Sample lineage error: {0}")] + SampleLineageError(#[from] rusqlite::Error), + #[error("Block group write error: {0}")] + BlockGroupError(#[from] BlockGroupError), + #[error("Edge write error: {0}")] + EdgeError(#[from] EdgeError), + #[error("Node creation error: {0}")] + NodeError(#[from] NodeError), + #[error("Path creation error: {0}")] + PathError(#[from] PathError), + #[error("Sequence save error: {0}")] + SequenceError(#[from] SequenceError), + #[error("Alignment parse error: {0}")] + ParseError(#[from] ParseError), + #[error("Alignment file reader did not reach EOF before checksum was requested")] + MissingChecksum, +} + +pub fn import_alignment_aln( + context: &DbContext, + alignment_path: &String, + collection_name: &str, + sample_override: Option<&str>, +) -> Result { + let conn = context.graph().conn(); + let progress_bar = get_handler(); + let mut session = start_operation(conn); + let path = PathBuf::from(alignment_path); + + let asset_uri = ::new(context.workspace(), alignment_path); + let file = asset_uri.reader(context.workspace())?; + let checksum_handle = file.checksum_handle(); + let reader_stream: Box = match path.extension().and_then(|ext| ext.to_str()) { + Some("gz") => Box::new(BufReader::new(MultiGzDecoder::new(file))), + Some("bgz") => Box::new(bgzf::io::Reader::new(file)), + _ => Box::new(BufReader::new(file)), + }; + + let collection = match Collection::create(conn, collection_name) { + Ok(collection) => collection, + Err(CollectionError::Duplicate(collection)) => collection, + Err(e) => return Err(AlignmentImportError::CollectionError(e)), + }; + + let _ = progress_bar.println("Parsing alignment"); + let bar = progress_bar.add(get_progress_bar(None)); + bar.set_message("Alignments Processed."); + let mut summary = String::new(); + for alignment in ClustalwParser::new(reader_stream) { + let alignment = alignment?; + persist_alignment(conn, &collection.name, &alignment, sample_override)?; + summary.push_str(&format!( + " {}: {} aligned sequences.\n", + alignment.base_name, + alignment.sequence_order.len() + )); + bar.inc(1); + } + bar.finish(); + + let checksum_override = checksum_handle + .checksum() + .ok_or(AlignmentImportError::MissingChecksum)?; + let bar = add_saving_operation_bar(&progress_bar); + let op = end_operation( + context, + &mut session, + &OperationInfo { + files: vec![ + OperationFile::new(alignment_path.to_string()) + .set_file_type(FileTypes::None) + .set_checksum_override(checksum_override), + ], + description: "alignment_import".to_string(), + }, + &summary, + None, + ) + .map_err(AlignmentImportError::OperationError); + bar.finish(); + op +} + +fn persist_alignment( + conn: &gen_models::db::GraphConnection, + collection_name: &str, + alignment: &ParsedAlignment, + sample_override: Option<&str>, +) -> Result<(), AlignmentImportError> { + let graph = build_import_graph(conn, collection_name, alignment)?; + let sequence_names = parse_alignment_sequence_names(alignment, sample_override); + persist_alignment_samples(conn, &sequence_names)?; + + let parent_sequence_name = sequence_names + .first() + .expect("should contain at least one alignment sequence"); + let parent_block_group = BlockGroup::create( + conn, + NewBlockGroup { + collection_name, + sample_name: &parent_sequence_name.sample_name, + name: &parent_sequence_name.block_group_name, + ..Default::default() + }, + )?; + persist_alignment_path( + conn, + &graph, + &parent_block_group.id, + parent_sequence_name, + 0, + )?; + + for (index, sequence_name) in sequence_names.iter().enumerate().skip(1) { + let block_group = BlockGroup::create( + conn, + NewBlockGroup { + collection_name, + sample_name: &sequence_name.sample_name, + name: &sequence_name.block_group_name, + parent_block_group_id: Some(&parent_block_group.id), + ..Default::default() + }, + )?; + block_group.copy_contents_from(conn, &parent_block_group)?; + delete_paths_for_block_group(conn, &block_group.id); + persist_alignment_path(conn, &graph, &block_group.id, sequence_name, index as i64)?; + } + + Ok(()) +} + +fn persist_alignment_path( + conn: &gen_models::db::GraphConnection, + graph: &gen_graph::GenGraph, + block_group_id: &HashId, + sequence_name: &AlignmentSequenceName, + chromosome_index: i64, +) -> Result<(), AlignmentImportError> { + let (edge_ids, block_group_edges) = + persist_path_edges(conn, graph, block_group_id, chromosome_index)?; + BlockGroupEdge::bulk_create(conn, &block_group_edges); + Path::create(conn, &sequence_name.sample_name, block_group_id, &edge_ids)?; + Ok(()) +} + +fn delete_paths_for_block_group(conn: &gen_models::db::GraphConnection, block_group_id: &HashId) { + for path in Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![block_group_id], + ) { + Path::delete(conn, &path.name, block_group_id); + } +} + +fn parse_alignment_sequence_names( + alignment: &ParsedAlignment, + sample_override: Option<&str>, +) -> Vec { + alignment + .sequence_order + .iter() + .map(|name| parse_alignment_sequence_name(name, sample_override)) + .collect() +} + +fn parse_alignment_sequence_name( + name: &str, + sample_override: Option<&str>, +) -> AlignmentSequenceName { + if let Some(sample_name) = sample_override { + return AlignmentSequenceName { + sample_name: sample_name.to_string(), + block_group_name: name.to_string(), + }; + } + + match name.split_once('.') { + Some((sample_name, block_group_name)) + if !sample_name.is_empty() && !block_group_name.is_empty() => + { + AlignmentSequenceName { + sample_name: sample_name.to_string(), + block_group_name: block_group_name.to_string(), + } + } + _ => AlignmentSequenceName { + sample_name: name.to_string(), + block_group_name: name.to_string(), + }, + } +} + +fn persist_alignment_samples( + conn: &gen_models::db::GraphConnection, + sequence_names: &[AlignmentSequenceName], +) -> Result<(), AlignmentImportError> { + let parent_sample = sequence_names + .first() + .expect("should contain at least one alignment sequence"); + + for name in sequence_names { + match Sample::get_or_create( + conn, + NewSample { + name: &name.sample_name, + ..Default::default() + }, + ) { + Ok(_) => {} + Err(SampleError::Duplicate(_)) => {} + Err(e) => return Err(AlignmentImportError::SampleError(e)), + } + } + + for child_sample in sequence_names.iter().skip(1) { + if parent_sample.sample_name != child_sample.sample_name { + SampleLineage::create(conn, &parent_sample.sample_name, &child_sample.sample_name)?; + } + } + + Ok(()) +} + +fn persist_path_edges( + conn: &gen_models::db::GraphConnection, + graph: &gen_graph::GenGraph, + block_group_id: &HashId, + chromosome_index: i64, +) -> Result<(Vec, Vec), AlignmentImportError> { + let mut edge_ids = Vec::new(); + let mut block_group_edges = Vec::new(); + let mut current = terminal_node(PATH_START_NODE_ID); + + while current.node_id != PATH_END_NODE_ID { + let Some((target, metadata)) = next_path_edge(graph, current, chromosome_index) else { + break; + }; + let stored_edge = Edge::create( + conn, + current.node_id, + current.sequence_end, + metadata.source_strand, + target.node_id, + target.sequence_start, + metadata.target_strand, + )?; + edge_ids.push(stored_edge.id); + block_group_edges.push(BlockGroupEdgeData { + block_group_id: *block_group_id, + edge_id: stored_edge.id, + chromosome_index, + phased: metadata.phased, + }); + current = target; + } + + Ok((edge_ids, block_group_edges)) +} + +fn next_path_edge( + graph: &gen_graph::GenGraph, + source: GraphNode, + chromosome_index: i64, +) -> Option<(GraphNode, gen_graph::GraphEdge)> { + graph.edge_references().find_map(|edge| { + if edge.source() != source { + return None; + } + edge.weight() + .iter() + .find(|metadata| metadata.chromosome_index == chromosome_index) + .copied() + .map(|metadata| (edge.target(), metadata)) + }) +} + +fn build_import_graph( + conn: &gen_models::db::GraphConnection, + collection_name: &str, + alignment: &ParsedAlignment, +) -> Result { + let base_aligned = alignment + .aligned_sequences + .get(&alignment.base_name) + .expect("should contain base sequence"); + let base_sequence = ungapped(base_aligned); + let base_sequence_record = Sequence::new() + .sequence_type("DNA") + .sequence(&base_sequence) + .save(conn)?; + let base_node_id = Node::create( + conn, + &base_sequence_record.hash, + &HashId::convert_str(&format!( + "alignment:{collection_name}:{}:{}", + alignment.base_name, base_sequence_record.hash + )), + )?; + + let mut graph = gen_graph::GenGraph::new(); + let mut previous_by_sequence = + vec![terminal_node(PATH_START_NODE_ID); alignment.sequence_order.len()]; + let end_node = terminal_node(PATH_END_NODE_ID); + let mut base_offset = 0_i64; + let runs = alignment_runs(alignment); + let mut created_variant_nodes = HashSet::new(); + + for run in runs { + let base_text = ungapped(&base_aligned[run.start..run.end]); + let base_node = (!base_text.is_empty()).then(|| { + let node = GraphNode { + node_id: base_node_id, + sequence_start: base_offset, + sequence_end: base_offset + base_text.len() as i64, + }; + base_offset += base_text.len() as i64; + node + }); + + for (sequence_index, name) in alignment.sequence_order.iter().enumerate() { + let aligned = alignment + .aligned_sequences + .get(name) + .expect("should contain aligned sequence"); + let text = ungapped(&aligned[run.start..run.end]); + let next = if run.kind == AlignmentRunKind::Invariant || sequence_index == 0 { + base_node + } else if text.is_empty() { + None + } else if text == base_text { + base_node + } else { + let node_id = HashId::convert_str(&format!( + "alignment:{collection_name}:{name}:{}:{}:{text}", + run.start, run.end + )); + if created_variant_nodes.insert(node_id) { + let sequence = Sequence::new() + .sequence_type("DNA") + .sequence(&text) + .save(conn)?; + Node::create(conn, &sequence.hash, &node_id)?; + } + Some(GraphNode { + node_id, + sequence_start: 0, + sequence_end: text.len() as i64, + }) + }; + + if let Some(next_node) = next { + add_import_edge( + &mut graph, + previous_by_sequence[sequence_index], + next_node, + sequence_index as i64, + ); + previous_by_sequence[sequence_index] = next_node; + } + } + } + + for (sequence_index, previous) in previous_by_sequence.iter().enumerate() { + add_import_edge(&mut graph, *previous, end_node, sequence_index as i64); + } + + Ok(graph) +} + +fn add_import_edge( + graph: &mut gen_graph::GenGraph, + source: GraphNode, + target: GraphNode, + sequence_index: i64, +) { + let edge = gen_graph::GraphEdge { + edge_id: HashId::convert_str(&format!( + "alignment:edge:{sequence_index}:{}:{}:{}:{}", + source.node_id, source.sequence_end, target.node_id, target.sequence_start + )), + source_strand: Strand::Forward, + target_strand: Strand::Forward, + chromosome_index: sequence_index, + phased: 0, + created_on: 0, + }; + + if let Some(edges) = graph.edge_weight_mut(source, target) { + edges.push(edge); + } else { + graph.add_edge(source, target, vec![edge]); + } +} + +fn alignment_runs(alignment: &ParsedAlignment) -> Vec { + let base_aligned = alignment + .aligned_sequences + .get(&alignment.base_name) + .expect("should contain base sequence"); + let mut runs = Vec::new(); + let mut start = 0; + let mut current_kind = None; + + for column in 0..base_aligned.len() { + let kind = alignment_column_kind(alignment, column); + if current_kind.is_some_and(|current| current != kind) { + runs.push(AlignmentRun { + start, + end: column, + kind: current_kind.expect("should have run kind"), + }); + start = column; + } + current_kind = Some(kind); + } + + if let Some(kind) = current_kind { + runs.push(AlignmentRun { + start, + end: base_aligned.len(), + kind, + }); + } + + runs +} + +fn alignment_column_kind(alignment: &ParsedAlignment, column: usize) -> AlignmentRunKind { + let base_aligned = alignment + .aligned_sequences + .get(&alignment.base_name) + .expect("should contain base sequence"); + let base = base_aligned.as_bytes()[column]; + if base == b'-' { + return AlignmentRunKind::Variant; + } + + let invariant = alignment.sequence_order.iter().all(|name| { + alignment + .aligned_sequences + .get(name) + .expect("should contain aligned sequence") + .as_bytes()[column] + == base + }); + + if invariant { + AlignmentRunKind::Invariant + } else { + AlignmentRunKind::Variant + } +} + +fn terminal_node(node_id: HashId) -> GraphNode { + GraphNode { + node_id, + sequence_start: 0, + sequence_end: 0, + } +} + +fn ungapped(sequence: &str) -> String { + sequence + .chars() + .filter(|character| *character != '-') + .collect() +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +struct AlignmentRun { + start: usize, + end: usize, + kind: AlignmentRunKind, +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum AlignmentRunKind { + Invariant, + Variant, +} + +#[derive(Debug, Eq, PartialEq)] +struct AlignmentSequenceName { + sample_name: String, + block_group_name: String, +} + +#[cfg(test)] +mod tests { + use std::path::PathBuf; + + use gen_models::{ + block_group::BlockGroup, path::Path, sample_lineage::SampleLineage, traits::Query as _, + }; + + use crate::{ + imports::alignment::import_alignment_aln, test_helpers::setup_gen, track_database, + }; + + #[test] + fn imports_aln_alignment_as_samples_with_base_lineage() { + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + track_database(conn, op_conn).unwrap(); + + let aln_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/alignments/simple.aln"); + + import_alignment_aln( + &context, + &aln_path + .to_str() + .expect("should have fixture path") + .to_string(), + "test", + None, + ) + .unwrap(); + + let base_block_group_id = BlockGroup::get_id("test", "SeqA", "chr1", None); + let base_paths = Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![base_block_group_id], + ); + assert_eq!( + base_paths + .iter() + .map(|path| path.name.as_str()) + .collect::>(), + vec!["SeqA"] + ); + assert_eq!(base_paths[0].sequence(conn).unwrap(), "ACGT"); + + let child_block_group_id = + BlockGroup::get_id("test", "SeqB", "chr1", Some(&base_block_group_id)); + let child_block_group = BlockGroup::get_by_id(conn, &child_block_group_id).unwrap(); + assert_eq!( + child_block_group.parent_block_group_id, + Some(base_block_group_id) + ); + let child_paths = Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![child_block_group_id], + ); + assert_eq!( + child_paths + .iter() + .map(|path| path.name.as_str()) + .collect::>(), + vec!["SeqB"] + ); + assert_eq!(child_paths[0].sequence(conn).unwrap(), "AGT"); + let second_child_block_group_id = + BlockGroup::get_id("test", "SeqC", "chr1", Some(&base_block_group_id)); + let second_child_block_group = + BlockGroup::get_by_id(conn, &second_child_block_group_id).unwrap(); + assert_eq!( + second_child_block_group.parent_block_group_id, + Some(base_block_group_id) + ); + let second_child_paths = Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![second_child_block_group_id], + ); + assert_eq!( + second_child_paths + .iter() + .map(|path| path.name.as_str()) + .collect::>(), + vec!["SeqC"] + ); + assert_eq!(second_child_paths[0].sequence(conn).unwrap(), "ACG"); + assert_eq!( + SampleLineage::get_children(conn, "SeqA"), + vec!["SeqB", "SeqC"] + ); + assert_eq!(SampleLineage::get_parents(conn, "SeqB"), vec!["SeqA"]); + assert_eq!(SampleLineage::get_parents(conn, "SeqC"), vec!["SeqA"]); + } + + #[test] + fn imports_plain_alignment_names_as_sample_and_block_group_names() { + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + track_database(conn, op_conn).unwrap(); + + let aln_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/alignments/plain_names.aln"); + + import_alignment_aln( + &context, + &aln_path + .to_str() + .expect("should have fixture path") + .to_string(), + "test", + None, + ) + .unwrap(); + + let base_block_group_id = BlockGroup::get_id("test", "Base", "Base", None); + let base_paths = Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![base_block_group_id], + ); + assert_eq!( + base_paths + .iter() + .map(|path| path.name.as_str()) + .collect::>(), + vec!["Base"] + ); + + let child_block_group_id = + BlockGroup::get_id("test", "Child", "Child", Some(&base_block_group_id)); + let child_block_group = BlockGroup::get_by_id(conn, &child_block_group_id).unwrap(); + assert_eq!( + child_block_group.parent_block_group_id, + Some(base_block_group_id) + ); + let child_paths = Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![child_block_group_id], + ); + assert_eq!( + child_paths + .iter() + .map(|path| path.name.as_str()) + .collect::>(), + vec!["Child"] + ); + assert_eq!(SampleLineage::get_children(conn, "Base"), vec!["Child"]); + } + + #[test] + fn imports_alignment_with_sample_override_and_raw_block_group_names() { + let context = setup_gen(); + let conn = context.graph().conn(); + let op_conn = context.operations().conn(); + track_database(conn, op_conn).unwrap(); + + let aln_path = + PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("fixtures/alignments/plain_names.aln"); + + import_alignment_aln( + &context, + &aln_path + .to_str() + .expect("should have fixture path") + .to_string(), + "test", + Some("override"), + ) + .unwrap(); + + let base_block_group_id = BlockGroup::get_id("test", "override", "Base", None); + let base_paths = Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![base_block_group_id], + ); + assert_eq!( + base_paths + .iter() + .map(|path| path.name.as_str()) + .collect::>(), + vec!["override"] + ); + + let child_block_group_id = + BlockGroup::get_id("test", "override", "Child", Some(&base_block_group_id)); + let child_block_group = BlockGroup::get_by_id(conn, &child_block_group_id).unwrap(); + assert_eq!( + child_block_group.parent_block_group_id, + Some(base_block_group_id) + ); + let child_paths = Path::query( + conn, + "select * from paths where block_group_id = ?1 order by name;", + rusqlite::params![child_block_group_id], + ); + assert_eq!( + child_paths + .iter() + .map(|path| path.name.as_str()) + .collect::>(), + vec!["override"] + ); + assert!(SampleLineage::get_children(conn, "override").is_empty()); + } +}