From 1d12a97e923323917f5d61ecdb3688562b4a8fb1 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Sun, 16 Aug 2026 13:29:06 -0700 Subject: [PATCH] fix(consensus)!: error on reads with absent base qualities MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A mapped read that reaches a quality-weighted consensus caller with no base qualities (BAM QUAL of '*', encoded as 0xFF per base) is structurally invalid input, not a biological filter: it cannot be weighted by the model. fgumi dropped such reads silently, folding them under ZeroLengthAfterTrimming, which can quietly degrade a molecule's depth when an upstream tool strips QUAL — with no signal to the user. create_source_read now returns Result>: Err on absent or length-mismatched qualities (the run aborts, naming the read), Ok(None) on the legitimate zero-length-after-trimming path, Ok(Some) otherwise. The vanilla and duplex production callers propagate the error; codec builds source reads through its own path and is unaffected. This matches fgbio, whose toSourceRead throws on missing base qualities. BREAKING CHANGE: inputs containing mapped reads with absent base qualities now abort consensus calling instead of silently dropping those reads. Closes #794 --- crates/fgumi-consensus/src/duplex_caller.rs | 38 +++- crates/fgumi-consensus/src/vanilla_caller.rs | 178 +++++++++++++++---- tests/integration/test_duplex_command.rs | 122 ------------- 3 files changed, 174 insertions(+), 164 deletions(-) diff --git a/crates/fgumi-consensus/src/duplex_caller.rs b/crates/fgumi-consensus/src/duplex_caller.rs index b4ec5e148..c45acefcc 100644 --- a/crates/fgumi-consensus/src/duplex_caller.rs +++ b/crates/fgumi-consensus/src/duplex_caller.rs @@ -2023,10 +2023,11 @@ impl DuplexConsensusCaller { // Keep references to original raw records for later tag extraction. // - // `create_source_read` returns `None` for reads with absent qualities or that trim to zero - // length. The single-strand path counts those as `ZeroLengthAfterTrimming` and writes them; - // building the source vectors with a bare `filter_map` would drop them from both `--stats` - // and `--rejects` (#792). Collect the dropped raws — tagged with their index in the X/Y + // `create_source_read` returns `None` for reads that trim to zero length (reads with absent + // qualities return `Err`, propagated below via `?`). The single-strand path counts the + // zero-length reads as `ZeroLengthAfterTrimming` and writes them; building the source + // vectors with a bare `filter_map` would drop them from both `--stats` and `--rejects` + // (#792). Collect the dropped raws — tagged with their index in the X/Y // partition so their group ordinal can be recovered — and hand them to the sub-caller so // they drain through the same statistics/rejects path as its alignment-filter rejections. let x_raws: Vec<&RawRecord> = ab_r1s.iter().chain(ba_r2s.iter()).copied().collect(); @@ -2035,7 +2036,7 @@ impl DuplexConsensusCaller { let mut x_sources: Vec = Vec::with_capacity(x_raws.len()); for (i, r) in x_raws.iter().enumerate() { let mate_clip = bam_fields::num_bases_extending_past_mate_raw(r); - match ss_caller.create_source_read(r, i, mate_clip) { + match ss_caller.create_source_read(r, i, mate_clip)? { Some(source) => x_sources.push(source), None => x_zero_length.push((i, *r)), } @@ -2046,7 +2047,7 @@ impl DuplexConsensusCaller { let mut y_sources: Vec = Vec::with_capacity(y_raws.len()); for (i, r) in y_raws.iter().enumerate() { let mate_clip = bam_fields::num_bases_extending_past_mate_raw(r); - match ss_caller.create_source_read(r, i, mate_clip) { + match ss_caller.create_source_read(r, i, mate_clip)? { Some(source) => y_sources.push(source), None => y_zero_length.push((i, *r)), } @@ -4239,6 +4240,31 @@ mod tests { Ok(()) } + /// #794: the duplex path must also abort on a read with absent base qualities. + /// + /// `create_source_read` is reached through the single-strand caller when building the X/Y + /// source vectors; the error must propagate out of `consensus_reads` rather than the read being + /// silently dropped, matching fgbio and the vanilla path. + #[test] + fn test_absent_base_qualities_abort_duplex_consensus() -> Result<()> { + let mut caller = rejection_accounting_caller(vec![1], true)?; + + let mut reads = duplex_molecule(2, false, false); + // One extra /A read with absent qualities (0xFF per base), which create_source_read rejects. + let mut b = SamBuilder::new(); + let cigar_10m = &[encode_op(0, 10)]; + let absent = &[0xFFu8; 10]; + reads.push(ab_r1(&mut b, b"absent", b"AAAAAAAAAA", absent, cigar_10m, b"foo/A", &[])); + reads.push(ab_r2(&mut b, b"absent", b"CCCCCCCCCC", absent, cigar_10m, b"foo/A", &[])); + + let result = caller.consensus_reads(reads); + assert!( + result.is_err(), + "a read with absent base qualities must abort the duplex path, not be silently dropped" + ); + Ok(()) + } + /// R2-UCC-01: a single MI group with *mixed* family composition — paired AB/BA /// reads plus stray unpaired fragments carrying the same MI. fgbio partitions /// the fragments out (`partition(_.paired)`) and rejects them as `NonPairedReads` diff --git a/crates/fgumi-consensus/src/vanilla_caller.rs b/crates/fgumi-consensus/src/vanilla_caller.rs index 36b04f9cb..bfbb61040 100644 --- a/crates/fgumi-consensus/src/vanilla_caller.rs +++ b/crates/fgumi-consensus/src/vanilla_caller.rs @@ -553,8 +553,8 @@ impl VanillaUmiConsensusCaller { /// Records `records` as rejected for [`RejectionReason::ZeroLengthAfterTrimming`]. /// - /// [`Self::create_source_read`] returns `None` for reads with absent qualities or that trim to - /// zero length; [`Self::process_subgroup`] counts and writes those. A composing caller (the + /// [`Self::create_source_read`] returns `None` for reads that trim to zero length; + /// [`Self::process_subgroup`] counts and writes those. A composing caller (the /// duplex path) that calls `create_source_read` directly hands the dropped records here so they /// receive the same accounting: counted on the reason breakdown and, when tracking is enabled, /// drained through [`Self::take_rejected_reads`] alongside this caller's other rejections. @@ -613,6 +613,13 @@ impl VanillaUmiConsensusCaller { return None; } + // Mirror `create_source_read`'s guard: absent qualities (0xFF per base) are not valid + // evidence, so drop such a read rather than treating it as high-quality. (`quals` is + // non-empty here, so `all` is not vacuously true.) + if quals.iter().all(|&q| q == 0xFF) { + return None; + } + if is_negative_strand { bases = reverse_complement(&bases); quals.reverse(); @@ -1021,13 +1028,18 @@ impl VanillaUmiConsensusCaller { /// 5. Trailing N removal /// 6. CIGAR transformation (reverse if negative strand, truncate to final length) /// - /// Returns None if the final length is 0 (`ZeroPostAfterTrimming`). + /// Returns `Ok(None)` if the final length is 0 (`ZeroPostAfterTrimming`), which is a legitimate + /// filtering outcome. Returns `Err` if the read has absent base qualities (BAM `QUAL` of `*`, + /// encoded as `0xFF` per base) or a quality string whose length does not match the sequence: + /// those are structurally invalid input, not a biological filter, and — matching fgbio's + /// `toSourceRead` (`UmiConsensusCaller.scala:267-270`) — abort the run rather than silently + /// dropping the read (which would mask an upstream pipeline error). pub(crate) fn create_source_read( &self, raw: &[u8], original_idx: usize, mate_overlap_clip: usize, - ) -> Option { + ) -> Result> { use fgumi_raw_bam as bam_fields; let view = RawRecordView::new(raw); @@ -1040,15 +1052,32 @@ impl VanillaUmiConsensusCaller { let mut quals = view.quality_scores().to_vec(); let read_len = bases.len(); + // A legal zero-length record (`SEQ=*`) has no bases and no qualities: it is not missing + // qualities, it simply has nothing to weight. Route it as `ZeroLengthAfterTrimming` + // (`Ok(None)`) before the quality checks below, which would otherwise treat the empty + // quality string as a length mismatch (and an empty `all(|&q| q == 0xFF)` is vacuously true). + if read_len == 0 { + return Ok(None); + } + if quals.is_empty() || quals.len() != read_len { - return None; + bail!( + "input read has invalid base qualities (length {} does not match sequence length \ + {read_len}): {}", + quals.len(), + String::from_utf8_lossy(view.read_name()), + ); } - // Per the BAM spec, absent quality scores are encoded as 0xFF per base. Reject - // such reads rather than treating 0xFF as genuine very-high-quality evidence - // (which would let unqualified reads dominate consensus). + // Per the BAM spec, absent quality scores are encoded as 0xFF per base. This is + // structurally invalid input for a quality-weighted consensus — treating 0xFF as genuine + // very-high-quality evidence would let unqualified reads dominate — so abort the run + // (matching fgbio) rather than silently dropping the read and quietly degrading depth. if quals.iter().all(|&q| q == 0xFF) { - return None; + bail!( + "input read is missing base qualities (BAM QUAL is '*'): {}", + String::from_utf8_lossy(view.read_name()), + ); } // If negative strand, reverse complement bases and reverse quals @@ -1084,7 +1113,7 @@ impl VanillaUmiConsensusCaller { } if final_len == 0 { - return None; + return Ok(None); } bases.truncate(final_len); @@ -1104,7 +1133,7 @@ impl VanillaUmiConsensusCaller { let rid = view.ref_id(); let astart = i64::from(view.pos()); - Some(SourceRead { + Ok(Some(SourceRead { original_idx, bases, quals, @@ -1114,7 +1143,7 @@ impl VanillaUmiConsensusCaller { alignment_start: astart, original_cigar, name_hash: self.source_read_name_rank(view.read_name()), - }) + })) } /// The downsampling rank for a `SourceRead`, or `0` when no cap is configured. @@ -1369,7 +1398,7 @@ impl VanillaUmiConsensusCaller { for (idx, (raw, &mate_clip)) in group_reads.iter().zip(mate_overlap_clips.iter()).enumerate() { - if let Some(sr) = self.create_source_read(raw.as_ref(), idx, mate_clip) { + if let Some(sr) = self.create_source_read(raw.as_ref(), idx, mate_clip)? { source_reads.push(sr); } else { zero_length_indices.push(idx); @@ -2239,9 +2268,11 @@ mod tests { let r2_raw = create_test_read_at(&name, false, 1, pos_r2, b"TTTTGGGGCCCC"); let sr1 = caller .create_source_read(r1_raw.as_ref(), i, 0) + .expect("valid qualities") .expect("R1 source read should be created"); let sr2 = caller .create_source_read(r2_raw.as_ref(), i, 0) + .expect("valid qualities") .expect("R2 source read should be created"); (sr1, sr2) }) @@ -2361,7 +2392,10 @@ mod tests { let srcs: Vec = (0..10) .map(|i| { let read = create_test_read(&format!("q{i}"), b"ACGT", b"####", false, false); - caller.create_source_read(read.as_ref(), i, 0).expect("source read") + caller + .create_source_read(read.as_ref(), i, 0) + .expect("valid qualities") + .expect("source read") }) .collect(); @@ -3570,7 +3604,7 @@ mod tests { b.build() }; - let source = caller.create_source_read(record.as_ref(), 0, 0); + let source = caller.create_source_read(record.as_ref(), 0, 0).expect("valid qualities"); assert!(source.is_some(), "Should produce a SourceRead"); let sr = source.expect("source read should be Some"); @@ -3608,8 +3642,10 @@ mod tests { VanillaUmiConsensusOptions { min_reads: 1, max_reads, ..Default::default() }, ); let read = create_test_read(name, b"ACGT", b"####", false, false); - let sr = - caller.create_source_read(read.as_ref(), 0, 0).expect("source read should be created"); + let sr = caller + .create_source_read(read.as_ref(), 0, 0) + .expect("valid qualities") + .expect("source read should be created"); assert_eq!(sr.name_hash, expected); } @@ -3634,7 +3670,7 @@ mod tests { b.build() }; - let source = caller.create_source_read(record.as_ref(), 0, 0); + let source = caller.create_source_read(record.as_ref(), 0, 0).expect("valid qualities"); assert!(source.is_some(), "Should produce a SourceRead"); let sr = source.expect("source read should be Some"); @@ -3664,7 +3700,7 @@ mod tests { b.build() }; - let source = caller.create_source_read(record.as_ref(), 0, 0); + let source = caller.create_source_read(record.as_ref(), 0, 0).expect("valid qualities"); assert!(source.is_some(), "Should produce a SourceRead"); let sr = source.expect("source read should be Some"); @@ -3696,7 +3732,7 @@ mod tests { b.build() }; - let source = caller.create_source_read(record.as_ref(), 0, 0); + let source = caller.create_source_read(record.as_ref(), 0, 0).expect("valid qualities"); assert!(source.is_some(), "Should produce a SourceRead"); let sr = source.expect("source read should be Some"); @@ -3731,7 +3767,7 @@ mod tests { b.build() }; - let source = caller.create_source_read(record.as_ref(), 0, 0); + let source = caller.create_source_read(record.as_ref(), 0, 0).expect("valid qualities"); // fgbio expected: None assert!(source.is_none(), "Should return None when all bases are masked or N"); } @@ -3768,7 +3804,7 @@ mod tests { let clip = num_bases_extending_past_mate_raw(r1.as_ref()); assert_eq!(clip, 0, "FF pair should not trigger mate overlap clipping"); - let source = caller.create_source_read(r1.as_ref(), 0, clip); + let source = caller.create_source_read(r1.as_ref(), 0, clip).expect("valid qualities"); assert!(source.is_some(), "Should produce a SourceRead"); let sr = source.expect("source read should be Some"); @@ -3805,7 +3841,7 @@ mod tests { }; let clip = num_bases_extending_past_mate_raw(r1.as_ref()); - let source = caller.create_source_read(r1.as_ref(), 0, clip); + let source = caller.create_source_read(r1.as_ref(), 0, clip).expect("valid qualities"); assert!(source.is_some(), "Should produce a SourceRead"); let sr = source.expect("source read should be Some"); @@ -3845,7 +3881,7 @@ mod tests { // R1 ends at 148, mate ends at 168, so R1 doesn't extend past mate assert_eq!(clip, 0, "R1 should not extend past mate"); - let source = caller.create_source_read(r1.as_ref(), 0, clip); + let source = caller.create_source_read(r1.as_ref(), 0, clip).expect("valid qualities"); assert!(source.is_some()); assert_eq!(source.expect("source should be Some").bases.len(), 50); } @@ -3881,7 +3917,7 @@ mod tests { // R1 ends at 148, mate ends at 128, so R1 extends 20 bases past mate assert_eq!(clip, 20, "R1 should extend 20 bases past mate"); - let source = caller.create_source_read(r1.as_ref(), 0, clip); + let source = caller.create_source_read(r1.as_ref(), 0, clip).expect("valid qualities"); assert!(source.is_some()); let sr = source.expect("source read should be Some"); assert_eq!(sr.bases.len(), 30, "Should be trimmed to 30 bases"); @@ -3997,7 +4033,8 @@ mod tests { // R1 extends 10 bases past mate's end assert_eq!(clip_r1, 10, "R1 should extend 10 bases past mate"); - let source_r1 = caller.create_source_read(r1.as_ref(), 0, clip_r1); + let source_r1 = + caller.create_source_read(r1.as_ref(), 0, clip_r1).expect("valid qualities"); assert!(source_r1.is_some(), "R1 should produce SourceRead"); let sr1 = source_r1.expect("failed to get sr1"); @@ -4033,7 +4070,8 @@ mod tests { // R2's first 10 bases (positions 1-10) extend before mate's start (11) assert_eq!(clip_r2, 10, "R2 should extend 10 bases before mate start"); - let source_r2 = caller.create_source_read(r2.as_ref(), 0, clip_r2); + let source_r2 = + caller.create_source_read(r2.as_ref(), 0, clip_r2).expect("valid qualities"); assert!(source_r2.is_some(), "R2 should produce SourceRead"); let sr2 = source_r2.expect("failed to get sr2"); @@ -4075,7 +4113,8 @@ mod tests { let clip_r1 = num_bases_extending_past_mate_raw(r1.as_ref()); - let source_r1 = caller.create_source_read(r1.as_ref(), 0, clip_r1); + let source_r1 = + caller.create_source_read(r1.as_ref(), 0, clip_r1).expect("valid qualities"); assert!(source_r1.is_some(), "R1 should produce SourceRead"); let sr1 = source_r1.expect("failed to get sr1"); @@ -4116,7 +4155,8 @@ mod tests { let clip_r1 = num_bases_extending_past_mate_raw(r1.as_ref()); - let source_r1 = caller.create_source_read(r1.as_ref(), 0, clip_r1); + let source_r1 = + caller.create_source_read(r1.as_ref(), 0, clip_r1).expect("valid qualities"); assert!(source_r1.is_some(), "R1 should produce SourceRead"); let sr1 = source_r1.expect("failed to get sr1"); @@ -4164,7 +4204,8 @@ mod tests { // R2 unclipped end is 57 (20+30-1+8=57) // R1 extends 2 bases past R2's end - let source_r1 = caller.create_source_read(r1.as_ref(), 0, clip_r1); + let source_r1 = + caller.create_source_read(r1.as_ref(), 0, clip_r1).expect("valid qualities"); assert!(source_r1.is_some(), "R1 should produce SourceRead"); let sr1 = source_r1.expect("failed to get sr1"); @@ -4666,7 +4707,7 @@ mod tests { let caller = VanillaUmiConsensusCaller::new("test".to_string(), "UMI1".to_string(), options); - let source = caller.create_source_read(record.as_ref(), 0, 0); + let source = caller.create_source_read(record.as_ref(), 0, 0).expect("valid qualities"); assert!(source.is_some(), "Should produce a source read"); let sr = source.expect("source read should be Some"); @@ -4680,12 +4721,13 @@ mod tests { // Test 40: Exception when reads lack base qualities // ========================================================================= - /// Port of fgbio test: "except when the reads do not have base qualities" - /// Tests that reads without quality scores are handled appropriately + /// Port of fgbio test: "except when the reads do not have base qualities" (#794). + /// fgbio's `toSourceRead` throws when a read has no base qualities; `create_source_read` must + /// likewise return an error rather than silently dropping the read. #[test] fn test_reads_without_base_qualities() { // The BAM spec represents absent quality scores as 0xFF per base. SamBuilder - // emits 0xFF when qualities are not set. create_source_read must reject such + // emits 0xFF when qualities are not set. create_source_read must error on such // records rather than treating 0xFF as usable very-high-quality evidence. let record = { let mut b = SamBuilder::new(); @@ -4705,8 +4747,72 @@ mod tests { let result = caller.create_source_read(record.as_ref(), 0, 0); assert!( - result.is_none(), - "reads without base qualities should be rejected, not treated as usable qualities" + result.is_err(), + "reads without base qualities must error, not be silently rejected" + ); + } + + /// #794: a read with absent base qualities must abort consensus calling, matching fgbio. + /// + /// fgbio's `toSourceRead` throws `IllegalArgumentException` when a read has no base qualities + /// (`UmiConsensusCaller.scala:267-270`). fgumi previously dropped such reads silently, which can + /// mask an upstream pipeline error (a tool that stripped `QUAL`) by quietly degrading a + /// molecule's depth. This drives the whole consensus path — not just `create_source_read` — so + /// it fails loudly regardless of which caller sees the read. + #[test] + fn test_absent_base_qualities_abort_consensus() { + // The BAM spec encodes absent qualities as 0xFF per base; SamBuilder emits that when + // qualities are not set. + let record = { + let mut b = SamBuilder::new(); + b.read_name(b"test") + .flags(flags::PAIRED | flags::FIRST_SEGMENT) + .ref_id(0) + .pos(0) + .sequence(b"AAAAAAAAAA") + .cigar_ops(&[encode_op(0, 10)]) + .add_string_tag(SamTag::MI, b"UMI1"); + b.build() + }; + + let options = VanillaUmiConsensusOptions { min_reads: 1, ..Default::default() }; + let mut caller = + VanillaUmiConsensusCaller::new("test".to_string(), "UMI1".to_string(), options); + + let result = caller.consensus_reads(vec![record]); + assert!( + result.is_err(), + "a read with absent base qualities must abort the run, not be silently dropped" + ); + } + + /// A legal zero-length record (`SEQ=*`) has no bases and no qualities. It must be routed as + /// `ZeroLengthAfterTrimming` (`Ok(None)`) rather than tripping the absent-qualities guard: an + /// empty quality string is not "missing qualities on a read with bases", and an empty + /// `all(|&q| q == 0xFF)` is vacuously true. `process_subgroup` calls `create_source_read` before + /// unmapped filtering, so erroring here would abort otherwise-valid `--allow-unmapped` runs. + #[test] + fn test_zero_length_record_is_not_treated_as_absent_qualities() { + let record = { + let mut b = SamBuilder::new(); + // SEQ=* / QUAL=*: no bases, no qualities. + b.read_name(b"empty") + .flags(flags::UNMAPPED) + .sequence(b"") + .add_string_tag(SamTag::MI, b"UMI1"); + b.build() + }; + + let options = VanillaUmiConsensusOptions::default(); + let caller = + VanillaUmiConsensusCaller::new("test".to_string(), "UMI1".to_string(), options); + + let source = caller + .create_source_read(record.as_ref(), 0, 0) + .expect("a zero-length record must not error, unlike an absent-quality read"); + assert!( + source.is_none(), + "a zero-length record must route as ZeroLengthAfterTrimming (None)" ); } diff --git a/tests/integration/test_duplex_command.rs b/tests/integration/test_duplex_command.rs index 010b94c07..55369a411 100644 --- a/tests/integration/test_duplex_command.rs +++ b/tests/integration/test_duplex_command.rs @@ -632,128 +632,6 @@ fn test_duplex_zero_length_rejects_reconcile_with_stats(#[case] threads: Option< ); } -/// #792: reads with absent qualities must reach both the rejects BAM and `--stats`. -/// -/// `create_source_read` returns `None` for two distinct reasons — a read that quality-trims to -/// zero length, and a read whose qualities are absent (all `0xFF`, the BAM spec's sentinel). -/// `test_duplex_zero_length_rejects_reconcile_with_stats` covers the first; this covers the -/// second, which the duplex path also counts as `ZeroLengthAfterTrimming` and writes to -/// `--rejects`. Absent qualities are rejected before any trimming, so `--trim` is deliberately -/// not passed here — the sentinel check alone must route the pair. -/// -/// The `absent_qual` template's two reads have all-`0xFF` qualities while the three majority -/// templates per strand still consense. Asserted single-threaded and threaded, since they drain -/// rejects and statistics through different code. -#[rstest] -#[case::single_threaded(None)] -#[case::threaded(Some("2"))] -fn test_duplex_absent_quality_rejects_reconcile_with_stats(#[case] threads: Option<&str>) { - let temp_dir = TempDir::new().unwrap(); - let input_bam = temp_dir.path().join("input.bam"); - let output_bam = temp_dir.path().join("output.bam"); - let rejects_bam = temp_dir.path().join("rejects.bam"); - let stats_path = temp_dir.path().join("stats.txt"); - - let mut molecule = create_duplex_molecule("1", "ACGTACGT", 30, 100, 3); - // One extra /A template whose qualities are all 0xFF — the BAM spec's "quality absent" - // sentinel — so create_source_read rejects it regardless of trimming. - molecule.push(create_duplex_read_pair_with_cigar( - "absent_qual", - "1/A", - "ACGTACGT", - "ACGTACGT", - 0xFF, - 100, - false, - &[8u32 << 4], - )); - create_duplex_bam(&input_bam, vec![molecule]); - - let mut args = vec![ - "duplex", - "--input", - input_bam.to_str().unwrap(), - "--output", - output_bam.to_str().unwrap(), - "--rejects", - rejects_bam.to_str().unwrap(), - "--stats", - stats_path.to_str().unwrap(), - "--min-reads", - "1", - "--compression-level", - "1", - ]; - if let Some(threads) = threads { - args.extend_from_slice(&["--threads", threads]); - } - Duplex::try_parse_from(args) - .expect("failed to parse duplex args") - .execute("fgumi duplex") - .expect("Duplex command failed"); - - let mut output_reader = bam::io::Reader::new(fs::File::open(&output_bam).unwrap()); - output_reader.read_header().expect("read output header"); - assert_eq!( - output_reader.records().count(), - 2, - "the majority-alignment reads must still emit a duplex R1/R2 pair" - ); - - // Key input records by (name, flags) so each reject is compared byte-for-byte against the - // exact bytes it came from. - let expected_rejects: HashMap<(Vec, u16), Vec> = read_raw_records(&input_bam) - .into_iter() - .filter(|record| fgumi_raw_bam::read_name(record) == b"absent_qual") - .map(|record| { - let view = RawRecordView::new(&record); - ((view.read_name().to_vec(), view.flags()), record) - }) - .collect(); - assert_eq!(expected_rejects.len(), 2, "the absent-quality template contributes two records"); - - let mut seen: Vec<(Vec, u16)> = Vec::new(); - for record in read_raw_records(&rejects_bam) { - let view = RawRecordView::new(&record); - let key = (view.read_name().to_vec(), view.flags()); - let expected = expected_rejects.get(&key).unwrap_or_else(|| { - panic!("unexpected reject record {}", String::from_utf8_lossy(&key.0)) - }); - assert_eq!( - &record, expected, - "each reject must be byte-for-byte identical to its input record" - ); - assert!(!seen.contains(&key), "a reject record was written more than once"); - seen.push(key); - } - - let stats = fs::read_to_string(&stats_path).expect("read stats"); - let stat = |key: &str| -> usize { - stats - .lines() - .find_map(|line| line.strip_prefix(&format!("{key}\t"))) - .and_then(|rest| rest.split('\t').next()) - .unwrap_or_else(|| panic!("stats must contain a {key} row")) - .parse() - .unwrap_or_else(|_| panic!("{key} must be an integer")) - }; - assert_eq!( - seen.len(), - expected_rejects.len(), - "the rejects BAM must hold both ends of the absent-quality template" - ); - assert_eq!( - stat("raw_reads_rejected_for_zero_bases_post_trimming"), - 2, - "both ends of the absent-quality template must be counted as zero-bases-post-trimming" - ); - assert_eq!( - seen.len(), - stat("raw_reads_rejected"), - "the rejects BAM record count must equal raw_reads_rejected" - ); -} - /// #792 (follow-up): zero-length rejects must reach the rejects BAM in input order. /// /// A read that trims to zero length is dropped by `create_source_read` while its X/Y source