diff --git a/crates/fgumi-consensus/src/duplex_caller.rs b/crates/fgumi-consensus/src/duplex_caller.rs index 2190ec93a..3797fa329 100644 --- a/crates/fgumi-consensus/src/duplex_caller.rs +++ b/crates/fgumi-consensus/src/duplex_caller.rs @@ -205,6 +205,7 @@ use crate::caller::{ ConsensusCaller, ConsensusCallingStats, RejectionReason, clamp_combined_error_to_fgbio_short, clamp_per_base_to_fgbio_short, consensus_read_name_too_long_context, write_consensus_read_name, }; +use crate::mate_clip::{MateClipper, missing_mate_error}; use crate::methylation::ConversionPattern; use crate::phred::MAX_PHRED; use crate::phred::{MIN_PHRED, PhredScore}; @@ -213,7 +214,7 @@ use crate::vanilla_caller::{ VanillaConsensusRead, VanillaUmiConsensusCaller, VanillaUmiConsensusOptions, }; use crate::{ReadType, SourceRead}; -use fgumi_raw_bam::{self as bam_fields, RawRecord, RawRecordView, UnmappedSamBuilder, flags}; +use fgumi_raw_bam::{RawRecord, RawRecordView, UnmappedSamBuilder, flags}; use fgumi_sam::SamTag; use std::collections::HashSet; @@ -2164,12 +2165,19 @@ impl DuplexConsensusCaller { // (#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. + // + // Read-through clipping looks up the mate of a read lacking `MC` across the whole group: + // its mate always lands in the other R1/R2 partition. It is set up only now, after the + // whole-group rejections above, so a read whose needed mate is missing fails the run only + // in a group that passed its minimum-reads and strand-orientation checks. + let mut mate_clipper = + MateClipper::new(a_records.iter().chain(b_records.iter()).map(AsRef::as_ref)); let x_raws: Vec<&RawRecord> = ab_r1s.iter().chain(ba_r2s.iter()).copied().collect(); let mut x_zero_length: Vec<(usize, &RawRecord)> = Vec::new(); 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); + let mate_clip = mate_clipper.clip(r).ok_or_else(|| missing_mate_error(r))?; match ss_caller.create_source_read(r, i, mate_clip)? { Some(source) => x_sources.push(source), None => x_zero_length.push((i, *r)), @@ -2180,7 +2188,7 @@ impl DuplexConsensusCaller { let mut y_zero_length: Vec<(usize, &RawRecord)> = Vec::new(); 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); + let mate_clip = mate_clipper.clip(r).ok_or_else(|| missing_mate_error(r))?; match ss_caller.create_source_read(r, i, mate_clip)? { Some(source) => y_sources.push(source), None => y_zero_length.push((i, *r)), @@ -2750,6 +2758,9 @@ impl ConsensusCaller for DuplexConsensusCaller { )] mod tests { use super::*; + use crate::mate_clip::read_through_fixture::{ + McTags, READ_THROUGH_ADAPTER, READ_THROUGH_INSERT, read_through_pair, + }; use fgumi_raw_bam::{ParsedBamRecord, SamBuilder, testutil::encode_op}; use noodles::sam::alignment::record_buf::data::field::Value; use rstest::rstest; @@ -4165,6 +4176,151 @@ mod tests { Ok(()) } + /// Duplex caller (`--min-reads min_reads`) for the mate-CIGAR backfill tests. + fn backfill_caller(min_reads: Vec) -> Result { + DuplexConsensusCaller::new( + "consensus".to_string(), + "A".to_string(), + min_reads, + 10, + false, + false, + None, + None, + false, + 45, + 40, + ) + } + + /// Port of fgbio test: "add the mate cigar when not present before consensus calling" + /// (`VanillaUmiConsensusCallerTest.scala:818`), applied to the duplex caller and extended to + /// assert the effect of the backfill. + /// + /// fgbio's `UmiConsensusCaller.updateMateCigars` (`UmiConsensusCaller.scala:328-343`, called + /// at `:356`) runs for `CallDuplexConsensusReads` too, filling each read's `MC` from its + /// mate before any read-through trim. Three AB and three BA read-through pairs (80 bp insert, + /// 20 bp adapter) must therefore yield 80 bp duplex consensus reads with depths + /// `aD=3 bD=3 cD=6`, whether or not the input carries `MC`. A pair that already carries `MC` + /// is used as-is, even when stale (fgumi does not mirror fgbio's overwrite of an existing + /// `MC`). + /// + /// The `forward_reads_missing_mc` case drops `MC` from every forward-strand read (AB-R1 and + /// BA-R2), which together form the duplex R1. Dropping it from only one of them would not + /// test the backfill: the duplex is trimmed to the shorter single-strand consensus, so the + /// trimmed partner would mask the untrimmed one. + #[rstest] + #[case::both_mc_present(McTags::Both, false)] + #[case::neither_mc_backfilled(McTags::Neither, false)] + #[case::forward_reads_missing_mc_backfilled(McTags::ForwardMissing, false)] + #[case::stale_mc_is_used_not_overwritten(McTags::StaleBoth, true)] + fn test_duplex_mate_cigar_backfilled_before_consensus_trims_read_through( + #[case] mc: McTags, + #[case] r1_keeps_adapter: bool, + ) -> Result<()> { + let mut reads = Vec::new(); + for i in 0..3 { + let (r1, r2) = read_through_pair(format!("ab{i}").as_bytes(), b"1/A", mc, true); + reads.push(r1); + reads.push(r2); + } + for i in 0..3 { + let (r1, r2) = read_through_pair(format!("ba{i}").as_bytes(), b"1/B", mc, false); + reads.push(r1); + reads.push(r2); + } + + let mut caller = backfill_caller(vec![1])?; + let output = caller.consensus_reads(reads)?; + + assert_eq!(output.count, 2, "one duplex consensus read per end"); + let consensus = ParsedBamRecord::parse_all(&output.data); + let r1 = consensus.iter().find(|r| r.flag & flags::FIRST_SEGMENT != 0).expect("R1"); + let r2 = consensus.iter().find(|r| r.flag & flags::LAST_SEGMENT != 0).expect("R2"); + + // R1 combines the forward reads (AB-R1 + BA-R2): the adapter is their 3' end. R2 combines + // the reverse reads and is emitted in sequencing orientation, so the insert appears + // reverse-complemented. With the stale `100M` mate CIGAR the forward reads see a mate that + // ends 20 bp past the insert and keep their adapter (as fgbio does), while the reverse + // reads' clip depends only on the mate start and is unaffected. + let mut expected_r1 = READ_THROUGH_INSERT.to_vec(); + if r1_keeps_adapter { + expected_r1.extend_from_slice(READ_THROUGH_ADAPTER); + } + assert_eq!( + String::from_utf8(r1.bases.clone())?, + String::from_utf8(expected_r1)?, + "R1 duplex consensus bases" + ); + assert_eq!( + String::from_utf8(r2.bases.clone())?, + String::from_utf8(fgumi_dna::dna::reverse_complement(READ_THROUGH_INSERT))?, + "R2 duplex consensus bases" + ); + for r in [r1, r2] { + assert_eq!(r.get_int_tag(SamTag::AD), Some(3)); + assert_eq!(r.get_int_tag(SamTag::BD), Some(3)); + assert_eq!(r.get_int_tag(SamTag::CD), Some(6)); + } + Ok(()) + } + + /// Port of the fgbio behavior at `UmiConsensusCaller.scala:328-343` (`updateMateCigars`) + /// exercised by `VanillaUmiConsensusCallerTest.scala:818`, for the duplex caller: a read with + /// no `MC` whose mate is not in the group cannot have its mate CIGAR backfilled. fgbio throws + /// (a `MatchError`); fgumi must stop with an error that names the read. + #[test] + fn test_duplex_missing_mate_without_mc_errors_naming_the_read() -> Result<()> { + let mut reads = Vec::new(); + let (orphan_r1, _orphan_r2) = + read_through_pair(b"ab_orphan", b"1/A", McTags::Neither, true); + reads.push(orphan_r1); + let (ba_r1, ba_r2) = read_through_pair(b"ba0", b"1/B", McTags::Neither, false); + reads.push(ba_r1); + reads.push(ba_r2); + + let mut caller = backfill_caller(vec![1])?; + let err = caller + .consensus_reads(reads) + .err() + .expect("a read lacking MC whose mate is absent must be an error"); + + assert_eq!( + format!("{err:#}"), + "Mate cigar (MC SAM tag) needed for read 'ab_orphan': the read has no MC tag and \ + its primary mate is not in the same group. Add MC tags (e.g. with fgumi zipper or \ + samtools fixmate), or keep both reads of each template in the same group." + ); + Ok(()) + } + + /// A read lacking `MC` whose mate is missing fails the run only in a group that passes the + /// duplex whole-group rejections. Here the AB R1 lacks `MC` and its R2 is absent, but the + /// group is rejected first, so it is counted as rejected rather than failing the run: + /// - `potential_collision`: the BA template is on the AB orientation (its R1 forward), so the + /// duplex R1 holds both a forward and a reverse read; + /// - `too_few_reads`: `--min-reads 2` needs two R1s on each strand, and each strand has one. + #[rstest] + #[case::potential_collision(vec![1], true, RejectionReason::PotentialCollision)] + #[case::too_few_reads(vec![2], false, RejectionReason::InsufficientReads)] + fn test_duplex_missing_mate_is_not_an_error_in_a_rejected_group( + #[case] min_reads: Vec, + #[case] ba_r1_forward: bool, + #[case] expected_reason: RejectionReason, + ) -> Result<()> { + let (orphan_r1, _) = read_through_pair(b"ab_orphan", b"1/A", McTags::Neither, true); + let (ba_r1, ba_r2) = read_through_pair(b"ba0", b"1/B", McTags::Neither, ba_r1_forward); + + let mut caller = backfill_caller(min_reads)?; + let output = caller.consensus_reads(vec![orphan_r1, ba_r1, ba_r2])?; + + assert_eq!(output.count, 0); + let rejections: Vec<(RejectionReason, usize)> = + caller.stats.rejection_reasons.iter().map(|(reason, n)| (*reason, *n)).collect(); + assert_eq!(rejections, vec![(expected_reason, 3)]); + Ok(()) + } + /// Builds one duplex molecule, optionally seeded with minority-alignment templates. /// /// Each strand gets `majority_templates` templates carrying `10M`. When `gapped_a` (or diff --git a/crates/fgumi-consensus/src/lib.rs b/crates/fgumi-consensus/src/lib.rs index 8197659df..4bb00235f 100644 --- a/crates/fgumi-consensus/src/lib.rs +++ b/crates/fgumi-consensus/src/lib.rs @@ -24,6 +24,9 @@ pub mod sequence; pub mod simple_umi; pub mod tags; +#[cfg(feature = "simplex")] +pub(crate) mod mate_clip; + #[cfg(feature = "simplex")] pub mod vanilla_caller; diff --git a/crates/fgumi-consensus/src/mate_clip.rs b/crates/fgumi-consensus/src/mate_clip.rs new file mode 100644 index 000000000..ead6f5746 --- /dev/null +++ b/crates/fgumi-consensus/src/mate_clip.rs @@ -0,0 +1,484 @@ +//! Read-through clipping for the simplex and duplex callers: the bases at a read's 3' end that +//! extend past its mate, with the mate CIGAR taken from the read's `MC` tag or, when that tag is +//! absent, from the mate record in the same group. + +use ahash::AHashMap; +use anyhow::anyhow; +use fgumi_raw_bam::{RawRecordView, flags}; + +/// Computes the read-through clip — the bases a read extends past its mate — for the reads of +/// one group, taking the mate CIGAR from the read's `MC` tag or, when that tag is absent, from +/// the mate record in the group. +/// +/// fgbio's `UmiConsensusCaller.updateMateCigars` (`UmiConsensusCaller.scala:328-343`, called at +/// `:356`) fills in `MC` from the in-group mate before any read is trimmed, so input without +/// `MC` is still clipped at the mate's end; without the backfill the clip would be 0 and +/// read-through adapter would reach the consensus. Deliberate differences from fgbio: +/// - a read that carries `MC` is measured against that tag, even when it disagrees with the +/// mate's actual CIGAR (fgbio overwrites both reads' `MC` when either read lacks it); +/// - a read lacking `MC` whose mate is not in the group gets `None` from [`Self::clip`], which +/// the callers turn into an error naming the read ([`missing_mate_error`]), and a name seen +/// more than twice uses the first mate found (fgbio fails both with a `MatchError`); +/// - a mate is needed by every primary read mapped to the same reference as its mate on the +/// opposite strand, RF (outie) pairs included, where fgbio only backfills, and so only fails +/// on, pairs it judges FR from `TLEN` — which misclassifies dovetailed FR pairs; +/// - the callers fail the run on a missing mate only for a read that can still reach a +/// consensus — simplex for a read in an emitted consensus, duplex for any read of a group that +/// passes its minimum-reads and strand-orientation checks — whereas fgbio backfills every +/// group before any filtering (`UmiConsensusCaller.scala:356`), so it also fails on groups it +/// then rejects; +/// - an in-group mate that is unmapped, or whose CIGAR is not a valid CIGAR for a mapped +/// record, gives a clip of 0, as an invalid `MC` tag does. +/// +/// Only primary reads are backfilled or indexed as mates; a secondary or supplementary read +/// keeps its own `MC` tag (or no clip). fgbio never backfills them either: its +/// `ConsensusCallingIterator` drops secondary and supplementary reads before any caller sees +/// the group (`ConsensusCallingIterator.scala:57`). +/// +/// The mate index is built on the first read that needs it, so input that carries `MC` — the +/// norm — pays only the one `MC` lookup per read that clipping always needed. +pub(crate) struct MateClipper<'a, I> { + /// Every read of the group: where `mates` is built from. + records: I, + /// Primary paired reads keyed by `(read name, is first segment)`; built on first use. + mates: Option>, + /// Reused buffer for the decoded CIGAR of the mate being clipped against. + mate_ops: Vec, +} + +/// Primary paired reads of a group keyed by `(read name, is first segment)`. +type MatesByName<'a> = AHashMap<(&'a [u8], bool), &'a [u8]>; + +impl<'a, I> MateClipper<'a, I> +where + I: Iterator + Clone, +{ + /// Creates a clipper over `records`, all reads of one group. Nothing is indexed yet. + pub(crate) fn new(records: I) -> Self { + Self { records, mates: None, mate_ops: Vec::new() } + } + + /// Returns the number of bases `rec` extends past its mate (0 unless an FR pair), or `None` + /// when `rec` has no `MC` tag, needs its mate's CIGAR, and that mate is not in the group. + pub(crate) fn clip(&mut self, rec: &[u8]) -> Option { + if let Some(clip) = fgumi_raw_bam::num_bases_extending_past_mate_from_mc_raw(rec) { + return Some(clip); + } + if !needs_mate(rec) { + return Some(0); + } + let records = &self.records; + let mates = self.mates.get_or_insert_with(|| index_mates(records.clone())); + let mate = RawRecordView::new(mates.get(&mate_key(rec))?); + if mate.flags() & flags::UNMAPPED != 0 { + return Some(0); + } + self.mate_ops.clear(); + self.mate_ops.extend(mate.cigar_ops_iter()); + Some(fgumi_raw_bam::num_bases_extending_past_mate_cigar_raw(rec, &self.mate_ops)) + } +} + +/// The error for a read whose [`MateClipper::clip`] is `None`: it needs its mate's CIGAR and has +/// neither an `MC` tag nor its mate in the group. +pub(crate) fn missing_mate_error(rec: &[u8]) -> anyhow::Error { + let name = String::from_utf8_lossy(RawRecordView::new(rec).read_name()); + anyhow!( + "Mate cigar (MC SAM tag) needed for read '{name}': the read has no MC tag and its primary \ + mate is not in the same group. Add MC tags (e.g. with fgumi zipper or samtools fixmate), \ + or keep both reads of each template in the same group." + ) +} + +/// Indexes the primary paired reads of `records` by name and segment, keeping the first read +/// seen for each key. +fn index_mates<'a>(records: impl Iterator) -> MatesByName<'a> { + let mut mates = MatesByName::new(); + for rec in records { + if RawRecordView::new(rec).flags() & flags::PAIRED != 0 && is_primary(rec) { + mates.entry(key(rec)).or_insert(rec); + } + } + mates +} + +fn is_primary(rec: &[u8]) -> bool { + RawRecordView::new(rec).flags() & (flags::SECONDARY | flags::SUPPLEMENTARY) == 0 +} + +/// True when `rec`, a read lacking `MC`, could still extend past its mate: a primary read and its +/// mate both mapped to the same reference on opposite strands. FR orientation is not required +/// here because, without the mate CIGAR, it can only be inferred from `TLEN`, which +/// misclassifies dovetailed pairs (htsjdk/samtools#1771). +fn needs_mate(rec: &[u8]) -> bool { + is_primary(rec) && fgumi_raw_bam::is_mapped_opposite_strand_pair_raw(rec) +} + +fn key(rec: &[u8]) -> (&[u8], bool) { + let view = RawRecordView::new(rec); + (view.read_name(), view.flags() & flags::FIRST_SEGMENT != 0) +} + +fn mate_key(rec: &[u8]) -> (&[u8], bool) { + let (name, is_first) = key(rec); + (name, !is_first) +} + +/// Read-through pair fixture shared by the mate-clip, simplex and duplex tests. +#[cfg(test)] +pub(crate) mod read_through_fixture { + use fgumi_dna::dna::reverse_complement; + use fgumi_raw_bam::{RawRecord, SamBuilder, encode_op, flags}; + use fgumi_sam::SamTag; + + /// 80 bp genomic insert of the read-through pair. + pub(crate) const READ_THROUGH_INSERT: &[u8] = + b"ACCTCTCCATCTGACCCAAGATTGTGCTTGTTCAATTCTTCTTAACGTGATAACAGAATCAAACCTGCCAGGCGGTCGTC"; + /// 20 bp adapter sequenced past the end of the insert. + pub(crate) const READ_THROUGH_ADAPTER: &[u8] = b"AGATCGGAAGAGCACACGTC"; + + /// Which reads of a read-through pair carry an `MC` tag. + #[derive(Clone, Copy, Debug)] + pub(crate) enum McTags { + /// Both reads carry the correct mate CIGAR. + Both, + /// Neither read carries an `MC` tag. + Neither, + /// Only the reverse-strand read carries `MC`; the forward-strand read lacks it. + ForwardMissing, + /// Both reads carry a (stale) `100M`, which is not the mate's real CIGAR. + StaleBoth, + } + + /// Builds a read-through FR pair `(R1, R2)`: a forward `80M20S` read and a reverse `20S80M` + /// read, both 100 bp starting at 1-based position 1001, so each reads 20 bp of adapter past + /// its mate. When `r1_forward` is set R1 is the forward read (a simplex or AB-strand + /// template), otherwise R2 is (a BA-strand template). + pub(crate) fn read_through_pair( + name: &[u8], + mi: &[u8], + mc: McTags, + r1_forward: bool, + ) -> (RawRecord, RawRecord) { + let mut fwd_seq = READ_THROUGH_INSERT.to_vec(); + fwd_seq.extend_from_slice(READ_THROUGH_ADAPTER); + let mut rev_seq = reverse_complement(READ_THROUGH_ADAPTER); + rev_seq.extend_from_slice(READ_THROUGH_INSERT); + // (forward MC, reverse MC): each read's MC names its mate's CIGAR. + let (fwd_mc, rev_mc): (Option<&[u8]>, Option<&[u8]>) = match mc { + McTags::Both => (Some(b"20S80M"), Some(b"80M20S")), + McTags::Neither => (None, None), + McTags::ForwardMissing => (None, Some(b"80M20S")), + McTags::StaleBoth => (Some(b"100M"), Some(b"100M")), + }; + let build = |segment: u16, reverse: bool, mc: Option<&[u8]>| -> RawRecord { + let (strand, tlen, seq, cigar) = if reverse { + (flags::REVERSE, -80, &rev_seq, [encode_op(4, 20), encode_op(0, 80)]) + } else { + (flags::MATE_REVERSE, 80, &fwd_seq, [encode_op(0, 80), encode_op(4, 20)]) + }; + let mut b = SamBuilder::new(); + b.read_name(name) + .ref_id(0) + .pos(1000) + .flags(flags::PAIRED | flags::PROPER_PAIR | segment | strand) + .mate_ref_id(0) + .mate_pos(1000) + .template_length(tlen) + .sequence(seq) + .qualities(&[40; 100]) + .cigar_ops(&cigar) + .add_string_tag(SamTag::MI, mi); + if let Some(mc) = mc { + b.add_string_tag(SamTag::MC, mc); + } + b.build() + }; + let (r1_mc, r2_mc) = if r1_forward { (fwd_mc, rev_mc) } else { (rev_mc, fwd_mc) }; + ( + build(flags::FIRST_SEGMENT, !r1_forward, r1_mc), + build(flags::LAST_SEGMENT, r1_forward, r2_mc), + ) + } +} + +#[cfg(test)] +mod tests { + use super::read_through_fixture::{McTags, read_through_pair}; + use super::*; + use fgumi_raw_bam::{RawRecord, SamBuilder, encode_op}; + use fgumi_sam::SamTag; + use rstest::rstest; + + /// Clips every read of `records` with one [`MateClipper`] over the group. + fn clips(records: &[RawRecord]) -> Vec> { + let mut clipper = MateClipper::new(records.iter().map(AsRef::as_ref)); + records.iter().map(|r| clipper.clip(r.as_ref())).collect() + } + + /// One read of a same-reference pair for the [`MateClipper`] tests, with `cigar` at 0-based + /// `pos` and its mate at 0-based `mate_pos` on the opposite strand (the read is forward unless + /// `reverse`), with `extra_flags` OR-ed in. + #[expect(clippy::too_many_arguments, reason = "test fixture builder")] + fn clipper_read( + name: &[u8], + first: bool, + reverse: bool, + pos: i32, + mate_pos: i32, + cigar: &[u32], + mc: Option<&[u8]>, + extra_flags: u16, + ) -> RawRecord { + let segment = if first { flags::FIRST_SEGMENT } else { flags::LAST_SEGMENT }; + let strand = if reverse { flags::REVERSE } else { flags::MATE_REVERSE }; + let len: usize = cigar + .iter() + .filter(|&&op| matches!(op & 0xf, 0 | 1 | 4)) + .map(|&op| (op >> 4) as usize) + .sum(); + let mut b = SamBuilder::new(); + b.read_name(name) + .ref_id(0) + .pos(pos) + .flags(flags::PAIRED | segment | strand | extra_flags) + .mate_ref_id(0) + .mate_pos(mate_pos) + .sequence(&vec![b'A'; len]) + .qualities(&vec![40; len]) + .cigar_ops(cigar); + if let Some(mc) = mc { + b.add_string_tag(SamTag::MC, mc); + } + b.build() + } + + /// `80M20S`, `20S80M` and `100M`, packed as `(len << 4) | op_code`. + const FWD: [u32; 2] = [80 << 4, (20 << 4) | 4]; + const REV: [u32; 2] = [(20 << 4) | 4, 80 << 4]; + const FULL: [u32; 1] = [100 << 4]; + + /// The groups the [`MateClipper`] clip tests run over. + #[derive(Clone, Copy, Debug)] + enum ClipperGroup { + /// A read-through pair with the given `MC` tags. + Pair(McTags), + /// A forward R1 lacking `MC`, then two R2s of the same name: the read-through mate + /// (`20S80M`) first and a `100M` one second. + DuplicateNameReadThroughMateFirst, + /// As [`ClipperGroup::DuplicateNameReadThroughMateFirst`], with the `100M` R2 first. + DuplicateNameFullMateFirst, + /// A read-through pair lacking `MC` whose reads carry no `SEQ`/`QUAL` (`*`). + EmptySeqWithoutMc, + /// An RF (outie) pair lacking `MC`: the forward read starts after its reverse mate ends. + RfPairWithoutMc, + /// A read-through pair lacking `MC`, plus a supplementary R1 of the same name that also + /// lacks `MC`. + SupplementaryWithoutMc, + } + + impl ClipperGroup { + fn records(self) -> Vec { + let r1 = || clipper_read(b"d", true, false, 1000, 1000, &FWD, None, 0); + match self { + Self::Pair(mc) => { + let (r1, r2) = read_through_pair(b"p", b"1", mc, true); + vec![r1, r2] + } + Self::DuplicateNameReadThroughMateFirst | Self::DuplicateNameFullMateFirst => { + let rt = clipper_read(b"d", false, true, 1000, 1000, &REV, Some(b"80M20S"), 0); + let fm = clipper_read(b"d", false, true, 1000, 1000, &FULL, Some(b"80M20S"), 0); + if matches!(self, Self::DuplicateNameReadThroughMateFirst) { + vec![r1(), rt, fm] + } else { + vec![r1(), fm, rt] + } + } + Self::EmptySeqWithoutMc => [(true, false, &FWD), (false, true, &REV)] + .into_iter() + .map(|(first, reverse, cigar)| { + let segment = + if first { flags::FIRST_SEGMENT } else { flags::LAST_SEGMENT }; + let strand = if reverse { flags::REVERSE } else { flags::MATE_REVERSE }; + let mut b = SamBuilder::new(); + b.read_name(b"e") + .ref_id(0) + .pos(1000) + .flags(flags::PAIRED | segment | strand) + .mate_ref_id(0) + .mate_pos(1000) + .cigar_ops(cigar); + b.build() + }) + .collect(), + Self::RfPairWithoutMc => vec![ + clipper_read(b"outie", true, false, 1200, 1000, &FULL, None, 0), + clipper_read(b"outie", false, true, 1000, 1200, &FULL, None, 0), + ], + Self::SupplementaryWithoutMc => { + let (r1, r2) = read_through_pair(b"p", b"1", McTags::Neither, true); + let supplementary = clipper_read( + b"p", + true, + false, + 1000, + 1000, + &FWD, + None, + flags::SUPPLEMENTARY, + ); + vec![r1, r2, supplementary] + } + } + } + } + + /// [`MateClipper::clip`] measures a read with `MC` against the tag, and a read without `MC` + /// against the CIGAR of its in-group mate. The read-through pair's reads each extend 20 bp + /// past their mate; a stale `100M` on the forward read hides its overhang. When a name has + /// more than two reads, the first mate in group order is used. The clip is measured on the + /// CIGARs alone, so a read with no `SEQ` is clipped like any other. An RF pair never + /// extends past its mate. A supplementary read is never backfilled, matching fgbio, which + /// drops it before backfilling. + #[rstest] + #[case::both_mc(ClipperGroup::Pair(McTags::Both), &[20, 20])] + #[case::neither_mc_backfilled(ClipperGroup::Pair(McTags::Neither), &[20, 20])] + #[case::forward_mc_backfilled(ClipperGroup::Pair(McTags::ForwardMissing), &[20, 20])] + #[case::stale_mc_used(ClipperGroup::Pair(McTags::StaleBoth), &[0, 20])] + #[case::duplicate_name_first_mate_used(ClipperGroup::DuplicateNameReadThroughMateFirst, &[20, 20, 0])] + #[case::duplicate_name_first_mate_used_swapped(ClipperGroup::DuplicateNameFullMateFirst, &[0, 0, 20])] + #[case::empty_seq_measured_on_cigars(ClipperGroup::EmptySeqWithoutMc, &[20, 20])] + #[case::rf_pair_not_clipped(ClipperGroup::RfPairWithoutMc, &[0, 0])] + #[case::supplementary_not_backfilled(ClipperGroup::SupplementaryWithoutMc, &[20, 20, 0])] + fn test_mate_clipper_clips_from_mc_or_in_group_mate( + #[case] group: ClipperGroup, + #[case] expected: &[usize], + ) { + let expected: Vec> = expected.iter().copied().map(Some).collect(); + assert_eq!(clips(&group.records()), expected); + } + + /// The in-group mate, at 1-based 1051, that the forward `80M20S` R1 at 1001 (lacking `MC`) + /// is clipped against. + #[derive(Clone, Copy, Debug)] + enum MateRecord { + /// A `30M` mate ending at 1080, where the R1's alignment ends: the usable control, past + /// which the R1's 20 bp soft clip extends. + Valid, + /// Unmapped with no CIGAR, while the R1's own `MATE_UNMAPPED` flag is (inconsistently) + /// clear. + Unmapped, + /// Unmapped but still carrying the `30M`, so only its flag rules it out. + UnmappedWithCigar, + /// Flagged mapped, with no CIGAR. + MappedWithoutCigar, + /// Fully soft-clipped (`100S`): no reference-consuming operation. + SoftClipOnly, + /// The `30M` mate truncated after its name, so its CIGAR is cut off. + Truncated, + } + + /// An in-group mate that is unmapped or whose CIGAR is not a valid CIGAR for a mapped record + /// is not clipped against, exactly as an `MC` tag of `*` or `100S` is not. Before that check + /// such a mate's empty reference span passed the FR gate and clipped the forward read at the + /// mate's start, taking 30 aligned bases with its soft clip (a clip of 50). + #[rstest] + #[case::valid_control(MateRecord::Valid, Some(20))] + #[case::unmapped(MateRecord::Unmapped, Some(0))] + #[case::unmapped_with_cigar(MateRecord::UnmappedWithCigar, Some(0))] + #[case::mapped_without_cigar(MateRecord::MappedWithoutCigar, Some(0))] + #[case::soft_clip_only(MateRecord::SoftClipOnly, Some(0))] + #[case::truncated(MateRecord::Truncated, Some(0))] + fn test_mate_clipper_does_not_clip_against_an_unusable_mate( + #[case] mate: MateRecord, + #[case] expected: Option, + ) { + let r1 = clipper_read(b"d", true, false, 1000, 1050, &FWD, None, 0); + let thirty_m = [encode_op(0, 30)]; + let (cigar, extra_flags): (&[u32], u16) = match mate { + MateRecord::Valid | MateRecord::Truncated => (&thirty_m, 0), + MateRecord::Unmapped => (&[], flags::UNMAPPED), + MateRecord::UnmappedWithCigar => (&thirty_m, flags::UNMAPPED), + MateRecord::MappedWithoutCigar => (&[], 0), + MateRecord::SoftClipOnly => (&[encode_op(4, 100)], 0), + }; + let r2 = clipper_read(b"d", false, true, 1050, 1000, cigar, None, extra_flags); + // Keep the 32 fixed bytes and the name (`d\0`), but cut off the CIGAR and the rest. + let r2_bytes = + if matches!(mate, MateRecord::Truncated) { &r2.as_ref()[..34] } else { r2.as_ref() }; + let mut clipper = MateClipper::new([r1.as_ref(), r2_bytes].into_iter()); + assert_eq!(clipper.clip(r1.as_ref()), expected); + } + + /// `clip` is `None` — the mate is required — for a primary read lacking `MC` whose mate is + /// mapped to its reference on the opposite strand but is absent from the group. Every other + /// case differs from `needs_mate` in exactly one guard, so it is that guard which makes the + /// absent mate unnecessary; `clip` then returns 0 without consulting the group. + #[rstest] + #[case::needs_mate(flags::PAIRED | flags::MATE_REVERSE, 0, None, None)] + #[case::read_carries_mc(flags::PAIRED | flags::MATE_REVERSE, 0, Some(&b"100M"[..]), Some(0))] + #[case::mate_unmapped(flags::PAIRED | flags::MATE_REVERSE | flags::MATE_UNMAPPED, 0, None, Some(0))] + #[case::read_unmapped(flags::PAIRED | flags::MATE_REVERSE | flags::UNMAPPED, 0, None, Some(0))] + #[case::mate_on_other_reference(flags::PAIRED | flags::MATE_REVERSE, 1, None, Some(0))] + #[case::mate_on_same_strand(flags::PAIRED, 0, None, Some(0))] + #[case::unpaired(flags::MATE_REVERSE, 0, None, Some(0))] + #[case::secondary(flags::PAIRED | flags::MATE_REVERSE | flags::SECONDARY, 0, None, Some(0))] + #[case::supplementary(flags::PAIRED | flags::MATE_REVERSE | flags::SUPPLEMENTARY, 0, None, Some(0))] + fn test_mate_clipper_requires_the_mate_only_when_its_cigar_is_needed( + #[case] flag: u16, + #[case] mate_ref_id: i32, + #[case] mc: Option<&[u8]>, + #[case] expected: Option, + ) { + let mut b = SamBuilder::new(); + b.read_name(b"lonely") + .ref_id(0) + .pos(1000) + .flags(flag | flags::FIRST_SEGMENT) + .mate_ref_id(mate_ref_id) + .mate_pos(1100) + .sequence(&[b'A'; 100]) + .qualities(&[40; 100]) + .cigar_ops(&FULL); + if let Some(mc) = mc { + b.add_string_tag(SamTag::MC, mc); + } + let read = b.build(); + assert_eq!(clips(std::slice::from_ref(&read)), vec![expected]); + } + + /// The mate must be a primary read of the group: an FR or RF pair whose mate is absent, or + /// present only as a secondary alignment, cannot be backfilled. An RF pair is included, + /// which fgbio would not backfill or fail on: without the mate CIGAR, orientation can only + /// be judged from `TLEN`, which misclassifies dovetailed FR pairs. + #[rstest] + #[case::fr_mate_absent(false, 1000, None)] + #[case::rf_mate_absent(false, 1200, None)] + #[case::mate_only_secondary(true, 1000, Some(flags::SECONDARY))] + fn test_mate_clipper_needs_a_primary_mate_in_the_group( + #[case] include_mate: bool, + #[case] r1_pos: i32, + #[case] mate_extra_flags: Option, + ) { + // R1 forward at `r1_pos`, its reverse mate at 1100: FR when R1 starts first, RF after. + let r1 = clipper_read(b"needy", true, false, r1_pos, 1100, &FULL, None, 0); + let mut records = vec![r1]; + if include_mate { + let flags = mate_extra_flags.unwrap_or(0); + records.push(clipper_read(b"needy", false, true, 1100, r1_pos, &FULL, None, flags)); + } + assert_eq!(clips(&records)[0], None); + } + + /// The missing-mate error names the read and says how to fix the input. + #[test] + fn test_missing_mate_error_names_the_read_and_the_remedy() { + let (r1, _) = read_through_pair(b"orphan_read", b"1", McTags::Neither, true); + assert_eq!( + format!("{:#}", missing_mate_error(r1.as_ref())), + "Mate cigar (MC SAM tag) needed for read 'orphan_read': the read has no MC tag and \ + its primary mate is not in the same group. Add MC tags (e.g. with fgumi zipper or \ + samtools fixmate), or keep both reads of each template in the same group." + ); + } +} diff --git a/crates/fgumi-consensus/src/vanilla_caller.rs b/crates/fgumi-consensus/src/vanilla_caller.rs index dcbf55263..161abece7 100644 --- a/crates/fgumi-consensus/src/vanilla_caller.rs +++ b/crates/fgumi-consensus/src/vanilla_caller.rs @@ -9,6 +9,7 @@ use crate::caller::{ ConsensusCaller, ConsensusCallingStats, ConsensusOutput, RejectionReason, consensus_read_name_too_long_context, select_lowest_ranking, write_consensus_read_name, }; +use crate::mate_clip::{MateClipper, missing_mate_error}; use crate::phred::{ MIN_PHRED, NO_CALL_BASE, NO_CALL_BASE_LOWER, PhredScore, ln_error_prob_two_trials, ln_prob_to_phred, phred_to_ln_error_prob, @@ -1390,14 +1391,21 @@ impl VanillaUmiConsensusCaller { // Sub-group by read type let (fragment_reads, r1_reads, r2_reads) = self.subgroup_reads(reads); + // Read-through clipping looks up the mate of a read lacking `MC` across both ends of the + // pair: its mate lands in the other subgroup. Fragments are unpaired, so they never need + // one. + let mut clipper = + MateClipper::new(r1_reads.iter().chain(&r2_reads).map(|(_, raw)| raw.as_ref())); + let mut output = ConsensusOutput::default(); // Process fragment subgroup - let (fragment_ok, _, _) = self.process_subgroup( + let (fragment_ok, _, _, _) = self.process_subgroup( &mut output, umi, ReadType::Fragment, - fragment_reads, + &fragment_reads, + &mut clipper, &mut group_rejects, )?; if fragment_ok { @@ -1406,20 +1414,36 @@ impl VanillaUmiConsensusCaller { // Process R1/R2 subgroups let mut r1r2_output = ConsensusOutput::default(); - let (r1_ok, r1_surviving_count, r1_surviving_reads) = self.process_subgroup( - &mut r1r2_output, - umi, - ReadType::R1, - r1_reads, - &mut group_rejects, - )?; - let (r2_ok, r2_surviving_count, r2_surviving_reads) = self.process_subgroup( - &mut r1r2_output, - umi, - ReadType::R2, - r2_reads, - &mut group_rejects, - )?; + let (r1_ok, r1_surviving_count, r1_surviving_reads, r1_missing_mate) = self + .process_subgroup( + &mut r1r2_output, + umi, + ReadType::R1, + &r1_reads, + &mut clipper, + &mut group_rejects, + )?; + let (r2_ok, r2_surviving_count, r2_surviving_reads, r2_missing_mate) = self + .process_subgroup( + &mut r1r2_output, + umi, + ReadType::R2, + &r2_reads, + &mut clipper, + &mut group_rejects, + )?; + + // A read whose read-through clip needs a mate missing from the group fails the run only + // when it reaches an emitted consensus; one the group rejects anyway (too few reads, + // a minority alignment, downsampled, or an orphan end) does not. + if r1_ok && r2_ok { + let missing = r1_missing_mate + .map(|idx| &r1_reads[idx].1) + .or_else(|| r2_missing_mate.map(|idx| &r2_reads[idx].1)); + if let Some(raw) = missing { + return Err(missing_mate_error(raw.as_ref())); + } + } match (r1_ok, r2_ok) { (true, true) => { @@ -1466,6 +1490,11 @@ impl VanillaUmiConsensusCaller { /// - `bool` — whether a consensus was produced /// - `usize` — count of reads that survived all internal filtering (not rejected) /// - `Vec>` — the surviving raw reads as bytes (only populated when `self.track_rejects`) + /// - `Option` — when a consensus was produced, the index in `group_reads` of the first + /// read in it whose read-through clip needs a mate that is not in the group (see + /// [`MateClipper::clip`]). Such a read is clipped by 0 until then, so that whether it + /// reaches the consensus is decided by the usual filters; the caller fails the run only if + /// the group's consensus is emitted. /// /// The return type's third element is `Vec<(usize, Vec)>` (not /// `Vec`) because orphan-consensus rejections are forwarded to the @@ -1475,43 +1504,42 @@ impl VanillaUmiConsensusCaller { /// restore input order — `group_rejects` collects every non-orphan rejection /// with that same position. #[allow(clippy::too_many_lines, clippy::type_complexity)] - fn process_subgroup( + fn process_subgroup<'a, I>( &mut self, output: &mut ConsensusOutput, umi: &str, read_type: ReadType, - group_reads: Vec<(usize, RawRecord)>, + group_reads: &[(usize, RawRecord)], + clipper: &mut MateClipper<'a, I>, group_rejects: &mut Vec<(usize, Vec)>, - ) -> Result<(bool, usize, Vec<(usize, Vec)>)> { - use fgumi_raw_bam as bam_fields; - + ) -> Result<(bool, usize, Vec<(usize, Vec)>, Option)> + where + I: Iterator + Clone, + { if group_reads.is_empty() { - return Ok((false, 0, Vec::new())); + return Ok((false, 0, Vec::new(), None)); } if group_reads.len() < self.options.min_reads { self.stats.record_rejection(RejectionReason::InsufficientReads, group_reads.len()); if self.track_rejects { - group_rejects - .extend(group_reads.into_iter().map(|(pos, raw)| (pos, raw.into_inner()))); + group_rejects.extend(group_reads.iter().map(|(pos, raw)| (*pos, raw.to_vec()))); } - return Ok((false, 0, Vec::new())); + return Ok((false, 0, Vec::new(), None)); } - // Calculate mate overlap clips from raw bytes - let mate_overlap_clips: Vec = group_reads - .iter() - .map(|(_, raw)| bam_fields::num_bases_extending_past_mate_raw(raw.as_ref())) - .collect(); - // Create SourceReads from raw bytes. `idx` indexes `group_reads` (this // subgroup); `group_reads[idx].0` recovers the group-input position. let mut source_reads: Vec = Vec::new(); let mut zero_length_indices: Vec = Vec::new(); - - for (idx, ((_, raw), &mate_clip)) in - group_reads.iter().zip(mate_overlap_clips.iter()).enumerate() - { + // Reads whose mate CIGAR is needed but unavailable; see the return value. + let mut missing_mate_indices: HashSet = HashSet::new(); + + for (idx, (_, raw)) in group_reads.iter().enumerate() { + let mate_clip = clipper.clip(raw.as_ref()).unwrap_or_else(|| { + missing_mate_indices.insert(idx); + 0 + }); if let Some(sr) = self.create_source_read(raw.as_ref(), idx, mate_clip)? { source_reads.push(sr); } else { @@ -1542,7 +1570,7 @@ impl VanillaUmiConsensusCaller { } } } - return Ok((false, 0, Vec::new())); + return Ok((false, 0, Vec::new(), None)); } // Filter by alignment @@ -1569,7 +1597,7 @@ impl VanillaUmiConsensusCaller { } } } - return Ok((false, 0, Vec::new())); + return Ok((false, 0, Vec::new(), None)); } // Apply the per-end `--max-reads` cap (fgumi#723). fgbio applies `--max-reads` @@ -1614,13 +1642,21 @@ impl VanillaUmiConsensusCaller { } } } - return Ok((false, 0, Vec::new())); + return Ok((false, 0, Vec::new(), None)); } // Capture surviving count and reads before building consensus. The // orphan-consensus path rejects survivors, so they too carry their // group-input position. let surviving_count = filtered_source_reads.len(); + let missing_mate = if missing_mate_indices.is_empty() { + None + } else { + filtered_source_reads + .iter() + .map(|sr| sr.original_idx) + .find(|idx| missing_mate_indices.contains(idx)) + }; let surviving_reads = if self.track_rejects { filtered_source_reads .iter() @@ -1665,7 +1701,7 @@ impl VanillaUmiConsensusCaller { methylation.as_ref(), )?; - Ok((true, surviving_count, surviving_reads)) + Ok((true, surviving_count, surviving_reads, missing_mate)) } /// Creates consensus from `SourceReads` (which already have transformed bases/quals). @@ -1965,6 +2001,9 @@ pub(crate) enum ReadType { )] mod tests { use super::*; + use crate::mate_clip::read_through_fixture::{ + McTags, READ_THROUGH_ADAPTER, READ_THROUGH_INSERT, read_through_pair, + }; use fgumi_raw_bam::{ ParsedBamRecord, SamBuilder, encode_op, num_bases_extending_past_mate_raw, }; @@ -5555,7 +5594,9 @@ mod tests { // ========================================================================= /// Port of fgbio test: "add the mate cigar when not present before consensus calling" - /// Tests that consensus calling works even when MC tag is missing + /// (`VanillaUmiConsensusCallerTest.scala:818`). Tests that consensus calling works even when + /// the MC tag is missing; the effect of the backfill on the read-through trim is asserted by + /// [`test_mate_cigar_backfilled_before_consensus_trims_read_through`]. #[test] fn test_mate_cigar_handling() { let len = 10; @@ -5849,6 +5890,181 @@ mod tests { ); } + /// Port of fgbio test: "add the mate cigar when not present before consensus calling" + /// (`VanillaUmiConsensusCallerTest.scala:818`), extended to assert the effect of the backfill. + /// + /// fgbio's `UmiConsensusCaller.updateMateCigars` (`UmiConsensusCaller.scala:328-343`, called + /// at `:356`) copies each mate's CIGAR into the other read's `MC` before any read is + /// converted to a `SourceRead`, so the read-through trim in `toSourceRead` always sees the + /// mate CIGAR. The scala test only counts the consensus reads; here the pair reads through + /// the 80 bp insert into 20 bp of adapter, so a backfilled mate CIGAR must trim both + /// consensus reads to the 80 bp insert. A pair that already carries `MC` is used as-is, even + /// when it is stale (fgumi deliberately does not mirror fgbio's overwrite of an existing + /// `MC`). + #[rstest] + #[case::both_mc_present(McTags::Both, false)] + #[case::neither_mc_backfilled(McTags::Neither, false)] + #[case::forward_r1_missing_mc_backfilled(McTags::ForwardMissing, false)] + #[case::stale_mc_is_used_not_overwritten(McTags::StaleBoth, true)] + fn test_mate_cigar_backfilled_before_consensus_trims_read_through( + #[case] mc: McTags, + #[case] r1_keeps_adapter: bool, + ) { + let mut records = Vec::new(); + for i in 0..3 { + let name = format!("t{i}"); + let (r1, r2) = read_through_pair(name.as_bytes(), b"1", mc, true); + records.push(r1); + records.push(r2); + } + + let options = VanillaUmiConsensusOptions { min_reads: 1, ..Default::default() }; + let mut caller = VanillaUmiConsensusCaller::new("c".to_string(), "1".to_string(), options); + let output = consensus_reads_from_raw(&mut caller, records) + .expect("read-through pairs must produce consensus reads"); + + assert_eq!(output.count, 2, "one consensus read per end"); + let consensus = ParsedBamRecord::parse_all(&output.data); + let r1 = consensus.iter().find(|r| r.flag & flags::FIRST_SEGMENT != 0).expect("R1"); + let r2 = consensus.iter().find(|r| r.flag & flags::LAST_SEGMENT != 0).expect("R2"); + + // R1 is the forward read: the adapter is its 3' end. R2 is emitted in sequencing + // orientation, so the insert appears reverse-complemented. With the stale `100M` mate + // CIGAR, R1 sees a mate that ends 20 bp past the insert and so keeps its adapter (as + // fgbio does), while R2's clip depends only on the mate start and is unaffected. + let mut expected_r1 = READ_THROUGH_INSERT.to_vec(); + if r1_keeps_adapter { + expected_r1.extend_from_slice(READ_THROUGH_ADAPTER); + } + assert_eq!( + String::from_utf8(r1.bases.clone()).expect("ASCII bases"), + String::from_utf8(expected_r1).expect("ASCII bases"), + "R1 consensus bases" + ); + assert_eq!( + String::from_utf8(r2.bases.clone()).expect("ASCII bases"), + String::from_utf8(reverse_complement(READ_THROUGH_INSERT)).expect("ASCII bases"), + "R2 consensus bases" + ); + assert_eq!(r1.get_int_tag(SamTag::CD), Some(3)); + assert_eq!(r2.get_int_tag(SamTag::CD), Some(3)); + } + + /// Calls simplex consensus (`--min-reads min_reads`) on `records`, returning the outcome and + /// the rejection counts, sorted by reason name so comparisons do not depend on hash-map + /// iteration order. + fn call_simplex( + records: Vec, + min_reads: usize, + ) -> (anyhow::Result, Vec<(RejectionReason, usize)>) { + let options = VanillaUmiConsensusOptions { min_reads, ..Default::default() }; + let mut caller = VanillaUmiConsensusCaller::new("c".to_string(), "1".to_string(), options); + let output = consensus_reads_from_raw(&mut caller, records); + let mut rejections: Vec<(RejectionReason, usize)> = + caller.stats.rejection_reasons.iter().map(|(reason, n)| (*reason, *n)).collect(); + rejections.sort_by_key(|&(reason, _)| format!("{reason:?}")); + (output, rejections) + } + + /// `count` read-through pairs `t0`, `t1`, ... with the given `MC` tags, R1 then R2. + fn read_through_family(count: usize, mc: McTags) -> Vec { + (0..count) + .flat_map(|i| { + let (r1, r2) = read_through_pair(format!("t{i}").as_bytes(), b"1", mc, true); + [r1, r2] + }) + .collect() + } + + /// Port of the fgbio behavior at `UmiConsensusCaller.scala:328-343` (`updateMateCigars`) + /// exercised by `VanillaUmiConsensusCallerTest.scala:818` ("add the mate cigar when not + /// present before consensus calling"): a read with no `MC` whose mate is not in the group + /// cannot have its mate CIGAR backfilled. fgbio throws (a `MatchError`); fgumi must stop + /// with an error that names the read rather than silently skipping the read-through trim, + /// here because the read reaches the emitted R1 consensus. + #[test] + fn test_missing_mate_without_mc_errors_naming_the_read() { + let mut records = read_through_family(2, McTags::Neither); + let (orphan, _) = read_through_pair(b"orphan_read", b"1", McTags::Neither, true); + records.push(orphan); + + let err = call_simplex(records, 1) + .0 + .err() + .expect("a read lacking MC whose mate is absent must be an error"); + assert_eq!( + format!("{err:#}"), + "Mate cigar (MC SAM tag) needed for read 'orphan_read': the read has no MC tag and \ + its primary mate is not in the same group. Add MC tags (e.g. with fgumi zipper or \ + samtools fixmate), or keep both reads of each template in the same group." + ); + } + + /// A read lacking `MC` whose mate is missing fails the run only if it would reach an emitted + /// consensus. Each group here is rejected or filtered by the usual checks, and must be handled + /// exactly as when every read carries `MC`: + /// - `end_below_min_reads`: with `--min-reads 2`, one pair plus a lone R1 give an R2 end of + /// one read, so no consensus pair can be built (the R1 consensus is an orphan). fgbio, which + /// backfills before any filtering, would fail here. + /// - `minority_alignment`: the lone R1 carries an insertion its family lacks, so it is + /// dropped as a minority alignment and the family's consensus pair is still emitted. + #[rstest] + #[case::end_below_min_reads( + 1, + 2, + false, + 0, + &[(RejectionReason::InsufficientReads, 1), (RejectionReason::OrphanConsensus, 2)] + )] + #[case::minority_alignment(2, 1, true, 2, &[(RejectionReason::MinorityAlignment, 1)])] + fn test_missing_mate_is_not_an_error_for_a_read_that_reaches_no_consensus( + #[case] pairs: usize, + #[case] min_reads: usize, + #[case] lone_read_is_minority_alignment: bool, + #[case] expected_count: usize, + #[case] expected_rejections: &[(RejectionReason, usize)], + ) { + // The lone R1 carries `MC` exactly when the rest of the group does. + let group = |mc: McTags| { + let mut records = read_through_family(pairs, mc); + let (mut lone, _) = read_through_pair(b"lone", b"1", mc, true); + if lone_read_is_minority_alignment { + let mut b = SamBuilder::new(); + b.read_name(b"lone") + .ref_id(0) + .pos(1000) + .flags(flags::PAIRED | flags::FIRST_SEGMENT | flags::MATE_REVERSE) + .mate_ref_id(0) + .mate_pos(1000) + .sequence(&[READ_THROUGH_INSERT, READ_THROUGH_ADAPTER].concat()) + .qualities(&[40; 100]) + .cigar_ops(&[ + encode_op(0, 40), + encode_op(1, 2), + encode_op(0, 38), + encode_op(4, 20), + ]) + .add_string_tag(SamTag::MI, b"1"); + if matches!(mc, McTags::Both) { + b.add_string_tag(SamTag::MC, b"20S80M"); + } + lone = b.build(); + } + records.push(lone); + records + }; + + let (with_mc, with_mc_rejections) = call_simplex(group(McTags::Both), min_reads); + let (without_mc, without_mc_rejections) = call_simplex(group(McTags::Neither), min_reads); + + let with_mc = with_mc.expect("the group with MC is valid"); + let without_mc = without_mc.expect("a read that reaches no consensus needs no mate"); + assert_eq!(without_mc.count, expected_count); + assert_eq!(without_mc.data, with_mc.data, "output must match the input with MC"); + assert_eq!(without_mc_rejections, expected_rejections); + assert_eq!(without_mc_rejections, with_mc_rejections); + } + // ========================================================================= // Test 42: Consensus UMI (RX) on filtered reads only // ========================================================================= diff --git a/crates/fgumi-raw-bam/src/lib.rs b/crates/fgumi-raw-bam/src/lib.rs index a376facdb..7676c9650 100644 --- a/crates/fgumi-raw-bam/src/lib.rs +++ b/crates/fgumi-raw-bam/src/lib.rs @@ -98,8 +98,9 @@ pub use cigar::{ // -- overlap -- pub use overlap::{ - bases_extending_past_mate_ops, is_primary_fr_pair_raw, num_bases_extending_past_mate_raw, - num_bases_extending_past_mate_vs_mate_raw, + bases_extending_past_mate_ops, is_mapped_opposite_strand_pair_raw, is_primary_fr_pair_raw, + num_bases_extending_past_mate_cigar_raw, num_bases_extending_past_mate_from_mc_raw, + num_bases_extending_past_mate_raw, num_bases_extending_past_mate_vs_mate_raw, }; // -- raw_bam_record -- diff --git a/crates/fgumi-raw-bam/src/overlap.rs b/crates/fgumi-raw-bam/src/overlap.rs index 3e9504d4e..dbe0acda8 100644 --- a/crates/fgumi-raw-bam/src/overlap.rs +++ b/crates/fgumi-raw-bam/src/overlap.rs @@ -5,6 +5,24 @@ use crate::cigar::{ use crate::fields::{RawRecordView, aux_data_slice, flags, mate_pos, mate_ref_id, template_length}; use crate::tags::RawTagsView; +/// Whether a read and its mate, judged from the read's own fields, are both mapped to the same +/// reference on opposite strands: the read is paired, neither it nor its mate is unmapped, its +/// mate reference is its own, and its strand differs from its mate's. +/// +/// These are the preconditions for the read to be one end of an FR or RF pair, before any +/// orientation test. They are shared by every per-record FR test here and by the consensus +/// callers' check for whether a read lacking `MC` needs its mate's CIGAR, so the guards cannot +/// drift apart. +#[must_use] +pub fn is_mapped_opposite_strand_pair_raw(bam: &[u8]) -> bool { + let v = RawRecordView::new(bam); + let flg = v.flags(); + flg & flags::PAIRED != 0 + && flg & (flags::UNMAPPED | flags::MATE_UNMAPPED) == 0 + && v.ref_id() == mate_ref_id(bam) + && (flg & flags::REVERSE == 0) != (flg & flags::MATE_REVERSE == 0) +} + /// Crate-internal per-record FR (forward-reverse) classification from raw BAM bytes. /// /// **Per-record, and not a gate for pair decisions.** The forward-strand arm derives the mate's @@ -25,32 +43,12 @@ use crate::tags::RawTagsView; /// a dovetail whose reverse read's aligned end equals the forward read's aligned start — is FR. #[must_use] pub(crate) fn is_fr_pair_raw(bam: &[u8]) -> bool { - let v = RawRecordView::new(bam); - let flg = v.flags(); - - // Must be paired - if flg & flags::PAIRED == 0 { - return false; - } - - // Both read and mate must be mapped - if flg & flags::UNMAPPED != 0 || flg & flags::MATE_UNMAPPED != 0 { - return false; - } - - // Must be on the same reference - let this_ref_id = v.ref_id(); - let m_ref_id = mate_ref_id(bam); - if this_ref_id != m_ref_id { - return false; - } - - // Must be on opposite strands for FR or RF - let is_reverse = flg & flags::REVERSE != 0; - let mate_is_reverse = flg & flags::MATE_REVERSE != 0; - if is_reverse == mate_is_reverse { + // Paired, both mapped, same reference, opposite strands (FR or RF). + if !is_mapped_opposite_strand_pair_raw(bam) { return false; } + let v = RawRecordView::new(bam); + let is_reverse = v.flags() & flags::REVERSE != 0; // Determine if FR or RF using htsjdk 5.0.0's logic (TLEN fallback for the forward arm): // positiveStrandFivePrimePos = readIsOnReverseStrand ? mateStart : alignmentStart @@ -150,27 +148,14 @@ pub fn is_primary_fr_pair_raw(a: &[u8], b: &[u8]) -> bool { /// pairs without needing the mate record in hand. #[must_use] fn is_fr_pair_with_mate_cigar_raw(bam: &[u8], mate_ref_len: i32) -> bool { - let v = RawRecordView::new(bam); - let flg = v.flags(); - - // Paired, both read and mate mapped, same reference, opposite strands — identical to the - // non-orientation guards in `is_fr_pair_raw`. - if flg & flags::PAIRED == 0 { - return false; - } - if flg & flags::UNMAPPED != 0 || flg & flags::MATE_UNMAPPED != 0 { - return false; - } - if v.ref_id() != mate_ref_id(bam) { - return false; - } - let is_reverse = flg & flags::REVERSE != 0; - let mate_is_reverse = flg & flags::MATE_REVERSE != 0; - if is_reverse == mate_is_reverse { + // Paired, both read and mate mapped, same reference, opposite strands — the same guards as + // `is_fr_pair_raw`. + if !is_mapped_opposite_strand_pair_raw(bam) { return false; } + let v = RawRecordView::new(bam); - if is_reverse { + if v.flags() & flags::REVERSE != 0 { // This read is the reverse record; its own arm is CIGAR-derived and TLEN-independent. return is_fr_pair_raw(bam); } @@ -200,36 +185,76 @@ fn is_fr_pair_with_mate_cigar_raw(bam: &[u8], mate_ref_len: i32) -> bool { /// [`crate::cigar::mate_unclipped_5prime`]), which follows samtools' soft+hard convention. /// /// Used by the simplex (vanilla) and duplex consensus callers, in the place of fgbio's base -/// `UmiConsensusCaller.numBasesExtendingPastMate`. The codec caller instead uses +/// `UmiConsensusCaller.numBasesExtendingPastMate`. For a read with no MC tag those callers take +/// the mate CIGAR from the mate record in the group instead, through +/// [`num_bases_extending_past_mate_cigar_raw`]. The codec caller instead uses /// [`num_bases_extending_past_mate_vs_mate_raw`], which reads the mate's alignment from the mate /// record in hand rather than its MC tag (see CODEC3-04). #[must_use] pub fn num_bases_extending_past_mate_raw(bam: &[u8]) -> usize { - // Need the MC tag for mate CIGAR information. Parsed before the FR gate because the gate - // now classifies orientation from the mate CIGAR (see below); with no usable MC there is - // nothing to compute the overhang from anyway, so absent/invalid MC fails closed to 0. + num_bases_extending_past_mate_from_mc_raw(bam).unwrap_or(0) +} + +/// As [`num_bases_extending_past_mate_raw`], but distinguishes a read with no MC tag: returns +/// `None` when the tag is absent, and otherwise the clip, which is `0` for non-FR pairs or an +/// invalid MC tag. +/// +/// Lets a caller that falls back to another source for the mate CIGAR (the simplex and duplex +/// callers' backfill from the in-group mate) look the tag up once per read. +#[must_use] +pub fn num_bases_extending_past_mate_from_mc_raw(bam: &[u8]) -> Option { + // Parsed before the FR gate because the gate classifies orientation from the mate CIGAR (see + // `num_bases_extending_past_mate_cigar_raw`); with no usable MC there is nothing to compute + // the overhang from anyway, so an invalid MC fails closed to 0. let aux = aux_data_slice(bam); - let Some(mc_bytes) = RawTagsView::new(aux).find_string(SamTag::MC) else { - return 0; - }; + let mc_bytes = RawTagsView::new(aux).find_string(SamTag::MC)?; let Ok(mc_cigar) = std::str::from_utf8(mc_bytes) else { - return 0; + return Some(0); }; - // A malformed MC CIGAR fails closed to 0 (no clip), honoring the doc contract above - // ("Returns 0 ... when the MC tag is absent/invalid"). + // A malformed MC CIGAR fails closed to 0 (no clip), honoring the doc contract of + // `num_bases_extending_past_mate_raw` ("Returns 0 ... when the MC tag is absent/invalid"). let Some(mate_ops) = parse_mc_cigar_ops(mc_cigar) else { - return 0; + return Some(0); }; + Some(clip_against_valid_mate_cigar(bam, &mate_ops)) +} +/// Number of bases a read extends past its mate for FR pairs, given the mate's CIGAR ops +/// (BAM-encoded, as from [`crate::cigar::get_cigar_ops`]) in place of the read's MC tag. Returns +/// `0` for non-FR pairs, and for mate ops that are not a structurally valid CIGAR for a mapped +/// record — the same validation [`num_bases_extending_past_mate_raw`] applies to an MC tag, so +/// an empty (`*`), fully clipped or otherwise degenerate mate CIGAR is never clipped against. +/// +/// This is the core of [`num_bases_extending_past_mate_raw`] once the MC tag has been parsed: +/// the mate's position still comes from the read's own mate-position field, exactly as when the +/// MC tag supplies the CIGAR. The simplex and duplex callers use it to backfill a missing MC tag +/// from the mate record in the same group, as fgbio's `UmiConsensusCaller.updateMateCigars` does. +/// +/// [`num_bases_extending_past_mate_vs_mate_raw`] also measures against an in-group mate, but it +/// is not interchangeable: it takes the mate's position from the mate record rather than the +/// read's mate-position field, and gates on both records' flags rather than the read's own. A +/// read lacking MC must be clipped exactly as it would be with a correct MC tag, so this path +/// keeps the MC path's geometry and gate and swaps in only the CIGAR. +#[must_use] +pub fn num_bases_extending_past_mate_cigar_raw(bam: &[u8], mate_ops: &[u32]) -> usize { + if !is_valid_mapped_cigar(mate_ops) { + return 0; + } + clip_against_valid_mate_cigar(bam, mate_ops) +} + +/// Shared body of the MC-tag and mate-CIGAR entry points, for mate ops already validated by +/// [`is_valid_mapped_cigar`]. +fn clip_against_valid_mate_cigar(bam: &[u8], mate_ops: &[u32]) -> usize { // Only applies to FR pairs. Classify per-pair from the read plus its mate CIGAR rather than // the per-record `is_fr_pair_raw`, whose forward-strand TLEN arm misclassifies dovetail FR // pairs (#839 / htsjdk/samtools#1771) and would zero the read-through clip for the forward // read, leaking adapter into the simplex/duplex consensus. - if !is_fr_pair_with_mate_cigar_raw(bam, saturating_reference_length(&mate_ops)) { + if !is_fr_pair_with_mate_cigar_raw(bam, saturating_reference_length(mate_ops)) { return 0; } - bases_extending_past_mate(bam, mate_pos(bam) + 1, &mate_ops) + bases_extending_past_mate(bam, mate_pos(bam) + 1, mate_ops) } /// Number of bases a read extends past its mate for FR pairs, taking the mate's alignment from @@ -450,19 +475,20 @@ const MAX_CIGAR_OP_LEN: u32 = (1 << 28) - 1; /// mate-record paths share one set of CIGAR walkers. /// /// Returns `None` for any CIGAR that is not a structurally valid SAM CIGAR for a mapped record, -/// so malformed MC-tag input fails closed rather than yielding a partial answer. Rejected: +/// so malformed MC-tag input fails closed rather than yielding a partial answer. Lexically +/// rejected here: /// - an unknown operator byte (e.g. the `f`/`o` in `10Mfoo5S`, or a non-ASCII byte); /// - an operator with no preceding run-length (a bare `MS`) or a zero run-length (`0M`, `10M0S`); /// - a run-length past [`MAX_CIGAR_OP_LEN`], which no BAM CIGAR operation can hold; -/// - a trailing run-length with no operator (`10M5`), or an empty CIGAR; -/// - a soft (`S`) or hard (`H`) clip that is not at the end of the CIGAR — `S` may only sit at the -/// ends (inside any `H`), and `H` only as the first/last operation (e.g. `10M5S10M` is invalid); -/// - a CIGAR with no reference-consuming operation (`10S`), which a mapped mate cannot have. +/// - a trailing run-length with no operator (`10M5`), or an empty CIGAR. +/// +/// The parsed ops are then held to [`is_valid_mapped_cigar`], the structural check shared with a +/// mate CIGAR taken from the mate record (misplaced clips, or no reference-consuming operation). fn parse_mc_cigar_ops(cigar: &str) -> Option> { - // Phase 1: tokenize into (run_length, operator) pairs. Accumulate the run-length from digit + // Tokenize and pack into (run_length << 4) | op_code. Accumulate the run-length from digit // bytes directly rather than slicing the string by byte offset — slicing on a non-ASCII byte - // (a malformed CIGAR) would panic on a char boundary. Reject lexical errors here. - let mut tokens: Vec<(u32, u8)> = Vec::new(); + // (a malformed CIGAR) would panic on a char boundary. + let mut ops: Vec = Vec::new(); let mut num = 0u32; let mut have_digits = false; for &c in cigar.as_bytes() { @@ -477,50 +503,59 @@ fn parse_mc_cigar_ops(cigar: &str) -> Option> { continue; } // Every operator needs a positive run-length, and must be a known CIGAR operator. - if !have_digits || num == 0 || cigar_op_code(c).is_none() { + if !have_digits || num == 0 { return None; } - tokens.push((num, c)); + // `num <= MAX_CIGAR_OP_LEN` (28 bits), so the shift cannot drop bits. + ops.push((num << 4) | cigar_op_code(c)?); num = 0; have_digits = false; } - // A dangling run-length with no operator (`10M5`), or an empty CIGAR, is malformed. - if have_digits || tokens.is_empty() { + // A dangling run-length with no operator (`10M5`) is malformed; an empty CIGAR is rejected + // by the structural check. + if have_digits { return None; } + is_valid_mapped_cigar(&ops).then_some(ops) +} - // Phase 2: validate structure and pack the ops in one pass. Soft clips must sit at the ends - // (inside any hard clip); hard clips only as the first/last operation. - let last = tokens.len() - 1; - let mut ops = Vec::with_capacity(tokens.len()); +/// Whether BAM-packed CIGAR ops (`(len << 4) | op_code`) are a structurally valid CIGAR for a +/// mapped record. The single structural check for a mate CIGAR, whether it came from an MC tag +/// ([`parse_mc_cigar_ops`]) or from the mate record itself +/// ([`num_bases_extending_past_mate_cigar_raw`]), so the two paths cannot drift. Rejected: +/// - no operations at all (an unmapped record's `*` CIGAR, or a truncated record); +/// - a zero-length operation, or an operation code that is not a CIGAR operator; +/// - a soft (`S`) or hard (`H`) clip that is not at the end of the CIGAR — `S` may only sit at the +/// ends (inside any `H`), and `H` only as the first/last operation (e.g. `10M5S10M` is invalid); +/// - no reference-consuming operation (`10S`), which a mapped record cannot have. +fn is_valid_mapped_cigar(ops: &[u32]) -> bool { + let Some(last) = ops.len().checked_sub(1) else { + return false; + }; + let is_hard = |op: &u32| op & 0xF == 5; let mut saw_ref_op = false; - for (i, &(len, op)) in tokens.iter().enumerate() { - // Phase 1 accepted every operator, so this lookup cannot fail. - let op_code = cigar_op_code(op)?; - match op { - b'M' | b'D' | b'N' | b'=' | b'X' => saw_ref_op = true, - b'I' | b'P' => {} // interior op: no placement constraint - b'S' => { - let leading = tokens[..i].iter().all(|&(_, o)| o == b'H'); - let trailing = tokens[i + 1..].iter().all(|&(_, o)| o == b'H'); + for (i, &op) in ops.iter().enumerate() { + if op >> 4 == 0 { + return false; + } + match op & 0xF { + 0 | 2 | 3 | 7 | 8 => saw_ref_op = true, // M, D, N, =, X + 1 | 6 => {} // I, P: interior op, no placement constraint + 4 => { + // S: only inside the hard clips at either end. + let leading = ops[..i].iter().all(is_hard); + let trailing = ops[i + 1..].iter().all(is_hard); if !leading && !trailing { - return None; // internal soft clip: invalid placement + return false; } } - b'H' if i == 0 || i == last => {} // hard clip only permitted at the ends - // An internal hard clip is invalid placement; the tokenizer already rejected unknown - // operators, so every other operator byte here also fails closed. - _ => return None, + 5 if i == 0 || i == last => {} // H: only as the first/last operation + // An internal hard clip, or a code past `X`, is invalid. + _ => return false, } - // `len <= MAX_CIGAR_OP_LEN` (28 bits), so the shift cannot drop bits. - ops.push((len << 4) | op_code); - } - // A mapped mate's CIGAR must consume reference; a fully-clipped CIGAR (`10S`) cannot. - if !saw_ref_op { - return None; } - - Some(ops) + // A mapped record's CIGAR must consume reference; a fully-clipped CIGAR (`10S`) cannot. + saw_ref_op } /// BAM operation code for a CIGAR operator character, in `BAM_CIGAR_STR` (`MIDNSHP=X`) order, or @@ -842,6 +877,97 @@ mod tests { assert_eq!(num_bases_extending_past_mate_raw(&rec), 0); } + /// A forward `100M` read at 1-based 1001 whose reverse mate starts at 1-based 1051, carrying + /// `mc` as its MC tag when given. + fn forward_read_with_mate_at_1051(mc: Option<&str>) -> Vec { + let mut aux = Vec::new(); + if let Some(mc) = mc { + aux.extend_from_slice(b"MCZ"); + aux.extend_from_slice(mc.as_bytes()); + aux.push(0); + } + make_bam_bytes_with_tlen( + 0, + 1000, + flags::PAIRED | flags::MATE_REVERSE, + b"rea", + &[encode_op(0, 100)], + 100, + 0, + 1050, + 90, + &aux, + ) + } + + /// A mate CIGAR taken from the mate record is held to the same structural check as an MC + /// tag: an empty CIGAR (an unmapped mate's `*`, or a truncated record), a CIGAR with no + /// reference-consuming operation, a misplaced clip, a zero-length operation or an unknown + /// operation code is never clipped against. Before that check an empty mate CIGAR gave the + /// mate a reference span of 0, which passed the FR gate and clipped the forward read at the + /// mate's start, removing 50 aligned bases. + #[rstest] + #[case::empty(&[])] + #[case::soft_clip_only(&[encode_op(4, 100)])] + #[case::internal_soft_clip(&[encode_op(0, 10), encode_op(4, 5), encode_op(0, 10)])] + #[case::internal_hard_clip(&[encode_op(0, 10), encode_op(5, 5), encode_op(0, 10)])] + #[case::zero_length_op(&[encode_op(0, 0), encode_op(0, 40)])] + #[case::unknown_op_code(&[encode_op(0, 40), 9])] + fn mate_cigar_path_rejects_what_the_mc_path_rejects(#[case] mate_ops: &[u32]) { + let rec = forward_read_with_mate_at_1051(None); + assert_eq!(num_bases_extending_past_mate_cigar_raw(&rec, mate_ops), 0); + } + + /// A valid mate CIGAR clips the same whether it comes from the MC tag or is passed in: the + /// `40M` mate ends at 1090, 10 bases before the read's end at 1100. + #[test] + fn mate_cigar_path_matches_the_mc_path_for_a_valid_cigar() { + let with_mc = forward_read_with_mate_at_1051(Some("40M")); + let without_mc = forward_read_with_mate_at_1051(None); + assert_eq!(num_bases_extending_past_mate_raw(&with_mc), 10); + assert_eq!(num_bases_extending_past_mate_cigar_raw(&without_mc, &[encode_op(0, 40)]), 10); + } + + /// `num_bases_extending_past_mate_from_mc_raw` returns `None` only when the MC tag is + /// absent; an MC tag that is present but invalid fails closed to `Some(0)`. + #[rstest] + #[case::absent(None, None)] + #[case::valid(Some("40M"), Some(10))] + #[case::unmapped_mate_star(Some("*"), Some(0))] + #[case::soft_clip_only(Some("100S"), Some(0))] + fn from_mc_distinguishes_an_absent_tag( + #[case] mc: Option<&str>, + #[case] expected: Option, + ) { + let rec = forward_read_with_mate_at_1051(mc); + assert_eq!(num_bases_extending_past_mate_from_mc_raw(&rec), expected); + assert_eq!(num_bases_extending_past_mate_raw(&rec), expected.unwrap_or(0)); + } + + // ======================================================================== + // is_mapped_opposite_strand_pair_raw tests + // ======================================================================== + + /// Each case differs from `fr_candidate` in exactly one guard, so each guard is what rejects + /// it. + #[rstest] + #[case::fr_candidate(flags::PAIRED | flags::MATE_REVERSE, 0, true)] + #[case::reverse_read(flags::PAIRED | flags::REVERSE, 0, true)] + #[case::not_paired(flags::MATE_REVERSE, 0, false)] + #[case::read_unmapped(flags::PAIRED | flags::MATE_REVERSE | flags::UNMAPPED, 0, false)] + #[case::mate_unmapped(flags::PAIRED | flags::MATE_REVERSE | flags::MATE_UNMAPPED, 0, false)] + #[case::mate_on_other_reference(flags::PAIRED | flags::MATE_REVERSE, 1, false)] + #[case::both_forward(flags::PAIRED, 0, false)] + #[case::both_reverse(flags::PAIRED | flags::REVERSE | flags::MATE_REVERSE, 0, false)] + fn test_is_mapped_opposite_strand_pair_raw( + #[case] flag: u16, + #[case] mate_tid: i32, + #[case] expected: bool, + ) { + let rec = make_bam_bytes(0, 100, flag, b"rea", &[encode_op(0, 10)], 10, mate_tid, 200, &[]); + assert_eq!(is_mapped_opposite_strand_pair_raw(&rec), expected); + } + /// A mate record whose trailing soft clips sum past `i32::MAX` must clamp the mate /// boundary outward (to `i32::MAX`), not collapse to `0`. Collapsing to `0` shrank /// `mate_unclipped_end`, making the read look like it extends further past the mate diff --git a/docs/src/guide/best-practices.md b/docs/src/guide/best-practices.md index d1233f373..de68b4717 100644 --- a/docs/src/guide/best-practices.md +++ b/docs/src/guide/best-practices.md @@ -55,9 +55,13 @@ missing `MC` is caught up front; the rest degrade silently: longer filtered on its mate's alignment quality. - `fgumi simplex` and `fgumi duplex` clip the bases at each read's 3' end that run past its mate — read-through past the fragment boundary, which is adapter rather than template — and - they locate that boundary from the read's own `MC` tag. A read with no `MC` tag gets no clip, - so those read-through bases are called into the consensus, silently and without a warning. - (`fgumi codec` is unaffected: it reads the boundary from the mate record in hand.) + they locate that boundary from the read's own `MC` tag. A read with no `MC` tag takes the mate + CIGAR from its mate in the same UMI group instead, as fgbio does; if that mate is mapped to + the same reference on the opposite strand but is not in the group, the run fails with an + error naming the read, once that read would reach a consensus (see the + [fgbio migration guide](migration-from-fgbio.md)). An `MC` tag that is present is trusted as + written, so a stale `MC` still mis-places the boundary silently. (`fgumi codec` reads the + boundary from the mate record in hand and never consults `MC`.) Correct mate tags are necessary but not sufficient: `fgumi group` and `fgumi dedup` also constrain the sort order, and `fgumi dedup` requires a tag that no mate-fixing tool writes. diff --git a/docs/src/guide/consensus-calling.md b/docs/src/guide/consensus-calling.md index fdd910be6..34118aafb 100644 --- a/docs/src/guide/consensus-calling.md +++ b/docs/src/guide/consensus-calling.md @@ -84,6 +84,7 @@ Any consensus base with quality below the minimum threshold is masked to `N`. - Each end of a pair is treated independently; overlapping bases within a pair are jointly called by default (disable with `--consensus-call-overlapping-bases false`) - Indel errors in the reads are not considered in the consensus model +- Bases at a read's 3' end that extend past its mate (read-through into adapter) are clipped before consensus calling. `simplex` and `duplex` find the mate's alignment from the read's `MC` (mate CIGAR) tag; a read without `MC` takes the CIGAR of its mate in the same UMI group instead, and a read lacking `MC` whose mate is mapped to the same reference on the opposite strand but is not in the group is an error naming the read. That error is raised only for a read that can still reach a consensus: in `simplex`, a read in an emitted consensus (not one rejected by `--min-reads`, the alignment filter or `--max-reads`, or left in an orphaned single-end consensus); in `duplex`, a read of a group that passes the `--min-reads` and strand-orientation checks. An `MC` tag that is present is used as written, so it must match the mate's alignment (`fgumi zipper` writes it correctly); an in-group mate that is unmapped, or whose CIGAR is not valid for a mapped read, is ignored just as an invalid `MC` tag is. Only primary reads are filled in. `codec` always uses the in-group mate. - `simplex` and `codec` do not accept a `--sort-order` flag; consensus reads are emitted as unmapped and should be sorted by the downstream pipeline (`fgumi zipper` + `fgumi sort`) - By default, reads that are unmapped and lack a mapped mate — along with all secondary and diff --git a/docs/src/guide/migration-from-fgbio.md b/docs/src/guide/migration-from-fgbio.md index 1b958dd3f..0783faa96 100644 --- a/docs/src/guide/migration-from-fgbio.md +++ b/docs/src/guide/migration-from-fgbio.md @@ -176,6 +176,30 @@ Use the `--metrics PREFIX` flag to write all three files in one step. fgumi's `simplex` and `duplex` stats output uses the same three-column key-value format as fgbio's `CallMolecularConsensusReads`, allowing direct comparison with `fgumi compare metrics`. +### Mate CIGAR (`MC`) Handling in Consensus Calling + +Like fgbio, `simplex` and `duplex` need each read's mate CIGAR to clip read-through bases past +the mate, and fill it in from the read's mate in the same UMI group when the `MC` tag is absent. +Four edge cases differ: + +- fgbio overwrites the `MC` of **both** reads of a template whenever either one lacks it; fgumi + keeps an `MC` that is present and fills in only the missing one. +- fgbio fails with a `MatchError` when a read name lacking `MC` does not appear exactly twice in + the group. fgumi fails only when a read lacking `MC` has no mate in the group, with an error + naming the read; when a name appears more than twice, the first mate found is used. +- fgbio backfills, and so can fail on, only pairs it judges FR, which it infers from `TLEN`. + fgumi requires the mate for any read lacking `MC` whose mate is mapped to the same reference + on the opposite strand, RF pairs included, because `TLEN` misclassifies dovetailed FR pairs. + An RF pair lacking `MC` whose mate is missing from the group is therefore an error in fgumi + but not in fgbio. +- fgbio backfills every group before any filtering, so a missing mate fails the run even when + the group would produce no consensus. fgumi fails only for a read that can still reach a + consensus. In `simplex` that is a read in an emitted consensus, so a read whose end falls + below `--min-reads`, that is dropped as a minority alignment or by `--max-reads`, or whose + single-end consensus is an orphan is rejected as usual. In `duplex` it is any read of a group + that passes the `--min-reads` checks and the strand-orientation (potential UMI collision) + check. + ### Sort Orders fgumi's `sort` command supports the same sort orders as fgbio: diff --git a/src/lib/commands/duplex.rs b/src/lib/commands/duplex.rs index 727f67067..fc0fd088c 100644 --- a/src/lib/commands/duplex.rs +++ b/src/lib/commands/duplex.rs @@ -46,6 +46,11 @@ Because of the nature of duplex sequencing, this tool does not support fragment input they are ignored. Similarly, read pairs for which consensus reads cannot be generated for one or other read (R1 or R2) are omitted from the output. +Bases at a read's 3' end that extend past its mate (read-through into adapter) are trimmed before consensus +calling. The mate's alignment comes from the read's MC tag or, when MC is absent, from its mate in the same group; +a read lacking MC whose mate is mapped to the same reference on the opposite strand but is not in the group is an +error when the read is in a group that passes the --min-reads and strand-orientation checks. + The consensus reads produced are unaligned, due to the difficulty and error-prone nature of inferring the consensus alignment. Consensus reads should therefore be aligned after, which should not be too expensive as likely there are far fewer consensus reads than input raw reads. diff --git a/src/lib/commands/simplex.rs b/src/lib/commands/simplex.rs index 8cb54eb5e..032352a12 100644 --- a/src/lib/commands/simplex.rs +++ b/src/lib/commands/simplex.rs @@ -64,6 +64,11 @@ overlap within a read pair are consensus called jointly before UMI consensus cal --consensus-call-overlapping-bases false to call each end of a pair independently, as fgbio does. Insertion or deletion errors in the reads are not considered in the consensus model. +Bases at a read's 3' end that extend past its mate (read-through into adapter) are trimmed before consensus +calling. The mate's alignment comes from the read's MC tag or, when MC is absent, from its mate in the same group; +a read lacking MC whose mate is mapped to the same reference on the opposite strand but is not in the group is an +error when the read is in an emitted consensus. + The consensus reads produced are unaligned, due to the difficulty and error-prone nature of inferring the consensus alignment. Consensus reads should therefore be aligned after, which should not be too expensive as likely there are far fewer consensus reads than input raw reads. diff --git a/tests/integration/test_duplex_command.rs b/tests/integration/test_duplex_command.rs index f4b394009..557c1fbf8 100644 --- a/tests/integration/test_duplex_command.rs +++ b/tests/integration/test_duplex_command.rs @@ -1789,6 +1789,8 @@ fn test_duplex_command_emits_unmapped_consensus_header(#[case] threads: Option<& /// [`create_duplex_read_pair_with_sequences`] so it lands in the same molecule. /// /// The R2 keeps its REVERSE flag so strand grouping is unaffected; only the mapping is removed. +/// The R1 carries `MATE_UNMAPPED` to match: the pre-group filter drops the unmapped R2, and a +/// mapped read with no `MC` that claims a mapped mate absent from its group is an error. fn create_half_mapped_a_strand_pair( name: &str, mi_tag: &str, @@ -1811,7 +1813,9 @@ fn create_half_mapped_a_strand_pair( b.read_name(name.as_bytes()) .sequence(r1_sequence.as_bytes()) .qualities(&vec![quality; read_len]) - .flags(flags::PAIRED | flags::FIRST_SEGMENT | flags::MATE_REVERSE) + .flags( + flags::PAIRED | flags::FIRST_SEGMENT | flags::MATE_REVERSE | flags::MATE_UNMAPPED, + ) .ref_id(0) .pos(ref_start - 1) .mapq(60) @@ -2487,3 +2491,165 @@ fn test_duplex_chain_ref_without_methylation_mode_rejected_same_as_oracle() { "--ref without --methylation-mode must be rejected with the same message on both paths", ); } + +// ───────────────────────────────────────────────────────────────────────────── +// Mate-CIGAR (`MC`) backfill for read-through clipping. +// ───────────────────────────────────────────────────────────────────────────── + +/// 80 bp insert of the read-through duplex templates, at 1-based 1001-1080. +const READ_THROUGH_INSERT: &str = + "ACCTCTCCATCTGACCCAAGATTGTGCTTGTTCAATTCTTCTTAACGTGATAACAGAATCAAACCTGCCAGGCGGTCGTC"; +/// 20 bp adapter each read sequences past the end of the insert. +const READ_THROUGH_ADAPTER: &str = "AGATCGGAAGAGCACACGTC"; + +/// One read-through duplex template `(R1, R2)`: a forward `80M20S` read and a reverse `20S80M` +/// read, both starting at 1-based 1001, so each reads 20 bp of adapter past its mate. On the AB +/// strand R1 is the forward read, on the BA strand R2 is. Each read carries its mate's CIGAR in +/// `MC` when `with_mc` is set. +fn read_through_duplex_template( + name: &str, + mi: &str, + is_b_strand: bool, + with_mc: bool, +) -> (RawRecord, RawRecord) { + let forward_seq = format!("{READ_THROUGH_INSERT}{READ_THROUGH_ADAPTER}"); + // BAM stores a reverse read's bases on the forward reference strand. + let reverse_seq = format!("{}{READ_THROUGH_INSERT}", reverse_complement(READ_THROUGH_ADAPTER)); + let build = |segment: u16, reverse: bool| { + let (strand, seq, cigar, mate_cigar, tlen) = if reverse { + (flags::REVERSE, &reverse_seq, [(20u32 << 4) | 4, 80 << 4], "80M20S", -80) + } else { + (flags::MATE_REVERSE, &forward_seq, [80u32 << 4, (20 << 4) | 4], "20S80M", 80) + }; + let mut b = SamBuilder::new(); + b.read_name(name.as_bytes()) + .sequence(seq.as_bytes()) + .qualities(&[40; 100]) + .flags(flags::PAIRED | flags::PROPER_PAIR | segment | strand) + .ref_id(0) + .pos(1000) + .mapq(60) + .cigar_ops(&cigar) + .mate_ref_id(0) + .mate_pos(1000) + .template_length(tlen) + .add_string_tag(SamTag::MI, mi.as_bytes()); + if with_mc { + b.add_string_tag(SamTag::MC, mate_cigar.as_bytes()); + } + b.build() + }; + (build(flags::FIRST_SEGMENT, is_b_strand), build(flags::LAST_SEGMENT, !is_b_strand)) +} + +/// One read-through duplex molecule (MI `1`) with three templates per strand. +fn read_through_duplex_molecule(with_mc: bool) -> Vec<(RawRecord, RawRecord)> { + let mut molecule = Vec::new(); + for (strand, is_b_strand) in [("A", false), ("B", true)] { + for i in 0..3 { + let name = format!("{}_{i}", strand.to_lowercase()); + molecule.push(read_through_duplex_template( + &name, + &format!("1/{strand}"), + is_b_strand, + with_mc, + )); + } + } + molecule +} + +/// The `duplex` arguments for the mate-CIGAR backfill tests, with `--threads` when given. +fn duplex_backfill_args<'a>( + input: &'a Path, + output: &'a Path, + threads: Option<&'a str>, +) -> Vec<&'a str> { + let mut args = vec![ + "duplex", + "--input", + input.to_str().unwrap(), + "--output", + output.to_str().unwrap(), + "--min-reads", + "1", + ]; + if let Some(threads) = threads { + args.extend(["--threads", threads]); + } + args +} + +/// A read with no `MC` tag takes its mate CIGAR from its mate in the same molecule, as fgbio's +/// `UmiConsensusCaller.updateMateCigars` (`UmiConsensusCaller.scala:328-343`, called at `:356`) +/// does for `CallDuplexConsensusReads` too. The `/A` and `/B` reads are merged into one group, so +/// each read's mate is found and the read-through molecule is trimmed exactly as when it carries +/// `MC`: both duplex consensus reads keep the 80 bp insert. Checked on the single-worker and the +/// multi-worker chain. +#[rstest] +#[case::single_threaded(None)] +#[case::threaded(Some("2"))] +fn test_duplex_backfills_a_missing_mate_cigar_from_the_mate(#[case] threads: Option<&str>) { + let temp_dir = TempDir::new().unwrap(); + let with_mc = temp_dir.path().join("with_mc.bam"); + let without_mc = temp_dir.path().join("without_mc.bam"); + create_duplex_bam(&with_mc, vec![read_through_duplex_molecule(true)]); + create_duplex_bam(&without_mc, vec![read_through_duplex_molecule(false)]); + + let with_mc_out = temp_dir.path().join("with_mc.out.bam"); + let without_mc_out = temp_dir.path().join("without_mc.out.bam"); + for (input, output) in [(&with_mc, &with_mc_out), (&without_mc, &without_mc_out)] { + Duplex::try_parse_from(duplex_backfill_args(input, output, threads)) + .expect("failed to parse duplex args") + .execute("fgumi duplex") + .expect("duplex run failed"); + } + + let (_, expected) = crate::helpers::read_bam_output(&with_mc_out); + let (_, actual) = crate::helpers::read_bam_output(&without_mc_out); + let lengths: Vec = actual.iter().map(|r| r.sequence().len()).collect(); + assert_eq!(lengths, vec![80, 80], "both ends trimmed to the insert"); + assert_eq!(actual, expected, "output without MC must match output with MC"); +} + +/// A read with no `MC` tag whose mate is missing from its molecule cannot be trimmed at the +/// mate's end, so `duplex` fails naming the molecule and the read rather than calling the +/// read-through bases into the consensus (fgbio fails too, in +/// `UmiConsensusCaller.updateMateCigars`, `UmiConsensusCaller.scala:328-343`). The error must +/// surface, with the duplex step's `MI` context, from the multi-worker chain too. +#[rstest] +#[case::single_threaded(None)] +#[case::threaded(Some("2"))] +fn test_duplex_fails_naming_a_read_whose_missing_mate_cigar_cannot_be_backfilled( + #[case] threads: Option<&str>, +) { + let temp_dir = TempDir::new().unwrap(); + let input = temp_dir.path().join("input.bam"); + let output = temp_dir.path().join("output.bam"); + let header = create_minimal_header("chr1", 10000); + let mut writer = bam::io::Writer::new(fs::File::create(&input).expect("create input BAM")); + writer.write_header(&header).expect("write header"); + // Drop the R2 of template `a_0`, leaving its R1 without MC or a mate. + for (i, (r1, r2)) in read_through_duplex_molecule(false).into_iter().enumerate() { + writer.write_alignment_record(&header, &to_record_buf(&r1)).expect("write R1"); + if i != 0 { + writer.write_alignment_record(&header, &to_record_buf(&r2)).expect("write R2"); + } + } + writer.try_finish().expect("finish input BAM"); + + let err = Duplex::try_parse_from(duplex_backfill_args(&input, &output, threads)) + .expect("failed to parse duplex args") + .execute("fgumi duplex") + .expect_err("a read lacking MC whose mate is absent must fail the run"); + + // Both paths run the consensus as a pipeline step, which flattens the caller's error into + // one message prefixed with the step and the molecule. + assert_eq!( + err.root_cause().to_string(), + "Pipeline::run: step \"DuplexConsensus\" failed: Duplex consensus error for MI 1\t: Mate \ + cigar (MC SAM tag) needed for read 'a_0': the read has no MC tag and its primary mate is \ + not in the same group. Add MC tags (e.g. with fgumi zipper or samtools fixmate), or keep \ + both reads of each template in the same group." + ); +} diff --git a/tests/integration/test_simplex_command.rs b/tests/integration/test_simplex_command.rs index c314f6723..241ad3d1c 100644 --- a/tests/integration/test_simplex_command.rs +++ b/tests/integration/test_simplex_command.rs @@ -925,6 +925,99 @@ fn test_simplex_indel_at_overlap_boundary_still_calls_a_consensus() { } } +/// [`indel_readthrough_family`] with every `MC` tag removed. +fn indel_readthrough_family_without_mc(depth: usize) -> Vec { + indel_readthrough_family(depth) + .into_iter() + .map(|mut record| { + fgumi_raw_bam::remove_tag(record.as_mut_vec(), SamTag::MC); + record + }) + .collect() +} + +/// The `simplex` arguments for the mate-CIGAR backfill tests, with `--threads` when given. +fn backfill_args<'a>(input: &'a Path, output: &'a Path, threads: Option<&'a str>) -> Vec<&'a str> { + let mut args = vec![ + "simplex", + "--input", + input.to_str().unwrap(), + "--output", + output.to_str().unwrap(), + "--min-reads", + "1", + ]; + if let Some(threads) = threads { + args.extend(["--threads", threads]); + } + args +} + +/// A read with no `MC` tag takes its mate CIGAR from its mate in the same group, as fgbio's +/// `UmiConsensusCaller.updateMateCigars` (`UmiConsensusCaller.scala:328-343`, called at `:356`) +/// does, so the read-through family is trimmed exactly as when it carries `MC`: both consensus +/// reads keep 127 of their 129 bases. Checked on the single-worker and the multi-worker chain. +#[rstest] +#[case::single_threaded(None)] +#[case::threaded(Some("2"))] +fn test_simplex_backfills_a_missing_mate_cigar_from_the_mate(#[case] threads: Option<&str>) { + let temp_dir = TempDir::new().unwrap(); + let with_mc = temp_dir.path().join("with_mc.bam"); + let without_mc = temp_dir.path().join("without_mc.bam"); + create_grouped_bam(&with_mc, vec![("1", indel_readthrough_family(3))]); + create_grouped_bam(&without_mc, vec![("1", indel_readthrough_family_without_mc(3))]); + + let with_mc_out = temp_dir.path().join("with_mc.out.bam"); + let without_mc_out = temp_dir.path().join("without_mc.out.bam"); + for (input, output) in [(&with_mc, &with_mc_out), (&without_mc, &without_mc_out)] { + Simplex::try_parse_from(backfill_args(input, output, threads)) + .expect("failed to parse simplex args") + .execute("fgumi simplex") + .expect("simplex run failed"); + } + + let (_, expected) = read_bam_output(&with_mc_out); + let (_, actual) = read_bam_output(&without_mc_out); + let lengths: Vec = actual.iter().map(|r| r.sequence().len()).collect(); + assert_eq!(lengths, vec![127, 127], "both ends trimmed by the 2 bases past the mate"); + assert_eq!(actual, expected, "output without MC must match output with MC"); +} + +/// A read with no `MC` tag whose mate is missing from its group cannot be trimmed at the mate's +/// end, so `simplex` fails naming the read rather than calling the read-through bases into the +/// consensus (fgbio fails too, in `UmiConsensusCaller.updateMateCigars`, +/// `UmiConsensusCaller.scala:328-343`). The error, naming the molecule and the read, must +/// surface from the multi-worker chain too. +#[rstest] +#[case::single_threaded(None)] +#[case::threaded(Some("2"))] +fn test_simplex_fails_naming_a_read_whose_missing_mate_cigar_cannot_be_backfilled( + #[case] threads: Option<&str>, +) { + let temp_dir = TempDir::new().unwrap(); + let input = temp_dir.path().join("input.bam"); + let output = temp_dir.path().join("output.bam"); + // Drop the R2 of template `rt_0`, the second record of the family. + let mut family = indel_readthrough_family_without_mc(3); + family.remove(1); + create_grouped_bam(&input, vec![("1", family)]); + + let err = Simplex::try_parse_from(backfill_args(&input, &output, threads)) + .expect("failed to parse simplex args") + .execute("fgumi simplex") + .expect_err("a read lacking MC whose mate is absent must fail the run"); + + // Both paths run the consensus as a pipeline step, which flattens the caller's error into + // one message prefixed with the step and the molecule. + assert_eq!( + err.root_cause().to_string(), + "Pipeline::run: step \"SimplexConsensus\" failed: Consensus error for MI 1\t: Mate cigar \ + (MC SAM tag) needed for read 'rt_0': the read has no MC tag and its primary mate is not \ + in the same group. Add MC tags (e.g. with fgumi zipper or samtools fixmate), or keep both \ + reads of each template in the same group." + ); +} + ////////////////////////////////////////////////////////////////////////////// // Chain-path worker-count determinism tests //