From 14a3f055624bd36375612185a12daed58510eb1f Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Fri, 23 Oct 2026 11:34:07 +0200 Subject: [PATCH 1/6] test(bam): pin the unmapped sweep to records with no reference On the indexed (parallel) path, `fetch(Unmapped)` seeks to the last linear-index window of the last reference and keeps every UNMAP-flagged record, so a placed unmapped mate in that window is read by both its reference's fetch and the sweep. htslib's HTS_IDX_NOCOOR returns only records with no reference (tid -1). - Unit tests: BAM (noodles-written) and two samtools-written CRAMs (per-reference slices, one multi-reference slice); the sweep must return exactly the two unplaced records. - Integration test: flagstat at --threads 1 and 4 must match samtools 1.22.1; at 4 the placed mate is counted twice (in total 7, not 6). Co-Authored-By: Claude Opus 5.5 --- src/rna/bam/io.rs | 150 +++++++++++++++++- tests/data/unmapped_sweep/README.md | 12 ++ .../data/unmapped_sweep/multi_ref_slice.cram | Bin 0 -> 931 bytes .../unmapped_sweep/multi_ref_slice.cram.crai | Bin 0 -> 58 bytes .../unmapped_sweep/single_ref_slices.cram | Bin 0 -> 1369 bytes .../single_ref_slices.cram.crai | Bin 0 -> 68 bytes tests/data/unmapped_sweep/sweep.fa | 4 + tests/data/unmapped_sweep/sweep.sam | 9 ++ tests/integration_test.rs | 69 ++++++++ 9 files changed, 243 insertions(+), 1 deletion(-) create mode 100644 tests/data/unmapped_sweep/README.md create mode 100644 tests/data/unmapped_sweep/multi_ref_slice.cram create mode 100644 tests/data/unmapped_sweep/multi_ref_slice.cram.crai create mode 100644 tests/data/unmapped_sweep/single_ref_slices.cram create mode 100644 tests/data/unmapped_sweep/single_ref_slices.cram.crai create mode 100644 tests/data/unmapped_sweep/sweep.fa create mode 100644 tests/data/unmapped_sweep/sweep.sam diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index 6ebc23a3..cea02fc4 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -613,7 +613,7 @@ mod tests { use noodles::sam::header::record::value::Map; use noodles::{cram, sam}; - use super::{Read, Reader}; + use super::{FetchTarget, IndexedReader, Read, Reader}; use crate::rna::bam::Record; static COUNTER: AtomicU64 = AtomicU64::new(0); @@ -693,4 +693,152 @@ mod tests { assert_eq!(got, expected.to_vec()); } + + // ------------------------------------------------------------------------ + // Unmapped sweep: htslib's HTS_IDX_NOCOOR returns exactly the records with + // no reference (tid -1). A placed unmapped mate (0x4 with a tid) belongs + // to its reference's fetch, so the sweep must not return it again. + // ------------------------------------------------------------------------ + + const CONTIG_LEN: usize = 100_000; + + /// Two mapped reads on chr1; on chr2, a mapped read whose unmapped mate is + /// placed at its position in chr2's last 16 kb linear-index window; two + /// unplaced reads (no reference) at the end. + fn sweep_world() -> (sam::Header, Vec) { + let header: sam::Header = format!( + "@HD\tVN:1.6\tSO:coordinate\n@SQ\tSN:chr1\tLN:{CONTIG_LEN}\n@SQ\tSN:chr2\tLN:{CONTIG_LEN}\n" + ) + .parse() + .unwrap(); + let rec = |name: &str, flags: u16, tid: Option, start: Option| { + let mut builder = RecordBuf::builder() + .set_name(name) + .set_flags(Flags::from(flags)) + .set_sequence(b"AAAA".to_vec().into()) + .set_quality_scores(vec![30u8; 4].into()); + if let (Some(tid), Some(start)) = (tid, start) { + let start = Position::new(start).unwrap(); + builder = builder + .set_reference_sequence_id(tid) + .set_alignment_start(start); + if flags & 0x1 != 0 { + builder = builder + .set_mate_reference_sequence_id(tid) + .set_mate_alignment_start(start); + } + if flags & 0x4 == 0 { + builder = builder + .set_mapping_quality( + sam::alignment::record::MappingQuality::new(60).unwrap(), + ) + .set_cigar([Op::new(Kind::Match, 4)].into_iter().collect()); + } + } + builder.build() + }; + let records = vec![ + rec("r0", 0x0, Some(0), Some(1_001)), + rec("r1", 0x0, Some(0), Some(1_201)), + rec("r2", 0x1 | 0x8 | 0x40, Some(1), Some(99_001)), + rec("r3", 0x1 | 0x4 | 0x80, Some(1), Some(99_001)), + rec("u1", 0x4, None, None), + rec("u2", 0x4, None, None), + ]; + (header, records) + } + + /// Fetch each reference, then the unmapped bucket; return (fetch, read name) + /// for every record returned. + fn fetch_everything( + path: &std::path::Path, + fasta: Option<&std::path::Path>, + ) -> Vec<(String, String)> { + let mut reader = IndexedReader::from_path(path).unwrap(); + if let Some(fasta) = fasta { + reader.set_reference(fasta).unwrap(); + } + let mut got = Vec::new(); + let mut record = Record::new(); + let targets = [ + ("tid 0", FetchTarget::Tid(0)), + ("tid 1", FetchTarget::Tid(1)), + ("unmapped", FetchTarget::Unmapped), + ]; + for (label, target) in targets { + reader.fetch(target).unwrap(); + while let Some(result) = reader.read(&mut record) { + result.unwrap(); + got.push(( + label.to_string(), + String::from_utf8(record.qname().to_vec()).unwrap(), + )); + } + } + got + } + + fn assert_sweep(got: &[(String, String)]) { + let swept: Vec<&str> = got + .iter() + .filter(|(label, _)| label == "unmapped") + .map(|(_, name)| name.as_str()) + .collect(); + assert_eq!( + swept, + ["u1", "u2"], + "unmapped sweep must return only the unplaced records; all fetches: {got:?}" + ); + } + + fn assert_each_record_once(got: &[(String, String)]) { + let mut names: Vec<&str> = got.iter().map(|(_, name)| name.as_str()).collect(); + names.sort_unstable(); + assert_eq!( + names, + ["r0", "r1", "r2", "r3", "u1", "u2"], + "each record must be fetched exactly once; all fetches: {got:?}" + ); + } + + #[test] + fn bam_unmapped_sweep_returns_only_unplaced_records() { + let (header, records) = sweep_world(); + let bam_path = temp_path("sweep.bam"); + { + let mut writer = noodles::bam::io::writer::Builder + .build_from_path(&bam_path) + .unwrap(); + writer.write_header(&header).unwrap(); + for record in &records { + writer.write_alignment_record(&header, record).unwrap(); + } + writer.try_finish().unwrap(); + } + let bai_path = std::path::PathBuf::from(format!("{}.bai", bam_path.display())); + let index = noodles::bam::fs::index(&bam_path).unwrap(); + noodles::bam::bai::fs::write(&bai_path, &index).unwrap(); + + let got = fetch_everything(&bam_path, None); + let _ = std::fs::remove_file(&bam_path); + let _ = std::fs::remove_file(&bai_path); + assert_sweep(&got); + assert_each_record_once(&got); + } + + /// samtools-written CRAMs of the same records (see the fixture README). + #[test] + fn cram_unmapped_sweep_returns_only_unplaced_records() { + let dir = std::path::Path::new("tests/data/unmapped_sweep"); + let fasta = dir.join("sweep.fa"); + + let got = fetch_everything(&dir.join("single_ref_slices.cram"), Some(&fasta)); + 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); + } } diff --git a/tests/data/unmapped_sweep/README.md b/tests/data/unmapped_sweep/README.md new file mode 100644 index 00000000..5f230af5 --- /dev/null +++ b/tests/data/unmapped_sweep/README.md @@ -0,0 +1,12 @@ +# Unmapped-sweep CRAM fixtures + +`sweep.sam` has two mapped reads on chr1, a mapped read on chr2 with its unmapped mate placed at the same position (flag 0x85, tid 1), and two unplaced reads (no reference). The CRAMs are that file written by samtools 1.22.1 against `sweep.fa`: + +```bash +samtools view -C -T sweep.fa -o single_ref_slices.cram sweep.sam +samtools view -C -T sweep.fa --output-fmt-option multi_seq_per_slice=1 -o multi_ref_slice.cram sweep.sam +samtools index single_ref_slices.cram +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. diff --git a/tests/data/unmapped_sweep/multi_ref_slice.cram b/tests/data/unmapped_sweep/multi_ref_slice.cram new file mode 100644 index 0000000000000000000000000000000000000000..9613ea845da5b9bd9f986b9dd51a814641eb9c74 GIT binary patch literal 931 zcmZ<`a`a_p%q`6+$<#|OO3Y<|gCmSECKD4wqngIvhv3^X3 zBCQ`cu304g;?6n#37JVwOr5h77Imy&vZZ{<$r7tAa~oxM`Z}uaV!WUC&hqCuOL_C= zgJ0_W?oCtjJ(luUWlC|}Y&SON+(m{UTp8eWs}VU&hRfaviU3x_rXThV|*vwSEG99l1d#?FDb{cU|-D#Nim0 z%*Dl}3ag*K|7&V}deW3LVoYX?wGGEs>|LGeVe(C5xA2*b?JoP0g=bDJQ)irQ*YMuS zM&sE=^NNpfoT#eAKC(`sPw2!F5<1%^$-KXwMCnWOGN3E}0|6fc0~><` zH)F$+y@~TB7#JBEW;V~RV%*YVr z1Eg6P85{%Hm>9T%+?hETLB{bjGWY_OiTVXIGK%;GFfxjV1OpWU%?$x*0Mc%L5Zcp^ zm5G7Z$q~qNWMpu1Vqs%s)M8{*^$lcX6b=GOat8tp;R*6&WK;}x2J-nB8S;Lg4QF6t z(9nQ8os)x`hhJDkR9q2kJ$filR}Y)Wz{pF literal 0 HcmV?d00001 diff --git a/tests/data/unmapped_sweep/multi_ref_slice.cram.crai b/tests/data/unmapped_sweep/multi_ref_slice.cram.crai new file mode 100644 index 0000000000000000000000000000000000000000..4eb5796b8217afee8ae087b6611da0fc4f4fce6e GIT binary patch literal 58 zcmb2|=3oE=VWS5oCI&1?!hXWu#eUAl!hXum%zDm#-U~EV%~+?noN>8gmw-b9Bg6Y8 JlR0gHS^#~t4zK_K literal 0 HcmV?d00001 diff --git a/tests/data/unmapped_sweep/single_ref_slices.cram b/tests/data/unmapped_sweep/single_ref_slices.cram new file mode 100644 index 0000000000000000000000000000000000000000..d9341523ff6306408e2ffd5272bd53b16aff3ba2 GIT binary patch literal 1369 zcmc)JeN2pT902g|_j$Rl?zq);#j>(C%Uds380m0#^;pgw?m}xqB}!gvHk>~miHuT~ zS<8eb{JpX*3eV^ak_w)Ncdb8F@V?pk! zWh=6yvhp)lg6FA3&l^R-X5XVXGz`h`{hJHJc1p!h-x|t{rpqybMK>Ipq^R}yy2x|(usN% zeQSQJEO&WS*_09c47_P$_$qAso9^hl%Ug>>?RV9MnT{4Gu%y&zpRc{Z$Rv{d7|2({ zgkjlnEF)*72cJo3gb!u1_3+tz4YP|InWU6b4 ZuB2I;X9^MNqztT literal 0 HcmV?d00001 diff --git a/tests/data/unmapped_sweep/single_ref_slices.cram.crai b/tests/data/unmapped_sweep/single_ref_slices.cram.crai new file mode 100644 index 0000000000000000000000000000000000000000..77f378f8def1b3e695734af47a0ecea93746e9ef GIT binary patch literal 68 zcmb2|=3oE=VVN#NJ_ZGjBhNmG3;P9r_%A=%q9h~OeVtmz-T<{y-#dMd>)w=1IMi0~ U?`GWONw$n^3|`?nB{o3K0QTb-^Z)<= literal 0 HcmV?d00001 diff --git a/tests/data/unmapped_sweep/sweep.fa b/tests/data/unmapped_sweep/sweep.fa new file mode 100644 index 00000000..a108fc96 --- /dev/null +++ b/tests/data/unmapped_sweep/sweep.fa @@ -0,0 +1,4 @@ +>chr1 +AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA +>chr2 +AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA diff --git a/tests/data/unmapped_sweep/sweep.sam b/tests/data/unmapped_sweep/sweep.sam new file mode 100644 index 00000000..2c86967e --- /dev/null +++ b/tests/data/unmapped_sweep/sweep.sam @@ -0,0 +1,9 @@ +@HD VN:1.6 SO:coordinate +@SQ SN:chr1 LN:1000 +@SQ SN:chr2 LN:1000 +r0 0 chr1 101 60 4M * 0 0 AAAA IIII +r1 0 chr1 201 60 4M * 0 0 AAAA IIII +r2 73 chr2 901 60 4M = 901 0 AAAA IIII +r3 133 chr2 901 0 * = 901 0 AAAA IIII +u1 4 * 0 0 * * 0 0 AAAA IIII +u2 4 * 0 0 * * 0 0 AAAA IIII diff --git a/tests/integration_test.rs b/tests/integration_test.rs index 45fc4509..f9d0415d 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -1026,3 +1026,72 @@ fn test_dup_check_parallel_uses_global_duplicate_state() { let _ = fs::remove_dir_all(root); } + +// =================================================================== +// Unmapped sweep on the parallel path +// =================================================================== + +/// A placed unmapped mate (0x4 with a reference) in the last reference's final +/// 16 kb, plus two unplaced reads. The parallel path's unmapped sweep must not +/// count the placed mate a second time. Expected text is `samtools flagstat` +/// 1.22.1 output for this input. +#[test] +fn test_parallel_unmapped_sweep_counts_placed_unmapped_mate_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 seq = "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC"; + let qual = "??????????????????????????????????????????????????"; + let lines: Vec = [ + "r0\t0\tchr1\t1001\t60\t50M\t*\t0\t0", + "r1\t0\tchr1\t1201\t60\t50M\t*\t0\t0", + "r2\t73\tchr2\t99001\t60\t50M\t=\t99001\t0", + "r3\t133\tchr2\t99001\t0\t*\t=\t99001\t0", + "u1\t4\t*\t0\t0\t*\t*\t0\t0", + "u2\t4\t*\t0\t0\t*\t*\t0\t0", + ] + .iter() + .map(|fields| format!("{fields}\t{seq}\t{qual}")) + .collect(); + let lines: Vec<&str> = lines.iter().map(String::as_str).collect(); + + for threads in [1, 4] { + let root = unique_test_dir(&format!("unmapped-sweep-t{threads}")); + let bam_path = root.join("s.bam"); + 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\t900\t3000\t.\t+\t.\tgene_id \"g1\"; gene_name \"g1\";\n\ + chr1\tt\texon\t900\t3000\t.\t+\t.\tgene_id \"g1\"; transcript_id \"t1\"; exon_number \"1\";\n", + ) + .unwrap(); + write_bam_fixture(&bam_path, &[("chr1", 100_000), ("chr2", 100_000)], &lines); + + let output = run_rustqc_custom(&outdir, &bam_path, >f_path, threads, true, &[]); + 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/s.flagstat")).unwrap(); + assert_eq!(flagstat, expected, "flagstat with --threads {threads}"); + let _ = fs::remove_dir_all(root); + } +} From 036b417409f8c4e0dc11e594f47154663e060a60 Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Fri, 23 Oct 2026 11:42:38 +0200 Subject: [PATCH 2/6] fix(bam): unmapped sweep returns only records with no reference `fetch(Unmapped)` now keeps a record only when it has no reference sequence ID, matching htslib's HTS_IDX_NOCOOR: - BAM: the sweep seeks to the last reference's final linear-index window, so it re-read placed records there and kept every UNMAP-flagged one; a placed unmapped mate in that window was counted twice on the parallel path (flagstat "in total" 7 instead of 6 in the new test). - CRAM: `query_unmapped` starts at the slice holding the unplaced reads, which in a multi-reference slice also holds placed records; the same filter applies. Co-Authored-By: Claude Opus 5.5 --- src/rna/bam/io.rs | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index cea02fc4..1f454a60 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -282,7 +282,7 @@ enum BamCursor { chunk_end: Option, reference_sequence_id: usize, }, - /// Streaming unmapped records from a seeked position to end of file. + /// Streaming records with no reference from a seeked position to end of file. Unmapped, } @@ -436,8 +436,10 @@ impl IndexedReader { let source_path = self.source_path.clone(); match &mut self.backend { IndexedBackend::Bam { reader, cursor, .. } => { - // Seek to the start of the unmapped read region and stream from - // there, filtering for unmapped records as we go. + // Seek to the last linear-index window of the last reference and + // stream from there, keeping only records with no reference + // (htslib's HTS_IDX_NOCOOR). Placed unmapped mates in that window + // belong to their reference's fetch. match reader.index().last_first_record_start_position() { Some(pos) => { reader @@ -448,7 +450,7 @@ impl IndexedReader { None => { // The index carries no metadata pseudo-bin, so there is // no recorded unmapped offset. Reopen and stream from the - // first record, relying on the per-record unmapped filter. + // first record, relying on the per-record no-reference filter. let mut fresh = noodles_bam::io::indexed_reader::Builder::default() .build_from_path(&source_path) .with_context(|| { @@ -470,9 +472,16 @@ impl IndexedReader { reader: heads.reader, header: heads.header, iter_builder: |reader, header| -> Result { + // The unplaced reads' slice can also hold placed records + // (multi-reference slice); keep only those with no reference. let query = reader .query_unmapped(header) - .context("failed to query unmapped CRAM reads")?; + .context("failed to query unmapped CRAM reads")? + .filter(|result| { + result.as_ref().map_or(true, |record| { + record.reference_sequence_id().is_none() + }) + }); Ok(Box::new(query)) }, } @@ -503,13 +512,13 @@ impl Read for IndexedReader { Ok(0) => return None, Ok(_) => { *cursor = BamCursor::Unmapped; - if scratch.flags().is_unmapped() { + if scratch.reference_sequence_id().is_none() { return Some( Record::from_bam(header.noodles_header(), scratch) .and_then(|loaded| store_record(record, loaded)), ); } - // Mapped record in the unmapped tail: skip and continue. + // Placed record before the unplaced tail: skip and continue. } Err(e) => return Some(Err(e.into())), }, From 2f813ef8f18d09447cf2642ebe15e2e40f0bc511 Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Sun, 1 Nov 2026 07:10:49 +0100 Subject: [PATCH 3/6] fix(bam): handle CRAM sweeps with no unplaced reads --- src/rna/bam/io.rs | 25 +++++++++ tests/data/unmapped_sweep/README.md | 7 +++ tests/data/unmapped_sweep/placed_only.cram | Bin 0 -> 962 bytes .../data/unmapped_sweep/placed_only.cram.crai | Bin 0 -> 58 bytes tests/data/unmapped_sweep/placed_only.sam | 7 +++ tests/integration_test.rs | 52 ++++++++++++++++++ 6 files changed, 91 insertions(+) create mode 100644 tests/data/unmapped_sweep/placed_only.cram create mode 100644 tests/data/unmapped_sweep/placed_only.cram.crai create mode 100644 tests/data/unmapped_sweep/placed_only.sam diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index 1f454a60..dfa7e1b9 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -472,6 +472,15 @@ impl IndexedReader { reader: heads.reader, header: heads.header, iter_builder: |reader, header| -> Result { + // With no unplaced CRAI entry, noodles 0.94 seeks EOF + // and reports a truncated container instead of an empty query. + if reader + .index() + .iter() + .all(|entry| entry.reference_sequence_id().is_some()) + { + return Ok(Box::new(std::iter::empty())); + } // The unplaced reads' slice can also hold placed records // (multi-reference slice); keep only those with no reference. let query = reader @@ -850,4 +859,20 @@ mod tests { let got = fetch_everything(&dir.join("multi_ref_slice.cram"), Some(&fasta)); assert_sweep(&got); } + #[test] + fn cram_unmapped_sweep_with_no_unplaced_reads_is_empty() { + let dir = std::path::Path::new("tests/data/unmapped_sweep"); + let got = fetch_everything(&dir.join("placed_only.cram"), Some(&dir.join("sweep.fa"))); + let swept: Vec<_> = got + .iter() + .filter(|(target, _)| target == "unmapped") + .collect(); + assert!( + swept.is_empty(), + "no unplaced records should be returned: {got:?}" + ); + let mut names: Vec<_> = got.iter().map(|(_, name)| name.as_str()).collect(); + names.sort_unstable(); + assert_eq!(names, ["r0", "r1", "r2", "r3"]); + } } diff --git a/tests/data/unmapped_sweep/README.md b/tests/data/unmapped_sweep/README.md index 5f230af5..f71d46e0 100644 --- a/tests/data/unmapped_sweep/README.md +++ b/tests/data/unmapped_sweep/README.md @@ -10,3 +10,10 @@ 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. + +`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 +samtools view --no-PG -C -T sweep.fa -o placed_only.cram placed_only.sam +samtools index placed_only.cram +``` diff --git a/tests/data/unmapped_sweep/placed_only.cram b/tests/data/unmapped_sweep/placed_only.cram new file mode 100644 index 0000000000000000000000000000000000000000..723129ebcc8abe77a6b5e2b515948941e6b1bdc8 GIT binary patch literal 962 zcmZ<`a`a_pEXYYrPECo=&&#RQOD;;xWdH&>Mi2=COiT<7JC$>-b^%2i)-{^ScXL1_ z{vYkz%XdgY!nO2BdFIa_f=t2e?CfhVx`jJ89OKIsT(NK8k(1^HI}h(wH*a4PzB+Bc z^Y*aGQ48KJd28ct?xXhG_wwAp>sh}$SY){v7+&nS8dRXCQhGo%^a6A6gDIh1kuN3E z`C|2@+Aa2+Vn3IBpwRPY=4ye}syl*~buo#yFlb6l6q>H39;4PgJ<~+7qTE0KJ@db` zyR+rRf0a+?f7G=z#}w$|(4;nQpmSqlAa14&xypleCl@HdQZtwsIM^5)Hs%#NpJ!lX zXjs;;Oo%1OkI^xJF~o(DA=p>C-FvchL%)!d8#4KsO2uwOD z|7Pj$Py&VxV{DYn4r6LxZAWaP&T8HkHf1Ba1O6c{or zjEsy7z))cY+QP^{AXprM${iUQkb{SjQ5ZdJfCj?CX7hc`ufVWTQiAw_jg^~+Usx0p zJRKA6KFEBzdTsK@S-kV)jg0nf*nt*2Mhr#93|Dj~SEC0HKQpsuL#x$MlQ#@3!iq21>V#46J)^`U5F8Hc(_UGQ3M>^acRjui&-- literal 0 HcmV?d00001 diff --git a/tests/data/unmapped_sweep/placed_only.cram.crai b/tests/data/unmapped_sweep/placed_only.cram.crai new file mode 100644 index 0000000000000000000000000000000000000000..6b2243b8db588d36da73cce1519eda54e2a46ced GIT binary patch literal 58 zcmb2|=3oE=VWS5oCI&1?!hXrg3b}dbG80dxp2*}BF?ncWY{JSlvrtNERskc61cU2S J!Bd(*Eda8n4_*KO literal 0 HcmV?d00001 diff --git a/tests/data/unmapped_sweep/placed_only.sam b/tests/data/unmapped_sweep/placed_only.sam new file mode 100644 index 00000000..2a401fca --- /dev/null +++ b/tests/data/unmapped_sweep/placed_only.sam @@ -0,0 +1,7 @@ +@HD VN:1.6 SO:coordinate +@SQ SN:chr1 LN:1000 +@SQ SN:chr2 LN:1000 +r0 0 chr1 101 60 4M * 0 0 AAAA IIII +r1 0 chr1 201 60 4M * 0 0 AAAA IIII +r2 73 chr2 901 60 4M = 901 0 AAAA IIII +r3 133 chr2 901 0 * = 901 0 AAAA IIII diff --git a/tests/integration_test.rs b/tests/integration_test.rs index f9d0415d..14d7ef40 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -1095,3 +1095,55 @@ fn test_parallel_unmapped_sweep_counts_placed_unmapped_mate_once() { let _ = fs::remove_dir_all(root); } } + +/// A CRAM with no unplaced records must succeed on the parallel sweep path too. +#[test] +fn test_parallel_cram_with_no_unplaced_reads() { + let expected = "\ +4 + 0 in total (QC-passed reads + QC-failed reads) +4 + 0 primary +0 + 0 secondary +0 + 0 supplementary +0 + 0 duplicates +0 + 0 primary duplicates +3 + 0 mapped (75.00% : N/A) +3 + 0 primary mapped (75.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-no-unplaced-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("placed_only.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/placed_only.flagstat")).unwrap(); + assert_eq!(flagstat, expected, "flagstat with --threads {threads}"); + fs::remove_dir_all(root).unwrap(); + } +} From fdcc7568150feef676964fdafcc991ea0e5b3e1d Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Sun, 1 Nov 2026 07:14:47 +0100 Subject: [PATCH 4/6] chore(lint): iterate corresponding count arrays Port the verified iterator-only delta from a9083c513071ef9c23d5b7e0282cd739ee672f5e. Preserve integer addition and channel order. Clears the three existing needless_range_loop warnings on the pinned native toolchain. --- src/rna/rseqc/accumulators.rs | 6 +++--- src/rna/rseqc/stats.rs | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/rna/rseqc/accumulators.rs b/src/rna/rseqc/accumulators.rs index a2ad5a13..6d0c4ad8 100644 --- a/src/rna/rseqc/accumulators.rs +++ b/src/rna/rseqc/accumulators.rs @@ -2326,9 +2326,9 @@ fn merge_vec_arrays(target: &mut Vec<[u64; N]>, source: Vec<[u64 if source.len() > target.len() { target.resize(source.len(), [0u64; N]); } - for (i, arr) in source.into_iter().enumerate() { - for j in 0..N { - target[i][j] += arr[j]; + for (target_row, source_row) in target.iter_mut().zip(source) { + for (target_count, source_count) in target_row.iter_mut().zip(source_row) { + *target_count += source_count; } } } diff --git a/src/rna/rseqc/stats.rs b/src/rna/rseqc/stats.rs index a583e006..c8af662c 100644 --- a/src/rna/rseqc/stats.rs +++ b/src/rna/rseqc/stats.rs @@ -525,13 +525,13 @@ fn write_gcc( for cycle in 0..max_cycles { let mut combined = [0u64; 6]; if cycle < fbc.len() { - for j in 0..6 { - combined[j] += fbc[cycle][j]; + for (count, &value) in combined.iter_mut().zip(&fbc[cycle]) { + *count += value; } } if cycle < lbc.len() { - for j in 0..6 { - combined[j] += lbc[cycle][j]; + for (count, &value) in combined.iter_mut().zip(&lbc[cycle]) { + *count += value; } } // Denominator is A+C+G+T only (first 4 elements), matching upstream From 74a31fffd29e666be210bd6898f6082bbfe9499a Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Sun, 1 Nov 2026 12:35:15 +0100 Subject: [PATCH 5/6] chore(bam): propagate exhausted chunks with question mark Use the equivalent Option propagation required by current Clippy. Preserve cursor exhaustion and seek-error behavior. --- src/rna/bam/io.rs | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/src/rna/bam/io.rs b/src/rna/bam/io.rs index dfa7e1b9..af86c8e2 100644 --- a/src/rna/bam/io.rs +++ b/src/rna/bam/io.rs @@ -536,22 +536,19 @@ impl Read for IndexedReader { chunk_end, reference_sequence_id, } => match chunk_end { - None => match chunks.next() { - Some(chunk) => { - if let Err(e) = - reader.get_mut().seek_to_virtual_position(chunk.start()) - { - return Some(Err(anyhow::Error::new(e) - .context("failed to seek to BAM index chunk"))); - } - *cursor = BamCursor::Region { - chunks, - chunk_end: Some(chunk.end()), - reference_sequence_id, - }; + None => { + let chunk = chunks.next()?; + if let Err(e) = reader.get_mut().seek_to_virtual_position(chunk.start()) + { + return Some(Err(anyhow::Error::new(e) + .context("failed to seek to BAM index chunk"))); } - None => return None, - }, + *cursor = BamCursor::Region { + chunks, + chunk_end: Some(chunk.end()), + reference_sequence_id, + }; + } Some(end) => { if reader.get_mut().virtual_position() >= end { // Reached the end of this chunk; advance to the next. From f45947dc3b597345748724d886e32c0300eeb7c7 Mon Sep 17 00:00:00 2001 From: Nico DV + Seqera Co-Scientist Date: Sun, 1 Nov 2026 12:35:17 +0100 Subject: [PATCH 6/6] fix(deps): update crossbeam-epoch security patch Pin crossbeam-epoch 0.9.20 using the exact version and checksum already validated in the combined candidate. Leave all other lock entries unchanged. --- Cargo.lock | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index c77b341b..ff396a4e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -375,9 +375,9 @@ dependencies = [ [[package]] name = "crossbeam-epoch" -version = "0.9.18" +version = "0.9.20" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e" +checksum = "2d6914041f254d6e9176c01941b21115dcfb7089e55135a35411081bd106ef3f" dependencies = [ "crossbeam-utils", ]