From 8d926ce6058a75308e0ff192c1a624fabd89d076 Mon Sep 17 00:00:00 2001 From: Logan Nye <87274608+logannye@users.noreply.github.com> Date: Tue, 2 Jun 2026 13:13:52 -0700 Subject: [PATCH] =?UTF-8?q?fix:=20real-genome=20correctness=20=E2=80=94=20?= =?UTF-8?q?multi-contig=20eval,=20contig-naming=20guard,=20IUPAC?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three ways the engine produced wrong results (or crashed) on a real multi-contig human reference while passing on the synthetic demo path: 1. Multi-contig eval (HIGH). eval-germline/eval-somatic loaded only the FIRST FASTA record and compared EVERY variant against it, so a real GIAB run either silently miscompared indels or crashed with an opaque out-of-bounds error. compare_callsets now takes a contig->sequence map and normalizes each variant against its OWN contig (read_fasta_map loads all records); an absent contig errors clearly (names the contig + the available ones), and the normalize out-of-bounds error names the contig. 2. Contig-naming guard (HIGH). A UCSC 'chr1' vs Ensembl '1' mismatch silently dropped every read and wrote an empty VCF with exit 0 — the worst failure mode. validate_contig_lengths now refuses up front when no BAM @SQ name resolves into the index, with an actionable message. 3. IUPAC ambiguity (medium). `rosalind index` aborted on R/Y/S/W/K/M/B/D/H/V, which appear in GRCh38's primary assembly and many references. sanitize_ reference maps them to N (matching bwa/bowtie); genuinely invalid bytes still error. Tests: multi-contig eval + naming-mismatch (integration), disjoint-naming guard + IUPAC-to-N (unit). compare_callsets call sites updated to the map API. Full suite green; rustc 0 warnings. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/genomics/eval/compare.rs | 31 +++++++-- src/genomics/eval/normalize.rs | 8 ++- src/genomics/fm_index.rs | 36 +++++++++-- src/io/bam.rs | 46 ++++++++++++++ src/main.rs | 28 +++++++- tests/clinical_eval_gates.rs | 7 +- tests/eval_multicontig.rs | 113 +++++++++++++++++++++++++++++++++ tests/germline_accuracy.rs | 6 +- 8 files changed, 258 insertions(+), 17 deletions(-) create mode 100644 tests/eval_multicontig.rs diff --git a/src/genomics/eval/compare.rs b/src/genomics/eval/compare.rs index 43c1c5e..0cb8b61 100644 --- a/src/genomics/eval/compare.rs +++ b/src/genomics/eval/compare.rs @@ -48,9 +48,13 @@ impl ComparisonReport { } } -/// Compare two VCF callsets against a reference sequence, optionally masked by BED regions. +/// Compare two VCF callsets against a per-contig reference map, optionally masked +/// by BED regions. Each variant is normalized against the sequence of ITS OWN +/// contig (`references[v.chrom]`) — loading only the first FASTA record and +/// applying it to every contig silently miscompares (or crashes) on any +/// multi-contig benchmark, e.g. a real GIAB run. pub fn compare_callsets( - reference: &[u8], + references: &BTreeMap>, calls: &[VcfVariant], truth: &[VcfVariant], bed: Option<&BedIndex>, @@ -64,7 +68,7 @@ pub fn compare_callsets( continue; } } - calls_set.insert(normalize_variant(reference, v)?); + calls_set.insert(normalize_variant(reference_for(references, v)?, v)?); } for v in truth { if let Some(bed) = bed { @@ -72,7 +76,7 @@ pub fn compare_callsets( continue; } } - truth_set.insert(normalize_variant(reference, v)?); + truth_set.insert(normalize_variant(reference_for(references, v)?, v)?); } let mut tp = 0usize; @@ -113,6 +117,25 @@ pub fn compare_callsets( }) } +/// The reference sequence for a variant's contig, or a clear error naming the +/// missing contig and the available ones (the common contig-naming mismatch). +fn reference_for<'a>( + references: &'a BTreeMap>, + v: &VcfVariant, +) -> Result<&'a [u8], anyhow::Error> { + references.get(&v.chrom).map(Vec::as_slice).ok_or_else(|| { + let available: Vec<&str> = references.keys().map(String::as_str).collect(); + anyhow::anyhow!( + "variant on contig '{}' (pos {}) has no matching sequence in the reference FASTA — \ + check the contig naming scheme (e.g. UCSC 'chr1' vs Ensembl '1'). \ + Reference contigs: [{}]", + v.chrom, + v.pos0 + 1, + available.join(", ") + ) + }) +} + fn variant_type(v: &NormalizedVariant) -> VariantType { if v.reference.len() == 1 && v.alternate.len() == 1 { VariantType::Snv diff --git a/src/genomics/eval/normalize.rs b/src/genomics/eval/normalize.rs index 2e1027b..52c405a 100644 --- a/src/genomics/eval/normalize.rs +++ b/src/genomics/eval/normalize.rs @@ -6,8 +6,13 @@ use super::VcfVariant; /// Errors that can occur during variant normalization. pub enum NormalizeError { /// The variant POS is outside the provided reference slice. - #[error("variant position {pos0} out of bounds for reference length {reference_len}")] + #[error( + "variant on contig '{chrom}' position {pos0} is out of bounds for that contig's \ + reference length {reference_len}" + )] OutOfBounds { + /// Contig name (so a multi-contig miscompare names the right contig). + chrom: String, /// 0-based position. pos0: u32, /// Reference length. @@ -55,6 +60,7 @@ pub fn normalize_variant( if pos0 as usize >= reference.len() { return Err(NormalizeError::OutOfBounds { + chrom: v.chrom.clone(), pos0, reference_len: reference.len(), }); diff --git a/src/genomics/fm_index.rs b/src/genomics/fm_index.rs index f7fe055..23881d7 100644 --- a/src/genomics/fm_index.rs +++ b/src/genomics/fm_index.rs @@ -426,10 +426,23 @@ fn sanitize_reference(reference: &[u8]) -> Result, FMIndexError> { clean.push(uppercase); } None => { - return Err(FMIndexError::UnsupportedCharacter { - ch: ch as char, - position: idx, - }); + // IUPAC ambiguity codes (R,Y,S,W,K,M,B,D,H,V) appear in real + // assemblies (GRCh38's primary assembly, many bacterial/viral + // references). Map them to N — matching bwa/bowtie — so the index + // step ingests a stock reference instead of aborting; the N-mask + // carries the ambiguity, and these positions never match an + // A/C/G/T query. Genuinely non-sequence bytes still error loudly. + if matches!( + ch.to_ascii_uppercase(), + b'R' | b'Y' | b'S' | b'W' | b'K' | b'M' | b'B' | b'D' | b'H' | b'V' + ) { + clean.push(b'N'); + } else { + return Err(FMIndexError::UnsupportedCharacter { + ch: ch as char, + position: idx, + }); + } } } } @@ -561,6 +574,21 @@ mod tests { bwt[..bounded].iter().filter(|&&ch| ch == base).count() as u32 } + #[test] + fn sanitize_maps_iupac_ambiguity_codes_to_n() { + // Real references carry IUPAC degeneracy codes; the index must ingest + // them (mapped to N) rather than aborting. Genuinely invalid bytes still + // error loudly. + let clean = sanitize_reference(b"ACGTRYSWKMryswkmBDHV").unwrap(); + assert_eq!(&clean, b"ACGTNNNNNNNNNNNNNNNN"); + // A stock-reference-shaped sequence with ambiguity codes builds an index. + let index = BlockedFMIndex::build(b"ACGTRYSWKMACGTACGTAC", 4) + .expect("index build must succeed over IUPAC codes"); + let _ = index; + // A non-sequence byte is still rejected. + assert!(sanitize_reference(b"ACGT@CGT").is_err()); + } + #[test] fn sa_at_recovers_reference_position() { let reference = b"ACGTACGT"; diff --git a/src/io/bam.rs b/src/io/bam.rs index 77f12b9..47a5a6b 100644 --- a/src/io/bam.rs +++ b/src/io/bam.rs @@ -84,10 +84,12 @@ pub(crate) fn validate_contig_lengths( header: &bam::HeaderView, contigs: &ContigSet, ) -> Result<(), CoreError> { + let mut matched = 0usize; for tid in 0..header.target_count() { let name = std::str::from_utf8(header.tid2name(tid)) .map_err(|_| CoreError::MalformedRecord("BAM reference name is not UTF-8".into()))?; if let Some(c) = contigs.by_name(name) { + matched += 1; let header_len = header.target_len(tid); if header_len != Some(c.length as u64) { return Err(CoreError::MalformedRecord(format!( @@ -101,6 +103,27 @@ pub(crate) fn validate_contig_lengths( } } } + // If the BAM names a reference but NONE of its @SQ contigs resolve into the + // index, the two use incompatible naming schemes (the classic UCSC 'chr1' vs + // Ensembl '1' / GRCh38-vs-hg38 mismatch). Left unchecked, every record is + // silently dropped and the run writes an empty VCF with exit 0 — the worst + // kind of failure. Refuse up front with an actionable message instead. + if header.target_count() > 0 && matched == 0 { + let bam_names: Vec = (0..header.target_count()) + .filter_map(|tid| { + std::str::from_utf8(header.tid2name(tid)) + .ok() + .map(str::to_string) + }) + .take(3) + .collect(); + return Err(CoreError::MalformedRecord(format!( + "no BAM @SQ contig name matches the index — the alignments and the index use \ + incompatible naming schemes (e.g. UCSC 'chr1' vs Ensembl '1'). BAM names start \ + with [{}]; check the reference/assembly used for alignment.", + bam_names.join(", ") + ))); + } Ok(()) } @@ -305,6 +328,29 @@ mod tests { std::fs::remove_dir_all(&dir).ok(); } + #[test] + fn streaming_source_rejects_disjoint_contig_naming() { + // BAM @SQ uses Ensembl '1'; the index uses UCSC 'chr1'. No name resolves, + // so opening must refuse up front rather than silently dropping every read + // into an empty VCF. + let bam_path = tmp("naming-mismatch"); + let header = test_header(&[("1", 1000)]); // Ensembl-style + { + let mut writer = bam::Writer::from_path(&bam_path, &header, bam::Format::Bam).unwrap(); + let hv = writer.header().clone(); + push_record(&mut writer, &hv, 0, 10, b"ACGT"); + } + let contigs = one_contig(); // UCSC-style 'chr1' + let err = StreamingBamSource::new(&bam_path, &contigs) + .expect_err("disjoint naming must be rejected"); + let msg = err.to_string(); + assert!( + msg.contains("no BAM @SQ contig name matches the index"), + "unexpected error: {msg}" + ); + std::fs::remove_file(&bam_path).ok(); + } + #[test] fn streaming_source_rejects_out_of_order() { let dir = std::env::temp_dir().join(format!( diff --git a/src/main.rs b/src/main.rs index 63a6eb9..59e5a41 100644 --- a/src/main.rs +++ b/src/main.rs @@ -484,9 +484,8 @@ fn run_eval( truth_path: PathBuf, regions_path: Option, ) -> Result<()> { - let fasta = read_fasta(&reference_path) + let references = read_fasta_map(&reference_path) .with_context(|| format!("failed to read reference from {}", reference_path.display()))?; - let reference = fasta.sequence; let calls_txt = std::fs::read_to_string(&calls_path) .with_context(|| format!("failed to read calls VCF {}", calls_path.display()))?; @@ -506,7 +505,7 @@ fn run_eval( None }; - let report = compare_callsets(&reference, &calls, &truth, bed.as_ref())?; + let report = compare_callsets(&references, &calls, &truth, bed.as_ref())?; println!("truth_total={}", report.total_truth); println!("calls_total={}", report.total_calls); println!("tp={}", report.true_positive); @@ -1877,6 +1876,29 @@ fn read_fasta(path: &PathBuf) -> Result { Ok(first) } +/// Read ALL FASTA records into a contig-name → sequence map. Used by `eval-*`, +/// which must normalize each variant against its OWN contig — loading only the +/// first record (the single-contig `read_fasta` policy) silently miscompares or +/// crashes on any multi-contig benchmark. +fn read_fasta_map(path: &PathBuf) -> Result>> { + let reader = open_input(path).with_context(|| format!("failed to open {}", path.display()))?; + let mut map = std::collections::BTreeMap::new(); + for rec in FastaReader::new(reader) { + let rec = rec.with_context(|| format!("failed to parse FASTA {}", path.display()))?; + if map.insert(rec.name.clone(), rec.sequence).is_some() { + bail!( + "FASTA {} has a duplicate contig name '{}'", + path.display(), + rec.name + ); + } + } + if map.is_empty() { + bail!("FASTA file {} is missing a record", path.display()); + } + Ok(map) +} + /// Read a FASTQ file (plain or gzip; `-` = stdin) into a vector of records. /// The streaming parser lives in `io::fastq`. fn read_fastq(path: &PathBuf) -> Result> { diff --git a/tests/clinical_eval_gates.rs b/tests/clinical_eval_gates.rs index 4fde515..ec1d835 100644 --- a/tests/clinical_eval_gates.rs +++ b/tests/clinical_eval_gates.rs @@ -1,9 +1,10 @@ use rosalind::genomics::{compare_callsets, read_vcf_variants, BedIndex}; +use std::collections::BTreeMap; #[test] fn eval_compare_smoke_and_thresholds() { - // Small deterministic reference. - let reference = vec![b'A'; 100]; + // Small deterministic reference (single contig 'chr1'). + let references = BTreeMap::from([("chr1".to_string(), vec![b'A'; 100])]); let calls_vcf = "\ ##fileformat=VCFv4.3 @@ -23,7 +24,7 @@ chr1\t31\t.\tA\tT\t50\tPASS\t. // Mask to only evaluate positions 0..40. let bed = BedIndex::from_str("chr1\t0\t40\n").unwrap(); - let rep = compare_callsets(&reference, &calls, &truth, Some(&bed)).unwrap(); + let rep = compare_callsets(&references, &calls, &truth, Some(&bed)).unwrap(); // TP: pos11 A>C // FP: pos21 A>G diff --git a/tests/eval_multicontig.rs b/tests/eval_multicontig.rs new file mode 100644 index 0000000..9b48932 --- /dev/null +++ b/tests/eval_multicontig.rs @@ -0,0 +1,113 @@ +//! `eval-germline` must normalize each variant against its OWN contig. Loading +//! only the first FASTA record silently miscompared (or crashed with an opaque +//! out-of-bounds error) on any multi-contig benchmark — the headline truth- +//! comparison surface, and exactly what a real GIAB run looks like. + +use std::path::PathBuf; +use std::process::Command; +use std::sync::atomic::{AtomicU64, Ordering}; + +fn bin() -> &'static str { + env!("CARGO_BIN_EXE_rosalind") +} + +fn unique_dir(prefix: &str) -> PathBuf { + static COUNTER: AtomicU64 = AtomicU64::new(0); + let n = COUNTER.fetch_add(1, Ordering::Relaxed); + let nanos = std::time::SystemTime::now() + .duration_since(std::time::UNIX_EPOCH) + .unwrap() + .as_nanos(); + let d = std::env::temp_dir().join(format!("{prefix}-{nanos}-{n}")); + std::fs::create_dir_all(&d).unwrap(); + d +} + +// chr1 = 20 bp, chr2 = 40 bp. +const REFERENCE: &str = + ">chr1\nACGTACGTACGTACGTACGT\n>chr2\nACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT\n"; + +#[test] +fn eval_germline_compares_each_variant_against_its_own_contig() { + let dir = unique_dir("rosalind-eval-mc"); + let fa = dir.join("ref.fa"); + std::fs::write(&fa, REFERENCE).unwrap(); + + // A chr1 SNV (pos 5) and a chr2 SNV at pos 30 — past chr1's length (20), so + // the old first-contig-only path crashed out-of-bounds on the chr2 variant. + let header = "##fileformat=VCFv4.3\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"; + let calls = dir.join("calls.vcf"); + std::fs::write( + &calls, + format!("{header}chr1\t5\t.\tA\tC\t50\tPASS\t.\nchr2\t30\t.\tG\tT\t50\tPASS\t.\n"), + ) + .unwrap(); + let truth = dir.join("truth.vcf"); + std::fs::write( + &truth, + format!("{header}chr1\t5\t.\tA\tC\t50\tPASS\t.\nchr2\t30\t.\tG\tT\t50\tPASS\t.\n"), + ) + .unwrap(); + + let out = Command::new(bin()) + .args(["eval-germline", "--reference"]) + .arg(&fa) + .arg("--calls") + .arg(&calls) + .arg("--truth") + .arg(&truth) + .output() + .unwrap(); + assert!( + out.status.success(), + "multi-contig eval must succeed: {}", + String::from_utf8_lossy(&out.stderr) + ); + let stdout = String::from_utf8_lossy(&out.stdout); + // Both contigs' variants match truth → 2 TP, 0 FP, 0 FN. + assert!( + stdout.contains("tp=2"), + "expected 2 true positives: {stdout}" + ); + assert!( + stdout.contains("fp=0"), + "expected 0 false positives: {stdout}" + ); + assert!( + stdout.contains("fn=0"), + "expected 0 false negatives: {stdout}" + ); + std::fs::remove_dir_all(&dir).ok(); +} + +#[test] +fn eval_germline_errors_clearly_on_a_contig_naming_mismatch() { + let dir = unique_dir("rosalind-eval-naming"); + let fa = dir.join("ref.fa"); + std::fs::write(&fa, REFERENCE).unwrap(); + + // A variant on a contig the FASTA does not contain (Ensembl '1' vs UCSC + // 'chr1') must error clearly, not silently miscompare. + let header = "##fileformat=VCFv4.3\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"; + let calls = dir.join("calls.vcf"); + std::fs::write(&calls, format!("{header}1\t5\t.\tA\tC\t50\tPASS\t.\n")).unwrap(); + let truth = dir.join("truth.vcf"); + std::fs::write(&truth, header).unwrap(); + + let out = Command::new(bin()) + .args(["eval-germline", "--reference"]) + .arg(&fa) + .arg("--calls") + .arg(&calls) + .arg("--truth") + .arg(&truth) + .output() + .unwrap(); + assert!(!out.status.success(), "naming mismatch must fail"); + let stderr = String::from_utf8_lossy(&out.stderr); + assert!( + stderr.contains("naming scheme") || stderr.contains("no matching sequence"), + "expected a clear naming-mismatch error: {stderr}" + ); + std::fs::remove_dir_all(&dir).ok(); +} diff --git a/tests/germline_accuracy.rs b/tests/germline_accuracy.rs index 613c748..a1aa846 100644 --- a/tests/germline_accuracy.rs +++ b/tests/germline_accuracy.rs @@ -251,8 +251,10 @@ fn run_accuracy(coverage: usize, error_rate: f64, cap: u32, seed: u64) -> Accura }) .collect(); - let report_all = compare_callsets(&reference, &calls_all, &truth_vcf, None).unwrap(); - let report_pass = compare_callsets(&reference, &calls_pass, &truth_vcf, None).unwrap(); + // Single-contig harness: all variants are on 'chr1'. + let references = BTreeMap::from([("chr1".to_string(), reference.clone())]); + let report_all = compare_callsets(&references, &calls_all, &truth_vcf, None).unwrap(); + let report_pass = compare_callsets(&references, &calls_pass, &truth_vcf, None).unwrap(); let deep_called = calls_pass.iter().any(|c| c.pos0 as usize == deep_site); std::fs::remove_dir_all(&dir).ok();