fix(bam): prevent duplicate reads from indexed CRAM queries

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<noodles_bam::io::Reader<bgzf::io::Read
     Ok(noodles_bam::io::Reader::new(file))
 }
 
+fn cram_container_query_index(path: &Path) -> Result<noodles_cram::crai::Index> {
+    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<CramRecordIter> {
+                            // 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, &region)
-                                .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() {

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

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('<i', data, offset)[0]
+    position = offset + 4
+    fields = []
+    for _ in range(4):
+        value, position = read_itf8(data, position)
+        fields.append(value)
+    for _ in range(2):
+        value, position = read_ltf8(data, position)
+        fields.append(value)
+    blocks, position = read_itf8(data, position)
+    count, position = read_itf8(data, position)
+    landmarks = []
+    for _ in range(count):
+        value, position = read_itf8(data, position)
+        landmarks.append(value)
+    assert struct.unpack_from('<I', data, position)[0] == zlib.crc32(data[offset:position])
+    position += 4
+    return fields, blocks, landmarks, data[position:position + length]
+
+
+def recounter_slice(data, counter):
+    method, content_type = data[:2]
+    assert content_type == 2
+    content_id, position = read_itf8(data, 2)
+    compressed_size, position = read_itf8(data, position)
+    raw_size, position = read_itf8(data, position)
+    end = position + compressed_size
+    assert struct.unpack_from('<I', data, end)[0] == zlib.crc32(data[:end])
+    payload = data[position:end]
+    if method == 1:
+        payload = gzip.decompress(payload)
+    else:
+        assert method == 0
+    assert len(payload) == raw_size
+    position = 0
+    for _ in range(4):
+        _, position = read_itf8(payload, position)
+    _, counter_end = read_ltf8(payload, position)
+    payload = payload[:position] + ltf8(counter) + payload[counter_end:]
+    block = bytes([0, content_type]) + itf8(content_id) + itf8(len(payload)) * 2 + payload
+    return block + struct.pack('<I', zlib.crc32(block)) + data[end + 4:]
+
+
+def main():
+    directory = Path(__file__).resolve().parent
+    source = directory / 'single_ref_slices.cram'
+    data = source.read_bytes()
+    assert data[:5] == b'CRAM\x03' and data[5] in (0, 1)
+    entries = [list(map(int, line.split('\t'))) for line in gzip.decompress(source.with_suffix('.cram.crai').read_bytes()).decode().splitlines()]
+    offset = next(entry[3] for entry in entries if entry[0] == 0)
+    unplaced = next(entry[3] for entry in entries if entry[0] == -1)
+    fields, blocks, landmarks, payload = container(data, offset)
+    assert len(landmarks) == 1 and fields[0] == 0 and fields[3] == 2 and fields[4] == 0
+    header_block = payload[:landmarks[0]]
+    first_slice = payload[landmarks[0]:]
+    second_slice = recounter_slice(first_slice, 2)
+    payload = header_block + first_slice + second_slice
+    fields[3] *= 2
+    fields[5] *= 2
+    landmarks.append(landmarks[0] + len(first_slice))
+    header = struct.pack('<i', len(payload)) + b''.join(itf8(n) for n in fields[:4])
+    header += b''.join(ltf8(n) for n in fields[4:])
+    header += itf8(1 + 2 * (blocks - 1)) + itf8(2) + b''.join(itf8(n) for n in landmarks)
+    output = data[:offset] + header + struct.pack('<I', zlib.crc32(header)) + payload + data[unplaced:]
+    (directory / 'same_reference_multi_slice.cram').write_bytes(output)
+
+
+if __name__ == '__main__':
+    main()

tests/integration_test.rs

@@ -1096,6 +1096,120 @@ fn test_parallel_unmapped_sweep_counts_placed_unmapped_mate_once() {
     }
 }
 
+/// 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(
+            &gtf_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(&gtf_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);
+    }
+}
+
+/// 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(&gtf_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(&gtf_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();
+    }
+}
+
 /// A CRAM with no unplaced records must succeed on the parallel sweep path too.
 #[test]
 fn test_parallel_cram_with_no_unplaced_reads() {