From 7c561a7cdd5aa5eb721bafb5062cc587dbb45c78 Mon Sep 17 00:00:00 2001 From: AndreaGuarracino Date: Fri, 24 Jul 2026 17:35:43 +0200 Subject: [PATCH 1/3] Accept AGC archives as sequence input -s/--seqs now accepts an AGC (.agc) archive in addition to FASTA/FASTQ. SeqIndex ingests every sample/contig via ragc-core into the same concatenated store and FM-index used for FASTA, so contig names round-trip (a PanSN archive reproduces sample#haplotype#contig verbatim). Ingestion of archives above 32 MB is parallelized across -t threads via disjoint writes to the memory-mapped sequence file; output is byte-identical to the serial path and to the equivalent FASTA input. --- Cargo.lock | 285 +++++++++++++++++++++++++++++++++++++++++++++++- Cargo.toml | 3 + src/main.rs | 2 +- src/seqindex.rs | 233 ++++++++++++++++++++++++++++++++++++++- 4 files changed, 518 insertions(+), 5 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 47bb040..3772347 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -8,6 +8,19 @@ version = "2.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" +[[package]] +name = "ahash" +version = "0.8.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a15f179cd60c4584b8a8c596927aadc462e27f2ca70c04e0071964a73ba7a75" +dependencies = [ + "cfg-if", + "getrandom 0.3.4", + "once_cell", + "version_check", + "zerocopy", +] + [[package]] name = "aho-corasick" version = "1.1.4" @@ -17,6 +30,12 @@ dependencies = [ "memchr", ] +[[package]] +name = "allocator-api2" +version = "0.2.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "683d7910e743518b0e34f1186f92494becacb047c7b6bf616c96772180fef923" + [[package]] name = "anes" version = "0.1.6" @@ -91,6 +110,15 @@ version = "1.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" +[[package]] +name = "bincode" +version = "1.3.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b1f45e9417d87227c7a56d22e471c6206462cba514c7590c09aff4cf6d1ddcad" +dependencies = [ + "serde", +] + [[package]] name = "bitflags" version = "2.13.1" @@ -109,6 +137,15 @@ dependencies = [ "wyz", ] +[[package]] +name = "block-buffer" +version = "0.10.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3078c7629b62d3f0439517fa394996acacc5cbc91c5a20d8c658e77abd503a71" +dependencies = [ + "generic-array", +] + [[package]] name = "block-pseudorand" version = "0.1.2" @@ -125,6 +162,12 @@ version = "3.20.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649" +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + [[package]] name = "cast" version = "0.3.0" @@ -157,6 +200,8 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c89588d05638b5b4594a3348a2d6c20277e43a7f5c5202b05cc56888475a47b8" dependencies = [ "find-msvc-tools", + "jobserver", + "libc", "shlex", ] @@ -248,6 +293,15 @@ version = "1.0.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1d07550c9036bf2ae0c684c4297d503f838287c83c53686d05370d0e139ae570" +[[package]] +name = "cpufeatures" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "59ed5838eebb26a2bb2e58f6d5b5316989ae9d08bab10e0e6d103e656d1b0280" +dependencies = [ + "libc", +] + [[package]] name = "crc32fast" version = "1.5.0" @@ -355,6 +409,40 @@ version = "0.2.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" +[[package]] +name = "crypto-common" +version = "0.1.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78c8292055d1c1df0cce5d180393dc8cce0abec0a7102adb6c7b1eef6016d60a" +dependencies = [ + "generic-array", + "typenum", +] + +[[package]] +name = "dashmap" +version = "6.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6361d5c062261c78a176addb82d4c821ae42bed6089de0e12603cd25de2059c" +dependencies = [ + "cfg-if", + "crossbeam-utils", + "hashbrown 0.14.5", + "lock_api", + "once_cell", + "parking_lot_core", +] + +[[package]] +name = "digest" +version = "0.10.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ed9a281f7bc9b7576e61468ba615a66a5c8cfdff42420a70aa82701a3b1e292" +dependencies = [ + "block-buffer", + "crypto-common", +] + [[package]] name = "either" version = "1.16.0" @@ -440,6 +528,28 @@ dependencies = [ "slab", ] +[[package]] +name = "generic-array" +version = "0.14.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85649ca51fd72272d7821adaf274ad91c288277713d9c18820d8499a7ff69e9a" +dependencies = [ + "typenum", + "version_check", +] + +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "libc", + "r-efi 5.3.0", + "wasip2", +] + [[package]] name = "getrandom" version = "0.4.3" @@ -448,7 +558,7 @@ checksum = "300e883d756b2e4ec94e02791f39b04b522276138852cfc41d9fb7e904106099" dependencies = [ "cfg-if", "libc", - "r-efi", + "r-efi 6.0.0", ] [[package]] @@ -462,6 +572,16 @@ dependencies = [ "zerocopy", ] +[[package]] +name = "hashbrown" +version = "0.14.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" +dependencies = [ + "ahash", + "allocator-api2", +] + [[package]] name = "hashbrown" version = "0.17.1" @@ -505,7 +625,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d466e9454f08e4a911e14806c24e16fba1b4c121d1ea474396f396069cf949d9" dependencies = [ "equivalent", - "hashbrown", + "hashbrown 0.17.1", ] [[package]] @@ -540,6 +660,16 @@ version = "1.0.18" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682" +[[package]] +name = "jobserver" +version = "0.1.35" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1c00acbd29eabad4a2392fa0e921c874934dbbf4194312ad20f04a0ed67a3cb3" +dependencies = [ + "getrandom 0.4.3", + "libc", +] + [[package]] name = "js-sys" version = "0.3.103" @@ -618,6 +748,16 @@ dependencies = [ "autocfg", ] +[[package]] +name = "num_cpus" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91df4bbde75afed763b708b7eee1e8e7651e02d97f6d5dd763e89367e957b23b" +dependencies = [ + "hermit-abi", + "libc", +] + [[package]] name = "once_cell" version = "1.21.4" @@ -671,6 +811,12 @@ version = "0.2.17" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a89322df9ebe1c1578d689c92318e070967d1042b512afbe49518723f4e6d5cd" +[[package]] +name = "pkg-config" +version = "0.3.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19f132c84eca552bf34cab8ec81f1c1dcc229b811638f9d283dceabe58c5569e" + [[package]] name = "plotters" version = "0.3.7" @@ -723,6 +869,12 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + [[package]] name = "r-efi" version = "6.0.0" @@ -735,6 +887,46 @@ version = "0.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "dc33ff2d4973d518d823d61aa239014831e521c75da58e3df4840d3f47749d09" +[[package]] +name = "ragc-common" +version = "0.1.1" +source = "git+https://github.com/AndreaGuarracino/ragc?rev=c2d7ac2f40b80bab7b17b3c038f755f061dd5593#c2d7ac2f40b80bab7b17b3c038f755f061dd5593" +dependencies = [ + "ahash", + "anyhow", + "bincode", + "byteorder", + "hashbrown 0.14.5", + "serde", + "thiserror", + "zstd", +] + +[[package]] +name = "ragc-core" +version = "0.1.1" +source = "git+https://github.com/AndreaGuarracino/ragc?rev=c2d7ac2f40b80bab7b17b3c038f755f061dd5593#c2d7ac2f40b80bab7b17b3c038f755f061dd5593" +dependencies = [ + "ahash", + "anyhow", + "bincode", + "byteorder", + "cc", + "crossbeam", + "dashmap", + "flate2", + "hashbrown 0.14.5", + "num_cpus", + "ragc-common", + "rayon", + "rdst", + "serde", + "sha2", + "thiserror", + "zstd", + "zstd-safe", +] + [[package]] name = "rayon" version = "1.12.0" @@ -865,6 +1057,7 @@ dependencies = [ "once_cell", "parking_lot", "portable-atomic", + "ragc-core", "rayon", "rdst", "rustc-hash", @@ -924,6 +1117,17 @@ dependencies = [ "serde", ] +[[package]] +name = "sha2" +version = "0.10.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a7507d819769d01a365ab707794a4084392c824f54a7a6a7862f8c3d0892b283" +dependencies = [ + "cfg-if", + "cpufeatures", + "digest", +] + [[package]] name = "shlex" version = "2.0.1" @@ -999,12 +1203,32 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "32497e9a4c7b38532efcdebeef879707aa9f794296a4f0244f6f69e9bc8574bd" dependencies = [ "fastrand", - "getrandom", + "getrandom 0.4.3", "once_cell", "rustix", "windows-sys", ] +[[package]] +name = "thiserror" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + [[package]] name = "tikv-jemalloc-sys" version = "0.5.4+5.3.0-patched" @@ -1076,6 +1300,12 @@ version = "0.1.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5d99f8c9a7727884afe522e9bd5edbfc91a3312b36a77b5fb8926e4c31a41801" +[[package]] +name = "typenum" +version = "1.20.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6f5e870be6c3b371b77fe0ee0bafb859fa4964b4404c27de1d380043c4dda20" + [[package]] name = "uf_rush" version = "0.1.1" @@ -1103,6 +1333,12 @@ dependencies = [ "serde", ] +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + [[package]] name = "voracious_radix_sort" version = "1.2.0" @@ -1122,6 +1358,15 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "wasip2" +version = "1.0.4+wasi-0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b67efb37e106e55ce722a510d6b5f9c17f083e5fc79afc2badeb12cc313d9487" +dependencies = [ + "wit-bindgen", +] + [[package]] name = "wasm-bindgen" version = "0.2.126" @@ -1210,6 +1455,12 @@ dependencies = [ "memchr", ] +[[package]] +name = "wit-bindgen" +version = "0.57.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e" + [[package]] name = "wyz" version = "0.5.1" @@ -1244,3 +1495,31 @@ name = "zmij" version = "1.0.23" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "29666d0abbfad1e3dc4dcf6144730dd3a3ab225bbbdac83319345b1b44ccfc1b" + +[[package]] +name = "zstd" +version = "0.13.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e91ee311a569c327171651566e07972200e76fcfe2242a4fa446149a3881c08a" +dependencies = [ + "zstd-safe", +] + +[[package]] +name = "zstd-safe" +version = "7.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f49c4d5f0abb602a93fb8736af2a4f4dd9512e36f7f570d66e65ff867ed3b9d" +dependencies = [ + "zstd-sys", +] + +[[package]] +name = "zstd-sys" +version = "2.0.16+zstd.1.5.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91e19ebc2adc8f83e43039e79776e3fda8ca919132d68a1fed6a5faca2683748" +dependencies = [ + "cc", + "pkg-config", +] diff --git a/Cargo.toml b/Cargo.toml index 903cbfd..1f71bc8 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -63,5 +63,8 @@ itoa = "1" # Fast integer-to-decimal formatting for GFA emi # CLI dependencies clap = { version = "4.5", features = ["derive"] } # Command-line argument parsing +# AGC (Assembled Genomes Compressor) input support: pure-Rust decompressor, same crate/rev impg and cigzip use +ragc-core = { git = "https://github.com/AndreaGuarracino/ragc", rev = "c2d7ac2f40b80bab7b17b3c038f755f061dd5593" } + [build-dependencies] cbindgen = "0.27" diff --git a/src/main.rs b/src/main.rs index 19697b8..0606c38 100644 --- a/src/main.rs +++ b/src/main.rs @@ -60,7 +60,7 @@ fn main() -> io::Result<()> { .short('s') .long("seqs") .value_name("FILE") - .help("The sequences used to generate the alignments (FASTA, FASTQ, .seq)") + .help("The sequences used to generate the alignments (FASTA, FASTQ, .seq, .agc)") .required(true)) .arg(Arg::new("gfa") .short('g') diff --git a/src/seqindex.rs b/src/seqindex.rs index e514774..4391333 100644 --- a/src/seqindex.rs +++ b/src/seqindex.rs @@ -1,9 +1,12 @@ use flate2::read::MultiGzDecoder; use fm_index::{FMIndexWithLocate, MatchWithLocate, Search, Text}; -use memmap2::Mmap; +use memmap2::{Mmap, MmapMut}; +use ragc_core::{Decompressor, DecompressorConfig}; +use rayon::prelude::*; use std::fs::{File, OpenOptions}; use std::io::{BufRead, BufReader, BufWriter, Write}; use std::path::PathBuf; +use std::sync::Mutex; /// Simple sparse bitvector for sequence boundaries /// Much faster than RsVec for sparse data (select is O(1) array access vs hierarchical search) @@ -110,6 +113,11 @@ impl SeqIndex { /// Time complexity: O(n log n) for suffix array construction /// Space complexity: O(n log σ) bits for CSA + O(m log(N/m)) bits for boundaries pub fn build_index(&mut self, filename: &str) -> Result<(), String> { + // AGC archives are read through ragc-core rather than the FASTA/FASTQ line parser. + if filename.ends_with(".agc") { + return self.build_index_agc(filename); + } + // Create a unique temp file for sequences using mkstemps so concurrent calls // (e.g. parallel tests or multi-partition builds) never share the same path. // This also honours the configured temp directory (seqwish::tempfile::set_dir or TMPDIR). @@ -252,6 +260,229 @@ impl SeqIndex { // Close sequence file drop(seq_out); + self.finalize( + name_text, + name_boundary_positions, + seq_boundary_positions, + seq_bytes_written, + ) + } + + /// Build the index from an AGC (Assembled Genomes Compressor) archive. + /// + /// Every contig of every sample in the archive is ingested exactly once (in the + /// archive's stored order) into the same concatenated sequence file + FM-index used + /// for FASTA/FASTQ input. Each contig is written under its own name verbatim (first + /// whitespace token, matching the FASTA header rule), so a PanSN archive reproduces + /// the original `sample#haplotype#contig` names 1:1. + fn build_index_agc(&mut self, filename: &str) -> Result<(), String> { + eprintln!("[seqindex] build_index_agc: opening AGC archive {filename}"); + + // Sequence temp file, identical scheme to the FASTA/FASTQ path. + let seq_file = crate::tempfile::create("seqwish", ".sqq") + .map_err(|e| format!("Failed to create sequence temp file: {e}"))?; + self.seq_filename = Some(seq_file.clone()); + + // --- Pass 1: enumerate contigs and lengths from AGC metadata only (get_contig_length + // reads segment descriptors, it does NOT decompress bases), and build the name text + + // boundary tables. This fixes the exact archive order and each contig's byte offset, so + // the concatenated sequence file is byte-identical whether we then fill it serially or + // in parallel. + let mut decompressor = Decompressor::open(filename, DecompressorConfig { verbosity: 0 }) + .map_err(|e| format!("Failed to open AGC archive {filename}: {e}"))?; + + let samples = decompressor.list_samples(); + eprintln!( + "[seqindex] build_index_agc: {} sample(s) in archive", + samples.len() + ); + + // (sample, contig, offset, len) for every non-empty contig, in archive order. + let mut jobs: Vec<(String, String, u64, usize)> = Vec::new(); + let mut name_text = String::new(); + let mut name_boundary_positions: Vec = Vec::new(); + let mut seq_boundary_positions: Vec = Vec::new(); + let mut total_bytes: u64 = 0; + let mut notified_empty_seqs = false; + + for sample in &samples { + let contigs = decompressor + .list_contigs_names_only(sample) + .map_err(|e| format!("Failed to list contigs for sample '{sample}': {e}"))?; + eprintln!( + "[seqindex] build_index_agc: sample '{sample}' has {} contig(s)", + contigs.len() + ); + + for contig in &contigs { + let len = decompressor + .get_contig_length(sample, contig) + .map_err(|e| { + format!("Failed to get length of '{contig}' in '{sample}': {e}") + })?; + + if len == 0 { + if !notified_empty_seqs { + notified_empty_seqs = true; + eprintln!("[seqindex] WARNING: AGC archive contains empty sequences, which will be ignored."); + } + continue; + } + + // Match the FASTA rule: the sequence name is the first whitespace token. + let seq_name = contig.split_whitespace().next().unwrap_or("").to_string(); + + // Record name boundary (position of '>' in concatenated name text) + name_boundary_positions.push(name_text.len() as u64); + name_text.push('>'); + name_text.push_str(&seq_name); + name_text.push(' '); + + // Record sequence boundary; the bases go in during pass 2. + seq_boundary_positions.push(total_bytes); + jobs.push((sample.clone(), contig.clone(), total_bytes, len)); + total_bytes += len as u64; + self.seq_count += 1; + } + } + + // Add final boundary for total length + seq_boundary_positions.push(total_bytes); + + // --- Pass 2: decompress the bases into the concatenated .sqq file at their fixed offsets. + // AGC decompression is the cost here and scales with total bp, so for large archives we + // fan out across threads, each writing a disjoint [offset, offset+len) range of a + // memory-mapped file. Small archives stay on the simple serial writer (no thread overhead, + // and it is the well-worn path for the common per-gene / per-chromosome pggb input). + const PARALLEL_MIN_BYTES: u64 = 32 * 1024 * 1024; + let threads = rayon::current_num_threads(); + + if threads > 1 && total_bytes >= PARALLEL_MIN_BYTES { + eprintln!( + "[seqindex] build_index_agc: parallel decompress of {total_bytes} bp across {threads} thread(s)" + ); + // Pre-size the file, then map it writable. + { + let f = OpenOptions::new() + .create(true) + .write(true) + .truncate(true) + .open(&seq_file) + .map_err(|e| format!("Failed to create sequence file: {e}"))?; + f.set_len(total_bytes) + .map_err(|e| format!("Failed to size sequence file: {e}"))?; + } + let file = OpenOptions::new() + .read(true) + .write(true) + .open(&seq_file) + .map_err(|e| format!("Failed to open sequence file for mmap: {e}"))?; + let mut mmap = unsafe { + MmapMut::map_mut(&file).map_err(|e| format!("Failed to mmap sequence file: {e}"))? + }; + // Base address as usize so it can cross the rayon closure (raw pointers are !Send). + // Each job writes a disjoint [offset, offset+len) slice, so the writes never alias. + let base = mmap.as_mut_ptr() as usize; + + // One decompressor per rayon thread (ragc's Decompressor is !Sync and its readers take + // &mut self); each thread only touches its own slot, so the Mutex never contends. + let n_slots = threads + 1; + let pool: Vec>> = + (0..n_slots).map(|_| Mutex::new(None)).collect(); + let path = filename.to_string(); + + jobs.par_iter().try_for_each( + |(sample, contig, offset, len)| -> Result<(), String> { + let tid = rayon::current_thread_index().unwrap_or(n_slots - 1); + let mut slot = pool[tid].lock().unwrap(); + let dec = slot.get_or_insert_with(|| { + Decompressor::open(&path, DecompressorConfig { verbosity: 0 }) + .unwrap_or_else(|e| { + panic!("Failed to open AGC '{path}' for thread: {e}") + }) + }); + let numeric = dec + .get_contig_range(sample, contig, 0, *len) + .map_err(|e| format!("Failed to read '{contig}' in '{sample}': {e}"))?; + // SAFETY: [offset, offset+len) is disjoint across jobs (prefix-sum offsets) and + // lies within the file sized to total_bytes, so no two threads write the same + // bytes and every byte is written exactly once. + let dst = unsafe { + std::slice::from_raw_parts_mut( + (base as *mut u8).add(*offset as usize), + *len, + ) + }; + for (d, &b) in dst.iter_mut().zip(numeric.iter()) { + *d = match b { + 0 => b'A', + 1 => b'C', + 2 => b'G', + 3 => b'T', + _ => b'N', + }; + } + Ok(()) + }, + )?; + + mmap.flush() + .map_err(|e| format!("Failed to flush sequence file: {e}"))?; + } else { + eprintln!("[seqindex] build_index_agc: serial decompress of {total_bytes} bp"); + let mut seq_out = BufWriter::with_capacity( + 1024 * 1024, + OpenOptions::new() + .create(true) + .write(true) + .truncate(true) + .open(&seq_file) + .map_err(|e| format!("Failed to create sequence file: {e}"))?, + ); + for (sample, contig, _offset, len) in &jobs { + let numeric = decompressor + .get_contig_range(sample, contig, 0, *len) + .map_err(|e| format!("Failed to read '{contig}' in '{sample}': {e}"))?; + let seq: Vec = numeric + .iter() + .map(|&b| match b { + 0 => b'A', + 1 => b'C', + 2 => b'G', + 3 => b'T', + _ => b'N', + }) + .collect(); + seq_out + .write_all(&seq) + .map_err(|e| format!("Failed to write sequence: {e}"))?; + } + drop(seq_out); + } + + eprintln!( + "[seqindex] build_index_agc: ingested {} sequence(s), {total_bytes} bp total", + self.seq_count + ); + + self.finalize( + name_text, + name_boundary_positions, + seq_boundary_positions, + total_bytes, + ) + } + + /// Finalize the index once every record has been written to the sequence temp file: + /// build the FM-index over the concatenated names, the sparse boundary bitvectors, and + /// memory-map the sequence file. Shared by the FASTA/FASTQ and AGC readers. + fn finalize( + &mut self, + name_text: String, + name_boundary_positions: Vec, + seq_boundary_positions: Vec, + seq_bytes_written: u64, + ) -> Result<(), String> { // Build FM-index (CSA) from name text // Space: O(n log σ) bits where n = name_text.len() // FM-index requires text to end with exactly one zero character From 6fce96e65e82d05fd4488dfa0af2f73cbdb9eabb Mon Sep 17 00:00:00 2001 From: AndreaGuarracino Date: Fri, 24 Jul 2026 18:10:49 +0200 Subject: [PATCH 2/3] Add AGC vs FASTA equivalence test Ingesting tests/data/agc/small.agc (built with 'agc create' from tests/data/agc/small.fa) must yield the same names, lengths, and bases as ingesting the FASTA. Runs under cargo test (ragc-core is a normal dep). --- tests/agc_input.rs | 47 +++++++++++++++++++++++++++++++++++++++ tests/data/agc/small.agc | Bin 0 -> 2108 bytes tests/data/agc/small.fa | 37 ++++++++++++++++++++++++++++++ 3 files changed, 84 insertions(+) create mode 100644 tests/agc_input.rs create mode 100644 tests/data/agc/small.agc create mode 100644 tests/data/agc/small.fa diff --git a/tests/agc_input.rs b/tests/agc_input.rs new file mode 100644 index 0000000..0816ea4 --- /dev/null +++ b/tests/agc_input.rs @@ -0,0 +1,47 @@ +// AGC input support: reading sequences from an AGC (.agc) archive must yield exactly the same +// sequence index as reading the FASTA the archive was created from. tests/data/agc/small.agc +// was produced with `agc create` from tests/data/agc/small.fa, so ingesting either must give +// identical names, lengths, and bases. +use seqwish::seqindex::SeqIndex; + +#[test] +fn agc_ingest_matches_fasta() { + let mut fa = SeqIndex::new(); + fa.build_index("tests/data/agc/small.fa") + .expect("build index from FASTA"); + + let mut agc = SeqIndex::new(); + agc.build_index("tests/data/agc/small.agc") + .expect("build index from AGC"); + + assert!(fa.n_seqs() > 0, "fixture should contain sequences"); + assert_eq!(fa.n_seqs(), agc.n_seqs(), "sequence count differs"); + + // Same names and lengths, in the same order. + for i in 1..=fa.n_seqs() { + assert_eq!( + fa.nth_name(i), + agc.nth_name(i), + "sequence name differs at index {i}" + ); + assert_eq!( + fa.nth_seq_length(i), + agc.nth_seq_length(i), + "sequence length differs at index {i}" + ); + } + + // Same concatenated bases, position by position. + assert_eq!( + fa.seq_length(), + agc.seq_length(), + "total concatenated length differs" + ); + for p in 0..fa.seq_length() { + assert_eq!( + fa.at(p), + agc.at(p), + "base differs at concatenated position {p}" + ); + } +} diff --git a/tests/data/agc/small.agc b/tests/data/agc/small.agc new file mode 100644 index 0000000000000000000000000000000000000000..ce48468e351ce9887430ca0970d575f574c9a43b GIT binary patch literal 2108 zcmb7@dr%Wc7{GTg5hH^oAwYOo5(F9_yb{C`MI-@4A}U0RLn~uLa%eQK<^l~`3W&8n zY8Bc-i)b~2;9F2Z%7_}RRPj-z>R27=h*m4w5fzY@R=T;|VZ^_BbIJVfxBGqj?d^Av zKcO+=Uje@mfM+&rm6+{2s%Exk#=;l?Fic90m6tIAs~|pyJNC+pX1-q0-tzL3E1qRt zv0Ypj!J^3e1nd5+HeVL=%l5e9(um!ST`;80H8K06!z$1ILZ0HGJnCERo`9{EZU4Me zIK#K1|8_`|VTMocg2%@PUfeI?h87=uzGnRXnJ|vELx;_4Vu2uD^P683wjm8Td7z=$CZ`@q(o7o0u$^U$VX zExD=XSqvbFz9CEqbbUI=tv<4VkYAkj{$~&PRqEbrp}|R#De|CQJF@#@?>`;4U_+F) z*fenzD?8aL@H~8YxowNf>1cwjOFSj(?2F}y3{~eEbc~j!R(4R+vDgCmtgtEf>Rx?D zy0`G;sz2Xb$B>Pzmcu?KJ9~Ld{#;ZNT`i7|0AeK`6F1hY9@y-qP9RH7RboJvk(r|Q zP=4;dF?42Z@bY8)mCw43ZspS_@#Mn5Jv??%Rf%lq-l03TitT+jE}S1Ycb7j+?U()y z|KZQp2v_&1zXbA^4Xk~*=Ex5q{*789vua8uiTmn2Ue!XigGZ!v_vyMC3p2^y4(Vh$ z6%9-fv$Ih#<{Gi>`mHCKF@Z+|qN>Y7!g@l^mGqoBWG=n5wp?DAS~mHPb;HX~zqskz z`nXv#ZdSu1{~G>e9oMuz%d7qJiMF_fE?L%GuVUe}PTNrK{&nZmS~W%WeX9nhHn(kB z-Mz@AukJv5%)Q&v=^vjtzM`>qTX$Se`1P`+w|*>I#x(b{4i*T@EAR5kJFS}^tY5hD z>YFXWVP{m%@7OuM#Qu^4$l?Q@B4iRuol%QRgG529r6v&=5D$XCWkaGt@LUXli@Teb z8z)%k=0!%hi0lyzbdis4N&wg#53Z+|H_yk{&)*@@>tUm5!2cvhir_u~>TO@ov~3!+ z|D(H1R35hJ79tD?R2SrL4<`c<+HZR_1Vuyt@8+X8PxM1&gAcn=4p#4#vD}>a0+w4v zKxg{Nv@xlbYq|@q?ZG{+feFGJ(URu;Q@7h=^%#hvd@T(2# z+FJ<6j?BP8x#=u)l%JCuqAp+LMy_J7ZR83(4<@-$ktKpZKaA`*$C!$sL3SEI!j=i??qW6&?wDVG^cAZmov z=#h*BMw3CEr@~D@xkNSmFdRx|2Cp)|IfG74=NS&wkqd~?p2(gXEkWBObV{ujzzCH= ztHo7jR9u+jn;`7Hh>hxoY^NcWL2uS9aSUO9K%vu+8aFF7TE`HB16Nui)no?7Mu=HN z6mLc>ruZwwQ507oo5GPT56!9#ITM#Fk$gdIlT9VbM5iGOGXXp>BEwf3H;t(xpMhYq94si_$IVt8i zL=Vif%ppZ8DKNMcO-UuqT*{O}Z_8ZL%%jW}m}!|un)#IJMN#ufljbm(J2GjcNvF&t zxZRRY8UsampleA#1#chr1 +AAGCCCAATAAACCACTCTGACTGGCCGAATAGGGATATAGGCAACGACATGTGCGGCGA +CCCTTGCGACAGTGACGCTTTCGCCGTTGCCTAAACCTATTTGAAGGAGTCTAGCAGCCG +CAGTAAGGCACAATACCTCGTCCGTGTTACCAGACCAAACAAGACGTCCTCTTCAATGTT +TAAATGACCCTCTCGTCATAAAACCTTTCTACTATGTGTTCCGCAAGAATCAACAACTAC +AATGGCGCGTCGTGAATAACGCGACGGCTGAGACGAACGGCGCGTGAATGAAGCGCTTAA +ACAGCTCAGGAGCCAGTCCCCTACGTCGCATATCCTGGCCACTGGAGGTGAAGCGAATGG +TATCGATACGTAGGAGGTGTGCCTTCGTAGGCTGTTTCTCAGGACGCCCAACTATTCTTT +CCAATCCTACATCTGTTTCTTGCGTCGTAGCGGGACCCTCCATTGTTACTTATTAGGTTC +TCGTTATGTCTCATAATCTCAGTGCTGGTGTGATAAGCAAACCACCCTACTGGCACGAAG +TTCACAGAAGTGAGATTATGTCTCGTTTGGCAGTCTT +>sampleA#2#chr1 +ATGCTCGGGGGACACTTCTTTAAGCTCGGTGTGGTGGGCACGACCCTGGACGCGCGACGA +AGCTAAGTTTGCAGTAATTAACCGACATCTTTGTGAACCGACCCACATTTGACGGTACGC +TACCGCAACGGTATGTGTTAATGGAACAGACTTGCTTATGTGGACGTTGTATAGGGATAT +TACGTTACGCGTTAACCGATACATACTGGTTTCTCTCCAGTGGAGGTCTTGGTTGCCTCT +AGTTTCTACGATATACTCATGGTAGTGTAACGCATAATCGAAGAGGGTCCTCCCATCTCC +TGTGATGCATGGTGTGCTTACTGGGATGAATGCGCCGCAAGTAGCAGGTCCCGGCGTGGA +TACCTGATAGATGGTGACTAGCATGTACAAGTAACCTTGTCTATTGAGCTTCGAGGATGC +AT +>sampleB#1#chr1 +CAAGCCCACCCGCAGCCGCAACAGCGACGACTAATTGATCAGTAATTTATTAAGCACGGT +GTTAACTTCTGTTTAGTGGGCTAAAATAGCAGATGTAGGGACCTCAGGAGCTAGACGGGG +ACCTACAACTTTGCGGGAACCAAGTTTTTGCAGTAGTGACTAACGCCGGGAATTCCTCGA +TATATAGTTTGATAGCTGATACTTATGGCGCAACGGCCACGCCCACTTTGGCTATTGGAG +AGTTAAGGAATTATCGTCATAGACACTTCGGGTTGAGAGATGGCGACGGTCAGTGC +>sampleB#1#chr2 +ATGAGGCCGTCCCCAGAAGCTCCCCTATGCTGTCCGTCGTTGTTCCCGATGAAGACGTCT +ACTGATATGCTAGCAGAGCCAGTCTTAAAGCCTAGCGAACTTAATACCGTAGCTCAGAAT +TATGGAGAGCAGCAGGCTTCCATAGCACAGGTTGACGGAGGAGTTTTGCTTGGATATCGG +AAGGGTTCTGTAGTGAATGCACTACACGGTACTGGTACGTGGCAACTTAGGTCGTCACAT +CTAGGAGGCCGCACCCTAGGTCAAGTTTTACGATTGCCCTAACGCCGCGGAGCGCGACCC +GAAAAGCTATGGTCTGTAACTTTTCGCGGGTCGAGCTAGTCCAAGTTCCGGCCTTTGTAA +TTCCGAAGTTGAATCGGTGATACGGATTGACATGGGCCTAAACGTTCCGGCTGGTGTAGG +ATGATGCATCTCCAACATGTCTCTTACCGTTGCTGGGTCCGGCGGCTGTGGGATTGCGAG +AGTGTCCGGCACCACCAATGTACACTTTCGGGAACACTCATTCGAAGAGGTTCTGCAGCT +GCAGGCCTTGATACCTGCAGTCTGGGAGGCAATGCTGAGGCCCTCTGTTCCAT From 5994b777d547596d6a84a3b8659f8ad087460dea Mon Sep 17 00:00:00 2001 From: AndreaGuarracino Date: Fri, 24 Jul 2026 19:46:44 +0200 Subject: [PATCH 3/3] Keep AGC records distinct when contig names repeat across samples AGC identifies a sequence by (sample, contig), seqwish by name alone. An archive built the canonical AGC way (agc create sampleA.fa sampleB.fa) repeats plain contig names across samples, which produced several sequences with the same name and made them unresolvable by name. A contig name that occurs in more than one sample is now qualified as contig@sample (AGC's query syntax, matching wfmash); already unique names, as in a PanSN archive, are unchanged. Duplicates that survive are a hard error rather than a silently broken index. --- src/seqindex.rs | 63 ++++++++++++++++++++++++++++------------- tests/agc_input.rs | 51 +++++++++++++++++++++++++++++++++ tests/data/agc/dup.agc | Bin 0 -> 1597 bytes tests/data/agc/dupA.fa | 12 ++++++++ tests/data/agc/dupB.fa | 13 +++++++++ 5 files changed, 120 insertions(+), 19 deletions(-) create mode 100644 tests/data/agc/dup.agc create mode 100644 tests/data/agc/dupA.fa create mode 100644 tests/data/agc/dupB.fa diff --git a/src/seqindex.rs b/src/seqindex.rs index 4391333..f3ab308 100644 --- a/src/seqindex.rs +++ b/src/seqindex.rs @@ -297,12 +297,10 @@ impl SeqIndex { samples.len() ); - // (sample, contig, offset, len) for every non-empty contig, in archive order. - let mut jobs: Vec<(String, String, u64, usize)> = Vec::new(); - let mut name_text = String::new(); - let mut name_boundary_positions: Vec = Vec::new(); - let mut seq_boundary_positions: Vec = Vec::new(); - let mut total_bytes: u64 = 0; + // (sample, contig, short name, len) for every non-empty contig, in archive order. + let mut records: Vec<(String, String, String, usize)> = Vec::new(); + let mut short_name_counts: std::collections::HashMap = + std::collections::HashMap::new(); let mut notified_empty_seqs = false; for sample in &samples { @@ -330,22 +328,49 @@ impl SeqIndex { } // Match the FASTA rule: the sequence name is the first whitespace token. - let seq_name = contig.split_whitespace().next().unwrap_or("").to_string(); - - // Record name boundary (position of '>' in concatenated name text) - name_boundary_positions.push(name_text.len() as u64); - name_text.push('>'); - name_text.push_str(&seq_name); - name_text.push(' '); - - // Record sequence boundary; the bases go in during pass 2. - seq_boundary_positions.push(total_bytes); - jobs.push((sample.clone(), contig.clone(), total_bytes, len)); - total_bytes += len as u64; - self.seq_count += 1; + let short_name = contig.split_whitespace().next().unwrap_or("").to_string(); + *short_name_counts.entry(short_name.clone()).or_insert(0) += 1; + records.push((sample.clone(), contig.clone(), short_name, len)); } } + // (sample, contig, offset, len) for every non-empty contig, in archive order. + let mut jobs: Vec<(String, String, u64, usize)> = Vec::new(); + let mut name_text = String::new(); + let mut name_boundary_positions: Vec = Vec::new(); + let mut seq_boundary_positions: Vec = Vec::new(); + let mut total_bytes: u64 = 0; + let mut seen_names: std::collections::HashSet = std::collections::HashSet::new(); + + for (sample, contig, short_name, len) in records { + // AGC identifies a sequence by (sample, contig), seqwish by name alone, so the + // name must be unique: a contig name that occurs in several samples is qualified + // as "contig@sample" (AGC's own query syntax, matching wfmash). Names that are + // already unique - as in a PanSN archive - are used verbatim. + let seq_name = if short_name_counts[&short_name] > 1 { + format!("{short_name}@{sample}") + } else { + short_name + }; + if !seen_names.insert(seq_name.clone()) { + return Err(format!( + "AGC archive {filename} contains duplicate sequence name '{seq_name}'" + )); + } + + // Record name boundary (position of '>' in concatenated name text) + name_boundary_positions.push(name_text.len() as u64); + name_text.push('>'); + name_text.push_str(&seq_name); + name_text.push(' '); + + // Record sequence boundary; the bases go in during pass 2. + seq_boundary_positions.push(total_bytes); + jobs.push((sample, contig, total_bytes, len)); + total_bytes += len as u64; + self.seq_count += 1; + } + // Add final boundary for total length seq_boundary_positions.push(total_bytes); diff --git a/tests/agc_input.rs b/tests/agc_input.rs index 0816ea4..61cf968 100644 --- a/tests/agc_input.rs +++ b/tests/agc_input.rs @@ -4,6 +4,22 @@ // identical names, lengths, and bases. use seqwish::seqindex::SeqIndex; +/// Read a FASTA into (name, sequence) pairs, applying the same "first whitespace token" +/// name rule the index uses. +fn read_fasta(path: &str) -> Vec<(String, String)> { + let text = std::fs::read_to_string(path).expect("read fixture FASTA"); + let mut out: Vec<(String, String)> = Vec::new(); + for line in text.lines() { + if let Some(header) = line.strip_prefix('>') { + let name = header.split_whitespace().next().unwrap_or("").to_string(); + out.push((name, String::new())); + } else if let Some(last) = out.last_mut() { + last.1.push_str(line.trim()); + } + } + out +} + #[test] fn agc_ingest_matches_fasta() { let mut fa = SeqIndex::new(); @@ -45,3 +61,38 @@ fn agc_ingest_matches_fasta() { ); } } + +/// AGC identifies a sequence by (sample, contig) while seqwish identifies it by name alone. +/// tests/data/agc/dup.agc holds two samples (dupA, dupB) that BOTH contain chr1 and chr2, the +/// canonical AGC layout (`agc create sampleA.fa sampleB.fa`). Such records must stay distinct: +/// each gets a "contig@sample" name and must carry its own sample's bases, not the bases of +/// whichever sample happens to be listed first. +#[test] +fn agc_duplicate_contig_names_stay_distinct() { + let mut idx = SeqIndex::new(); + idx.build_index("tests/data/agc/dup.agc") + .expect("build index from AGC with duplicate contig names"); + + // Expected content, taken from the FASTAs the archive was built from. + let mut expected: Vec<(String, String)> = Vec::new(); + for (sample, path) in [ + ("dupA", "tests/data/agc/dupA.fa"), + ("dupB", "tests/data/agc/dupB.fa"), + ] { + for (name, seq) in read_fasta(path) { + expected.push((format!("{name}@{sample}"), seq)); + } + } + + assert_eq!(idx.n_seqs(), expected.len(), "sequence count differs"); + + for (name, seq) in &expected { + let id = idx + .rank_of_seq_named(name) + .unwrap_or_else(|| panic!("sequence '{name}' not found (names must be unique)")); + let len = idx.nth_seq_length(id).expect("length of sequence"); + assert_eq!(len as usize, seq.len(), "length differs for '{name}'"); + let got = idx.subseq_by_id(id, 0, len).expect("bases of sequence"); + assert_eq!(&got, seq, "bases differ for '{name}' (wrong sample?)"); + } +} diff --git a/tests/data/agc/dup.agc b/tests/data/agc/dup.agc new file mode 100644 index 0000000000000000000000000000000000000000..adaed9f1445beb7d65998f60af0e1ca918092285 GIT binary patch literal 1597 zcma)*eM}o=9LJyMx}+G~Y-uZPz;R>&i-C?SW0PfWyf|13oiRA(WE1*wuu9u2cWo^? zbTevXbDIJ}nYhlyGNf$D0ADsGjsmXn)-U04)c{QQV6RCfrzQxKdal<100xseo56~Kx}pQ=FEP)wcU<*(e2d*@ zpZ%@Rzd7&4I8B7AG)b9JIDbcZawg~@1pUkzGLtGc`zVW1F*uhPLhx=eLhL^=ze5ctLxq>IO7cburLA#W`m2-n-4U^lp5z+{$GCo4exrsAa zH<1W~kUK|Y0U(h^$>a(M!UEtH@1g7gpr?LYPi+FbvrxnC?=JP!8p06}qDD{@24h&Q zGzM%JflsBX$X~_*NCd#AYmqB}94G*B)m(WkxZraS_@Hf5_dlSn_5}D)0ndgf7>t{Qw zPIPtm_WrJ4*r^-XJAI3wraH%)=l_}+TJP|Gbi+1MpA>UXy98o6CVjLi<` zXS7^wY^zrWYYW-$e)azUed3U}9KBlfmQJ8sTXm=Pf#dqV)@ZywH}_In@*@Y^zJ2S% z8?sGzuh%FF9{w>OiUkkBB1N%+cAsBSNE0#$k!&476Z{lN6Y>a}1Ze^bv9e{FU_%ms zVl!p8lMWZi%FDskSqwwkO*Ybk=aCM|PBM57Wp~mf!%%b*?$O|B8m)%Fi%o8;&4O!F zwHpXRo30`DfikO&Eb)*uW2GD=cH;{a4boPKRtFMOfs>{zZZk;(W4Sqe2#0c+!9SV* z>mei;@`OX}TmdbBk=!6h?m?quaXYgB+Oh#$XQpg6((FPFrZ7gk(?&7?4LY`J7nCnE znJI_MS{^dNC>XoUWFcKftBp6MLWeii!exMd$askmV+#0+SIct)#Jp)dHzD4@b2H*} zp8be7@_ZWcCY}R`ALqFhaR$#p#92I_L!8a?dBiz9cOia~=g$yt=lLSyT%NlT=h1Bc z0IIxs7Bmdr9o!klz}?;*oS+K{Zseik1bs+&2MKyk$PWoYsQ2b`LP1EFM;-;N;55=k zd+5%vOHjf%ZB~~H%?WpN*dW9ivYZ{jlogg4_%=6Ft|}*4Vs(^JT!82h9O`HPqGchr1 +TTCCCCCAGTATCTCGTCCTCGAATGTAGATCGATCTAGCCCTCCAAACTTATACGATGC +TCACTGTCGTACCTAAACGCCTCCGTCGAGCAGAAGCTTGTTTGACAGTTCGCGAGCCTG +TATTGAGGTCGTGTCGTTCTGCGAGGCCGTCTGTAACAGCTCGTTGCAAGAAATGGGCAT +TCGTGCCTTTCGGCGTTCTTCACTAAGTAGAGAGTGCCATAAGCGCTATCCATCGGGACC +GAACCCAATCGCCGAGTCGCACTACTACTACCTGGGATAGCCAGGGATTAAGCGGGATGC +>chr2 +CAGTTCCTACCGCATAACCCGAGCTCAATTACCGTCTAGCCTCATCCCTCCGCACCCTTA +CGGGTCAAATAGCAGCGGAGCGACGTCGGCTAAGCCATAGGTCCAAACCCGGAATATAGT +TTGGGAGCAACGGAGGCCTTAGATAAGTTTGCAGTCAGTTGAAAACTCTAAAGCCAATGC +CGCTAGCTTAAGGGCGAAAGTTCGTCGTGGCTGCCGCCGACTCAGTAATATTAAACACGC +ACGGCCAAGA diff --git a/tests/data/agc/dupB.fa b/tests/data/agc/dupB.fa new file mode 100644 index 0000000..6437bea --- /dev/null +++ b/tests/data/agc/dupB.fa @@ -0,0 +1,13 @@ +>chr1 +CGTGATGTGCGCAGCCTTCGCGAAGCCTAAGGGCGAGAGTCCACATTTTGTGACGCAGAA +GTGCAGGTCCCGAGCCTCGGTTAGATCTATGCACTAGCGAAACATCTTACACCTCCTGAA +GAACTGTCGGCAGGCTGCGTATTGGTACAGCCACCGGGTGTCGGCACTTGTGCGCATCTC +ATGTGGCAACCTCCTCGACAAATTCGTGACTTGGCCTGCTTTCTTGTTATATCAGTATTG +ACCATAACCCCGAGAACACCCCACTGTGGGGCGATACACTGGGCTGGTTATAGAGAATGA +CGGGGGGGTCCGAATAACAG +>chr2 +TGCCACGTCAGGCTGTCCAACAGGTCACGGGCGCACGTAACACGTTGCCCTTGTCTCGAC +GTCAACACGTAAACACAAAAAAGGGCTTATCTCTGAGCGGACCAGCTAGTCAACTACCCC +TAGACCATTGGTCTGTCACAACCCTGAGTCGAACGGCCTGTCGGCCGGCATCTGGTCCGA +TGAGATCGACTCGCCAATATGTTAATGCGGAAACCGCCTAGCTTGACCAATCTCCTGAAA +GTCCAATCGTCGTGCTCCCTTTGGTGAACTATTACATTTA