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..f3ab308 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,254 @@ 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, 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 { + 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 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); + + // --- 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 diff --git a/tests/agc_input.rs b/tests/agc_input.rs new file mode 100644 index 0000000..61cf968 --- /dev/null +++ b/tests/agc_input.rs @@ -0,0 +1,98 @@ +// 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; + +/// 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(); + 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}" + ); + } +} + +/// 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 0000000..adaed9f Binary files /dev/null and b/tests/data/agc/dup.agc differ diff --git a/tests/data/agc/dupA.fa b/tests/data/agc/dupA.fa new file mode 100644 index 0000000..879afbe --- /dev/null +++ b/tests/data/agc/dupA.fa @@ -0,0 +1,12 @@ +>chr1 +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 diff --git a/tests/data/agc/small.agc b/tests/data/agc/small.agc new file mode 100644 index 0000000..ce48468 Binary files /dev/null and b/tests/data/agc/small.agc differ diff --git a/tests/data/agc/small.fa b/tests/data/agc/small.fa new file mode 100644 index 0000000..812007b --- /dev/null +++ b/tests/data/agc/small.fa @@ -0,0 +1,37 @@ +>sampleA#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