From 91b44793739be2afc74d1266560854481d013411 Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Fri, 23 Oct 2026 11:55:18 +0200 Subject: [PATCH 1/3] test(bam): per-reference CRAM fetch must return only that reference's records In a multi-reference CRAM slice, `fetch(tid)` returns every record in the slice whose position overlaps the (whole-reference) interval, whatever its reference: noodles' query filters on the interval only. On the parallel path each such record is counted once per reference fetch (flagstat "in total" 10 instead of 6 for the samtools-written fixture at --threads 4). Co-Authored-By: Claude Opus 5.5 --- src/rna/bam/io.rs | 3 +- tests/integration_test.rs | 60 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 2 deletions(-) diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index 1f454a6..ffcb77c 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -845,9 +845,8 @@ mod tests { assert_sweep(&got); assert_each_record_once(&got); - // Only the sweep is checked here: a per-reference fetch on a - // multi-reference slice also returns the other references' records. let got = fetch_everything(&dir.join("multi_ref_slice.cram"), Some(&fasta)); assert_sweep(&got); + assert_each_record_once(&got); } } diff --git a/tests/integration_test.rs b/tests/integration_test.rs index f9d0415..afdc3e3 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -1095,3 +1095,63 @@ fn test_parallel_unmapped_sweep_counts_placed_unmapped_mate_once() { let _ = fs::remove_dir_all(root); } } + +/// The same records as a samtools-written CRAM holding all of them in one +/// multi-reference slice: each per-reference fetch on the parallel path must +/// return only its own reference's records. Expected text is `samtools +/// flagstat` 1.22.1 output for this file. +#[test] +fn test_parallel_cram_multi_reference_slice_counts_each_record_once() { + let expected = "\ +6 + 0 in total (QC-passed reads + QC-failed reads) +6 + 0 primary +0 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +3 + 0 mapped (50.00% : N/A) +3 + 0 primary mapped (50.00% : N/A) +2 + 0 paired in sequencing +1 + 0 read1 +1 + 0 read2 +0 + 0 properly paired (0.00% : N/A) +0 + 0 with itself and mate mapped +1 + 0 singletons (50.00% : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + let data = Path::new("tests/data/unmapped_sweep"); + for threads in [1, 4] { + let root = unique_test_dir(&format!("cram-multi-ref-t{threads}")); + let gtf_path = root.join("a.gtf"); + let outdir = root.join("out"); + fs::create_dir_all(&outdir).unwrap(); + fs::write( + >f_path, + "chr1\tt\tgene\t1\t1000\t.\t+\t.\tgene_id \"g1\"; gene_name \"g1\";\n\ + chr1\tt\texon\t1\t1000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\"; exon_number \"1\";\n", + ) + .unwrap(); + + let output = Command::new(rustqc_binary()) + .args(["rna", "--skip-dup-check", "--threads", &threads.to_string()]) + .arg("--reference") + .arg(data.join("sweep.fa")) + .arg("--gtf") + .arg(>f_path) + .arg("--outdir") + .arg(&outdir) + .arg(data.join("multi_ref_slice.cram")) + .output() + .expect("Failed to execute rustqc"); + assert!( + output.status.success(), + "rustqc failed with --threads {threads}:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + let flagstat = + fs::read_to_string(outdir.join("samtools/multi_ref_slice.flagstat")).unwrap(); + assert_eq!(flagstat, expected, "flagstat with --threads {threads}"); + let _ = fs::remove_dir_all(root); + } +} From 82fdd56f0505230c9e72f1be6237e09ece6eaa8e Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Fri, 23 Oct 2026 12:01:04 +0200 Subject: [PATCH 2/3] fix(bam): per-reference CRAM fetch keeps only that reference's records noodles' CRAM `query` filters records on the interval only, so in a multi-reference slice `fetch(tid)` also returned other references' records whose positions fall inside the whole-reference interval. Keep only records whose reference sequence ID is `tid`, as htslib does. Co-Authored-By: Claude Opus 5.5 --- src/rna/bam/io.rs | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index ffcb77c..fd087a5 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -419,9 +419,17 @@ impl IndexedReader { reader: heads.reader, header: heads.header, iter_builder: |reader, header| -> Result { + // noodles filters on the interval only, so a + // multi-reference slice also yields other references' + // records; keep this reference's. let query = reader .query(header, ®ion) - .with_context(|| format!("failed to query CRAM tid {tid}"))?; + .with_context(|| format!("failed to query CRAM tid {tid}"))? + .filter(move |result| { + result.as_ref().map_or(true, |record| { + record.reference_sequence_id() == Some(tid as usize) + }) + }); Ok(Box::new(query)) }, } From 41c84c7c2c2e91a1a9ae8ffb7538c5f29482d80f Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Sun, 1 Nov 2026 07:17:35 +0100 Subject: [PATCH 3/3] fix(bam): read each CRAM reference container once --- src/rna/bam/io.rs | 36 ++++++ tests/data/unmapped_sweep/README.md | 11 ++ tests/data/unmapped_sweep/make_multi_slice.py | 117 ++++++++++++++++++ .../same_reference_multi_slice.cram | Bin 0 -> 1151 bytes .../same_reference_multi_slice.cram.crai | Bin 0 -> 65 bytes tests/integration_test.rs | 54 ++++++++ 6 files changed, 218 insertions(+) create mode 100644 tests/data/unmapped_sweep/make_multi_slice.py create mode 100644 tests/data/unmapped_sweep/same_reference_multi_slice.cram create mode 100644 tests/data/unmapped_sweep/same_reference_multi_slice.cram.crai diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index fd087a5..2e4ac7d 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -1,5 +1,6 @@ //! BAM/SAM/CRAM readers with a rust-htslib-compatible API. +use std::collections::HashSet; use std::fs::File; use std::io::BufRead; use std::path::{Path, PathBuf}; @@ -75,6 +76,23 @@ fn open_bam_reader(path: &Path) -> Result Result { + let mut index_path = path.as_os_str().to_os_string(); + index_path.push(".crai"); + let mut index = noodles_cram::crai::fs::read(&index_path).with_context(|| { + format!( + "failed to read CRAM index {}", + Path::new(&index_path).display() + ) + })?; + // noodles 0.94 decodes the whole container for each CRAI slice entry. + // Keep one entry per reference/container, preserving other references' + // queries and all physical records, including repeated read names. + let mut seen = HashSet::new(); + index.retain(|entry| seen.insert((entry.reference_sequence_id(), entry.offset()))); + Ok(index) +} + /// Owned iterator over decoded CRAM records. Boxed so the `query` and /// `query_unmapped` iterators (which have distinct concrete types) can share one /// field on the indexed cursor. @@ -313,6 +331,7 @@ impl IndexedReader { let path = path.as_ref().to_path_buf(); let backend = if path.to_string_lossy().ends_with(".cram") { let mut reader = noodles_cram::io::indexed_reader::Builder::default() + .set_index(cram_container_query_index(&path)?) .build_from_path(&path) .with_context(|| format!("failed to open indexed CRAM {}", path.display()))?; let sam_header = reader @@ -361,6 +380,7 @@ impl IndexedReader { let repository = load_fasta_repository(&path)?; let mut reader = noodles_cram::io::indexed_reader::Builder::default() .set_reference_sequence_repository(repository) + .set_index(cram_container_query_index(&cram_path)?) .build_from_path(&cram_path)?; let sam_header = reader.read_header()?; *header = Header::from_noodles(sam_header.clone()); @@ -857,4 +877,20 @@ mod tests { assert_sweep(&got); assert_each_record_once(&got); } + #[test] + fn cram_same_reference_multi_slice_container_is_read_once() { + let dir = std::path::Path::new("tests/data/unmapped_sweep"); + let got = fetch_everything( + &dir.join("same_reference_multi_slice.cram"), + Some(&dir.join("sweep.fa")), + ); + let mut names: Vec<_> = got.iter().map(|(_, name)| name.as_str()).collect(); + names.sort_unstable(); + assert_eq!( + names, + ["r0", "r0", "r1", "r1", "u1", "u2"], + "each physical record, including repeated names, must be fetched once: {got:?}" + ); + assert_sweep(&got); + } } diff --git a/tests/data/unmapped_sweep/README.md b/tests/data/unmapped_sweep/README.md index 5f230af..ecac64a 100644 --- a/tests/data/unmapped_sweep/README.md +++ b/tests/data/unmapped_sweep/README.md @@ -10,3 +10,14 @@ samtools index multi_ref_slice.cram ``` `single_ref_slices.cram` has one slice per reference plus one for the unplaced reads; `multi_ref_slice.cram` holds all six records in one multi-reference slice, so its index points the unplaced-read entry at a slice that also holds placed records. + +`same_reference_multi_slice.cram` contains one chr1 container with two slices plus the unplaced-read container. Both chr1 slices hold r0 and r1, deliberately retaining repeated names as distinct physical records. Its CRAI has two chr1 entries at the same container offset. Regenerate with: + +```bash +python3 make_multi_slice.py +samtools index same_reference_multi_slice.cram +samtools quickcheck -v same_reference_multi_slice.cram +samtools view -T sweep.fa same_reference_multi_slice.cram +``` + +The generator reuses the original fixture's encoded slices, updates the second slice's record counter and the container's counts, landmarks, size and CRC. It omits the chr2 container, keeping the existing unplaced container's record counter correct. Samtools 1.22.1 decodes six records: r0, r1, r0, r1, u1, u2. Record-level or read-name deduplication would incorrectly discard two legitimate records. diff --git a/tests/data/unmapped_sweep/make_multi_slice.py b/tests/data/unmapped_sweep/make_multi_slice.py new file mode 100644 index 0000000..579d284 --- /dev/null +++ b/tests/data/unmapped_sweep/make_multi_slice.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +"""Make a two-slice CRAM container from the checked-in single-slice fixture.""" +import gzip +from pathlib import Path +import struct +import zlib + + +def read_itf8(data, position): + first = data[position] + size = next((n for n in range(1, 5) if first < (0x100 - (0x80 >> (n - 1)))), 5) + if size == 5: + value = ((first & 15) << 28) | (data[position + 1] << 20) | (data[position + 2] << 12) | (data[position + 3] << 4) | (data[position + 4] & 15) + else: + value = first & (0xff >> size) + for byte in data[position + 1:position + size]: + value = (value << 8) | byte + return value, position + size + + +def itf8(value): + value &= 0xffffffff + for size in range(1, 5): + if value < (1 << (7 * size)): + encoded = bytearray(value.to_bytes(size, 'big')) + encoded[0] |= (0xff << (9 - size)) & 0xff + return bytes(encoded) + return bytes([0xf0 | (value >> 28), (value >> 20) & 255, (value >> 12) & 255, (value >> 4) & 255, value & 15]) + + +def read_ltf8(data, position): + first = data[position] + leading = 0 + while leading < 8 and first & (0x80 >> leading): + leading += 1 + value = first & (0x7f >> leading) + for byte in data[position + 1:position + leading + 1]: + value = (value << 8) | byte + return value, position + leading + 1 + + +def ltf8(value): + if value < 128: + return bytes([value]) + raise ValueError('This fixture only needs small record counters') + + +def container(data, offset): + length = struct.unpack_from('uwE=tU0fPoo|FcuRNLnG(%)3z)?k%oVbP4e9w5UKxL`+fZm z8;G==zGUmFH^Jul^4>S0OWl~32yrcZ6Z?DVR*~BW{VuJuu-0AF#Jxef*Uk8T&Aqc} ztJ}NwT(|L^wxsgeCCl|&M48WppXJWX((AWVo$cKaW?q;jJ6rs}5`s%>Mf5)5ol`1No-X7zQfO1x#% zjQELd|E}KOzFEifWYytB@eAx5qU2)LofZ{18O~dHh*PN0!TN~XjJ*t7f*I`F9y%YB zQ}0#3k`UH2?>DqG zw7>#rAn7+V*4ih7qAE3mg+YLusbOOySN!ec=?shv4a*vq39$tEF**h?hPW^?1p7+2 zdry{b=ofNwV+KYhBZD)DW@cmv@&VF75uidw!60{L4#xjLz|Y9w3svYBj1cn+fM^a0 zhSDI-K*QYpAhf3+D-#2+lM@RYBcm20qpEKpP}~M+2}qnfFqn~%C&&|`G}sv^#K*`m zOZ4DP1||l1d7v+t7+9FOc?7`@;N)LDcl$->7Zc95%Wb%9*pzhX=zRuOpl9{j|9oU% z;bvniGGHh&WH6ISHD+buVdYrB!KlFyZ`mHfz#<4twkiK+>Fl>z^`st4Gnk(1sN1Lm=9RJdnD`Fvg3bJQs&3A(iu68H@Md{o!g*q`}wKLFGX>U Sc8hfk%nT>yYr5M2bpin3Mi*fK literal 0 HcmV?d00001 diff --git a/tests/integration_test.rs b/tests/integration_test.rs index afdc3e3..51352f7 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -1155,3 +1155,57 @@ fn test_parallel_cram_multi_reference_slice_counts_each_record_once() { let _ = fs::remove_dir_all(root); } } + +/// Multiple CRAI slice entries for one reference/container must not duplicate reads. +#[test] +fn test_parallel_cram_same_reference_multi_slice_counts_each_record_once() { + let expected = "\ +6 + 0 in total (QC-passed reads + QC-failed reads) +6 + 0 primary +0 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +4 + 0 mapped (66.67% : N/A) +4 + 0 primary mapped (66.67% : N/A) +0 + 0 paired in sequencing +0 + 0 read1 +0 + 0 read2 +0 + 0 properly paired (N/A : N/A) +0 + 0 with itself and mate mapped +0 + 0 singletons (N/A : N/A) +0 + 0 with mate mapped to a different chr +0 + 0 with mate mapped to a different chr (mapQ>=5) +"; + let data = Path::new("tests/data/unmapped_sweep"); + for threads in [1, 4] { + let root = unique_test_dir(&format!("cram-same-reference-multi-slice-t{threads}")); + let gtf_path = root.join("a.gtf"); + let outdir = root.join("out"); + fs::write(>f_path, + "chr1\tt\tgene\t1\t1000\t.\t+\t.\tgene_id \"g1\"; gene_name \"g1\";\n\ + chr1\tt\texon\t1\t1000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\"; exon_number \"1\";\n" + ).unwrap(); + let output = Command::new(env!("CARGO_BIN_EXE_rustqc")) + .args(["rna", "--skip-dup-check", "--threads", &threads.to_string()]) + .arg("--reference") + .arg(data.join("sweep.fa")) + .arg("--gtf") + .arg(>f_path) + .arg("--outdir") + .arg(&outdir) + .arg(data.join("same_reference_multi_slice.cram")) + .output() + .expect("Failed to execute rustqc"); + assert!( + output.status.success(), + "rustqc failed with --threads {threads}:\n{}", + String::from_utf8_lossy(&output.stderr) + ); + let flagstat = + fs::read_to_string(outdir.join("samtools/same_reference_multi_slice.flagstat")) + .unwrap(); + assert_eq!(flagstat, expected, "flagstat with --threads {threads}"); + fs::remove_dir_all(root).unwrap(); + } +}