fix: real-genome correctness — multi-contig eval, contig-naming guard, IUPAC (Act-1 PR-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<String, Vec<u8>>,
     calls: &[VcfVariant],
     truth: &[VcfVariant],
     bed: Option<&BedIndex>,
@@ -64,15 +68,15 @@ 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 {
             if !bed.contains(&v.chrom, v.pos0) {
                 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<String, Vec<u8>>,
+    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

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(),
         });

src/genomics/fm_index.rs

@@ -426,10 +426,23 @@ fn sanitize_reference(reference: &[u8]) -> Result<Vec<u8>, 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";

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<String> = (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!(

src/main.rs

@@ -484,9 +484,8 @@ fn run_eval(
     truth_path: PathBuf,
     regions_path: Option<PathBuf>,
 ) -> 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<FastaRecord> {
     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<std::collections::BTreeMap<String, Vec<u8>>> {
+    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<Vec<FastqRecord>> {

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

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();
+}