diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index af86c8e..4dc5e84 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()); @@ -419,9 +439,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)) }, } @@ -851,10 +879,25 @@ 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); + } + #[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); } #[test] fn cram_unmapped_sweep_with_no_unplaced_reads_is_empty() { diff --git a/tests/data/unmapped_sweep/README.md b/tests/data/unmapped_sweep/README.md index f71d46e..e63fdf6 100644 --- a/tests/data/unmapped_sweep/README.md +++ b/tests/data/unmapped_sweep/README.md @@ -11,6 +11,17 @@ 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. + `placed_only.sam` drops the two unplaced reads from `sweep.sam`. It pins the empty-unplaced-sweep case; regenerate with samtools 1.22.1: ```bash 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('