From 7d72248de16dacb8e67ceae7f37620040ffd6012 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Sat, 8 Aug 2026 12:28:38 +0200 Subject: [PATCH 1/4] feat: add swap, liftover and convert subcommands - swap: UCSC chainSwap-compatible target/query swap, metadata preserved - liftover: BED3/6/12 mapping via rust-lapper, --min-match, --multiple, --unmapped - convert: chains + sequences to VCF/MAF/BAM, BAM behind new optional bam feature - lib: SequenceResolver::sequences(), validate_block_spans, validate_distinct_paths --- Cargo.lock | 321 +++++++- Cargo.toml | 11 +- src/cli/convert.rs | 1845 +++++++++++++++++++++++++++++++++++++++++++ src/cli/liftover.rs | 1676 +++++++++++++++++++++++++++++++++++++++ src/cli/mod.rs | 127 ++- src/cli/swap.rs | 438 ++++++++++ src/seq/sequence.rs | 30 + 7 files changed, 4441 insertions(+), 7 deletions(-) create mode 100644 src/cli/convert.rs create mode 100644 src/cli/liftover.rs create mode 100644 src/cli/swap.rs diff --git a/Cargo.lock b/Cargo.lock index e93f94a..7c1dd94 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -58,15 +58,83 @@ dependencies = [ "windows-sys", ] +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + +[[package]] +name = "bit-vec" +version = "0.10.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5727b15fa97d4f4fee0a3b7c3d550ed0269f54329207b86388de918604e31269" +dependencies = [ + "borsh", + "serde", +] + +[[package]] +name = "bitflags" +version = "2.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b588b76d00fde79687d7646a9b5bdf3cc0f655e0bbd080335a95d7e96f3587da" + +[[package]] +name = "borsh" +version = "1.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a88b7ea17d208c4193f2c1e6de3c35fe71f98c96982d5ced308bdcc749ff6e1f" +dependencies = [ + "borsh-derive", + "bytes", + "cfg_aliases", +] + +[[package]] +name = "borsh-derive" +version = "1.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d8f347189c62a579b8cd5f80714efa178f52e461dc2e6d701d264f5ff22e566c" +dependencies = [ + "once_cell", + "proc-macro-crate", + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "bstr" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1f7dc094d718f2e1c1559ad110e27eeaae14a5465d3d56dd6dbd793079fbd530" +dependencies = [ + "memchr", + "serde_core", +] + +[[package]] +name = "bytes" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc652a48c352aef3ea3aed32080501cf3ef6ed5da78602a020c991775b0aff04" + [[package]] name = "cfg-if" version = "1.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" +[[package]] +name = "cfg_aliases" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f079e83a288787bcd14a6aea84cee5c87a67c5a3e660c30f557a3d24761b3527" + [[package]] name = "chaintools" -version = "0.0.9" +version = "0.0.10" dependencies = [ "clap", "flate2", @@ -74,8 +142,12 @@ dependencies = [ "log", "memchr", "memmap2", + "noodles-bam", + "noodles-core", + "noodles-sam", "num_cpus", "rayon", + "rust-lapper", "simple_logger", "twobit", ] @@ -144,6 +216,15 @@ dependencies = [ "cfg-if", ] +[[package]] +name = "crossbeam-channel" +version = "0.5.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d85363c37faeca707aef026efa9f3b34d077bce547e48f770770625c6013679e" +dependencies = [ + "crossbeam-utils", +] + [[package]] name = "crossbeam-deque" version = "0.8.6" @@ -184,6 +265,12 @@ version = "1.15.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" +[[package]] +name = "equivalent" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f" + [[package]] name = "flate2" version = "1.1.5" @@ -204,6 +291,12 @@ dependencies = [ "memchr", ] +[[package]] +name = "hashbrown" +version = "0.17.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed5909b6e89a2db4456e54cd5f673791d7eca6732202bbf2a9cc504fe2f9b84a" + [[package]] name = "heck" version = "0.5.0" @@ -216,6 +309,16 @@ version = "0.5.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "fc0fef456e4baa96da950455cd02c081ca953b141298e41db3fc7e36b1da849c" +[[package]] +name = "indexmap" +version = "2.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d466e9454f08e4a911e14806c24e16fba1b4c121d1ea474396f396069cf949d9" +dependencies = [ + "equivalent", + "hashbrown", +] + [[package]] name = "is_terminal_polyfill" version = "1.70.2" @@ -228,6 +331,63 @@ version = "1.0.18" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682" +[[package]] +name = "lexical-core" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7d8d125a277f807e55a77304455eb7b1cb52f2b18c143b60e766c120bd64a594" +dependencies = [ + "lexical-parse-float", + "lexical-parse-integer", + "lexical-util", + "lexical-write-float", + "lexical-write-integer", +] + +[[package]] +name = "lexical-parse-float" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "52a9f232fbd6f550bc0137dcb5f99ab674071ac2d690ac69704593cb4abbea56" +dependencies = [ + "lexical-parse-integer", + "lexical-util", +] + +[[package]] +name = "lexical-parse-integer" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a7a039f8fb9c19c996cd7b2fcce303c1b2874fe1aca544edc85c4a5f8489b34" +dependencies = [ + "lexical-util", +] + +[[package]] +name = "lexical-util" +version = "1.0.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2604dd126bb14f13fb5d1bd6a66155079cb9fa655b37f875b3a742c705dbed17" + +[[package]] +name = "lexical-write-float" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "50c438c87c013188d415fbabbb1dceb44249ab81664efbd31b14ae55dabb6361" +dependencies = [ + "lexical-util", + "lexical-write-integer", +] + +[[package]] +name = "lexical-write-integer" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "409851a618475d2d5796377cad353802345cba92c867d9fbcde9cf4eac4e14df" +dependencies = [ + "lexical-util", +] + [[package]] name = "libc" version = "0.2.177" @@ -240,7 +400,7 @@ version = "0.5.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c10501e7805cee23da17c7790e59df2870c0d4043ec6d03f67d31e2b53e77415" dependencies = [ - "zlib-rs", + "zlib-rs 0.5.5", ] [[package]] @@ -274,12 +434,86 @@ dependencies = [ "simd-adler32", ] +[[package]] +name = "noodles-bam" +version = "0.94.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e018289a9da86771349be387224f6634ad56ec48eb16a3a29a7670def5d206c2" +dependencies = [ + "bstr", + "indexmap", + "memchr", + "noodles-bgzf", + "noodles-core", + "noodles-csi", + "noodles-sam", +] + +[[package]] +name = "noodles-bgzf" +version = "0.51.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "280d401d0358893ad706fb6094a18e7f9e8f3002630f4960889cb702b6e2714e" +dependencies = [ + "bytes", + "crossbeam-channel", + "rayon", + "zlib-rs 0.6.7", +] + +[[package]] +name = "noodles-core" +version = "0.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c8dbac7c5f9a7de9fe45590f198a09697df631cd13d2060b4742cc48144555b0" +dependencies = [ + "bstr", +] + +[[package]] +name = "noodles-csi" +version = "0.60.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32789eeba9c4001d50de853a4c6a8f0551d586a91d6919c0008dd462231feacb" +dependencies = [ + "bit-vec", + "bstr", + "indexmap", + "noodles-bgzf", + "noodles-core", +] + +[[package]] +name = "noodles-sam" +version = "0.89.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8063dd9e92825db46ae21813525589d41cff9ce2dcb2a0d25d99af86b274a1c5" +dependencies = [ + "bitflags", + "bstr", + "indexmap", + "lexical-core", + "memchr", + "noodles-bgzf", + "noodles-core", + "noodles-csi", +] + [[package]] name = "num-conv" version = "0.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c6673768db2d862beb9b39a78fdcb1a69439615d5794a1be50caa9bc92c81967" +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + [[package]] name = "num_cpus" version = "1.17.0" @@ -299,6 +533,12 @@ dependencies = [ "libc", ] +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + [[package]] name = "once_cell_polyfill" version = "1.70.2" @@ -311,6 +551,15 @@ version = "0.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "439ee305def115ba05938db6eb1644ff94165c5ab5e9420d1c1bcedbba909391" +[[package]] +name = "proc-macro-crate" +version = "3.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e67ba7e9b2b56446f1d419b1d807906278ffa1a658a8a5d8a39dcb1f5a78614f" +dependencies = [ + "toml_edit", +] + [[package]] name = "proc-macro2" version = "1.0.106" @@ -331,9 +580,9 @@ dependencies = [ [[package]] name = "rayon" -version = "1.11.0" +version = "1.12.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "368f01d005bf8fd9b1206fb6fa653e6c4a81ceb1466406b81792d87c5677a58f" +checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d" dependencies = [ "either", "rayon-core", @@ -349,6 +598,25 @@ dependencies = [ "crossbeam-utils", ] +[[package]] +name = "rust-lapper" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5aafed57460ca80eb88c33989af3a31348dbde8e68f86f4fe2107c1a5fcd55b8" +dependencies = [ + "num-traits", +] + +[[package]] +name = "serde" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9a8e94ea7f378bd32cbbd37198a4a91436180c5bb472411e48b5ec2e2124ae9e" +dependencies = [ + "serde_core", + "serde_derive", +] + [[package]] name = "serde_core" version = "1.0.228" @@ -437,6 +705,36 @@ dependencies = [ "time-core", ] +[[package]] +name = "toml_datetime" +version = "1.1.1+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3165f65f62e28e0115a00b2ebdd37eb6f3b641855f9d636d3cd4103767159ad7" +dependencies = [ + "serde_core", +] + +[[package]] +name = "toml_edit" +version = "0.25.13+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6975367e4d2ef766d86af01ffad14b622fecc8d4357a998fbc4deb6e9bacaf9b" +dependencies = [ + "indexmap", + "toml_datetime", + "toml_parser", + "winnow", +] + +[[package]] +name = "toml_parser" +version = "1.1.3+spec-1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d38ac1cf9b95face32296c0a3ede1fdc270627c9d9c02a7274dd6d960dc4d56" +dependencies = [ + "winnow", +] + [[package]] name = "twobit" version = "0.2.2" @@ -470,8 +768,23 @@ dependencies = [ "windows-link", ] +[[package]] +name = "winnow" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23b97319f7b8343df12cc98938e5c3eb436064524c8d2b4e30a1d3a36eecdf81" +dependencies = [ + "memchr", +] + [[package]] name = "zlib-rs" version = "0.5.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "40990edd51aae2c2b6907af74ffb635029d5788228222c4bb811e9351c0caad3" + +[[package]] +name = "zlib-rs" +version = "0.6.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "34b31d188d9d685a4f9c7b46d6e36631b07058d2cfe190267adce54dc230bf12" diff --git a/Cargo.toml b/Cargo.toml index 857fb47..3c8bab2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "chaintools" -version = "0.0.9" +version = "0.0.10" edition = "2024" authors = ["alejandrogzi "] license = "MIT" @@ -13,13 +13,16 @@ default-run = "chaintools" [features] default = ["mmap"] -cli = ["dep:clap", "dep:genepred", "dep:log", "dep:simple_logger", "sequence", "write", "dep:num_cpus", "parallel"] +cli = ["dep:clap", "dep:genepred", "dep:log", "dep:simple_logger", "sequence", "write", "dep:num_cpus", "dep:rust-lapper", "parallel"] mmap = ["memmap2"] gzip = ["dep:flate2", "flate2/zlib-rs"] parallel = ["rayon"] index = [] write = [] sequence = ["dep:twobit"] +# BAM output for `convert --to bam`. Optional so library and non-BAM CLI builds +# do not pull the BAM writer; VCF and MAF conversion never touch it. +bam = ["dep:noodles-bam", "dep:noodles-core", "dep:noodles-sam"] [dependencies] clap = { version = "4", features = ["derive"], optional = true } @@ -27,7 +30,11 @@ flate2 = { version = "1", optional = true, default-features = false } genepred = { version = "0.0.16", optional = true } memchr = "2" memmap2 = { version = "0.9", optional = true } +noodles-bam = { version = "0.94", optional = true } +noodles-core = { version = "0.20", optional = true } +noodles-sam = { version = "0.89", optional = true } rayon = { version = "1", optional = true } +rust-lapper = { version = "1", optional = true } log = { version = "0.4.29", optional = true } simple_logger = { version = "5.2.0", optional = true, features = ["stderr"] } twobit = { version = "0.2.2", optional = true } diff --git a/src/cli/convert.rs b/src/cli/convert.rs new file mode 100644 index 0000000..4085df4 --- /dev/null +++ b/src/cli/convert.rs @@ -0,0 +1,1845 @@ +// Copyright (c) 2026 Alejandro Gonzales-Irribarren +// Distributed under the terms of the Apache License, Version 2.0. + +use std::fs::File; +use std::io::{BufRead, BufWriter, Write}; +use std::path::PathBuf; + +use chaintools::seq::revcomp::reverse_complement_in_place; +use chaintools::seq::sequence::{SequenceCache, SequenceResolver}; +use chaintools::{ChainError, OwnedChain, Strand, StreamingReader}; +use clap::{Args, ValueEnum}; +#[cfg(feature = "gzip")] +use flate2::{Compression, write::GzEncoder}; + +use super::CliError; + +const OUTPUT_BUFFER_CAPACITY: usize = 1024 * 1024; + +/// Command-line arguments for the convert subcommand. +/// +/// Converts chain alignments into a sequence-aware format. A chain records +/// alignment geometry but no bases, so both genomes are required. +/// +/// # Examples +/// +/// ```bash +/// chaintools convert --to vcf -r reference.2bit -q query.2bit -c input.chain -o out.vcf +/// ``` +#[derive(Debug, Args)] +pub struct ConvertArgs { + #[arg( + long = "to", + value_name = "FORMAT", + value_enum, + help = "Output format. The generic --output option carries the result." + )] + to: ConvertTarget, + + #[arg( + short = 'r', + long = "reference", + value_name = "PATH", + help = "Path to the reference (target) sequence file (.2bit, .fa, .fasta, .fna, and gzip variants)." + )] + reference: PathBuf, + + #[arg( + short = 'q', + long = "query", + value_name = "PATH", + help = "Path to the query sequence file (.2bit, .fa, .fasta, .fna, and gzip variants)." + )] + query: PathBuf, + + #[arg( + short = 'c', + long = "chain", + value_name = "PATH", + help = "Path to the input .chain file. If not provided, chain data is read from standard input." + )] + chain: Option, + + #[arg( + short = 'o', + long = "output", + value_name = "PATH", + help = "Path for the converted output. If not provided, output is written to standard output." + )] + output: Option, + + #[arg( + short = 'G', + long = "gzip", + help = "Compress convert output with gzip. Requires the `gzip` feature." + )] + gzip: bool, +} + +/// Conversion targets supported by `chaintools convert`. +#[derive(Debug, Clone, Copy, PartialEq, Eq, ValueEnum)] +pub enum ConvertTarget { + Vcf, + Maf, + Bam, +} + +impl std::fmt::Display for ConvertTarget { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + ConvertTarget::Vcf => f.write_str("vcf"), + ConvertTarget::Maf => f.write_str("maf"), + ConvertTarget::Bam => f.write_str("bam"), + } + } +} + +/// Runs the convert subcommand. +/// +/// # Arguments +/// +/// * `args` - Convert arguments +/// * `stdin` - Input stream (used if no --chain provided) +/// * `stdout` - Output stream (used if no --output provided) +/// * `_stderr` - Error/logging output +/// +/// # Output +/// +/// Returns `Ok(())` on success or `Err(CliError)` on failure +pub fn run( + args: ConvertArgs, + stdin: &mut R, + stdout: &mut W, + _stderr: &mut E, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, + E: Write, +{ + validate_output_args(&args)?; + super::ensure_inputs_exist( + &[("reference", &args.reference), ("query", &args.query)], + &[("input chain", args.chain.as_deref())], + )?; + if let Some(output) = &args.output { + super::validate_distinct_paths("--output", output, args.chain.as_deref())?; + } + + let input_desc = args + .chain + .as_deref() + .map_or_else(|| "".to_owned(), |path| path.display().to_string()); + let output_desc = args + .output + .as_deref() + .map_or_else(|| "".to_owned(), |path| path.display().to_string()); + log::info!( + "convert: to={}, reference={}, query={}, input={input_desc}, output={output_desc}", + args.to, + args.reference.display(), + args.query.display() + ); + + if let Some(path) = &args.output { + let file = File::create(path)?; + let writer = BufWriter::with_capacity(OUTPUT_BUFFER_CAPACITY, file); + if args.gzip { + run_gzip_output(&args, stdin, writer)?; + } else { + let mut writer = writer; + run_to_writer(&args, stdin, &mut writer)?; + writer.flush()?; + } + } else if args.gzip { + run_gzip_output(&args, stdin, stdout)?; + } else { + run_to_writer(&args, stdin, stdout)?; + stdout.flush()?; + } + + Ok(()) +} + +#[cfg(feature = "gzip")] +fn run_gzip_output(args: &ConvertArgs, stdin: &mut R, writer: W) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + let mut writer = GzEncoder::new(writer, Compression::fast()); + run_to_writer(args, stdin, &mut writer)?; + writer.try_finish()?; + writer.get_mut().flush()?; + Ok(()) +} + +#[cfg(not(feature = "gzip"))] +fn run_gzip_output(_args: &ConvertArgs, _stdin: &mut R, _writer: W) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + Err(CliError::Message( + "--gzip requires chaintools to be built with the `gzip` feature".to_owned(), + )) +} + +fn validate_output_args(args: &ConvertArgs) -> Result<(), CliError> { + if args.gzip && args.to == ConvertTarget::Bam { + // BAM is block-compressed by its own writer; wrapping it in gzip would + // produce a file no BAM reader accepts. + return Err(CliError::Message( + "--gzip is not supported with --to bam: BAM output is already block-compressed" + .to_owned(), + )); + } + + #[cfg(not(feature = "gzip"))] + if args.gzip { + return Err(CliError::Message( + "--gzip requires chaintools to be built with the `gzip` feature".to_owned(), + )); + } + Ok(()) +} + +fn run_to_writer(args: &ConvertArgs, stdin: &mut R, writer: &mut W) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + let sequences = Sequences { + reference: SequenceResolver::new(&args.reference)?, + query: SequenceResolver::new(&args.query)?, + cache: SequenceCache::default(), + }; + + match args.to { + ConvertTarget::Vcf => convert_vcf(args, stdin, writer, sequences), + ConvertTarget::Maf => convert_maf(args, stdin, writer, sequences), + ConvertTarget::Bam => convert_bam(args, stdin, writer, sequences), + } +} + +/// Streams chains and writes BAM assembly-alignment records. +/// +/// `@SQ` entries are taken from the reference genome rather than from a pre-scan +/// of the chains: a BAM header must be complete before the first record, and the +/// chain stream may be an unrewindable pipe. Every chain's declared target size +/// is checked against that header, so contradictory sizes are rejected instead +/// of silently accepted. +#[cfg(feature = "bam")] +fn convert_bam( + args: &ConvertArgs, + stdin: &mut R, + writer: &mut W, + mut sequences: Sequences, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + use noodles_bam as bam; + use noodles_sam::alignment::io::Write as _; + + let header = bam_header(&sequences.reference)?; + let mut bam_writer = bam::io::Writer::new(writer); + bam_writer.write_header(&header)?; + + let mut records = 0u64; + let chains = for_each_chain(args, stdin, |chain, offset| { + let reference_id = reference_sequence_id(&header, chain, offset)?; + let (reference_seq, query_seq) = sequences.fetch_chain(chain)?; + + // One chain's segments are buffered so the QNAME suffix can say whether + // the chain was split at all; nothing beyond one chain is held. + let mut segments = Vec::new(); + walk_segments(chain, offset, &reference_seq, &query_seq, |segment| { + segments.push(segment); + Ok(1) + })?; + + let split = segments.len() > 1; + for (index, segment) in segments.iter().enumerate() { + let record = bam_record(chain, reference_id, segment, split.then_some(index + 1))?; + bam_writer.write_alignment_record(&header, &record)?; + records += 1; + } + Ok(()) + })?; + + bam_writer.try_finish()?; + super::log_summary("convert", &[("chains", chains), ("records", records)]); + Ok(()) +} + +/// Builds a `@HD`/`@SQ`-only header from the reference genome's sequences. +#[cfg(feature = "bam")] +fn bam_header(reference: &SequenceResolver) -> Result { + use std::num::NonZero; + + use noodles_sam::{ + self as sam, + header::record::value::{Map, map::ReferenceSequence}, + }; + + let mut builder = sam::Header::builder().set_header(Default::default()); + for (name, length) in reference.sequences()? { + let length = NonZero::new(length as usize).ok_or_else(|| { + CliError::Message(format!( + "reference sequence {} is empty and cannot be a BAM @SQ entry", + String::from_utf8_lossy(name) + )) + })?; + builder = + builder.add_reference_sequence(name.to_vec(), Map::::new(length)); + } + Ok(builder.build()) +} + +/// Resolves a chain's target contig to its `@SQ` index, checking its length. +#[cfg(feature = "bam")] +fn reference_sequence_id( + header: &noodles_sam::Header, + chain: &OwnedChain, + offset: usize, +) -> Result { + let (index, _, reference_sequence) = header + .reference_sequences() + .get_full(chain.reference_name.as_slice()) + .ok_or_else(|| { + CliError::Chain(ChainError::MissingSequence { + name: String::from_utf8_lossy(&chain.reference_name) + .into_owned() + .into(), + }) + })?; + + if usize::from(reference_sequence.length()) != chain.reference_size as usize { + return Err(CliError::Chain(format_error( + offset, + format!( + "chain declares target {} as {} bases but the reference has {}", + String::from_utf8_lossy(&chain.reference_name), + chain.reference_size, + reference_sequence.length() + ), + ))); + } + Ok(index) +} + +/// Builds one BAM record from an aligned segment. +/// +/// The record is an assembly-alignment segment, not a sequencing read: MAPQ is +/// 0, qualities are absent, and no read-group or platform metadata is invented. +/// A chain split into several segments yields several records with +/// deterministically suffixed names, so no supplementary bookkeeping is needed. +#[cfg(feature = "bam")] +fn bam_record( + chain: &OwnedChain, + reference_id: usize, + segment: &AlignedSegment, + segment_index: Option, +) -> Result { + use noodles_core::Position; + use noodles_sam::{ + alignment::RecordBuf, + alignment::record::{Flags, MappingQuality}, + }; + + let mut name = chain.query_name.clone(); + write!(&mut name, ".{}", chain.id).expect("writing to Vec cannot fail"); + if let Some(index) = segment_index { + write!(&mut name, ".{index}").expect("writing to Vec cannot fail"); + } + + let start = Position::new(segment.reference_start as usize + 1).ok_or_else(|| { + CliError::Message("BAM alignment start overflows a 1-based position".to_owned()) + })?; + let flags = match chain.query_strand { + Strand::Plus => Flags::empty(), + Strand::Minus => Flags::REVERSE_COMPLEMENTED, + }; + // BAM encodes SEQ in 4 bits per base, which has no lowercase form. + let bases: Vec = segment + .query_text + .iter() + .filter(|base| **base != b'-') + .map(u8::to_ascii_uppercase) + .collect(); + + Ok(RecordBuf::builder() + .set_name(name) + .set_flags(flags) + .set_reference_sequence_id(reference_id) + .set_alignment_start(start) + .set_mapping_quality(MappingQuality::new(0).expect("0 is a valid mapping quality")) + .set_cigar(segment_cigar(segment).into()) + .set_sequence(bases.into()) + .build()) +} + +/// Derives an `=`/`X`/`I`/`D` CIGAR from a segment's two aligned rows. +/// +/// The exact bases are known, so matches and mismatches are distinguished +/// instead of being collapsed into `M`. One-sided chain gaps appear as `-` in one +/// row and become `I`/`D`; no secondary alignment is computed. +#[cfg(feature = "bam")] +fn segment_cigar(segment: &AlignedSegment) -> Vec { + use noodles_sam::alignment::record::cigar::{Op, op::Kind}; + + let mut ops: Vec = Vec::new(); + for (&reference_base, &query_base) in segment.reference_text.iter().zip(&segment.query_text) { + let kind = if reference_base == b'-' { + Kind::Insertion + } else if query_base == b'-' { + Kind::Deletion + } else if reference_base.eq_ignore_ascii_case(&query_base) { + Kind::SequenceMatch + } else { + Kind::SequenceMismatch + }; + + match ops.last_mut() { + Some(last) if last.kind() == kind => *last = Op::new(kind, last.len() + 1), + _ => ops.push(Op::new(kind, 1)), + } + } + ops +} + +#[cfg(not(feature = "bam"))] +fn convert_bam( + _args: &ConvertArgs, + _stdin: &mut R, + _writer: &mut W, + _sequences: Sequences, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + Err(CliError::Message( + "--to bam requires chaintools to be built with the `bam` feature".to_owned(), + )) +} + +/// Reference and query sequence access for one conversion run. +struct Sequences { + reference: SequenceResolver, + query: SequenceResolver, + cache: SequenceCache, +} + +impl Sequences { + /// Fetches a chain's target span and its strand-oriented query span. + /// + /// The target is always read on `+`; a `-` query is fetched from its forward + /// coordinates and reverse-complemented, so the returned query bases run in + /// the same direction as the chain's own query coordinates. This mirrors + /// [`chaintools::seq::score::chainscore::ChainScorer::score_chain`]. + fn fetch_chain(&mut self, chain: &OwnedChain) -> Result<(Vec, Vec), ChainError> { + let reference_len = span(chain.reference_start, chain.reference_end, "target")?; + let query_len = span(chain.query_start, chain.query_end, "query")?; + + let reference_seq = self.reference.fetch( + &mut self.cache, + &chain.reference_name, + chain.reference_start, + reference_len, + )?; + let query_seq = match chain.query_strand { + Strand::Plus => self.query.fetch( + &mut self.cache, + &chain.query_name, + chain.query_start, + query_len, + )?, + Strand::Minus => { + let start = chain + .query_size + .checked_sub(chain.query_end) + .ok_or_else(|| convert_error("query minus-strand fetch underflows"))?; + let mut seq = + self.query + .fetch(&mut self.cache, &chain.query_name, start, query_len)?; + reverse_complement_in_place(&mut seq); + seq + } + }; + + Ok((reference_seq, query_seq)) + } +} + +/// Streams chains and writes one VCF describing the query as the ALT genome. +fn convert_vcf( + args: &ConvertArgs, + stdin: &mut R, + writer: &mut W, + mut sequences: Sequences, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + write_vcf_header(writer)?; + + let mut records = 0u64; + let chains = for_each_chain(args, stdin, |chain, offset| { + let (reference_seq, query_seq) = sequences.fetch_chain(chain)?; + records += write_chain_variants(writer, chain, offset, &reference_seq, &query_seq)?; + Ok(()) + })?; + + super::log_summary("convert", &[("chains", chains), ("records", records)]); + Ok(()) +} + +/// Streams chains from the file input or standard input, returning the count. +/// +/// One chain is held at a time: it is read, converted, and dropped before the +/// next header is parsed. +fn for_each_chain(args: &ConvertArgs, stdin: &mut R, mut emit: F) -> Result +where + R: BufRead, + F: FnMut(&OwnedChain, usize) -> Result<(), CliError>, +{ + if let Some(path) = &args.chain { + drive_chains(&mut StreamingReader::from_path(path)?, &mut emit) + } else { + drive_chains(&mut StreamingReader::new(stdin), &mut emit) + } +} + +fn drive_chains(reader: &mut StreamingReader, emit: &mut F) -> Result +where + R: BufRead, + F: FnMut(&OwnedChain, usize) -> Result<(), CliError>, +{ + let mut chains = 0u64; + while let Some(header) = reader.next_header()? { + let offset = header.offset; + let blocks = reader.read_blocks(offset)?; + emit(&header.into_chain(blocks), offset)?; + chains += 1; + } + Ok(chains) +} + +fn write_vcf_header(writer: &mut W) -> Result<(), CliError> { + writer.write_all(b"##fileformat=VCFv4.2\n")?; + writeln!(writer, "##source=chaintools {}", env!("CARGO_PKG_VERSION"))?; + writer.write_all( + b"##INFO=\n\ + ##INFO=\n\ + ##INFO=\n\ + ##INFO=\n\ + #CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n", + )?; + Ok(()) +} + +/// Emits the VCF records for one chain, returning how many were written. +/// +/// Walks the chain block by block against the already-fetched target and +/// oriented query spans: aligned bases yield one record per mismatching base, +/// and each gap between blocks yields one indel or replacement record. Records +/// follow chain traversal order and are never buffered for sorting. +fn write_chain_variants( + writer: &mut W, + chain: &OwnedChain, + offset: usize, + reference_seq: &[u8], + query_seq: &[u8], +) -> Result { + let mut records = 0u64; + let mut reference_pos = chain.reference_start; + let mut query_pos = chain.query_start; + let last = chain.blocks.len().saturating_sub(1); + + for (index, block) in chain.blocks.iter().enumerate() { + let size = block.size as usize; + let reference_off = (reference_pos - chain.reference_start) as usize; + let query_off = (query_pos - chain.query_start) as usize; + let reference_block = slice(reference_seq, reference_off, size, offset, "target")?; + let query_block = slice(query_seq, query_off, size, offset, "query")?; + + for (step, (&reference_base, &query_base)) in + reference_block.iter().zip(query_block).enumerate() + { + let reference_base = reference_base.to_ascii_uppercase(); + let query_base = query_base.to_ascii_uppercase(); + // Soft-masked bases differ only in case, and a substitution against + // an ambiguous base is not a callable variant. + if reference_base == query_base || !is_acgt(reference_base) || !is_acgt(query_base) { + continue; + } + let step = step as u32; + write_vcf_record( + writer, + chain, + u64::from(reference_pos + step) + 1, + &[reference_base], + &[query_base], + query_pos + step, + )?; + records += 1; + } + + if index < last { + records += write_gap_variant( + writer, + chain, + offset, + GapContext { + reference_seq, + query_seq, + reference_off: reference_off + size, + query_off: query_off + size, + reference_pos: reference_pos + block.size, + query_pos: query_pos + block.size, + reference_gap: block.gap_reference as usize, + query_gap: block.gap_query as usize, + }, + )?; + } + + reference_pos += block.size + block.gap_reference; + query_pos += block.size + block.gap_query; + } + + verify_walk_end(chain, offset, reference_pos, query_pos)?; + Ok(records) +} + +/// Checks that a completed block walk landed exactly on both header ends. +fn verify_walk_end( + chain: &OwnedChain, + offset: usize, + reference_pos: u32, + query_pos: u32, +) -> Result<(), CliError> { + if reference_pos != chain.reference_end || query_pos != chain.query_end { + return Err(CliError::Chain(format_error( + offset, + "chain blocks do not sum to the header span", + ))); + } + Ok(()) +} + +/// One inter-block gap, resolved against the fetched chain spans. +struct GapContext<'a> { + reference_seq: &'a [u8], + query_seq: &'a [u8], + reference_off: usize, + query_off: usize, + reference_pos: u32, + query_pos: u32, + reference_gap: usize, + query_gap: usize, +} + +/// Writes the record for one inter-block gap, returning how many were written. +/// +/// A target-only gap is a deletion, a query-only gap is an insertion, and a gap +/// on both sides is emitted as a whole-interval replacement: a chain does not +/// define a base-level alignment inside a dual-sided gap, so none is invented. +fn write_gap_variant( + writer: &mut W, + chain: &OwnedChain, + offset: usize, + gap: GapContext<'_>, +) -> Result { + if gap.reference_gap == 0 && gap.query_gap == 0 { + return Ok(0); + } + + let deleted = slice( + gap.reference_seq, + gap.reference_off, + gap.reference_gap, + offset, + "target gap", + )? + .to_ascii_uppercase(); + let inserted = slice( + gap.query_seq, + gap.query_off, + gap.query_gap, + offset, + "query gap", + )? + .to_ascii_uppercase(); + + let (position, mut reference_allele, mut alternate_allele) = + if gap.reference_gap > 0 && gap.query_gap > 0 { + // Both alleles are non-empty, so no anchor base is needed. + (u64::from(gap.reference_pos) + 1, deleted, inserted) + } else if gap.reference_off > 0 { + // Anchor on the preceding target base (the last base of the block + // just walked), which is always inside the fetched span. + let anchor = gap.reference_seq[gap.reference_off - 1].to_ascii_uppercase(); + let mut reference_allele = vec![anchor]; + let mut alternate_allele = vec![anchor]; + reference_allele.extend_from_slice(&deleted); + alternate_allele.extend_from_slice(&inserted); + ( + u64::from(gap.reference_pos), + reference_allele, + alternate_allele, + ) + } else { + // The gap starts at the first base of the chain's target span, so + // there is no preceding base to anchor on: use the following one + // rather than underflowing the coordinate. + let anchor = slice( + gap.reference_seq, + gap.reference_off + gap.reference_gap, + 1, + offset, + "target anchor", + )?[0] + .to_ascii_uppercase(); + let mut reference_allele = deleted; + let mut alternate_allele = inserted; + reference_allele.push(anchor); + alternate_allele.push(anchor); + ( + u64::from(gap.reference_pos) + 1, + reference_allele, + alternate_allele, + ) + }; + + let position = trim_alleles(position, &mut reference_allele, &mut alternate_allele); + if reference_allele == alternate_allele { + return Ok(0); + } + + write_vcf_record( + writer, + chain, + position, + &reference_allele, + &alternate_allele, + gap.query_pos, + )?; + Ok(1) +} + +/// Trims redundant shared suffix then prefix bases, keeping both alleles non-empty. +/// +/// Returns the (possibly advanced) 1-based reference position. This is the +/// cheap, always-safe part of normalization; repeat-aware left alignment is +/// deliberately left to dedicated VCF tooling such as `bcftools norm`. +fn trim_alleles(position: u64, reference: &mut Vec, alternate: &mut Vec) -> u64 { + while reference.len() > 1 && alternate.len() > 1 && reference.last() == alternate.last() { + reference.pop(); + alternate.pop(); + } + + let mut trimmed = 0usize; + while reference.len() > trimmed + 1 + && alternate.len() > trimmed + 1 + && reference[trimmed] == alternate[trimmed] + { + trimmed += 1; + } + reference.drain(..trimmed); + alternate.drain(..trimmed); + position + trimmed as u64 +} + +/// Writes one VCF data line. `position` is 1-based; `query_pos` is 0-based. +fn write_vcf_record( + writer: &mut W, + chain: &OwnedChain, + position: u64, + reference_allele: &[u8], + alternate_allele: &[u8], + query_pos: u32, +) -> Result<(), CliError> { + writer.write_all(&chain.reference_name)?; + write!(writer, "\t{position}\t.\t")?; + writer.write_all(reference_allele)?; + writer.write_all(b"\t")?; + writer.write_all(alternate_allele)?; + write!(writer, "\t.\tPASS\tCHAIN_ID={};QUERY_CHROM=", chain.id)?; + writer.write_all(&chain.query_name)?; + writeln!( + writer, + ";QUERY_POS={};STRAND={}", + u64::from(query_pos) + 1, + strand_symbol(chain.query_strand) + )?; + Ok(()) +} + +/// Streams chains and writes a pairwise MAF. +fn convert_maf( + args: &ConvertArgs, + stdin: &mut R, + writer: &mut W, + mut sequences: Sequences, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + writeln!( + writer, + "##maf version=1 scoring=chain\n# chaintools {}\n", + env!("CARGO_PKG_VERSION") + )?; + + let mut stanzas = 0u64; + let chains = for_each_chain(args, stdin, |chain, offset| { + let (reference_seq, query_seq) = sequences.fetch_chain(chain)?; + stanzas += walk_segments(chain, offset, &reference_seq, &query_seq, |segment| { + flush_stanza(writer, chain, &segment) + })?; + Ok(()) + })?; + + super::log_summary("convert", &[("chains", chains), ("records", stanzas)]); + Ok(()) +} + +/// One gap-aware pairwise alignment segment of a chain. +/// +/// This is the unit both MAF and BAM emit: a run of the chain that has a +/// continuous base-level alignment, ending wherever the chain stops defining one. +/// Text is accumulated verbatim from the fetched sequences (soft-masking +/// included) so that stripping `-` reproduces the declared source interval +/// exactly. `start`/`size` are in each row's declared orientation: for a `-` +/// query that is the chain's own query coordinate space, which is already what +/// MAF asks for. +struct AlignedSegment { + reference_start: u32, + query_start: u32, + reference_size: u32, + query_size: u32, + reference_text: Vec, + query_text: Vec, +} + +impl AlignedSegment { + fn new(reference_start: u32, query_start: u32) -> Self { + AlignedSegment { + reference_start, + query_start, + reference_size: 0, + query_size: 0, + reference_text: Vec::new(), + query_text: Vec::new(), + } + } + + fn push_aligned(&mut self, reference: &[u8], query: &[u8]) { + self.reference_text.extend_from_slice(reference); + self.query_text.extend_from_slice(query); + self.reference_size += reference.len() as u32; + self.query_size += query.len() as u32; + } + + /// Adds unaligned reference bases against `-` on the query row. + fn push_reference_gap(&mut self, reference: &[u8]) { + self.reference_text.extend_from_slice(reference); + self.query_text + .extend(std::iter::repeat_n(b'-', reference.len())); + self.reference_size += reference.len() as u32; + } + + /// Adds unaligned query bases against `-` on the reference row. + fn push_query_gap(&mut self, query: &[u8]) { + self.query_text.extend_from_slice(query); + self.reference_text + .extend(std::iter::repeat_n(b'-', query.len())); + self.query_size += query.len() as u32; + } + + fn is_empty(&self) -> bool { + self.reference_text.is_empty() && self.query_text.is_empty() + } +} + +/// Walks a chain's blocks and hands each aligned segment to `emit`. +/// +/// Aligned blocks and one-sided gaps extend the current segment; a gap on both +/// sides ends it, because a chain records no alignment path through such a gap. +/// Returns the sum of the counts `emit` reports. MAF and BAM split identically, +/// so they share this walk and differ only in what they do with a segment. +fn walk_segments( + chain: &OwnedChain, + offset: usize, + reference_seq: &[u8], + query_seq: &[u8], + mut emit: F, +) -> Result +where + F: FnMut(AlignedSegment) -> Result, +{ + let mut emitted = 0u64; + let mut segment = AlignedSegment::new(chain.reference_start, chain.query_start); + let mut reference_pos = chain.reference_start; + let mut query_pos = chain.query_start; + let last = chain.blocks.len().saturating_sub(1); + + for (index, block) in chain.blocks.iter().enumerate() { + let size = block.size as usize; + let reference_off = (reference_pos - chain.reference_start) as usize; + let query_off = (query_pos - chain.query_start) as usize; + segment.push_aligned( + slice(reference_seq, reference_off, size, offset, "target")?, + slice(query_seq, query_off, size, offset, "query")?, + ); + + let reference_gap = block.gap_reference as usize; + let query_gap = block.gap_query as usize; + reference_pos += block.size + block.gap_reference; + query_pos += block.size + block.gap_query; + + if index == last { + break; + } + if reference_gap > 0 && query_gap > 0 { + // ponytail: chain does not define an alignment path through a + // dual-sided gap; split here rather than inventing one. Upgrade path + // is a local aligner over the two unaligned intervals, which is a + // different tool's job. + let finished = + std::mem::replace(&mut segment, AlignedSegment::new(reference_pos, query_pos)); + emitted += emit_non_empty(&mut emit, finished)?; + } else if reference_gap > 0 { + segment.push_reference_gap(slice( + reference_seq, + reference_off + size, + reference_gap, + offset, + "target gap", + )?); + } else if query_gap > 0 { + segment.push_query_gap(slice( + query_seq, + query_off + size, + query_gap, + offset, + "query gap", + )?); + } + } + + verify_walk_end(chain, offset, reference_pos, query_pos)?; + emitted += emit_non_empty(&mut emit, segment)?; + Ok(emitted) +} + +/// Hands a segment to `emit` unless it carries no alignment columns. +fn emit_non_empty(emit: &mut F, segment: AlignedSegment) -> Result +where + F: FnMut(AlignedSegment) -> Result, +{ + if segment.is_empty() { + return Ok(0); + } + emit(segment) +} + +/// Writes one MAF stanza. +fn flush_stanza( + writer: &mut W, + chain: &OwnedChain, + stanza: &AlignedSegment, +) -> Result { + writeln!(writer, "a score={} chain_id={}", chain.score, chain.id)?; + write_maf_row( + writer, + &chain.reference_name, + stanza.reference_start, + stanza.reference_size, + chain.reference_strand, + chain.reference_size, + &stanza.reference_text, + )?; + write_maf_row( + writer, + &chain.query_name, + stanza.query_start, + stanza.query_size, + chain.query_strand, + chain.query_size, + &stanza.query_text, + )?; + writer.write_all(b"\n")?; + Ok(1) +} + +fn write_maf_row( + writer: &mut W, + src: &[u8], + start: u32, + size: u32, + strand: Strand, + src_size: u32, + text: &[u8], +) -> Result<(), CliError> { + writer.write_all(b"s ")?; + writer.write_all(src)?; + write!( + writer, + " {start} {size} {} {src_size} ", + strand_symbol(strand) + )?; + writer.write_all(text)?; + writer.write_all(b"\n")?; + Ok(()) +} + +fn strand_symbol(strand: Strand) -> char { + match strand { + Strand::Plus => '+', + Strand::Minus => '-', + } +} + +fn is_acgt(base: u8) -> bool { + matches!(base, b'A' | b'C' | b'G' | b'T') +} + +/// Borrows `length` bases at `offset`, erroring instead of panicking or wrapping. +fn slice<'a>( + sequence: &'a [u8], + offset: usize, + length: usize, + chain_offset: usize, + label: &str, +) -> Result<&'a [u8], CliError> { + let end = offset.checked_add(length).ok_or_else(|| { + CliError::Chain(format_error( + chain_offset, + format!("{label} range overflows"), + )) + })?; + sequence.get(offset..end).ok_or_else(|| { + CliError::Chain(format_error( + chain_offset, + format!("{label} range exceeds the fetched chain sequence"), + )) + }) +} + +fn span(start: u32, end: u32, label: &str) -> Result { + end.checked_sub(start) + .ok_or_else(|| convert_error(format!("{label} span is inverted"))) +} + +fn format_error(offset: usize, message: impl Into) -> ChainError { + ChainError::Format { + offset, + msg: message.into().into(), + } +} + +fn convert_error(message: impl Into) -> ChainError { + ChainError::Unsupported { + msg: message.into().into(), + } +} + +#[cfg(test)] +mod tests { + use super::*; + use clap::Parser; + use std::io::Cursor; + use std::path::Path; + use std::sync::atomic::{AtomicUsize, Ordering}; + + #[derive(Debug, Parser)] + struct ConvertHarness { + #[command(flatten)] + args: ConvertArgs, + } + + static NEXT_TEMP_ID: AtomicUsize = AtomicUsize::new(0); + + /// Two tiny FASTA genomes in a self-cleaning temp directory. + struct Genomes { + dir: PathBuf, + reference: PathBuf, + query: PathBuf, + } + + impl Drop for Genomes { + fn drop(&mut self) { + let _ = std::fs::remove_dir_all(&self.dir); + } + } + + /// Writes `reference` and `query` as single-record FASTA genomes. + /// + /// The reference record is named `chr1` and the query record `qry1`. + fn genomes(reference: &str, query: &str) -> Genomes { + let id = NEXT_TEMP_ID.fetch_add(1, Ordering::Relaxed); + let dir = std::env::temp_dir().join(format!( + "chaintools-convert-test-{}-{id}", + std::process::id() + )); + std::fs::create_dir_all(&dir).expect("create temp dir"); + let reference_path = dir.join("reference.fa"); + let query_path = dir.join("query.fa"); + std::fs::write(&reference_path, format!(">chr1\n{reference}\n")).expect("write reference"); + std::fs::write(&query_path, format!(">qry1\n{query}\n")).expect("write query"); + Genomes { + dir, + reference: reference_path, + query: query_path, + } + } + + fn convert_bytes(genomes: &Genomes, to: ConvertTarget, chain: &str) -> Vec { + let mut stdin = Cursor::new(chain.as_bytes()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + run( + ConvertArgs { + to, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: None, + output: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect("convert run"); + stdout + } + + fn convert(genomes: &Genomes, to: ConvertTarget, chain: &str) -> String { + String::from_utf8(convert_bytes(genomes, to, chain)).expect("utf8 output") + } + + /// Runs a VCF conversion and returns only the data lines. + fn vcf_records(genomes: &Genomes, chain: &str) -> Vec { + convert(genomes, ConvertTarget::Vcf, chain) + .lines() + .filter(|line| !line.starts_with('#')) + .map(str::to_owned) + .collect() + } + + /// Asserts every record's REF equals the reference genome slice at POS. + /// + /// This is the high-value VCF invariant: a record whose REF does not match + /// the genome it claims to describe is unusable downstream. + fn assert_ref_matches_genome(reference: &str, records: &[String]) { + for record in records { + let fields: Vec<&str> = record.split('\t').collect(); + let position: usize = fields[1].parse().expect("numeric POS"); + let reference_allele = fields[3]; + let start = position - 1; + assert_eq!( + reference[start..start + reference_allele.len()].to_ascii_uppercase(), + reference_allele, + "REF mismatch in {record}" + ); + } + } + + #[test] + fn parses_minimal_args() { + let cli = ConvertHarness::try_parse_from([ + "chaintools", + "--to", + "vcf", + "--reference", + "reference.2bit", + "--query", + "query.2bit", + ]) + .expect("convert arguments should parse"); + assert_eq!(cli.args.to, ConvertTarget::Vcf); + assert!(cli.args.chain.is_none()); + assert!(cli.args.output.is_none()); + } + + #[test] + fn requires_both_genomes() { + for missing in ["--reference", "--query"] { + let args: Vec<&str> = ["chaintools", "--to", "vcf", missing, "genome.2bit"].to_vec(); + let err = ConvertHarness::try_parse_from(args) + .expect_err("a single genome should not be enough"); + assert_eq!(err.kind(), clap::error::ErrorKind::MissingRequiredArgument); + } + } + + #[test] + fn vcf_header_is_emitted_and_identical_alignment_has_no_records() { + let genomes = genomes("ACGTACGT", "ACGTACGT"); + let output = convert( + &genomes, + ConvertTarget::Vcf, + "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n", + ); + assert!(output.starts_with("##fileformat=VCFv4.2\n")); + assert!(output.contains("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n")); + assert!(!output.lines().any(|line| !line.starts_with('#'))); + } + + #[test] + fn vcf_emits_one_record_per_substitution() { + let genomes = genomes("ACGTACGT", "ACGTTCGT"); + let records = vcf_records(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n"); + assert_eq!( + records, + vec![ + "chr1\t5\t.\tA\tT\t.\tPASS\tCHAIN_ID=1;QUERY_CHROM=qry1;QUERY_POS=5;STRAND=+" + .to_owned() + ] + ); + assert_ref_matches_genome("ACGTACGT", &records); + } + + #[test] + fn vcf_emits_adjacent_substitutions_separately() { + let genomes = genomes("ACGTACGT", "AGTTACGT"); + let records = vcf_records(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n"); + assert_eq!(records.len(), 2, "{records:?}"); + assert!(records[0].starts_with("chr1\t2\t.\tC\tG\t")); + assert!(records[1].starts_with("chr1\t3\t.\tG\tT\t")); + assert_ref_matches_genome("ACGTACGT", &records); + } + + #[test] + fn vcf_ignores_soft_masking_and_ambiguous_bases() { + // Case differs on every base and the query has an N: neither is a + // callable substitution. + let genomes = genomes("ACGTACGT", "acgtNcgt"); + let records = vcf_records(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n"); + assert!(records.is_empty(), "{records:?}"); + } + + #[test] + fn vcf_anchors_insertions_on_the_preceding_base() { + let genomes = genomes("ACGTACGT", "ACGTGACGT"); + let records = vcf_records( + &genomes, + "chain 100 chr1 8 + 0 8 qry1 9 + 0 9 1\n4\t0\t1\n4\n\n", + ); + assert_eq!( + records, + vec![ + "chr1\t4\t.\tT\tTG\t.\tPASS\tCHAIN_ID=1;QUERY_CHROM=qry1;QUERY_POS=5;STRAND=+" + .to_owned() + ] + ); + assert_ref_matches_genome("ACGTACGT", &records); + } + + #[test] + fn vcf_encodes_multi_base_insertions() { + let genomes = genomes("ACGTACGT", "ACGTGGGACGT"); + let records = vcf_records( + &genomes, + "chain 100 chr1 8 + 0 8 qry1 11 + 0 11 1\n4\t0\t3\n4\n\n", + ); + assert!( + records[0].starts_with("chr1\t4\t.\tT\tTGGG\t"), + "{records:?}" + ); + assert_ref_matches_genome("ACGTACGT", &records); + } + + #[test] + fn vcf_encodes_one_base_deletions() { + let genomes = genomes("ACGTGACGT", "ACGTACGT"); + let records = vcf_records( + &genomes, + "chain 100 chr1 9 + 0 9 qry1 8 + 0 8 1\n4\t1\t0\n4\n\n", + ); + assert_eq!( + records, + vec![ + "chr1\t4\t.\tTG\tT\t.\tPASS\tCHAIN_ID=1;QUERY_CHROM=qry1;QUERY_POS=5;STRAND=+" + .to_owned() + ] + ); + assert_ref_matches_genome("ACGTGACGT", &records); + } + + #[test] + fn vcf_encodes_multi_base_deletions() { + let genomes = genomes("ACGTGGGACGT", "ACGTACGT"); + let records = vcf_records( + &genomes, + "chain 100 chr1 11 + 0 11 qry1 8 + 0 8 1\n4\t3\t0\n4\n\n", + ); + assert!( + records[0].starts_with("chr1\t4\t.\tTGGG\tT\t"), + "{records:?}" + ); + assert_ref_matches_genome("ACGTGGGACGT", &records); + } + + #[test] + fn vcf_encodes_dual_sided_gaps_as_replacements() { + let genomes = genomes("ACGTGGACGT", "ACGTTTTACGT"); + let records = vcf_records( + &genomes, + "chain 100 chr1 10 + 0 10 qry1 11 + 0 11 1\n4\t2\t3\n4\n\n", + ); + assert_eq!( + records, + vec![ + "chr1\t5\t.\tGG\tTTT\t.\tPASS\tCHAIN_ID=1;QUERY_CHROM=qry1;QUERY_POS=5;STRAND=+" + .to_owned() + ] + ); + assert_ref_matches_genome("ACGTGGACGT", &records); + } + + #[test] + fn vcf_trims_redundant_bases_from_replacements() { + // Reference gap "GA" against query gap "TA": the shared trailing A is + // redundant and is trimmed away. + let genomes = genomes("ACGTGAACGT", "ACGTTAACGT"); + let records = vcf_records( + &genomes, + "chain 100 chr1 10 + 0 10 qry1 10 + 0 10 1\n4\t2\t2\n4\n\n", + ); + assert_eq!(records.len(), 1, "{records:?}"); + assert!(records[0].starts_with("chr1\t5\t.\tG\tT\t"), "{records:?}"); + assert_ref_matches_genome("ACGTGAACGT", &records); + } + + #[test] + fn vcf_handles_negative_query_strand() { + // Query forward "TTGACTAA": the chain's query interval [2,6) read on the + // minus strand is revcomp("GACT") = "AGTC", which differs from the + // reference "AGTG" at its last base. + let genomes = genomes("AGTG", "TTGACTAA"); + let records = vcf_records(&genomes, "chain 100 chr1 4 + 0 4 qry1 8 - 2 6 5\n4\n\n"); + assert_eq!( + records, + vec![ + "chr1\t4\t.\tG\tC\t.\tPASS\tCHAIN_ID=5;QUERY_CHROM=qry1;QUERY_POS=6;STRAND=-" + .to_owned() + ] + ); + assert_ref_matches_genome("AGTG", &records); + } + + #[test] + fn vcf_handles_a_variant_at_reference_position_zero() { + let genomes = genomes("ACGT", "TCGT"); + let records = vcf_records(&genomes, "chain 100 chr1 4 + 0 4 qry1 4 + 0 4 1\n4\n\n"); + assert!(records[0].starts_with("chr1\t1\t.\tA\tT\t"), "{records:?}"); + assert_ref_matches_genome("ACGT", &records); + } + + #[test] + fn vcf_uses_a_following_anchor_at_the_start_of_the_target_span() { + // A leading zero-length block puts the gap at target offset 0, where + // there is no preceding base to anchor on. + let genomes = genomes("ACGT", "GACGT"); + let records = vcf_records( + &genomes, + "chain 100 chr1 4 + 0 4 qry1 5 + 0 5 1\n0\t0\t1\n4\n\n", + ); + assert_eq!(records.len(), 1, "{records:?}"); + assert!(records[0].starts_with("chr1\t1\t.\tA\tGA\t"), "{records:?}"); + assert_ref_matches_genome("ACGT", &records); + } + + #[test] + fn vcf_emits_a_variant_adjacent_to_the_final_block() { + let genomes = genomes("ACGTGT", "ACGTA"); + let records = vcf_records( + &genomes, + "chain 100 chr1 6 + 0 6 qry1 5 + 0 5 1\n4\t1\t0\n1\n\n", + ); + assert_eq!(records.len(), 2, "{records:?}"); + assert!(records[0].starts_with("chr1\t4\t.\tTG\tT\t"), "{records:?}"); + assert!(records[1].starts_with("chr1\t6\t.\tT\tA\t"), "{records:?}"); + assert_ref_matches_genome("ACGTGT", &records); + } + + /// Runs a MAF conversion and returns its stanzas as line groups. + fn maf_stanzas(genomes: &Genomes, chain: &str) -> Vec> { + convert(genomes, ConvertTarget::Maf, chain) + .lines() + .skip_while(|line| line.starts_with('#') || line.is_empty()) + .collect::>() + .split(|line| line.is_empty()) + .filter(|stanza| !stanza.is_empty()) + .map(|stanza| stanza.iter().map(|line| (*line).to_owned()).collect()) + .collect() + } + + /// Splits an `s` row into its MAF fields. + fn maf_row(line: &str) -> (String, u32, u32, char, u32, String) { + let fields: Vec<&str> = line.split(' ').collect(); + assert_eq!(fields[0], "s", "not an s row: {line}"); + ( + fields[1].to_owned(), + fields[2].parse().expect("numeric start"), + fields[3].parse().expect("numeric size"), + fields[4].chars().next().expect("strand"), + fields[5].parse().expect("numeric srcSize"), + fields[6].to_owned(), + ) + } + + /// Asserts a row's gap-free text reproduces its declared source interval. + /// + /// `source` is the forward sequence; a `-` row is checked against the + /// reverse complement, as MAF declares coordinates in the row's own + /// orientation. + fn assert_row_matches_source(source: &str, line: &str) { + let (_, start, size, strand, src_size, text) = maf_row(line); + assert_eq!(src_size as usize, source.len(), "srcSize mismatch: {line}"); + let bases: String = text.chars().filter(|base| *base != '-').collect(); + assert_eq!(bases.len(), size as usize, "size mismatch: {line}"); + + let oriented = match strand { + '+' => source.as_bytes().to_vec(), + _ => { + let mut reversed = source.as_bytes().to_vec(); + reverse_complement_in_place(&mut reversed); + reversed + } + }; + let start = start as usize; + assert_eq!( + bases.as_bytes(), + &oriented[start..start + size as usize], + "text does not reproduce the declared interval: {line}" + ); + } + + #[test] + fn maf_emits_a_pairwise_stanza_for_an_identical_alignment() { + let genomes = genomes("ACGTACGT", "ACGTACGT"); + let output = convert( + &genomes, + ConvertTarget::Maf, + "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 3\n8\n\n", + ); + assert!(output.starts_with("##maf version=1 scoring=chain\n")); + let stanzas = maf_stanzas(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 3\n8\n\n"); + assert_eq!( + stanzas, + vec![vec![ + "a score=100 chain_id=3".to_owned(), + "s chr1 0 8 + 8 ACGTACGT".to_owned(), + "s qry1 0 8 + 8 ACGTACGT".to_owned(), + ]] + ); + } + + #[test] + fn maf_carries_substitutions_in_the_aligned_text() { + let genomes = genomes("ACGTACGT", "ACGTTCGT"); + let stanzas = maf_stanzas(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n"); + assert_eq!(stanzas[0][1], "s chr1 0 8 + 8 ACGTACGT"); + assert_eq!(stanzas[0][2], "s qry1 0 8 + 8 ACGTTCGT"); + assert_row_matches_source("ACGTACGT", &stanzas[0][1]); + assert_row_matches_source("ACGTTCGT", &stanzas[0][2]); + } + + #[test] + fn maf_emits_reference_only_gaps_as_query_dashes() { + let genomes = genomes("ACGTGACGT", "ACGTACGT"); + let stanzas = maf_stanzas( + &genomes, + "chain 100 chr1 9 + 0 9 qry1 8 + 0 8 1\n4\t1\t0\n4\n\n", + ); + assert_eq!(stanzas.len(), 1, "{stanzas:?}"); + assert_eq!(stanzas[0][1], "s chr1 0 9 + 9 ACGTGACGT"); + assert_eq!(stanzas[0][2], "s qry1 0 8 + 8 ACGT-ACGT"); + assert_row_matches_source("ACGTGACGT", &stanzas[0][1]); + assert_row_matches_source("ACGTACGT", &stanzas[0][2]); + } + + #[test] + fn maf_emits_query_only_gaps_as_reference_dashes() { + let genomes = genomes("ACGTACGT", "ACGTGACGT"); + let stanzas = maf_stanzas( + &genomes, + "chain 100 chr1 8 + 0 8 qry1 9 + 0 9 1\n4\t0\t1\n4\n\n", + ); + assert_eq!(stanzas[0][1], "s chr1 0 8 + 8 ACGT-ACGT"); + assert_eq!(stanzas[0][2], "s qry1 0 9 + 9 ACGTGACGT"); + assert_row_matches_source("ACGTACGT", &stanzas[0][1]); + assert_row_matches_source("ACGTGACGT", &stanzas[0][2]); + } + + #[test] + fn maf_keeps_several_one_sided_gaps_in_one_stanza() { + let genomes = genomes("ACGTGACGTACGT", "ACGTACGTTTACGT"); + let stanzas = maf_stanzas( + &genomes, + "chain 100 chr1 13 + 0 13 qry1 14 + 0 14 1\n4\t1\t0\n4\t0\t2\n4\n\n", + ); + assert_eq!(stanzas.len(), 1, "{stanzas:?}"); + assert_eq!(stanzas[0][1], "s chr1 0 13 + 13 ACGTGACGT--ACGT"); + assert_eq!(stanzas[0][2], "s qry1 0 14 + 14 ACGT-ACGTTTACGT"); + assert_row_matches_source("ACGTGACGTACGT", &stanzas[0][1]); + assert_row_matches_source("ACGTACGTTTACGT", &stanzas[0][2]); + } + + #[test] + fn maf_splits_at_dual_sided_gaps() { + let genomes = genomes("ACGTGGACGT", "ACGTTTTACGT"); + let stanzas = maf_stanzas( + &genomes, + "chain 100 chr1 10 + 0 10 qry1 11 + 0 11 1\n4\t2\t3\n4\n\n", + ); + assert_eq!(stanzas.len(), 2, "{stanzas:?}"); + assert_eq!(stanzas[0][1], "s chr1 0 4 + 10 ACGT"); + assert_eq!(stanzas[0][2], "s qry1 0 4 + 11 ACGT"); + assert_eq!(stanzas[1][1], "s chr1 6 4 + 10 ACGT"); + assert_eq!(stanzas[1][2], "s qry1 7 4 + 11 ACGT"); + for stanza in &stanzas { + assert_row_matches_source("ACGTGGACGT", &stanza[1]); + assert_row_matches_source("ACGTTTTACGT", &stanza[2]); + } + } + + #[test] + fn maf_keeps_minus_strand_query_coordinates_and_orientation() { + // Query forward "TTGACTAA"; the chain's query interval [2,6) on the + // minus strand is revcomp("GACT") = "AGTC". The row must declare the + // minus-strand coordinates, not forward ones relabelled `+`. + let genomes = genomes("AGTC", "TTGACTAA"); + let stanzas = maf_stanzas(&genomes, "chain 100 chr1 4 + 0 4 qry1 8 - 2 6 9\n4\n\n"); + assert_eq!( + stanzas, + vec![vec![ + "a score=100 chain_id=9".to_owned(), + "s chr1 0 4 + 4 AGTC".to_owned(), + "s qry1 2 4 - 8 AGTC".to_owned(), + ]] + ); + assert_row_matches_source("TTGACTAA", &stanzas[0][2]); + } + + #[test] + fn maf_preserves_soft_masking() { + let genomes = genomes("acgtACGT", "acgtACGT"); + let stanzas = maf_stanzas(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n"); + assert_eq!(stanzas[0][1], "s chr1 0 8 + 8 acgtACGT"); + } + + /// Renders each BAM record as `name flags pos cigar seq` for comparison. + #[cfg(feature = "bam")] + fn bam_records(genomes: &Genomes, chain: &str) -> (Vec, Vec<(String, u32)>) { + use noodles_bam as bam; + use noodles_sam::alignment::record::cigar::op::Kind; + + let bytes = convert_bytes(genomes, ConvertTarget::Bam, chain); + let mut reader = bam::io::Reader::new(Cursor::new(bytes)); + let header = reader.read_header().expect("read BAM header"); + let references = header + .reference_sequences() + .iter() + .map(|(name, map)| { + ( + String::from_utf8_lossy(name).into_owned(), + map.length().get() as u32, + ) + }) + .collect(); + + let mut rendered = Vec::new(); + for record in reader.records() { + let record = record.expect("read BAM record"); + let cigar: String = record + .cigar() + .iter() + .map(|op| { + let op = op.expect("valid CIGAR op"); + let kind = match op.kind() { + Kind::SequenceMatch => '=', + Kind::SequenceMismatch => 'X', + Kind::Insertion => 'I', + Kind::Deletion => 'D', + other => panic!("unexpected CIGAR kind {other:?}"), + }; + format!("{}{kind}", op.len()) + }) + .collect(); + rendered.push(format!( + "{} {} {} {cigar} {} {}", + String::from_utf8_lossy(record.name().expect("record name")), + record.flags().bits(), + usize::from( + record + .alignment_start() + .expect("mapped record") + .expect("valid position") + ), + String::from_utf8_lossy(&record.sequence().iter().collect::>()), + record + .mapping_quality() + .map_or(255, |quality| quality.get()), + )); + } + (rendered, references) + } + + #[cfg(feature = "bam")] + #[test] + fn bam_derives_sq_entries_and_emits_one_record_per_chain() { + let genomes = genomes("ACGTACGT", "ACGTACGT"); + let (records, references) = + bam_records(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 4\n8\n\n"); + assert_eq!(references, vec![("chr1".to_owned(), 8)]); + assert_eq!(records, vec!["qry1.4 0 1 8= ACGTACGT 0".to_owned()]); + } + + #[cfg(feature = "bam")] + #[test] + fn bam_distinguishes_matches_from_mismatches() { + let genomes = genomes("ACGTACGT", "ACGTTCGT"); + let (records, _) = bam_records(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n"); + assert_eq!(records, vec!["qry1.1 0 1 4=1X3= ACGTTCGT 0".to_owned()]); + } + + #[cfg(feature = "bam")] + #[test] + fn bam_encodes_one_sided_gaps_in_the_cigar() { + let deletion = genomes("ACGTGACGT", "ACGTACGT"); + let (records, _) = bam_records( + &deletion, + "chain 100 chr1 9 + 0 9 qry1 8 + 0 8 1\n4\t1\t0\n4\n\n", + ); + assert_eq!(records, vec!["qry1.1 0 1 4=1D4= ACGTACGT 0".to_owned()]); + + let insertion = genomes("ACGTACGT", "ACGTGACGT"); + let (records, _) = bam_records( + &insertion, + "chain 100 chr1 8 + 0 8 qry1 9 + 0 9 1\n4\t0\t1\n4\n\n", + ); + assert_eq!(records, vec!["qry1.1 0 1 4=1I4= ACGTGACGT 0".to_owned()]); + } + + #[cfg(feature = "bam")] + #[test] + fn bam_splits_dual_sided_gaps_into_suffixed_records() { + let genomes = genomes("ACGTGGACGT", "ACGTTTTACGT"); + let (records, _) = bam_records( + &genomes, + "chain 100 chr1 10 + 0 10 qry1 11 + 0 11 2\n4\t2\t3\n4\n\n", + ); + assert_eq!( + records, + vec![ + "qry1.2.1 0 1 4= ACGT 0".to_owned(), + "qry1.2.2 0 7 4= ACGT 0".to_owned(), + ] + ); + } + + #[cfg(feature = "bam")] + #[test] + fn bam_flags_and_orients_a_minus_strand_query() { + let genomes = genomes("AGTC", "TTGACTAA"); + let (records, _) = bam_records(&genomes, "chain 100 chr1 4 + 0 4 qry1 8 - 2 6 6\n4\n\n"); + // FLAG 16 is REVERSE; SEQ is the oriented query, revcomp("GACT"). + assert_eq!(records, vec!["qry1.6 16 1 4= AGTC 0".to_owned()]); + } + + #[cfg(feature = "bam")] + #[test] + fn bam_uppercases_soft_masked_query_bases() { + let genomes = genomes("acgtACGT", "acgtACGT"); + let (records, _) = bam_records(&genomes, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n"); + assert_eq!(records, vec!["qry1.1 0 1 8= ACGTACGT 0".to_owned()]); + } + + #[cfg(feature = "bam")] + #[test] + fn bam_rejects_a_chain_target_size_that_contradicts_the_reference() { + let genomes = genomes("ACGTACGT", "ACGTACGT"); + let mut stdin = Cursor::new(b"chain 100 chr1 9 + 0 8 qry1 8 + 0 8 1\n8\n\n".to_vec()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + ConvertArgs { + to: ConvertTarget::Bam, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: None, + output: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("contradictory target size should be rejected"); + assert!(err.to_string().contains("but the reference has 8"), "{err}"); + } + + #[test] + fn rejects_gzip_for_bam_output() { + let genomes = genomes("ACGT", "ACGT"); + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + ConvertArgs { + to: ConvertTarget::Bam, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: None, + output: None, + gzip: true, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("gzip-wrapped BAM should be rejected"); + assert!( + err.to_string().contains("already block-compressed"), + "{err}" + ); + } + + #[test] + fn rejects_blocks_that_do_not_sum_to_the_header_span() { + let genomes = genomes("ACGTACGT", "ACGTACGT"); + let mut stdin = Cursor::new(b"chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n4\n\n".to_vec()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + ConvertArgs { + to: ConvertTarget::Vcf, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: None, + output: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("truncated block list should be rejected"); + assert!(err.to_string().contains("do not sum to the header span")); + } + + #[test] + fn rejects_coordinates_outside_the_reference() { + let genomes = genomes("ACGT", "ACGT"); + let mut stdin = Cursor::new(b"chain 100 chr1 8 + 0 8 qry1 4 + 0 4 1\n8\n\n".to_vec()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + ConvertArgs { + to: ConvertTarget::Vcf, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: None, + output: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("out-of-bounds chain should be rejected"); + assert!(err.to_string().contains("exceeds"), "{err}"); + } + + #[test] + fn rejects_a_missing_query_sequence() { + let genomes = genomes("ACGT", "ACGT"); + let mut stdin = Cursor::new(b"chain 100 chr1 4 + 0 4 absent 4 + 0 4 1\n4\n\n".to_vec()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + ConvertArgs { + to: ConvertTarget::Vcf, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: None, + output: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("missing query sequence should be rejected"); + assert!(err.to_string().contains("missing sequence"), "{err}"); + } + + #[test] + fn rejects_output_equal_to_the_input_chain() { + let genomes = genomes("ACGT", "ACGT"); + let chain_path = genomes.dir.join("in.chain"); + std::fs::write(&chain_path, "chain 100 chr1 4 + 0 4 qry1 4 + 0 4 1\n4\n\n") + .expect("write chain"); + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + ConvertArgs { + to: ConvertTarget::Vcf, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: Some(chain_path.clone()), + output: Some(chain_path), + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("overwriting the input chain should be rejected"); + assert!( + err.to_string().contains("must not be the same path"), + "{err}" + ); + } + + #[test] + fn reads_chains_from_a_file_input() { + let genomes = genomes("ACGTACGT", "ACGTTCGT"); + let chain_path = genomes.dir.join("in.chain"); + std::fs::write(&chain_path, "chain 100 chr1 8 + 0 8 qry1 8 + 0 8 1\n8\n\n") + .expect("write chain"); + let output_path = genomes.dir.join("out.vcf"); + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + run( + ConvertArgs { + to: ConvertTarget::Vcf, + reference: genomes.reference.clone(), + query: genomes.query.clone(), + chain: Some(chain_path), + output: Some(output_path.clone()), + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect("convert run"); + + let written = std::fs::read_to_string(&output_path).expect("read output"); + assert!(written.contains("chr1\t5\t.\tA\tT\t"), "{written}"); + assert!(stdout.is_empty(), "file output should not touch stdout"); + let _ = Path::new(&output_path); + } +} diff --git a/src/cli/liftover.rs b/src/cli/liftover.rs new file mode 100644 index 0000000..96fc56d --- /dev/null +++ b/src/cli/liftover.rs @@ -0,0 +1,1676 @@ +// Copyright (c) 2026 Alejandro Gonzales-Irribarren +// Distributed under the terms of the Apache License, Version 2.0. + +use std::collections::HashMap; +use std::fmt; +use std::fs::File; +use std::io::{BufRead, BufWriter, Write}; +use std::path::{Path, PathBuf}; + +use chaintools::{ChainError, OwnedChain, Strand, StreamingReader}; +use clap::{Args, ValueEnum}; +use genepred::writer::TargetFormat; +use genepred::{Bed3, Bed6, Bed12, BedFormat, GenePred, Reader, Strand as BedStrand, Writer}; +use rust_lapper::{Interval, Lapper}; + +use super::CliError; + +const OUTPUT_BUFFER_CAPACITY: usize = 1024 * 1024; +const DEFAULT_MIN_MATCH: f64 = 0.95; + +/// Command-line arguments for the liftover subcommand. +/// +/// Maps intervals given on a chain's reference/target assembly onto its query +/// assembly. No genome sequence is needed: a chain carries all the geometry. +/// For the opposite direction, swap the chain first (`chaintools swap`). +/// +/// # Examples +/// +/// ```bash +/// chaintools liftover --chain map.chain --bed input.bed --type 6 --output lifted.bed +/// ``` +#[derive(Debug, Args)] +pub struct LiftoverArgs { + #[arg( + short = 'c', + long = "chain", + value_name = "PATH", + help = "Path to the .chain file mapping reference/target coordinates onto query coordinates." + )] + chain: PathBuf, + + #[arg( + short = 'b', + long = "bed", + value_name = "PATH", + help = "Path to the input BED file, in reference/target coordinates." + )] + bed: PathBuf, + + #[arg( + long = "type", + value_name = "WIDTH", + value_enum, + help = "BED width of the input, also used for the output. Required so no column is silently ignored." + )] + bed_type: BedWidth, + + #[arg( + short = 'o', + long = "output", + value_name = "PATH", + help = "Path for the lifted BED output. If not provided, output is written to standard output." + )] + output: Option, + + #[arg( + short = 'u', + long = "unmapped", + value_name = "PATH", + help = "Path for records that could not be lifted, each preceded by a reason comment." + )] + unmapped: Option, + + #[arg( + short = 'm', + long = "min-match", + value_name = "FRACTION", + default_value_t = DEFAULT_MIN_MATCH, + help = "Minimum fraction of input bases that must fall inside aligned chain blocks (0.0-1.0)." + )] + min_match: f64, + + #[arg( + long = "multiple", + help = "Emit every qualifying mapping instead of rejecting records that several chains can place." + )] + multiple: bool, +} + +/// BED widths supported by liftover. +#[derive(Debug, Clone, Copy, PartialEq, Eq, ValueEnum)] +enum BedWidth { + #[value(name = "3")] + Three, + #[value(name = "6")] + Six, + #[value(name = "12")] + Twelve, +} + +impl fmt::Display for BedWidth { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + BedWidth::Three => f.write_str("3"), + BedWidth::Six => f.write_str("6"), + BedWidth::Twelve => f.write_str("12"), + } + } +} + +/// Runs the liftover subcommand. +/// +/// # Arguments +/// +/// * `args` - Liftover arguments +/// * `_stdin` - Unused: both inputs are files +/// * `stdout` - Output stream (used if no --output provided) +/// * `_stderr` - Error/logging output +/// +/// # Output +/// +/// Returns `Ok(())` on success or `Err(CliError)` on failure +pub fn run( + args: LiftoverArgs, + _stdin: &mut R, + stdout: &mut W, + _stderr: &mut E, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, + E: Write, +{ + validate_min_match(args.min_match)?; + super::ensure_inputs_exist(&[("chain", &args.chain), ("BED", &args.bed)], &[])?; + validate_output_paths(&args)?; + + log::info!( + "liftover: chain={}, bed={}, type={}, min_match={}, multiple={}", + args.chain.display(), + args.bed.display(), + args.bed_type, + args.min_match, + args.multiple + ); + + let index = LiftoverIndex::build(&args.chain)?; + log::info!( + "indexed {} chains over {} reference sequence(s)", + index.chains.len(), + index.by_reference.len() + ); + if index.chains.is_empty() { + log::warn!("chain file contains no chains: every record will be unmapped"); + } + + let mut unmapped = match &args.unmapped { + Some(path) => Some(BufWriter::new(File::create(path)?)), + None => None, + }; + + let stats = if let Some(path) = &args.output { + let mut writer = BufWriter::with_capacity(OUTPUT_BUFFER_CAPACITY, File::create(path)?); + let stats = lift_stream(&args, &index, &mut writer, &mut unmapped)?; + writer.flush()?; + stats + } else { + let stats = lift_stream(&args, &index, stdout, &mut unmapped)?; + stdout.flush()?; + stats + }; + + if let Some(writer) = &mut unmapped { + writer.flush()?; + } + + super::log_summary( + "liftover", + &[ + ("chains", index.chains.len() as u64), + ("records", stats.read), + ("mapped", stats.mapped), + ("unmapped", stats.unmapped), + ("mappings", stats.mappings), + ], + ); + if stats.multiple > 0 { + log::info!( + "{} record(s) produced more than one mapping (--multiple)", + stats.multiple + ); + } + if stats.unmapped > 0 && args.unmapped.is_none() { + log::warn!( + "{} record(s) were dropped without being written; pass --unmapped PATH to keep them", + stats.unmapped + ); + } + Ok(()) +} + +/// Rejects a `--min-match` outside `0.0..=1.0`. +fn validate_min_match(min_match: f64) -> Result<(), CliError> { + if !(0.0..=1.0).contains(&min_match) { + return Err(CliError::Message(format!( + "--min-match must be between 0.0 and 1.0, got {min_match}" + ))); + } + Ok(()) +} + +/// Rejects output paths that would clobber an input or each other. +fn validate_output_paths(args: &LiftoverArgs) -> Result<(), CliError> { + for (label, path) in [ + ("--output", args.output.as_deref()), + ("--unmapped", args.unmapped.as_deref()), + ] { + if let Some(path) = path { + super::validate_distinct_paths(label, path, Some(&args.chain))?; + super::validate_distinct_paths(label, path, Some(&args.bed))?; + } + } + if let (Some(output), Some(unmapped)) = (&args.output, &args.unmapped) + && output == unmapped + { + return Err(CliError::Message( + "--output and --unmapped must not be the same path".to_owned(), + )); + } + Ok(()) +} + +/// Per-reference interval index over whole chain target spans. +/// +/// ponytail: index whole chain spans, then inspect candidate blocks. If +/// profiling shows candidate traversal dominates on pathological high-overlap +/// chain sets, upgrade the index to aligned-block spans. +#[derive(Debug)] +struct LiftoverIndex { + chains: Vec, + by_reference: HashMap, Lapper>, +} + +impl LiftoverIndex { + /// Parses the chain file once, validating each chain and indexing its span. + /// + /// Chains are kept parsed in memory because every BED record may query any + /// of them; only the target span plus a `Vec` index goes into the interval + /// index, so chains are never cloned. + fn build(path: &Path) -> Result { + let mut reader = StreamingReader::from_path(path)?; + let mut chains: Vec = Vec::new(); + let mut intervals: HashMap, Vec>> = HashMap::new(); + let mut reference_sizes: HashMap, u32> = HashMap::new(); + let mut query_sizes: HashMap, u32> = HashMap::new(); + + while let Some(header) = reader.next_header()? { + let offset = header.offset; + let blocks = reader.read_blocks(offset)?; + let chain = header.into_chain(blocks); + validate_chain(&chain, offset)?; + // Reference and query dictionaries stay separate: the same + // chromosome name legitimately has different sizes in two + // assemblies. + record_size( + &mut reference_sizes, + &chain.reference_name, + chain.reference_size, + offset, + "reference", + )?; + record_size( + &mut query_sizes, + &chain.query_name, + chain.query_size, + offset, + "query", + )?; + + intervals + .entry(chain.reference_name.clone()) + .or_default() + .push(Interval { + start: u64::from(chain.reference_start), + stop: u64::from(chain.reference_end), + val: chains.len(), + }); + chains.push(chain); + } + + let by_reference = intervals + .into_iter() + .map(|(name, spans)| (name, Lapper::new(spans))) + .collect(); + Ok(LiftoverIndex { + chains, + by_reference, + }) + } + + /// Returns the indices of chains whose target span overlaps `[start, end)`. + /// + /// Order follows the interval index, which sorts by span, so it never + /// depends on `HashMap` iteration order. + fn candidates(&self, chrom: &[u8], start: u64, end: u64) -> Vec { + match self.by_reference.get(chrom) { + Some(lapper) => lapper.find(start, end).map(|span| span.val).collect(), + None => Vec::new(), + } + } +} + +/// Validates the invariants liftover relies on, at index time. +fn validate_chain(chain: &OwnedChain, offset: usize) -> Result<(), CliError> { + if chain.reference_strand != Strand::Plus { + return Err(CliError::Chain(format_error( + offset, + "liftover requires chains whose target strand is +; run `chaintools swap` to canonicalize", + ))); + } + if chain.reference_start > chain.reference_end { + return Err(CliError::Chain(format_error( + offset, + "chain target span is inverted", + ))); + } + if chain.query_start > chain.query_end { + return Err(CliError::Chain(format_error( + offset, + "chain query span is inverted", + ))); + } + if chain.reference_end > chain.reference_size { + return Err(CliError::Chain(format_error( + offset, + "chain target span exceeds its sequence size", + ))); + } + if chain.query_end > chain.query_size { + return Err(CliError::Chain(format_error( + offset, + "chain query span exceeds its sequence size", + ))); + } + super::validate_block_spans(chain, offset)?; + Ok(()) +} + +/// Records a chromosome size, rejecting a conflicting one for the same name. +fn record_size( + sizes: &mut HashMap, u32>, + name: &[u8], + size: u32, + offset: usize, + side: &str, +) -> Result<(), CliError> { + match sizes.get(name) { + Some(known) if *known != size => Err(CliError::Chain(format_error( + offset, + format!( + "{side} sequence {} is declared as both {known} and {size} bases", + String::from_utf8_lossy(name) + ), + ))), + _ => { + if !sizes.contains_key(name) { + sizes.insert(name.to_vec(), size); + } + Ok(()) + } + } +} + +/// One contiguous piece of a source interval that a chain places on the query. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +struct MappedSegment { + source_start: u64, + source_end: u64, + query_start: u64, + query_end: u64, +} + +/// A source interval's mapping through one candidate chain. +struct CandidateMapping { + chain_index: usize, + mapped_bases: u64, + /// Mapped segments per source span, parallel to the requested spans. + spans: Vec>, +} + +/// Why a record could not be lifted. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +enum UnmappedReason { + Deleted, + Partial, + Split, + Duplicated, +} + +impl fmt::Display for UnmappedReason { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + UnmappedReason::Deleted => f.write_str("Deleted in new"), + UnmappedReason::Partial => f.write_str("Partially deleted in new"), + UnmappedReason::Split => f.write_str("Split in new"), + UnmappedReason::Duplicated => f.write_str("Duplicated in new"), + } + } +} + +/// Maps one source interval through one chain's aligned blocks. +/// +/// Returns one segment per aligned block the interval overlaps, in source +/// order. Query coordinates are returned on the forward query sequence, so a `-` +/// chain yields segments whose query coordinates descend as the source ascends. +/// +/// # Arguments +/// +/// * `chain` - Candidate chain, target side on `+` +/// * `start` - Source interval start, 0-based +/// * `end` - Source interval end, exclusive +/// +/// # Output +/// +/// Returns `Ok(Vec)`, empty when no aligned block overlaps +fn map_interval_on_chain( + chain: &OwnedChain, + start: u64, + end: u64, +) -> Result, CliError> { + let query_size = u64::from(chain.query_size); + let mut segments = Vec::new(); + let mut reference = u64::from(chain.reference_start); + let mut query = u64::from(chain.query_start); + + for block in &chain.blocks { + let size = u64::from(block.size); + let reference_end = reference + size; + if reference >= end { + break; + } + + let overlap_start = start.max(reference); + let overlap_end = end.min(reference_end); + if overlap_start < overlap_end { + let oriented_start = query + (overlap_start - reference); + let oriented_end = query + (overlap_end - reference); + let (query_start, query_end) = match chain.query_strand { + Strand::Plus => (oriented_start, oriented_end), + // The chain's own query coordinates run along the reverse + // strand, so flipping them around the chromosome size yields + // forward coordinates. + Strand::Minus => ( + checked_flip(query_size, oriented_end)?, + checked_flip(query_size, oriented_start)?, + ), + }; + segments.push(MappedSegment { + source_start: overlap_start, + source_end: overlap_end, + query_start, + query_end, + }); + } + + reference = reference_end + u64::from(block.gap_reference); + query = query + size + u64::from(block.gap_query); + } + + Ok(segments) +} + +/// Converts an oriented query coordinate to its forward-strand equivalent. +fn checked_flip(query_size: u64, coordinate: u64) -> Result { + query_size.checked_sub(coordinate).ok_or_else(|| { + CliError::Message("chain query coordinate exceeds its sequence size".to_owned()) + }) +} + +/// Evaluates every candidate chain for a record's source spans. +/// +/// Coverage counts only source bases that fall inside aligned blocks, never +/// bounding-span overlap, and every span must contribute at least one segment +/// (for BED12 that is the "every block must map" rule). +/// +/// # Output +/// +/// Returns `(qualifying, touching)`, where `touching` counts candidates that +/// mapped at least one base whether or not they qualified. +fn evaluate_candidates( + index: &LiftoverIndex, + chrom: &[u8], + spans: &[(u64, u64)], + lookup: (u64, u64), + min_match: f64, +) -> Result<(Vec, usize), CliError> { + let source_bases: u64 = spans.iter().map(|(start, end)| end - start).sum(); + let mut qualifying = Vec::new(); + let mut touching = 0usize; + + for chain_index in index.candidates(chrom, lookup.0, lookup.1) { + let chain = &index.chains[chain_index]; + let mut mapped_spans = Vec::with_capacity(spans.len()); + let mut mapped_bases = 0u64; + for (start, end) in spans { + let segments = map_interval_on_chain(chain, *start, *end)?; + mapped_bases += segments + .iter() + .map(|segment| segment.source_end - segment.source_start) + .sum::(); + mapped_spans.push(segments); + } + + if mapped_bases == 0 { + continue; + } + touching += 1; + if mapped_spans.iter().any(Vec::is_empty) { + continue; + } + // Cross-multiplied so no division is needed at the threshold boundary. + if (mapped_bases as f64) < min_match * source_bases as f64 { + continue; + } + qualifying.push(CandidateMapping { + chain_index, + mapped_bases, + spans: mapped_spans, + }); + } + + Ok((qualifying, touching)) +} + +/// Orders qualifying mappings deterministically for `--multiple` output. +fn sort_mappings(mappings: &mut [CandidateMapping], chains: &[OwnedChain]) { + mappings.sort_by(|a, b| { + let left = &chains[a.chain_index]; + let right = &chains[b.chain_index]; + b.mapped_bases + .cmp(&a.mapped_bases) + .then_with(|| right.score.cmp(&left.score)) + .then_with(|| left.query_name.cmp(&right.query_name)) + .then_with(|| left.query_start.cmp(&right.query_start)) + .then_with(|| left.id.cmp(&right.id)) + // Chain file order is the final tie-break, so duplicate ids cannot + // make the output nondeterministic. + .then_with(|| a.chain_index.cmp(&b.chain_index)) + }); +} + +/// Running counts for the end-of-run summary. +#[derive(Default)] +struct LiftoverStats { + read: u64, + mapped: u64, + unmapped: u64, + multiple: u64, + mappings: u64, +} + +/// Streams the BED input through the index at the requested width. +fn lift_stream( + args: &LiftoverArgs, + index: &LiftoverIndex, + output: &mut W, + unmapped: &mut Option>, +) -> Result { + match args.bed_type { + BedWidth::Three => lift_records::(args, index, output, unmapped), + BedWidth::Six => lift_records::(args, index, output, unmapped), + BedWidth::Twelve => lift_records::(args, index, output, unmapped), + } +} + +/// Lifts every BED record of one width, preserving input order. +/// +/// One record is held at a time: read, lifted, written, dropped. +fn lift_records( + args: &LiftoverArgs, + index: &LiftoverIndex, + output: &mut W, + unmapped: &mut Option>, +) -> Result +where + F: BedFormat + Into + TargetFormat, + W: Write, +{ + let mut reader = Reader::::from_path(&args.bed).map_err(bed_error)?; + let mut stats = LiftoverStats::default(); + let mut warned_about_extras = false; + + for record in reader.records() { + let record = record.map_err(bed_error)?; + stats.read += 1; + if !warned_about_extras && !record.extras().is_empty() { + // Trailing columns are copied through verbatim, which is right for + // BED6+N/BED12+N annotations but wrong if the user understated the + // width and a coordinate or strand column is sitting in there. + log::warn!( + "BED input has columns beyond --type {}; they are copied through unlifted", + args.bed_type + ); + warned_about_extras = true; + } + + match lift_record(index, &record, args)? { + Ok(lifted) => { + stats.mapped += 1; + stats.mappings += lifted.len() as u64; + if lifted.len() > 1 { + stats.multiple += 1; + } + for record in &lifted { + write_record::(record, output)?; + } + } + Err(reason) => { + stats.unmapped += 1; + if let Some(writer) = unmapped { + writeln!(writer, "# {reason}")?; + write_record::(&record, writer)?; + } + } + } + } + + Ok(stats) +} + +/// Lifts one BED record, returning either its mappings or why it failed. +fn lift_record( + index: &LiftoverIndex, + record: &GenePred, + args: &LiftoverArgs, +) -> Result, UnmappedReason>, CliError> { + let spans = source_spans(record, args.bed_type)?; + let lookup = (record.start, record.end); + let (mut qualifying, touching) = + evaluate_candidates(index, &record.chrom, &spans, lookup, args.min_match)?; + + if qualifying.is_empty() { + return Ok(Err(match touching { + 0 => UnmappedReason::Deleted, + 1 => UnmappedReason::Partial, + _ => UnmappedReason::Split, + })); + } + if qualifying.len() > 1 && !args.multiple { + return Ok(Err(UnmappedReason::Duplicated)); + } + + sort_mappings(&mut qualifying, &index.chains); + let mut lifted = Vec::with_capacity(qualifying.len()); + for mapping in &qualifying { + let chain = &index.chains[mapping.chain_index]; + match args.bed_type { + BedWidth::Three | BedWidth::Six => { + lifted.push(lift_interval_record(record, chain, &mapping.spans[0])?); + } + BedWidth::Twelve => match lift_block_record(record, chain, &mapping.spans)? { + Some(record) => lifted.push(record), + // The thick interval could not be placed consistently; the + // record is rejected rather than clamped. + None => return Ok(Err(UnmappedReason::Partial)), + }, + } + } + Ok(Ok(lifted)) +} + +/// Returns the source spans whose coverage decides a record's mapping. +/// +/// BED3/BED6 contribute their single interval; BED12 contributes one span per +/// block, so intronic bases never count toward coverage. +fn source_spans(record: &GenePred, bed_type: BedWidth) -> Result, CliError> { + if record.start >= record.end { + return Err(CliError::Message(format!( + "BED record {}:{}-{} is empty or inverted; liftover requires start < end", + String::from_utf8_lossy(&record.chrom), + record.start, + record.end + ))); + } + + if bed_type != BedWidth::Twelve { + return Ok(vec![(record.start, record.end)]); + } + + let (Some(starts), Some(ends)) = (record.block_starts(), record.block_ends()) else { + return Err(CliError::Message(format!( + "BED12 record {}:{}-{} has no blocks", + String::from_utf8_lossy(&record.chrom), + record.start, + record.end + ))); + }; + if starts.len() != ends.len() || starts.is_empty() { + return Err(CliError::Message(format!( + "BED12 record {}:{}-{} has inconsistent block arrays", + String::from_utf8_lossy(&record.chrom), + record.start, + record.end + ))); + } + + let mut spans = Vec::with_capacity(starts.len()); + for (start, end) in starts.iter().zip(ends) { + if start >= end { + return Err(CliError::Message(format!( + "BED12 record {}:{}-{} has an empty block", + String::from_utf8_lossy(&record.chrom), + record.start, + record.end + ))); + } + spans.push((*start, *end)); + } + Ok(spans) +} + +/// Builds a lifted BED3/BED6 record spanning the mapped extent. +/// +/// A source interval crossing a chain gap still yields one output interval: +/// chain blocks are not exposed as separate BED records. +fn lift_interval_record( + record: &GenePred, + chain: &OwnedChain, + segments: &[MappedSegment], +) -> Result { + let (start, end) = mapped_extent(segments)?; + let mut lifted = record.clone(); + lifted.chrom = chain.query_name.clone(); + lifted.start = start; + lifted.end = end; + lifted.strand = lift_strand(record.strand, chain.query_strand); + Ok(lifted) +} + +/// Builds a lifted BED12 record, projecting every block and the thick interval. +/// +/// Returns `Ok(None)` when the thick interval cannot be placed inside the lifted +/// record, which rejects the record instead of clamping it. +fn lift_block_record( + record: &GenePred, + chain: &OwnedChain, + spans: &[Vec], +) -> Result, CliError> { + let mut blocks = Vec::with_capacity(spans.len()); + for segments in spans { + blocks.push(mapped_extent(segments)?); + } + // Transcript order reverses through a `-` chain while BED block arrays stay + // ascending by genomic coordinate. + blocks.sort_unstable(); + + for pair in blocks.windows(2) { + if pair[0].1 > pair[1].0 { + return Err(CliError::Message(format!( + "lifted blocks of {} overlap on {}", + String::from_utf8_lossy(record.name().unwrap_or(b".")), + String::from_utf8_lossy(&chain.query_name) + ))); + } + } + + let start = blocks[0].0; + let end = blocks[blocks.len() - 1].1; + let (thick_start, thick_end) = match lift_thick(record, chain, start)? { + Some(thick) => thick, + None => return Ok(None), + }; + if thick_start < start || thick_end > end { + return Ok(None); + } + + let mut lifted = record.clone(); + lifted.chrom = chain.query_name.clone(); + lifted.start = start; + lifted.end = end; + lifted.strand = lift_strand(record.strand, chain.query_strand); + lifted.thick_start = Some(thick_start); + lifted.thick_end = Some(thick_end); + lifted.block_count = Some(u32::try_from(blocks.len()).map_err(|_| { + CliError::Message("lifted record has too many blocks for BED12".to_owned()) + })?); + // The writer derives blockSizes and blockStarts from these absolute + // coordinates, so source blockStarts are never reused. + lifted.block_starts = Some(blocks.iter().map(|(start, _)| *start).collect()); + lifted.block_ends = Some(blocks.iter().map(|(_, end)| *end).collect()); + Ok(Some(lifted)) +} + +/// Maps a BED12 thick interval through the same chain as its blocks. +/// +/// Returns `Ok(None)` when a non-empty thick interval has no mapped segment. +fn lift_thick( + record: &GenePred, + chain: &OwnedChain, + lifted_start: u64, +) -> Result, CliError> { + let thick_start = record.thick_start.unwrap_or(record.start); + let thick_end = record.thick_end.unwrap_or(record.end); + + if thick_start >= thick_end { + // ponytail: a zero-width thick interval (the noncoding convention, + // thickStart == thickEnd == chromStart) collapses to the lifted record + // start. If interior zero-width thick regions ever need their exact + // lifted boundary, map the adjacent base instead. + return Ok(Some((lifted_start, lifted_start))); + } + + let segments = map_interval_on_chain(chain, thick_start, thick_end)?; + if segments.is_empty() { + return Ok(None); + } + Ok(Some(mapped_extent(&segments)?)) +} + +/// Returns the extent covered by mapped segments on the forward query strand. +fn mapped_extent(segments: &[MappedSegment]) -> Result<(u64, u64), CliError> { + let start = segments + .iter() + .map(|segment| segment.query_start) + .min() + .ok_or_else(|| CliError::Message("no mapped segments to place".to_owned()))?; + let end = segments + .iter() + .map(|segment| segment.query_end) + .max() + .expect("non-empty segments have a maximum"); + Ok((start, end)) +} + +/// Transforms a BED strand through a chain's query strand. +/// +/// A `+` chain preserves the input strand; a `-` chain flips it. An absent or +/// unknown strand is never given a direction. +fn lift_strand(input: Option, chain_query: Strand) -> Option { + match (input, chain_query) { + (Some(BedStrand::Forward), Strand::Minus) => Some(BedStrand::Reverse), + (Some(BedStrand::Reverse), Strand::Minus) => Some(BedStrand::Forward), + (other, _) => other, + } +} + +fn write_record(record: &GenePred, writer: &mut W) -> Result<(), CliError> +where + F: TargetFormat, + W: Write, +{ + Writer::::from_record(record, writer).map_err(|err| CliError::Message(err.to_string())) +} + +fn bed_error(err: genepred::reader::ReaderError) -> CliError { + CliError::Message(format!("invalid BED input: {err}")) +} + +fn format_error(offset: usize, message: impl Into) -> ChainError { + ChainError::Format { + offset, + msg: message.into().into(), + } +} + +#[cfg(test)] +mod tests { + use super::*; + use clap::Parser; + use std::io::Cursor; + use std::sync::atomic::{AtomicUsize, Ordering}; + + #[derive(Debug, Parser)] + struct LiftoverHarness { + #[command(flatten)] + args: LiftoverArgs, + } + + static NEXT_TEMP_ID: AtomicUsize = AtomicUsize::new(0); + + /// One 100 bp block, `+/+`: chr1:100-200 maps to qry1:500-600. + const CHAIN_PLUS: &str = "chain 100 chr1 1000 + 100 200 qry1 2000 + 500 600 1\n100\n\n"; + /// The same block with a `-` query: oriented 500-600 is forward 1400-1500. + const CHAIN_MINUS: &str = "chain 200 chr1 1000 + 100 200 qry1 2000 - 500 600 2\n100\n\n"; + /// Two blocks with a gap on both sides: chr1:100-150 and chr1:200-250. + const CHAIN_GAP: &str = + "chain 300 chr1 1000 + 100 250 qry1 2000 + 500 640 3\n50\t50\t40\n50\n\n"; + /// A target-only gap: chr1:100-150 and chr1:200-250 stay contiguous on qry1. + const CHAIN_REF_GAP: &str = + "chain 310 chr1 1000 + 100 250 qry1 2000 + 500 600 5\n50\t50\t0\n50\n\n"; + /// A query-only gap: chr1:100-200 is contiguous but qry1 jumps 40 bases. + const CHAIN_QUERY_GAP: &str = + "chain 320 chr1 1000 + 100 200 qry1 2000 + 500 640 6\n50\t0\t40\n50\n\n"; + /// A chain on another reference sequence. + const CHAIN_CHR2: &str = "chain 400 chr2 800 + 0 100 qry2 900 + 0 100 4\n100\n\n"; + + /// A 1:1 chain over chr1:0-1000, offsetting coordinates by +2000. + const CHAIN_WIDE_PLUS: &str = + "chain 1000 chr1 10000 + 0 1000 qry1 10000 + 2000 3000 10\n1000\n\n"; + /// The same region with a `-` query: source [a,b) becomes [8000-b, 8000-a). + const CHAIN_WIDE_MINUS: &str = + "chain 1000 chr1 10000 + 0 1000 qry1 10000 - 2000 3000 11\n1000\n\n"; + + /// A transcript with two 100 bp exons and a coding region of 150-450. + const BED12_TWO_EXONS: &str = + "chr1\t100\t500\ttx1\t742\t+\t150\t450\t255,0,0\t2\t100,100\t0,300\n"; + + struct Fixture { + dir: PathBuf, + chain: PathBuf, + bed: PathBuf, + } + + impl Drop for Fixture { + fn drop(&mut self) { + let _ = std::fs::remove_dir_all(&self.dir); + } + } + + fn fixture(chain: &str, bed: &str) -> Fixture { + let id = NEXT_TEMP_ID.fetch_add(1, Ordering::Relaxed); + let dir = std::env::temp_dir().join(format!( + "chaintools-liftover-test-{}-{id}", + std::process::id() + )); + std::fs::create_dir_all(&dir).expect("create temp dir"); + let chain_path = dir.join("map.chain"); + let bed_path = dir.join("input.bed"); + std::fs::write(&chain_path, chain).expect("write chain"); + std::fs::write(&bed_path, bed).expect("write bed"); + Fixture { + dir, + chain: chain_path, + bed: bed_path, + } + } + + fn args_for( + fixture: &Fixture, + bed_type: BedWidth, + min_match: f64, + multiple: bool, + ) -> LiftoverArgs { + LiftoverArgs { + chain: fixture.chain.clone(), + bed: fixture.bed.clone(), + bed_type, + output: None, + unmapped: Some(fixture.dir.join("unmapped.bed")), + min_match, + multiple, + } + } + + /// Lifts the fixture and returns `(output lines, unmapped file lines)`. + fn lift( + fixture: &Fixture, + bed_type: BedWidth, + min_match: f64, + multiple: bool, + ) -> (Vec, Vec) { + let args = args_for(fixture, bed_type, min_match, multiple); + let unmapped_path = args.unmapped.clone().expect("unmapped path"); + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + run(args, &mut stdin, &mut stdout, &mut stderr).expect("liftover run"); + + let output = String::from_utf8(stdout).expect("utf8 output"); + let unmapped = std::fs::read_to_string(&unmapped_path).unwrap_or_default(); + (lines(&output), lines(&unmapped)) + } + + fn lines(text: &str) -> Vec { + text.lines().map(str::to_owned).collect() + } + + fn lift_error(fixture: &Fixture, bed_type: BedWidth, min_match: f64) -> String { + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + run( + args_for(fixture, bed_type, min_match, false), + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("liftover should fail") + .to_string() + } + + fn parse_chain(text: &str) -> OwnedChain { + StreamingReader::new(Cursor::new(text.as_bytes())) + .next_chain() + .expect("chain parses") + .expect("one chain") + } + + fn map(chain: &str, start: u64, end: u64) -> Vec { + map_interval_on_chain(&parse_chain(chain), start, end).expect("mapping succeeds") + } + + #[test] + fn parses_minimal_args() { + let cli = LiftoverHarness::try_parse_from([ + "chaintools", + "--chain", + "map.chain", + "--bed", + "in.bed", + "--type", + "6", + ]) + .expect("liftover arguments should parse"); + assert_eq!(cli.args.bed_type, BedWidth::Six); + assert_eq!(cli.args.min_match, DEFAULT_MIN_MATCH); + assert!(!cli.args.multiple); + assert!(cli.args.unmapped.is_none()); + } + + #[test] + fn requires_chain_bed_and_type() { + for missing in [ + vec!["chaintools", "--bed", "in.bed", "--type", "3"], + vec!["chaintools", "--chain", "map.chain", "--type", "3"], + vec!["chaintools", "--chain", "map.chain", "--bed", "in.bed"], + ] { + let err = LiftoverHarness::try_parse_from(missing) + .expect_err("incomplete arguments should be rejected"); + assert_eq!(err.kind(), clap::error::ErrorKind::MissingRequiredArgument); + } + } + + #[test] + fn rejects_unsupported_bed_width() { + let err = LiftoverHarness::try_parse_from([ + "chaintools", + "--chain", + "map.chain", + "--bed", + "in.bed", + "--type", + "9", + ]) + .expect_err("BED9 is not supported"); + assert_eq!(err.kind(), clap::error::ErrorKind::InvalidValue); + } + + // --- Milestone 1: index ------------------------------------------------- + + fn index_of(chain: &str) -> LiftoverIndex { + let fixture = fixture(chain, ""); + LiftoverIndex::build(&fixture.chain).expect("index builds") + } + + #[test] + fn index_returns_only_chains_on_the_queried_reference() { + let index = index_of(&format!("{CHAIN_PLUS}{CHAIN_CHR2}")); + assert_eq!(index.chains.len(), 2); + assert_eq!(index.candidates(b"chr1", 100, 200), vec![0]); + assert_eq!(index.candidates(b"chr2", 0, 100), vec![1]); + assert!(index.candidates(b"chr3", 0, 100).is_empty()); + } + + #[test] + fn index_returns_no_candidates_outside_any_chain_span() { + let index = index_of(CHAIN_PLUS); + assert!(index.candidates(b"chr1", 900, 950).is_empty()); + // Half-open spans: touching the boundary is not an overlap. + assert!(index.candidates(b"chr1", 200, 300).is_empty()); + assert!(index.candidates(b"chr1", 50, 100).is_empty()); + } + + #[test] + fn index_returns_every_overlapping_chain() { + let index = index_of(&format!("{CHAIN_PLUS}{CHAIN_GAP}")); + let mut candidates = index.candidates(b"chr1", 120, 130); + candidates.sort_unstable(); + assert_eq!(candidates, vec![0, 1]); + } + + #[test] + fn index_rejects_conflicting_chromosome_sizes() { + let conflicting = + format!("{CHAIN_PLUS}chain 100 chr1 999 + 100 200 qry1 2000 + 700 800 9\n100\n\n"); + let fixture = fixture(&conflicting, ""); + let err = LiftoverIndex::build(&fixture.chain) + .expect_err("conflicting reference sizes should be rejected") + .to_string(); + assert!(err.contains("declared as both 1000 and 999"), "{err}"); + } + + #[test] + fn index_allows_one_name_to_differ_between_assemblies() { + // chr1 is 1000 bases on the reference and 2000 on the query: two + // assemblies, not a conflict. + let index = index_of("chain 100 chr1 1000 + 100 200 chr1 2000 + 500 600 1\n100\n\n"); + assert_eq!(index.chains.len(), 1); + } + + #[test] + fn index_rejects_a_minus_target_strand() { + let fixture = fixture( + "chain 100 chr1 1000 - 100 200 qry1 2000 + 500 600 1\n100\n\n", + "", + ); + let err = LiftoverIndex::build(&fixture.chain) + .expect_err("minus target strand should be rejected") + .to_string(); + assert!(err.contains("target strand is +"), "{err}"); + } + + #[test] + fn index_rejects_blocks_that_disagree_with_the_header() { + let fixture = fixture( + "chain 100 chr1 1000 + 100 200 qry1 2000 + 500 600 1\n50\n\n", + "", + ); + let err = LiftoverIndex::build(&fixture.chain) + .expect_err("short block list should be rejected") + .to_string(); + assert!(err.contains("do not sum to the target span"), "{err}"); + } + + // --- Milestone 2: interval mapper -------------------------------------- + + #[test] + fn maps_a_full_block_on_the_plus_strand() { + assert_eq!( + map(CHAIN_PLUS, 100, 200), + vec![MappedSegment { + source_start: 100, + source_end: 200, + query_start: 500, + query_end: 600 + }] + ); + } + + #[test] + fn maps_a_partial_block_on_the_plus_strand() { + assert_eq!( + map(CHAIN_PLUS, 120, 150), + vec![MappedSegment { + source_start: 120, + source_end: 150, + query_start: 520, + query_end: 550 + }] + ); + } + + #[test] + fn maps_a_full_block_on_the_minus_strand() { + assert_eq!( + map(CHAIN_MINUS, 100, 200), + vec![MappedSegment { + source_start: 100, + source_end: 200, + query_start: 1400, + query_end: 1500 + }] + ); + } + + #[test] + fn maps_a_partial_block_on_the_minus_strand() { + // Oriented [520, 550) flips to forward [2000-550, 2000-520). + assert_eq!( + map(CHAIN_MINUS, 120, 150), + vec![MappedSegment { + source_start: 120, + source_end: 150, + query_start: 1450, + query_end: 1480 + }] + ); + } + + #[test] + fn an_interval_inside_one_block_keeps_its_length_on_both_strands() { + for chain in [CHAIN_PLUS, CHAIN_MINUS] { + for (start, end) in [(100, 200), (120, 150), (199, 200), (100, 101)] { + let segments = map(chain, start, end); + assert_eq!(segments.len(), 1, "{chain}"); + let segment = segments[0]; + assert_eq!( + segment.query_end - segment.query_start, + end - start, + "length changed for {start}..{end}" + ); + } + } + } + + #[test] + fn negative_strand_mapping_satisfies_the_flip_invariant() { + let chain = parse_chain(CHAIN_MINUS); + let query_size = u64::from(chain.query_size); + for (start, end) in [(100, 200), (120, 150), (100, 101), (150, 200)] { + let segments = map(CHAIN_MINUS, start, end); + let segment = segments[0]; + // The oriented interval is the chain's own query coordinate space. + let oriented_start = u64::from(chain.query_start) + (start - 100); + let oriented_end = u64::from(chain.query_start) + (end - 100); + assert_eq!(segment.query_start, query_size - oriented_end); + assert_eq!(segment.query_end, query_size - oriented_start); + assert_eq!(segment.query_end - segment.query_start, end - start); + } + } + + #[test] + fn maps_across_a_target_only_gap() { + let segments = map(CHAIN_REF_GAP, 100, 250); + assert_eq!(segments.len(), 2); + assert_eq!((segments[0].query_start, segments[0].query_end), (500, 550)); + assert_eq!((segments[1].query_start, segments[1].query_end), (550, 600)); + } + + #[test] + fn maps_across_a_query_only_gap() { + let segments = map(CHAIN_QUERY_GAP, 100, 200); + assert_eq!(segments.len(), 2); + assert_eq!((segments[0].query_start, segments[0].query_end), (500, 550)); + assert_eq!((segments[1].query_start, segments[1].query_end), (590, 640)); + } + + #[test] + fn maps_across_a_dual_sided_gap_without_covering_it() { + let segments = map(CHAIN_GAP, 100, 250); + assert_eq!(segments.len(), 2); + assert_eq!( + (segments[0].source_start, segments[0].source_end), + (100, 150) + ); + assert_eq!( + (segments[1].source_start, segments[1].source_end), + (200, 250) + ); + let mapped: u64 = segments + .iter() + .map(|segment| segment.source_end - segment.source_start) + .sum(); + // The 50 unaligned target bases in the gap are not mapped. + assert_eq!(mapped, 100); + } + + #[test] + fn an_interval_inside_a_gap_maps_nothing() { + assert!(map(CHAIN_GAP, 150, 200).is_empty()); + } + + #[test] + fn block_boundaries_are_half_open() { + assert_eq!(map(CHAIN_GAP, 100, 150).len(), 1); + assert_eq!(map(CHAIN_GAP, 200, 250).len(), 1); + // Ends exactly where the first block ends, starts where the second does. + assert_eq!(map(CHAIN_GAP, 140, 150)[0].query_end, 550); + assert_eq!(map(CHAIN_GAP, 200, 210)[0].query_start, 590); + } + + // --- Milestone 2: BED3 ------------------------------------------------- + + #[test] + fn lifts_a_bed3_interval_through_a_plus_chain() { + let fixture = fixture(CHAIN_PLUS, "chr1\t120\t150\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert_eq!(output, vec!["qry1\t520\t550"]); + assert!(unmapped.is_empty()); + } + + #[test] + fn lifts_a_bed3_interval_through_a_minus_chain() { + let fixture = fixture(CHAIN_MINUS, "chr1\t120\t150\n"); + let (output, _) = lift(&fixture, BedWidth::Three, 0.95, false); + assert_eq!(output, vec!["qry1\t1450\t1480"]); + } + + #[test] + fn a_bed3_interval_crossing_a_gap_stays_one_record() { + let fixture = fixture(CHAIN_GAP, "chr1\t100\t250\n"); + let (output, _) = lift(&fixture, BedWidth::Three, 0.6, false); + // One output interval spanning the mapped extent, not two records. + assert_eq!(output, vec!["qry1\t500\t640"]); + } + + #[test] + fn reports_an_interval_with_no_chain_as_deleted() { + let fixture = fixture(CHAIN_PLUS, "chr9\t100\t150\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert!(output.is_empty()); + assert_eq!(unmapped, vec!["# Deleted in new", "chr9\t100\t150"]); + } + + #[test] + fn reports_an_interval_inside_a_gap_as_deleted() { + let fixture = fixture(CHAIN_GAP, "chr1\t150\t200\n"); + let (_, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert_eq!(unmapped[0], "# Deleted in new"); + } + + #[test] + fn reports_a_single_insufficient_chain_as_partially_deleted() { + let fixture = fixture(CHAIN_GAP, "chr1\t100\t250\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert!(output.is_empty()); + assert_eq!(unmapped[0], "# Partially deleted in new"); + } + + #[test] + fn reports_several_insufficient_chains_as_split() { + // Two chains each cover half of chr1:100-300, neither on its own. + let chains = concat!( + "chain 100 chr1 1000 + 100 200 qry1 2000 + 500 600 1\n100\n\n", + "chain 100 chr1 1000 + 200 300 qry2 2000 + 700 800 2\n100\n\n" + ); + let fixture = fixture(chains, "chr1\t100\t300\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert!(output.is_empty()); + assert_eq!(unmapped[0], "# Split in new"); + } + + #[test] + fn never_stitches_two_chains_into_one_mapping() { + // The same two half-covering chains must not combine into 100%. + let chains = concat!( + "chain 100 chr1 1000 + 100 200 qry1 2000 + 500 600 1\n100\n\n", + "chain 100 chr1 1000 + 200 300 qry1 2000 + 600 700 2\n100\n\n" + ); + let fixture = fixture(chains, "chr1\t100\t300\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert!(output.is_empty(), "{output:?}"); + assert_eq!(unmapped[0], "# Split in new"); + // Each chain alone covers only half, so a lower threshold reports them + // as competing candidates rather than a single stitched mapping. + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.5, false); + assert!(output.is_empty(), "{output:?}"); + assert_eq!(unmapped[0], "# Duplicated in new"); + } + + // --- Milestone 3: min-match ------------------------------------------- + + /// A chain covering `covered` of the 100 bases at chr1:100-200. + fn partial_chain(covered: u32) -> String { + format!( + "chain 100 chr1 1000 + 100 {} qry1 2000 + 500 {} 1\n{covered}\n\n", + 100 + covered, + 500 + covered + ) + } + + #[test] + fn min_match_accepts_full_coverage() { + let fixture = fixture(CHAIN_PLUS, "chr1\t100\t200\n"); + let (output, _) = lift(&fixture, BedWidth::Three, 1.0, false); + assert_eq!(output, vec!["qry1\t500\t600"]); + } + + #[test] + fn min_match_boundary_is_inclusive() { + let fixture = fixture(&partial_chain(95), "chr1\t100\t200\n"); + let (output, _) = lift(&fixture, BedWidth::Three, 0.95, false); + assert_eq!(output, vec!["qry1\t500\t595"], "95% must pass at 0.95"); + } + + #[test] + fn min_match_rejects_one_base_below_the_threshold() { + let fixture = fixture(&partial_chain(94), "chr1\t100\t200\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert!(output.is_empty(), "94% must fail at 0.95"); + assert_eq!(unmapped[0], "# Partially deleted in new"); + // The same chain passes once the threshold drops to 0.94. + let (output, _) = lift(&fixture, BedWidth::Three, 0.94, false); + assert_eq!(output, vec!["qry1\t500\t594"]); + } + + #[test] + fn min_match_counts_aligned_blocks_not_the_chain_span() { + // CHAIN_GAP spans chr1:100-250 but aligns only 100 of those 150 bases. + let fixture = fixture(CHAIN_GAP, "chr1\t100\t250\n"); + let (output, _) = lift(&fixture, BedWidth::Three, 0.66, false); + assert_eq!(output, vec!["qry1\t500\t640"]); + let (output, _) = lift(&fixture, BedWidth::Three, 0.67, false); + assert!(output.is_empty(), "bounding-span overlap must not count"); + } + + #[test] + fn min_match_zero_still_requires_one_mapped_base() { + let fixture = fixture(CHAIN_GAP, "chr1\t150\t200\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.0, false); + assert!(output.is_empty()); + assert_eq!(unmapped[0], "# Deleted in new"); + } + + #[test] + fn rejects_an_out_of_range_min_match() { + let fixture = fixture(CHAIN_PLUS, "chr1\t100\t200\n"); + for value in [-0.1, 1.1] { + let err = lift_error(&fixture, BedWidth::Three, value); + assert!( + err.contains("--min-match must be between 0.0 and 1.0"), + "{err}" + ); + } + } + + // --- Milestone 3: BED6 ------------------------------------------------ + + #[test] + fn lifts_bed6_preserving_name_and_flipping_strand() { + let fixture = fixture(CHAIN_MINUS, "chr1\t120\t150\tfeature\t742\t+\n"); + let (output, _) = lift(&fixture, BedWidth::Six, 0.95, false); + // Score is not preserved: genepred's GenePred does not carry it. + assert_eq!(output, vec!["qry1\t1450\t1480\tfeature\t0\t-"]); + } + + #[test] + fn lifts_every_strand_combination() { + for (chain, input, expected) in [ + (CHAIN_PLUS, "+", "+"), + (CHAIN_PLUS, "-", "-"), + (CHAIN_MINUS, "+", "-"), + (CHAIN_MINUS, "-", "+"), + (CHAIN_PLUS, ".", "."), + (CHAIN_MINUS, ".", "."), + ] { + let fixture = fixture(chain, &format!("chr1\t120\t150\tf\t0\t{input}\n")); + let (output, _) = lift(&fixture, BedWidth::Six, 0.95, false); + let strand = output[0].split('\t').nth(5).expect("strand column"); + assert_eq!(strand, expected, "{input} through {chain}"); + } + } + + // --- Milestone 4: multiple mappings ----------------------------------- + + /// Two chains that both place chr1:100-200 in full, on different targets. + const CHAIN_DUPLICATES: &str = concat!( + "chain 500 chr1 1000 + 100 200 qry1 2000 + 500 600 1\n100\n\n", + "chain 900 chr1 1000 + 100 200 qry2 2000 + 300 400 2\n100\n\n" + ); + + #[test] + fn duplicate_mappings_are_rejected_by_default() { + let fixture = fixture(CHAIN_DUPLICATES, "chr1\t100\t200\n"); + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, false); + assert!(output.is_empty()); + assert_eq!(unmapped, vec!["# Duplicated in new", "chr1\t100\t200"]); + } + + #[test] + fn multiple_emits_every_mapping_deterministically() { + let fixture = fixture(CHAIN_DUPLICATES, "chr1\t100\t200\n"); + let expected = vec!["qry2\t300\t400".to_owned(), "qry1\t500\t600".to_owned()]; + // Equal mapped_bases, so the higher chain score wins the tie-break. + for _ in 0..5 { + let (output, unmapped) = lift(&fixture, BedWidth::Three, 0.95, true); + assert_eq!(output, expected, "output order must not vary"); + assert!(unmapped.is_empty()); + } + } + + #[test] + fn multiple_orders_by_mapped_bases_first() { + // The lower-scoring chain maps more bases, so it must come first. + let chains = concat!( + "chain 100 chr1 1000 + 100 200 qry1 2000 + 500 600 1\n100\n\n", + "chain 900 chr1 1000 + 100 190 qry2 2000 + 300 390 2\n90\n\n" + ); + let fixture = fixture(chains, "chr1\t100\t200\n"); + let (output, _) = lift(&fixture, BedWidth::Three, 0.9, true); + assert_eq!(output, vec!["qry1\t500\t600", "qry2\t300\t390"]); + } + + #[test] + fn input_order_is_preserved() { + let fixture = fixture( + CHAIN_PLUS, + "chr1\t150\t160\nchr1\t120\t130\nchr1\t100\t110\n", + ); + let (output, _) = lift(&fixture, BedWidth::Three, 0.95, false); + assert_eq!( + output, + vec!["qry1\t550\t560", "qry1\t520\t530", "qry1\t500\t510"] + ); + } + + // --- Milestone 5: BED12 ---------------------------------------------- + + #[test] + fn lifts_bed12_through_a_plus_chain() { + let fixture = fixture(CHAIN_WIDE_PLUS, BED12_TWO_EXONS); + let (output, _) = lift(&fixture, BedWidth::Twelve, 0.95, false); + assert_eq!( + output, + vec!["qry1\t2100\t2500\ttx1\t0\t+\t2150\t2450\t0,0,0\t2\t100,100,\t0,300,"] + ); + } + + #[test] + fn lifts_bed12_through_a_minus_chain_reversing_blocks() { + let fixture = fixture(CHAIN_WIDE_MINUS, BED12_TWO_EXONS); + let (output, _) = lift(&fixture, BedWidth::Twelve, 0.95, false); + // Exons swap genomic order, blockStarts are recomputed, strand flips. + assert_eq!( + output, + vec!["qry1\t7500\t7900\ttx1\t0\t-\t7550\t7850\t0,0,0\t2\t100,100,\t0,300,"] + ); + assert_bed12_invariants(&output[0]); + } + + #[test] + fn a_query_insertion_inside_an_exon_grows_that_block() { + // 5 query bases are inserted 150 bases into the transcript's first exon. + let chain = "chain 1 chr1 10000 + 0 1000 qry1 10000 + 2000 3005 12\n150\t0\t5\n850\n\n"; + let fixture = fixture(chain, BED12_TWO_EXONS); + let (output, _) = lift(&fixture, BedWidth::Twelve, 0.95, false); + assert_eq!( + output, + vec!["qry1\t2100\t2505\ttx1\t0\t+\t2155\t2455\t0,0,0\t2\t105,100,\t0,305,"] + ); + assert_bed12_invariants(&output[0]); + } + + #[test] + fn a_chain_gap_inside_an_intron_does_not_reduce_coverage() { + // The chain drops chr1:250-350, which is entirely intronic, so both + // exons still map in full and the record passes at min_match 1.0. + let chain = "chain 1 chr1 10000 + 0 1000 qry1 10000 + 2000 2900 13\n250\t100\t0\n650\n\n"; + let fixture = fixture(chain, BED12_TWO_EXONS); + let (output, _) = lift(&fixture, BedWidth::Twelve, 1.0, false); + assert_eq!( + output, + vec!["qry1\t2100\t2400\ttx1\t0\t+\t2150\t2350\t0,0,0\t2\t100,100,\t0,200,"] + ); + } + + #[test] + fn a_completely_unmapped_exon_rejects_the_record() { + // The chain ends before the second exon. + let chain = "chain 1 chr1 10000 + 0 250 qry1 10000 + 2000 2250 14\n250\n\n"; + let fixture = fixture(chain, BED12_TWO_EXONS); + let (output, unmapped) = lift(&fixture, BedWidth::Twelve, 0.5, false); + assert!(output.is_empty(), "every block must contribute: {output:?}"); + assert_eq!(unmapped[0], "# Partially deleted in new"); + } + + #[test] + fn lifts_a_three_exon_transcript() { + let bed = "chr1\t100\t700\ttx3\t0\t+\t100\t700\t0,0,0\t3\t100,100,100\t0,300,500\n"; + let fixture = fixture(CHAIN_WIDE_PLUS, bed); + let (output, _) = lift(&fixture, BedWidth::Twelve, 1.0, false); + assert_eq!( + output, + vec!["qry1\t2100\t2700\ttx3\t0\t+\t2100\t2700\t0,0,0\t3\t100,100,100,\t0,300,500,"] + ); + assert_bed12_invariants(&output[0]); + } + + #[test] + fn preserves_a_zero_width_thick_interval() { + // The noncoding convention: thickStart == thickEnd == chromStart. + let bed = "chr1\t100\t500\tnc1\t0\t+\t100\t100\t0,0,0\t2\t100,100\t0,300\n"; + let fixture = fixture(CHAIN_WIDE_PLUS, bed); + let (output, _) = lift(&fixture, BedWidth::Twelve, 0.95, false); + let fields: Vec<&str> = output[0].split('\t').collect(); + assert_eq!((fields[6], fields[7]), ("2100", "2100")); + } + + #[test] + fn rejects_a_record_whose_thick_interval_cannot_be_placed() { + // Both exons map, but the coding region falls entirely in a chain gap. + let chain = "chain 1 chr1 10000 + 0 1000 qry1 10000 + 2000 2800 15\n200\t200\t0\n600\n\n"; + let bed = "chr1\t100\t500\ttx4\t0\t+\t250\t350\t0,0,0\t2\t100,100\t0,300\n"; + let fixture = fixture(chain, bed); + let (output, unmapped) = lift(&fixture, BedWidth::Twelve, 0.95, false); + assert!(output.is_empty(), "{output:?}"); + assert_eq!(unmapped[0], "# Partially deleted in new"); + } + + /// Checks the structural BED12 invariants on one output line. + fn assert_bed12_invariants(line: &str) { + let fields: Vec<&str> = line.split('\t').collect(); + let start: u64 = fields[1].parse().expect("chromStart"); + let end: u64 = fields[2].parse().expect("chromEnd"); + let count: usize = fields[9].parse().expect("blockCount"); + let sizes: Vec = parse_list(fields[10]); + let starts: Vec = parse_list(fields[11]); + + assert_eq!(count, sizes.len(), "blockCount vs blockSizes: {line}"); + assert_eq!(count, starts.len(), "blockCount vs blockStarts: {line}"); + let mut previous_end = start; + for (offset, size) in starts.iter().zip(&sizes) { + let block_start = start + offset; + let block_end = block_start + size; + assert!(block_start >= start, "block before chromStart: {line}"); + assert!(block_end <= end, "block past chromEnd: {line}"); + assert!(block_start >= previous_end, "blocks out of order: {line}"); + previous_end = block_end; + } + } + + fn parse_list(field: &str) -> Vec { + field + .split(',') + .filter(|item| !item.is_empty()) + .map(|item| item.parse().expect("numeric list entry")) + .collect() + } + + // --- Error handling --------------------------------------------------- + + #[test] + fn rejects_a_zero_width_bed_interval() { + let fixture = fixture(CHAIN_PLUS, "chr1\t150\t150\n"); + let err = lift_error(&fixture, BedWidth::Three, 0.95); + assert!(err.contains("liftover requires start < end"), "{err}"); + } + + #[test] + fn rejects_output_paths_that_collide_with_inputs() { + let fixture = fixture(CHAIN_PLUS, "chr1\t100\t200\n"); + for output in [fixture.chain.clone(), fixture.bed.clone()] { + let mut args = args_for(&fixture, BedWidth::Three, 0.95, false); + args.output = Some(output); + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run(args, &mut stdin, &mut stdout, &mut stderr) + .expect_err("clobbering an input should be rejected") + .to_string(); + assert!(err.contains("must not be the same path"), "{err}"); + } + } + + #[test] + fn rejects_the_same_path_for_output_and_unmapped() { + let fixture = fixture(CHAIN_PLUS, "chr1\t100\t200\n"); + let mut args = args_for(&fixture, BedWidth::Three, 0.95, false); + let shared = fixture.dir.join("both.bed"); + args.output = Some(shared.clone()); + args.unmapped = Some(shared); + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run(args, &mut stdin, &mut stdout, &mut stderr) + .expect_err("one file cannot be both outputs") + .to_string(); + assert!(err.contains("--output and --unmapped"), "{err}"); + } + + #[test] + fn missing_inputs_are_rejected_up_front() { + let fixture = fixture(CHAIN_PLUS, "chr1\t100\t200\n"); + let mut args = args_for(&fixture, BedWidth::Three, 0.95, false); + args.bed = fixture.dir.join("absent.bed"); + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run(args, &mut stdin, &mut stdout, &mut stderr) + .expect_err("a missing BED should be rejected") + .to_string(); + assert!(err.contains("BED file does not exist"), "{err}"); + } + + #[test] + fn unmapped_records_are_dropped_only_when_no_path_is_given() { + let fixture = fixture(CHAIN_PLUS, "chr9\t100\t150\n"); + let mut args = args_for(&fixture, BedWidth::Three, 0.95, false); + args.unmapped = None; + let mut stdin = Cursor::new(Vec::::new()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + run(args, &mut stdin, &mut stdout, &mut stderr).expect("liftover run"); + assert!(stdout.is_empty()); + } +} diff --git a/src/cli/mod.rs b/src/cli/mod.rs index 777391e..716d2fd 100644 --- a/src/cli/mod.rs +++ b/src/cli/mod.rs @@ -3,16 +3,19 @@ pub mod anti_repeat; pub mod bed; +pub mod convert; pub mod filter; +pub mod liftover; pub mod merge; pub mod score; pub mod sort; mod sort_core; pub mod split; +pub mod swap; use std::fmt; use std::io::{self, BufRead, Write}; -use std::path::Path; +use std::path::{Path, PathBuf}; use chaintools::ChainError; use clap::{Parser, Subcommand}; @@ -62,8 +65,12 @@ enum Command { AntiRepeat(anti_repeat::AntiRepeatArgs), #[command(about = "Convert chain files to BED")] Bed(bed::BedArgs), + #[command(about = "Convert chain files to VCF, MAF, or BAM")] + Convert(convert::ConvertArgs), #[command(about = "Filter chain files")] Filter(filter::FilterArgs), + #[command(about = "Map BED coordinates through a chain (reference/target -> query)")] + Liftover(liftover::LiftoverArgs), #[command(about = "Merge chain files")] Merge(merge::MergeArgs), #[command(about = "Recompute chain scores from sequence (UCSC chainScore-compatible)")] @@ -72,6 +79,8 @@ enum Command { Split(split::SplitArgs), #[command(about = "Sort chain files")] Sort(sort::SortArgs), + #[command(about = "Swap the target and query sides of chains (UCSC chainSwap)")] + Swap(swap::SwapArgs), } impl std::fmt::Display for Command { @@ -79,11 +88,14 @@ impl std::fmt::Display for Command { match self { Command::AntiRepeat(_) => f.write_str("antirepeat"), Command::Bed(_) => f.write_str("bed"), + Command::Convert(_) => f.write_str("convert"), Command::Filter(_) => f.write_str("filter"), + Command::Liftover(_) => f.write_str("liftover"), Command::Merge(_) => f.write_str("merge"), Command::Score(_) => f.write_str("score"), Command::Split(_) => f.write_str("split"), Command::Sort(_) => f.write_str("sort"), + Command::Swap(_) => f.write_str("swap"), } } } @@ -184,11 +196,14 @@ where let result = match cli.command { Command::AntiRepeat(args) => anti_repeat::run(args, stdin, stdout, stderr), Command::Bed(args) => bed::run(args, stdin, stdout, stderr), + Command::Convert(args) => convert::run(args, stdin, stdout, stderr), Command::Filter(args) => filter::run(args, stdin, stdout, stderr), + Command::Liftover(args) => liftover::run(args, stdin, stdout, stderr), Command::Merge(args) => merge::run(args, stdin, stdout, stderr), Command::Score(args) => score::run(args, stdin, stdout, stderr), Command::Split(args) => split::run(args, stdin, stdout, stderr), Command::Sort(args) => sort::run(args, stdin, stdout, stderr), + Command::Swap(args) => swap::run(args, stdin, stdout, stderr), }; log::info!("Execution time: {:?}", start.elapsed()); @@ -283,6 +298,116 @@ pub(crate) fn ensure_inputs_exist( Ok(()) } +/// Checks that a chain's block list accounts for exactly both header spans. +/// +/// `sum(sizes) + sum(inner target gaps)` must equal the target span, and the +/// same must hold on the query side. Widened to `u64` so a malformed chain is +/// reported instead of overflowing. Shared by `swap` and `liftover`, which both +/// have to trust these spans before transforming coordinates. +/// +/// # Arguments +/// +/// * `chain` - The chain to check +/// * `offset` - Byte offset of the chain header, for error reporting +/// +/// # Output +/// +/// Returns `Ok(())` when both sums match, or `Err(ChainError::Format)` +pub(crate) fn validate_block_spans( + chain: &chaintools::OwnedChain, + offset: usize, +) -> Result<(), ChainError> { + let mut reference = 0u64; + let mut query = 0u64; + let last = chain.blocks.len().saturating_sub(1); + for (index, block) in chain.blocks.iter().enumerate() { + reference += u64::from(block.size); + query += u64::from(block.size); + if index < last { + reference += u64::from(block.gap_reference); + query += u64::from(block.gap_query); + } + } + + let reference_span = span_len(chain.reference_start, chain.reference_end, offset, "target")?; + let query_span = span_len(chain.query_start, chain.query_end, offset, "query")?; + if reference != u64::from(reference_span) { + return Err(span_error( + offset, + "chain blocks do not sum to the target span", + )); + } + if query != u64::from(query_span) { + return Err(span_error( + offset, + "chain blocks do not sum to the query span", + )); + } + Ok(()) +} + +fn span_len(start: u32, end: u32, offset: usize, label: &str) -> Result { + end.checked_sub(start) + .ok_or_else(|| span_error(offset, format!("{label} span is inverted"))) +} + +fn span_error(offset: usize, message: impl Into) -> ChainError { + ChainError::Format { + offset, + msg: message.into().into(), + } +} + +/// Rejects an output path that resolves to the same file as the input chain. +/// +/// Writing the output would truncate the input before it is read, so this is a +/// pre-flight check rather than a recoverable error. Nothing to do when the +/// input is standard input. +/// +/// # Arguments +/// +/// * `label` - Name of the output flag, used in the error message +/// * `output` - The requested output path +/// * `input` - The input chain path, or `None` for standard input +/// +/// # Output +/// +/// Returns `Ok(())` when the paths differ, or `Err(CliError::Message)` +pub(crate) fn validate_distinct_paths( + label: &str, + output: &Path, + input: Option<&Path>, +) -> Result<(), CliError> { + let Some(input) = input else { + return Ok(()); + }; + if resolve_path(output)? == resolve_path(input)? { + return Err(CliError::Message(format!( + "{label} must not be the same path as input chain {}", + input.display() + ))); + } + Ok(()) +} + +/// Resolves a path for comparison, tolerating a not-yet-created output file. +fn resolve_path(path: &Path) -> Result { + if path.exists() { + return Ok(std::fs::canonicalize(path)?); + } + + let parent = path.parent().unwrap_or_else(|| Path::new("")); + let parent = if parent.as_os_str().is_empty() { + std::env::current_dir()? + } else { + std::fs::canonicalize(parent)? + }; + let file_name = path + .file_name() + .ok_or_else(|| CliError::Message(format!("invalid path {}", path.display())))?; + Ok(parent.join(file_name)) +} + /// Checks that a single input file exists, mapping the outcome to a `CliError`. fn check_input_exists(label: &str, path: &Path) -> Result<(), CliError> { match path.try_exists() { diff --git a/src/cli/swap.rs b/src/cli/swap.rs new file mode 100644 index 0000000..fe83555 --- /dev/null +++ b/src/cli/swap.rs @@ -0,0 +1,438 @@ +// Copyright (c) 2026 Alejandro Gonzales-Irribarren +// Distributed under the terms of the Apache License, Version 2.0. + +use std::fs::File; +use std::io::{BufRead, BufWriter, Write}; +use std::path::PathBuf; + +use chaintools::io::writer::write_chain_dense; +use chaintools::{ChainError, OwnedChain, Strand, StreamItem, StreamingReader}; +use clap::Args; +#[cfg(feature = "gzip")] +use flate2::{Compression, write::GzEncoder}; + +use super::CliError; + +const OUTPUT_BUFFER_CAPACITY: usize = 1024 * 1024; + +/// Command-line arguments for the swap subcommand. +/// +/// Swaps the target and query sides of every chain, equivalent to UCSC +/// `chainSwap`. The swapped chain describes the same alignment in the opposite +/// direction; scores and chain IDs are preserved and output order matches input +/// order (pipe into `chaintools sort` if sorted output is needed). +/// +/// # Examples +/// +/// ```bash +/// chaintools swap --chain input.chain --out-chain swapped.chain +/// cat input.chain | chaintools swap > swapped.chain +/// ``` +#[derive(Debug, Args)] +pub struct SwapArgs { + #[arg( + short = 'c', + long = "chain", + value_name = "PATH", + help = "Path to the input .chain file. If not provided, chain data is read from standard input." + )] + chain: Option, + + #[arg( + short = 'o', + long = "out-chain", + value_name = "PATH", + help = "Path for the swapped chain output. If not provided, output is written to standard output." + )] + out_chain: Option, + + #[arg( + short = 'G', + long = "gzip", + help = "Compress swap output with gzip. Requires the `gzip` feature." + )] + gzip: bool, +} + +/// Runs the swap subcommand. +/// +/// Streams chains, swaps each chain's target and query sides, and writes the +/// result. Metadata/comment lines are passed through unchanged. +/// +/// # Arguments +/// +/// * `args` - Swap arguments +/// * `stdin` - Input stream (used if no --chain provided) +/// * `stdout` - Output stream (used if no --out-chain provided) +/// * `_stderr` - Error/logging output +/// +/// # Output +/// +/// Returns `Ok(())` on success or `Err(CliError)` on failure +pub fn run( + args: SwapArgs, + stdin: &mut R, + stdout: &mut W, + _stderr: &mut E, +) -> Result<(), CliError> +where + R: BufRead, + W: Write, + E: Write, +{ + validate_output_args(&args)?; + super::ensure_inputs_exist(&[], &[("input chain", args.chain.as_deref())])?; + if let Some(out_chain) = &args.out_chain { + super::validate_distinct_paths("--out-chain", out_chain, args.chain.as_deref())?; + } + + let input_desc = args + .chain + .as_deref() + .map_or_else(|| "".to_owned(), |path| path.display().to_string()); + let output_desc = args + .out_chain + .as_deref() + .map_or_else(|| "".to_owned(), |path| path.display().to_string()); + log::info!("swap: reading chains from {input_desc}, writing to {output_desc}"); + + if let Some(path) = &args.out_chain { + let file = File::create(path)?; + let writer = BufWriter::with_capacity(OUTPUT_BUFFER_CAPACITY, file); + if args.gzip { + run_gzip_output(&args, stdin, writer)?; + } else { + let mut writer = writer; + run_to_writer(&args, stdin, &mut writer)?; + writer.flush()?; + } + } else if args.gzip { + run_gzip_output(&args, stdin, stdout)?; + } else { + run_to_writer(&args, stdin, stdout)?; + stdout.flush()?; + } + + Ok(()) +} + +#[cfg(feature = "gzip")] +fn run_gzip_output(args: &SwapArgs, stdin: &mut R, writer: W) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + let mut writer = GzEncoder::new(writer, Compression::fast()); + run_to_writer(args, stdin, &mut writer)?; + writer.try_finish()?; + writer.get_mut().flush()?; + Ok(()) +} + +#[cfg(not(feature = "gzip"))] +fn run_gzip_output(_args: &SwapArgs, _stdin: &mut R, _writer: W) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + Err(CliError::Message( + "--gzip requires chaintools to be built with the `gzip` feature".to_owned(), + )) +} + +#[cfg(feature = "gzip")] +fn validate_output_args(_args: &SwapArgs) -> Result<(), CliError> { + Ok(()) +} + +#[cfg(not(feature = "gzip"))] +fn validate_output_args(args: &SwapArgs) -> Result<(), CliError> { + if args.gzip { + return Err(CliError::Message( + "--gzip requires chaintools to be built with the `gzip` feature".to_owned(), + )); + } + Ok(()) +} + +fn run_to_writer(args: &SwapArgs, stdin: &mut R, writer: &mut W) -> Result<(), CliError> +where + R: BufRead, + W: Write, +{ + let written = if let Some(path) = &args.chain { + let mut reader = StreamingReader::from_path(path)?; + process_stream(&mut reader, writer)? + } else { + let mut reader = StreamingReader::new(stdin); + process_stream(&mut reader, writer)? + }; + + super::log_summary("swap", &[("chains", written), ("written", written)]); + Ok(()) +} + +/// Streams chains through the swap transform, returning the number written. +fn process_stream( + reader: &mut StreamingReader, + writer: &mut W, +) -> Result { + let mut written = 0u64; + while let Some(item) = reader.next_item()? { + match item { + StreamItem::MetaLine(line) => { + writer.write_all(&line)?; + writer.write_all(b"\n")?; + } + StreamItem::Header(header) => { + let offset = header.offset; + let blocks = reader.read_blocks(offset)?; + let mut chain = header.into_chain(blocks); + swap_chain(&mut chain, offset)?; + write_chain_dense(writer, &chain)?; + written += 1; + } + } + } + Ok(written) +} + +/// Swaps a chain's target and query sides in place (UCSC `chainSwap`). +/// +/// The swapped chain describes the same alignment in the opposite direction: +/// the old query becomes the new target (always on `+`) and the old target +/// becomes the new query, keeping the original query strand. Score and chain ID +/// are untouched. +/// +/// For a `+` query both sides walk forward, so block order is preserved and each +/// block's `dt`/`dq` are exchanged. For a `-` query the new coordinates run +/// opposite to the original traversal, so block order reverses; because gaps sit +/// *between* blocks, the gap following a reversed block is the one recorded on +/// its successor in the reversed list, not the one it carried itself. +/// +/// # Arguments +/// +/// * `chain` - Chain to swap, mutated in place +/// * `offset` - Byte offset of the chain header, for error reporting +/// +/// # Output +/// +/// Returns `Ok(())` on success or `Err(ChainError)` for a malformed chain +fn swap_chain(chain: &mut OwnedChain, offset: usize) -> Result<(), ChainError> { + if chain.reference_strand != Strand::Plus { + return Err(format_error( + offset, + "cannot swap a chain whose target strand is not +", + )); + } + super::validate_block_spans(chain, offset)?; + + // Derive both new header intervals from the chromosome sizes before any + // field is moved, so no step depends on a half-swapped chain. + let (reference_start, reference_end, query_start, query_end) = match chain.query_strand { + Strand::Plus => ( + chain.query_start, + chain.query_end, + chain.reference_start, + chain.reference_end, + ), + Strand::Minus => ( + flip(chain.query_size, chain.query_end, offset)?, + flip(chain.query_size, chain.query_start, offset)?, + flip(chain.reference_size, chain.reference_end, offset)?, + flip(chain.reference_size, chain.reference_start, offset)?, + ), + }; + + std::mem::swap(&mut chain.reference_name, &mut chain.query_name); + std::mem::swap(&mut chain.reference_size, &mut chain.query_size); + chain.reference_start = reference_start; + chain.reference_end = reference_end; + chain.query_start = query_start; + chain.query_end = query_end; + // The new target is the old query read on `+`, and the new query is the old + // target read in the old query's direction: the strand letter carries over. + + match chain.query_strand { + Strand::Plus => { + for block in &mut chain.blocks { + std::mem::swap(&mut block.gap_reference, &mut block.gap_query); + } + } + Strand::Minus => { + chain.blocks.reverse(); + for index in 0..chain.blocks.len().saturating_sub(1) { + // `index + 1` still holds its original gaps here: this loop only + // ever reads ahead of what it has already written. + let successor = chain.blocks[index + 1]; + chain.blocks[index].gap_reference = successor.gap_query; + chain.blocks[index].gap_query = successor.gap_reference; + } + if let Some(last) = chain.blocks.last_mut() { + last.gap_reference = 0; + last.gap_query = 0; + } + } + } + + Ok(()) +} + +/// Converts a coordinate to the opposite strand of a sequence of `size` bases. +fn flip(size: u32, coord: u32, offset: usize) -> Result { + size.checked_sub(coord) + .ok_or_else(|| format_error(offset, "chain coordinate exceeds its sequence size")) +} + +fn format_error(offset: usize, message: impl Into) -> ChainError { + ChainError::Format { + offset, + msg: message.into().into(), + } +} + +#[cfg(test)] +mod tests { + use super::*; + use clap::Parser; + use std::io::Cursor; + + #[derive(Debug, Parser)] + struct SwapHarness { + #[command(flatten)] + args: SwapArgs, + } + + /// One `+/+` block. + const PLUS_ONE_BLOCK: &str = "chain 100 chr1 1000 + 10 30 qry1 500 + 40 60 7\n20\n\n"; + /// Several `+/+` blocks with unequal target/query gaps. + const PLUS_MANY_BLOCKS: &str = + "chain 4200 chr1 1000 + 10 63 qry1 500 + 40 93 9\n20\t3\t1\n10\t5\t7\n15\n\n"; + /// One `+/-` block. + const MINUS_ONE_BLOCK: &str = "chain 250 chr2 100 + 0 4 qry2 20 - 3 7 11\n4\n\n"; + /// Several `+/-` blocks with unequal gaps on the `-` query. + const MINUS_MANY_BLOCKS: &str = + "chain 777 chr2 100 + 0 15 qry2 60 - 5 19 13\n4\t1\t2\n3\t2\t0\n5\n\n"; + + fn swap_text(input: &str) -> String { + let mut stdin = Cursor::new(input.as_bytes()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + run( + SwapArgs { + chain: None, + out_chain: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect("swap run"); + String::from_utf8(stdout).expect("utf8 output") + } + + #[test] + fn parses_minimal_args() { + let cli = SwapHarness::try_parse_from(["chaintools", "--chain", "in.chain"]) + .expect("swap arguments should parse"); + assert_eq!(cli.args.chain, Some(PathBuf::from("in.chain"))); + assert!(cli.args.out_chain.is_none()); + assert!(!cli.args.gzip); + } + + #[test] + fn swaps_single_plus_block() { + assert_eq!( + swap_text(PLUS_ONE_BLOCK), + "chain 100 qry1 500 + 40 60 chr1 1000 + 10 30 7\n20\n\n" + ); + } + + #[test] + fn swaps_many_plus_blocks_exchanging_gaps() { + // Block order is preserved; each block's dt/dq are exchanged. + assert_eq!( + swap_text(PLUS_MANY_BLOCKS), + "chain 4200 qry1 500 + 40 93 chr1 1000 + 10 63 9\n20\t1\t3\n10\t7\t5\n15\n\n" + ); + } + + #[test] + fn swaps_single_minus_block_from_sizes() { + // New target = query on +: 20 - 7 = 13 .. 20 - 3 = 17. + // New query = target on -: 100 - 4 = 96 .. 100 - 0 = 100. + assert_eq!( + swap_text(MINUS_ONE_BLOCK), + "chain 250 qry2 20 + 13 17 chr2 100 - 96 100 11\n4\n\n" + ); + } + + #[test] + fn swaps_many_minus_blocks_reversing_block_order_and_gaps() { + // Blocks reverse to 5, 3, 4; the gap after a reversed block is the one + // recorded on its successor in the reversed list, with dt/dq exchanged. + assert_eq!( + swap_text(MINUS_MANY_BLOCKS), + "chain 777 qry2 60 + 41 55 chr2 100 - 85 100 13\n5\t0\t2\n3\t2\t1\n4\n\n" + ); + } + + #[test] + fn double_swap_is_the_identity() { + for fixture in [ + PLUS_ONE_BLOCK, + PLUS_MANY_BLOCKS, + MINUS_ONE_BLOCK, + MINUS_MANY_BLOCKS, + ] { + assert_eq!(swap_text(&swap_text(fixture)), fixture); + } + } + + #[test] + fn preserves_metadata_lines_and_chain_order() { + let input = format!("#meta\n{PLUS_ONE_BLOCK}{MINUS_ONE_BLOCK}"); + let swapped = swap_text(&input); + assert!(swapped.starts_with("#meta\n")); + assert_eq!(swap_text(&swapped), input); + } + + #[test] + fn rejects_blocks_that_do_not_sum_to_the_span() { + let mut stdin = Cursor::new(b"chain 1 chr1 100 + 0 10 qry1 100 + 0 10 1\n4\n\n".to_vec()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + SwapArgs { + chain: None, + out_chain: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("short block list should be rejected"); + assert!(err.to_string().contains("do not sum to the target span")); + } + + #[test] + fn rejects_minus_target_strand() { + let mut stdin = Cursor::new(b"chain 1 chr1 100 - 0 4 qry1 100 + 0 4 1\n4\n\n".to_vec()); + let mut stdout = Vec::new(); + let mut stderr = Vec::new(); + let err = run( + SwapArgs { + chain: None, + out_chain: None, + gzip: false, + }, + &mut stdin, + &mut stdout, + &mut stderr, + ) + .expect_err("minus target strand should be rejected"); + assert!(err.to_string().contains("target strand is not +")); + } +} diff --git a/src/seq/sequence.rs b/src/seq/sequence.rs index 8d14947..e6b7b9c 100644 --- a/src/seq/sequence.rs +++ b/src/seq/sequence.rs @@ -245,6 +245,36 @@ impl SequenceResolver { } } + /// Returns every preloaded sequence name with its length, sorted by name. + /// + /// Formats that must declare their sequences before any record is written + /// (a BAM `@SQ` header, for instance) need the whole set up front, while + /// chains are still being streamed. The sort makes that declaration + /// deterministic across runs, which a `HashMap` iteration order is not. + /// + /// # Output + /// + /// Returns `Ok(Vec<(name, length)>)`, or `Err(ChainError)` if a length + /// exceeds the `u32` coordinate space or the source is not preloaded. + pub fn sequences(&self) -> Result, ChainError> { + match &self.source { + SequenceSource::Loaded { sequences, .. } => { + let mut entries = Vec::with_capacity(sequences.len()); + for (name, sequence) in sequences.iter() { + let len = u32::try_from(sequence.len()).map_err(|_| { + sequence_error("sequence length exceeds u32 coordinate space") + })?; + entries.push((name.as_slice(), len)); + } + entries.sort_unstable_by_key(|(name, _)| *name); + Ok(entries) + } + SequenceSource::TwoBit(_) => { + Err(sequence_error("sequences() requires preloaded sequences")) + } + } + } + /// Returns the full decoded sequence length. pub fn sequence_len(&self, seq_name: &[u8]) -> Result { let len = self.sequence(seq_name)?.len(); From 22f9eca0bf8a40712d1e8895ed1b6cf5be69168d Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Sat, 8 Aug 2026 12:28:38 +0200 Subject: [PATCH 2/4] docs: add convert, liftover and swap tool docs --- assets/tools/convert.md | 88 ++++++++++++++++++++++++++++ assets/tools/liftover.md | 123 +++++++++++++++++++++++++++++++++++++++ assets/tools/swap.md | 71 ++++++++++++++++++++++ 3 files changed, 282 insertions(+) create mode 100644 assets/tools/convert.md create mode 100644 assets/tools/liftover.md create mode 100644 assets/tools/swap.md diff --git a/assets/tools/convert.md b/assets/tools/convert.md new file mode 100644 index 0000000..2c12bb8 --- /dev/null +++ b/assets/tools/convert.md @@ -0,0 +1,88 @@ +

+

+ +

+ + +

+ chaintools +

+ + +

+ + Version Badge + + + Crates.io Version + + + GitHub License + + + Crates.io Total Downloads + +

+

+ + +# chaintools convert + +Converts chain alignments to VCF, MAF, or BAM. + +A chain records alignment geometry but no bases, so both genomes are always required. + +## Usage + +```bash +chaintools convert \ + --to vcf|maf|bam \ + --reference target.2bit-or-fasta \ + --query query.2bit-or-fasta \ + [--chain input.chain] \ + [--output out.vcf] \ + [--gzip] +``` + +If `--chain` is omitted, input is read from standard input. If `--output` is omitted, output is written to standard output. The output type is decided by `--to`, so there is one generic `--output` rather than a per-format flag. + +`--reference` and `--query` accept `.2bit`, `.fa`, `.fasta`, `.fna`, and gzipped FASTA variants. + +`--to bam` requires a build with the `bam` feature (included in `--all-features`), and rejects `--gzip` because BAM output is already block-compressed. + +## VCF + +The chain target is the REF genome and the chain query is the ALT genome. + +- Each mismatching aligned base becomes one record. Soft-masking is ignored (case-insensitive comparison) and a difference involving a non-`ACGT` base is not called. +- `dt > 0, dq == 0` becomes a deletion, `dt == 0, dq > 0` an insertion, both anchored on the preceding reference base — or on the following one when the event sits at the start of the chain's target span, so no coordinate underflows. +- `dt > 0 && dq > 0` becomes one replacement record carrying the whole reference gap against the whole query gap: a chain defines no base-level alignment inside a dual-sided gap and none is invented. +- Redundant shared suffix then prefix bases are trimmed while keeping both alleles non-empty. Repeat-aware left alignment is *not* performed; run `bcftools norm` if you need it. +- `INFO` carries `CHAIN_ID`, `QUERY_CHROM`, `QUERY_POS`, and `STRAND`. `QUAL` is `.` and `FILTER` is `PASS`. There are no sample or genotype columns. +- Records follow chain traversal order; nothing is buffered to coordinate-sort. + +## MAF + +Pairwise MAF, one stanza per continuously aligned segment. + +- Aligned blocks and one-sided gaps extend the current stanza; a one-sided gap is emitted as sequence against `-` on the other row. +- A dual-sided gap ends the stanza and the next aligned block starts a new one, for the same reason as above. +- Sequence text is verbatim, so soft-masking is preserved and stripping `-` reproduces the declared source interval exactly. +- The target row is `+`; the query row keeps the chain's strand, and its `start`/`size` are in that orientation, as MAF requires. + +## BAM + +Each emitted record is an assembly-alignment segment, not a sequencing read. + +- `@SQ` entries come from the reference genome, sorted by name: a BAM header must be complete before the first record, and the chain stream may be an unrewindable pipe. Every chain's declared target size is checked against that header, so contradictory sizes are rejected rather than silently accepted. +- `QNAME` is `query.chain_id`, suffixed `query.chain_id.1`, `.2`, … when a chain is split at a dual-sided gap. Distinct names mean no supplementary-alignment bookkeeping is needed. +- `CIGAR` uses `=`/`X` (the exact bases are known) plus `I`/`D` for one-sided chain gaps. +- `FLAG` carries the reverse bit for a `-` query and `SEQ` is the oriented query sequence. `MAPQ` is 0 and `QUAL` is absent. No `@RG`, sample names, read groups, or platform metadata are fabricated. + +## Errors + +Conversion fails clearly, rather than emitting partial or wrong records, on a missing reference or query sequence, chain coordinates outside a sequence, blocks that do not sum to the header span, an output path equal to the input chain, and — for BAM — an unavailable `bam` feature or a chromosome length that contradicts the reference. + +## Non-goals + +No liftover, PAF, SAM, AXT, PSL, BCF, tabix, bgzip-specific VCF, genotypes, full variant normalization, or local realignment inside chain gaps. diff --git a/assets/tools/liftover.md b/assets/tools/liftover.md new file mode 100644 index 0000000..e9e3db3 --- /dev/null +++ b/assets/tools/liftover.md @@ -0,0 +1,123 @@ +

+

+ +

+ + +

+ chaintools +

+ + +

+ + Version Badge + + + Crates.io Version + + + GitHub License + + + Crates.io Total Downloads + +

+

+ + +# chaintools liftover + +Maps BED intervals through a chain, in one direction: + +```text +chain reference/target → chain query +``` + +No genome sequence is required: a chain carries all the geometry. + +## Usage + +```bash +chaintools liftover \ + --chain map.chain \ + --bed input.bed \ + --type 3|6|12 \ + [--output lifted.bed] \ + [--unmapped unmapped.bed] \ + [--min-match 0.95] \ + [--multiple] +``` + +If `--output` is omitted, output goes to standard output. `--type` is required so no input column is silently ignored; columns beyond the declared width are copied through unlifted and a warning says so. + +### Positive-strand example + +```bash +$ cat plus.chain +chain 100 chr1 1000 + 100 200 qry1 2000 + 500 600 1 +100 + +$ printf 'chr1\t120\t150\n' | tee in.bed +$ chaintools liftover --chain plus.chain --bed in.bed --type 3 +qry1 520 550 +``` + +### Reverse-strand example + +The chain's query is on `-`, so coordinates flip around the query chromosome size and a BED6 strand flips with them: + +```bash +$ cat minus.chain +chain 200 chr1 1000 + 100 200 qry1 2000 - 500 600 2 +100 + +$ printf 'chr1\t120\t150\tfeature\t742\t+\n' > in6.bed +$ chaintools liftover --chain minus.chain --bed in6.bed --type 6 +qry1 1450 1480 feature 0 - +``` + +### Reverse direction + +There is no `--reverse`. Swap the chain first, which keeps one mapping direction in the engine: + +```bash +chaintools swap --chain aToB.chain --out-chain bToA.chain +chaintools liftover --chain bToA.chain --bed on_b.bed --type 3 --output on_a.bed +``` + +## Supported BED widths + +`BED3`, `BED6`, and `BED12`. Any other width is rejected. + +## Mapping rules + +- **`--min-match`** (default `0.95`) is the fraction of *input bases that land inside aligned chain blocks*. A chain's bounding span never counts, and at least one mapped base is always required — even at `--min-match 0`. +- An interval may cross chain gaps. BED3/BED6 still produce **one** output interval covering the mapped extent, not one record per chain block. +- **Mappings are never stitched across chains.** If chain A places 60% and chain B the other 40%, that is not a 100% mapping; each chain is judged on its own. +- **BED12** projects every block through the *same* chain, requires each block to contribute mapped sequence, and measures coverage over block bases only — intronic bases never count. Output blocks are re-sorted ascending by genomic coordinate (needed for `-` chains), `blockStarts` are recomputed, and `thickStart`/`thickEnd` are mapped through the same chain rather than copied. + +## Multiple mappings + +By default a record that several chains can independently place is **not** silently assigned to one of them; it is reported unmapped as `Duplicated in new`. `--multiple` emits every qualifying mapping instead, ordered deterministically by mapped bases, then chain score, query chromosome, query start, chain ID, and finally chain file order. + +## Unmapped records + +`--unmapped PATH` writes each rejected record preceded by its reason: + +```text +Deleted in new no chain aligns any base of the interval +Partially deleted in new one relevant chain, below --min-match +Split in new several chains touch it, none sufficient alone +Duplicated in new several chains qualify and --multiple is off +``` + +Without `--unmapped`, rejected records are dropped — never silently: the summary always reports `unmapped=N` and a warning recommends the flag when `N > 0`. + +## Known limitation + +BED `score` and `itemRgb` are not preserved: the pinned `genepred` model does not carry them, so lifted records emit `0` and `0,0,0`. Coordinates, name, strand, thick region, and blocks all round-trip. This matches the existing `chaintools bed` output behavior. + +## Not implemented + +GFF/GTF/genePred/VCF/PSL input, position/interbase coordinates, `--min-blocks`, `--fudgeThick`, sequence validation, and automatic direction inversion. diff --git a/assets/tools/swap.md b/assets/tools/swap.md new file mode 100644 index 0000000..5c7d963 --- /dev/null +++ b/assets/tools/swap.md @@ -0,0 +1,71 @@ +

+

+ +

+ + +

+ chaintools +

+ + +

+ + Version Badge + + + Crates.io Version + + + GitHub License + + + Crates.io Total Downloads + +

+

+ + +# chaintools swap + +Swaps the target and query sides of every chain, so the same alignment is described in the opposite direction. + +## Usage + +```bash +chaintools swap \ + [--chain input.chain] \ + [--out-chain output.chain] \ + [--gzip] +``` + +If `--chain` is omitted, input is read from standard input. If `--out-chain` is omitted, output is written to standard output. + +```bash +cat input.chain | chaintools swap > swapped.chain +``` + +## Behavior + +- The old query becomes the new target, always on `+`; the old target becomes the new query and keeps the original query strand. +- For a `+` query, block order is preserved and each block's `dt`/`dq` are exchanged. +- For a `-` query, both new header intervals are derived from the chromosome sizes and block order is reversed. Because gaps sit *between* blocks, the gap following a reversed block is the one recorded on its successor in the reversed list. +- Score and chain ID are preserved. +- Output order matches input order; nothing is sorted implicitly. Pipe into `chaintools sort` if you need sorted output. +- `#` metadata lines are preserved in stream order. +- Chains are streamed one at a time, so memory does not grow with input size. + +## Validation + +A chain is rejected, rather than silently mistranslated, when: + +- its target strand is not `+`; +- its block sizes and inner gaps do not sum to the target span; +- its block sizes and inner gaps do not sum to the query span; +- a header coordinate exceeds its declared sequence size. + +## Differences from UCSC `chainSwap` + +- Uses descriptive named arguments instead of positional UCSC arguments. +- Supports stdin/stdout and optional gzip-compressed output. +- Validates block/span consistency instead of assuming it. From 82fd52331c1d6f8364be78f219999962a4de2023 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Sat, 8 Aug 2026 12:28:38 +0200 Subject: [PATCH 3/4] bench: add swap and liftover benchmarks, --only flag, BED workload --- bench/benchmark.rs | 114 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 106 insertions(+), 8 deletions(-) diff --git a/bench/benchmark.rs b/bench/benchmark.rs index bbd3c35..257b217 100644 --- a/bench/benchmark.rs +++ b/bench/benchmark.rs @@ -18,6 +18,32 @@ const MERGE_SORT: [&str; 2] = [ "target/release/chaintools split -c {chain} --files 10 -o . && ./bench/chainMergeSort chains/* > out.chain && rm -rf chains", ]; +const SWAP: [&str; 2] = [ + "target/release/chaintools swap -c {chain} -o out.chain", + "./bench/chainSwap {chain} out.chain", +]; + +// Both tools default to a 0.95 minimum match and to rejecting records that +// several chains can place, so no flags are needed to line them up. `bench.bed` +// is built once by `prepare_liftover_bed`; regenerating it per run would time +// the setup instead of the liftover. +const LIFTOVER: [&str; 2] = [ + "target/release/chaintools liftover -c {chain} -b bench.bed --type 3 -o out.bed -u unmapped.bed", + "./bench/liftOver bench.bed {chain} out.bed unmapped.bed", +]; + +/// BED workload for the liftover benchmark. +const LIFTOVER_BED: &str = "bench.bed"; + +/// Records which chain the BED workload was built from. +/// +/// A BED built from one chain but benchmarked against another would silently +/// measure nothing useful, so the source is written next to it and checked. +const LIFTOVER_BED_SOURCE: &str = "bench.bed.source"; + +/// Keep every Nth aligned chain block when building the liftover workload. +const LIFTOVER_BED_STRIDE: u64 = 10; + /// Command line arguments for the benchmark utility. /// /// This struct defines configuration options for running performance benchmarks @@ -35,6 +61,18 @@ pub struct Args { help = "Extra arguments to pass to hyperfine" )] hyperfine_args: Vec, + + /// Run only the named benchmarks instead of all of them. + /// + /// Each benchmark needs its UCSC baseline binary present in `bench/`, so + /// this allows running the ones whose baseline is installed. + #[clap( + long = "only", + value_delimiter = ',', + num_args = 1.., + help = "Only run the named benchmarks (sorting, filter, merge_sort, swap, liftover)" + )] + only: Vec, } /// Configuration for hyperfine benchmark execution. @@ -164,26 +202,40 @@ fn benchmark() { let args = Args::parse(); std::fs::create_dir_all("runs").unwrap_or_else(|e| panic!("{}", e)); + prepare_liftover_bed(&args.chain); + + // (input, commands, run name, warmup, min runs, max runs). `swap` and + // `liftover` read every chain of a multi-hundred-megabyte input on each run, + // so they get fewer repetitions than the lighter benchmarks. let triplets = vec![ - (&args.chain, SORTING.to_vec(), "sorting"), - (&args.chain, FILTER.to_vec(), "filter"), - (&args.chain, MERGE_SORT.to_vec(), "merge_sort"), + (&args.chain, SORTING.to_vec(), "sorting", 3, 3, 10), + (&args.chain, FILTER.to_vec(), "filter", 3, 3, 10), + (&args.chain, MERGE_SORT.to_vec(), "merge_sort", 3, 3, 10), + (&args.chain, SWAP.to_vec(), "swap", 3, 3, 10), + (&args.chain, LIFTOVER.to_vec(), "liftover", 3, 3, 10), ]; - for (chain, tools, run_name) in triplets { + for (chain, tools, run_name, warmup, min_runs, max_runs) in triplets { + if !args.only.is_empty() && !args.only.iter().any(|name| name == run_name) { + continue; + } let csv = format!("bench_{}.csv", run_name); let md = format!("bench_{}.md", run_name); #[allow(clippy::needless_update)] let code = HyperfineCall { - warmup: 3, - min_runs: 3, - max_runs: Some(10), + warmup, + min_runs, + max_runs: Some(max_runs), export_csv: Some(format!("runs/{}", csv).to_string()), export_markdown: Some(format!("runs/{}", md).to_string()), parameters: vec![("chain".to_string(), vec![chain.to_string()])], setup: Some("cargo build --release --all-features".to_string()), - cleanup: Some(format!("rm -rf output {} chains", STDOUT)), + // Deliberately does not remove the liftover BED: that is an input. + cleanup: Some(format!( + "rm -rf output {} chains out.bed unmapped.bed", + STDOUT + )), commands: tools .iter() .map(|cmd| cmd.to_string()) @@ -205,6 +257,52 @@ fn benchmark() { } } +/// Builds the BED workload for the liftover benchmark, once. +/// +/// The intervals are the chain's own aligned blocks, sampled every +/// [`LIFTOVER_BED_STRIDE`] blocks, so the records genuinely map and the +/// benchmark measures candidate lookup plus mapping rather than the unmapped +/// path. Whole chain spans would be the wrong workload: a multi-megabase span +/// cannot meet a 0.95 minimum match. Skipped when the file already exists, so +/// repeated benchmark runs reuse the same input. +fn prepare_liftover_bed(chain: &str) { + let built_from_this_chain = std::fs::read_to_string(LIFTOVER_BED_SOURCE) + .is_ok_and(|source| source.trim() == chain.trim()); + if std::path::Path::new(LIFTOVER_BED).exists() && built_from_this_chain { + eprintln!("reusing existing {LIFTOVER_BED}"); + return; + } + + build_release(); + eprintln!("building {LIFTOVER_BED} from {chain}"); + let script = format!( + "target/release/chaintools --level off bed -c {chain} --side reference --type 12 \ + | awk 'BEGIN{{OFS=\"\\t\"}} \ + {{n=split($11,sz,\",\"); split($12,st,\",\"); \ + for(i=1;i<=n;i++) if(sz[i]!=\"\") {{k++; if(k%{LIFTOVER_BED_STRIDE}==0) \ + print $1, $2+st[i], $2+st[i]+sz[i]}}}}' > {LIFTOVER_BED}" + ); + let status = Command::new("sh") + .arg("-c") + .arg(&script) + .status() + .expect("Failed to build the liftover BED workload"); + assert!(status.success(), "building {LIFTOVER_BED} failed"); + std::fs::write(LIFTOVER_BED_SOURCE, chain).expect("Failed to record the BED source chain"); +} + +/// Builds the release binaries the benchmark commands invoke. +fn build_release() { + let status = Command::new("cargo") + .args(["build", "--release", "--all-features"]) + .status() + .expect("Failed to run cargo build"); + assert!( + status.success(), + "cargo build --release --all-features failed" + ); +} + /// Main entry point for the benchmark utility. /// /// This function parses command line arguments, runs the benchmark suite, From 1ad221b17003eddd0b0d2d743346a47b6ca2c053 Mon Sep 17 00:00:00 2001 From: alejandrogzi Date: Sat, 8 Aug 2026 12:28:38 +0200 Subject: [PATCH 4/4] chore: v0.0.10 changelog, README bump and gitignore --- .gitignore | 12 +++++++++ README.md | 18 +++++++++---- assets/changelog/Changelog.md | 48 +++++++++++++++++++++++++++++++++++ 3 files changed, 73 insertions(+), 5 deletions(-) diff --git a/.gitignore b/.gitignore index ea8c4bf..fa5842c 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,13 @@ /target + +# UCSC kent binaries used as benchmark baselines (see bench/benchmark.rs) +/bench/chain* +/bench/liftOver + +# hyperfine exports and benchmark scratch files +/runs +/bench.bed +/out.chain +/out.bed +/unmapped.bed +/bench.bed.source diff --git a/README.md b/README.md index 50268e9..e701d17 100644 --- a/README.md +++ b/README.md @@ -5,13 +5,13 @@

- chaintools + chaintools

- - Version Badge + + Version Badge Crates.io Version @@ -29,7 +29,7 @@ work with .chain files in Rust

- docs . + docs . usage . tools . chains @@ -60,7 +60,7 @@ Add this to your `Cargo.toml`: ```toml [dependencies] -chaintools = { version = "0.0.9", features = ["mmap", "gzip", "parallel"] } +chaintools = { version = "0.0.10", features = ["mmap", "gzip", "parallel"] } ``` ## Benchmarks @@ -76,7 +76,15 @@ chaintools = { version = "0.0.9", features = ["mmap", "gzip", "parallel"] } || | `chaintools merge` | 11.995 ± 0.045 | 11.943 | 12.024 | 1.00 | | `UCSC chainMergeSort` | 17.616 ± 0.071 | 17.557 | 17.696 | 1.47 ± 0.01 | +|| +| `chaintools swap` | 8.371 ± 0.163 | 8.183 | 8.475 | 1.00 | +| `UCSC chainSwap` | 14.152 ± 0.119 | 14.021 | 14.254 | 1.69 ± 0.04 | +|| +| `chaintools liftover` | 4.380 ± 0.009 | 4.372 | 4.390 | 1.00 | +| `UCSC liftOver` | 16.540 ± 0.121 | 16.424 | 16.666 | 3.78 ± 0.03 | *ran using hyperfine 1.18.0 on an AMD Ryzen 7 5700X with 128 GB of RAM and 16 cores with a 262 MB chain.gz file as input + +**`swap` and `liftover` were measured on the same machine with the same methodology, on their own inputs: `swap` on a 344 MB chain.gz (hg38 vs HLrhiRex1, 3,049,608 chains) and `liftover` on a 6.5 MB `.over.chain.gz` (16,278 chains) lifting 244,296 BED3 intervals sampled from its aligned blocks. Both commands produced byte-identical output to their UCSC counterpart on every input tested. `bed`, `split` and `convert` have no equivalent UCSC tool to compare against, and `score`/`antirepeat`/`convert` require genome sequence, so they are not listed. diff --git a/assets/changelog/Changelog.md b/assets/changelog/Changelog.md index 649b2b8..cb7269c 100644 --- a/assets/changelog/Changelog.md +++ b/assets/changelog/Changelog.md @@ -5,6 +5,53 @@ All notable changes to **chaintools** are documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.0.10] - 2026-08-08 + +### Added + +- **`swap` subcommand** (UCSC `chainSwap`-compatible) — swaps the target and + query sides of every chain. Metadata/comment lines pass through unchanged, + gzip output is supported (`-G`), and chains can be read from standard input + with output to standard output. Output verified byte-identical to UCSC + `chainSwap` on a 344 MB hg38 vs HLrhiRex1 input. + +- **`liftover` subcommand** (UCSC `liftOver`-compatible) — maps BED3/BED6/BED12 + intervals through a chain file, handling both query strands, gaps, and + multiple chains. Options include `--min-match` (default 0.95), `--multiple` + (emit every qualifying mapping instead of rejecting multi-chain records), + and `--unmapped` (unlifted records with a reason comment). Lookup uses the + new `rust-lapper` dependency over preloaded blocks. Output verified + byte-identical to UCSC `liftOver` on `.over.chain.gz` inputs. + +- **`convert` subcommand** — converts chains plus reference/query sequence + (`.2bit`, `.fa`, gzip variants) into VCF, MAF, or BAM (`--to vcf|maf|bam`). + BAM output is gated behind a new optional `bam` cargo feature + (`noodles-bam`/`noodles-core`/`noodles-sam`), so library and non-BAM CLI + builds do not pull the BAM writer. + +- **`SequenceResolver::sequences()`** — returns every preloaded sequence name + with its length, sorted by name, so formats that must declare sequences up + front (a BAM `@SQ` header) can while chains are still streamed. + +- **Shared CLI helpers** — `validate_block_spans` (a chain's blocks must sum + exactly to both header spans, before any coordinate transformation) and + `validate_distinct_paths` (output must not resolve to the input chain path), + shared by `swap` and `liftover`. + +### Changed + +- **Benchmark suite** — adds `swap` and `liftover` benchmarks against their + UCSC baselines (1.69x and 3.78x faster respectively), a `--only ` + flag to run a subset, and a once-built BED workload for `liftover` sampled + from the chain's own aligned blocks (with a `.source` stamp so it is not + reused against a different chain). + +### Documentation + +- Tool docs for the three new subcommands in `assets/tools/` + (`convert.md`, `liftover.md`, `swap.md`). +- README bumped to 0.0.10 with benchmark tables and methodology notes. + ## [0.0.9] - 2026-07-26 ### Added @@ -249,6 +296,7 @@ First stable release. - Automatic gzip (`.chain.gz`) detection and decompression (`gzip` feature). - Test suite, benchmark binary, CI workflows, Docker image, and rustdoc documentation. +[0.0.10]: https://github.com/alejandrogzi/chaintools/releases/tag/v0.0.10 [0.0.9]: https://github.com/alejandrogzi/chaintools/releases/tag/v0.0.9 [0.0.8]: https://github.com/alejandrogzi/chaintools/releases/tag/v0.0.8 [0.0.7]: https://github.com/alejandrogzi/chaintools/releases/tag/v0.0.7