diff --git a/README.md b/README.md index b8c8efbfc..fa3354929 100644 --- a/README.md +++ b/README.md @@ -33,6 +33,13 @@ ft --help [Help page for fibertools](https://fiberseq.github.io/fibertools/help.html#ft) +## Input BAM requirements + +Fiber-seq tags (`MM`/`ML`, `ns`/`nl`/`as`/`al`, `Ma`) describe the full read. Aligners that hard-clip supplementary alignments copy those tags unchanged onto the clipped record. When the copied `Ma` or legacy `ns`/`nl`/`as`/`al` tags describe the full read, `ft` keeps the nucleosome, MSP and FIRE calls and lifts them with the hard-clip offset, but the `MM`/`ML` m6A cannot be recovered and is dropped, so those reads count as NotCallable. Tags that match neither the record nor the full read are dropped. `ft` warns in both cases. To keep m6A on supplementary alignments, align with soft clipping: + +- PacBio: `pbmm2 align` (it never hard-clips). +- ONT: `dorado aligner --mm2-opts "-Y" ...`, or `samtools fastq -T '*' in.bam | minimap2 -Y -y -ax map-ont ref.fa -`. minimap2 also writes SEQ-less secondary alignments by default; their tags are dropped too, so add `--secondary=no` or filter with `-F 256`. + # Highlighted subcommands for `fibertools-rs` ### `ft predict-m6a` diff --git a/molecular-annotation/CHANGELOG.md b/molecular-annotation/CHANGELOG.md index 5f2d44272..c82701bca 100644 --- a/molecular-annotation/CHANGELOG.md +++ b/molecular-annotation/CHANGELOG.md @@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- `AlignedBlocks::with_query_offset` / `query_offset`, `MolecularAnnotations::query_offset` / `set_query_offset`, and `liftover::{hard_clips, query_span, full_read_query_offset}`: lift full-read MA annotations on hard-clipped records ([#136](https://github.com/fiberseq/fibertools-rs/issues/136)). `from_record` (Rust and pyMA) no longer parses MM/ML in that frame, and `project_query` is SEQ-relative. pyMA `set_aligned_blocks` gains a `query_offset` keyword. + ## [0.0.3](https://github.com/fiberseq/fibertools-rs/compare/molecular-annotation-v0.0.2...molecular-annotation-v0.0.3) - 2026-08-13 ### Fixed diff --git a/molecular-annotation/python/molecular_annotation/pysam_utils.py b/molecular-annotation/python/molecular_annotation/pysam_utils.py index d8e775929..a1f60e400 100644 --- a/molecular-annotation/python/molecular_annotation/pysam_utils.py +++ b/molecular-annotation/python/molecular_annotation/pysam_utils.py @@ -99,6 +99,38 @@ def _parse_mm_ml_into( annot.parse_mm_ml(mm, ml, forward_seq) +_BAM_CHARD_CLIP = 5 + + +def _hard_clips(record: "pysam.AlignedSegment") -> tuple[int, int]: + """(leading, trailing) hard clip lengths, in BAM/CIGAR orientation. + + The CIGAR is in reference orientation, so "leading" is the same end on + both strands: the bases before the first base of SEQ. + """ + ct = record.cigartuples or [] + lead = ct[0][1] if ct and ct[0][0] == _BAM_CHARD_CLIP else 0 + trail = ct[-1][1] if len(ct) > 1 and ct[-1][0] == _BAM_CHARD_CLIP else 0 + return lead, trail + + +def _query_span(record: "pysam.AlignedSegment") -> int: + """Query bases the record holds: len(SEQ), or the CIGAR query span + (hard clips excluded) when SEQ is absent.""" + if record.query_sequence is not None: + return len(record.query_sequence) + return record.infer_query_length() or 0 + + +def _full_read_query_offset(read_length: int, record: "pysam.AlignedSegment"): + """Mirror of the Rust `full_read_query_offset`: the leading hard clip + when `read_length` spans the record's hard-clipped bases, else None.""" + lead, trail = _hard_clips(record) + if lead + trail > 0 and read_length == _query_span(record) + lead + trail: + return lead + return None + + def from_record( record: "pysam.AlignedSegment", parse_tags: bool = True ) -> MolecularAnnotations: @@ -117,12 +149,18 @@ def from_record( Returns: MolecularAnnotations object with aligned blocks set for liftover support. + On a hard-clipped record whose MA read length spans the hard-clipped + bases (an aligner copied the full read's tags onto the clipped + record), the MA annotations are kept in the full read's frame, the + liftover carries the leading hard clip as `query_offset`, and MM/ML + are not parsed (they are SEQ-relative and cannot be recovered). Raises: KeyError: If parse_tags=True and MA tag is missing ValueError: If tag format is invalid """ is_reverse = record.is_reverse + query_offset = 0 if parse_tags: # Resolve the Ma/Aq/An family atomically (required when parsing @@ -135,8 +173,17 @@ def from_record( annot = MolecularAnnotations.from_tags(ma, aq=aq, an=an) annot.is_reverse_aligned = is_reverse - # Base modifications (m6A/5mC/…) live in MM/ML, not the MA tag set. - _parse_mm_ml_into(annot, record, is_reverse) + # A hard-clipping aligner copies the full read's tags onto the + # clipped record. The MA family is then still correct (molecular + # coordinates of the full read); only the lift needs to know SEQ + # starts H_lead bases in. MM/ML are SEQ-relative and cannot be + # recovered, so they are not parsed in that frame. + offset = _full_read_query_offset(annot.read_length, record) + if offset is not None: + query_offset = offset + else: + # Base modifications (m6A/5mC/…) live in MM/ML, not the MA tag set. + _parse_mm_ml_into(annot, record, is_reverse) else: # Create empty annotations with just read length annot = MolecularAnnotations(record.query_length) @@ -148,7 +195,9 @@ def from_record( if not record.is_unmapped and record.cigartuples: aligned_blocks = _extract_aligned_blocks(record) if aligned_blocks: - annot.set_aligned_blocks(aligned_blocks, is_reverse=is_reverse) + annot.set_aligned_blocks( + aligned_blocks, is_reverse=is_reverse, query_offset=query_offset + ) return annot diff --git a/molecular-annotation/python/src/lib.rs b/molecular-annotation/python/src/lib.rs index 975d217c3..bc569c4a6 100644 --- a/molecular-annotation/python/src/lib.rs +++ b/molecular-annotation/python/src/lib.rs @@ -40,6 +40,7 @@ use pyo3::prelude::*; // Use fully qualified path to avoid name collision with the pymodule use ::molecular_annotation::{ + AlignedBlocks as RustAlignedBlocks, Annotation as RustAnnotation, Encoding as RustEncoding, MolecularAnnotations as RustMolecularAnnotations, @@ -590,26 +591,42 @@ impl MolecularAnnotations { /// Set aligned blocks for liftover calculations. /// - /// Accepts 0-based half-open [start, end) intervals. + /// Accepts 0-based half-open [start, end) intervals. Block query + /// coordinates are SEQ-relative (as pysam reports them). /// /// Args: /// blocks: List of ((query_start, query_end), (ref_start, ref_end)) tuples. /// is_reverse: Whether the read is reverse-aligned. + /// query_offset: Where SEQ starts in the annotation frame. Pass the + /// leading hard clip when the annotations describe the full + /// read and this record is a hard-clipped part of it; 0 otherwise. /// /// Example: /// >>> # Query [0, 500) aligns to reference [1000, 1500) /// >>> annot.set_aligned_blocks([((0, 500), (1000, 1500))], is_reverse=False) - #[pyo3(signature = (blocks, is_reverse=false))] + #[pyo3(signature = (blocks, is_reverse=false, query_offset=0))] pub fn set_aligned_blocks( &mut self, blocks: Vec<((u32, u32), (u32, u32))>, is_reverse: bool, + query_offset: u32, ) { let block_pairs: Vec<([u32; 2], [u32; 2])> = blocks .into_iter() .map(|((qs, qe), (rs, re))| ([qs, qe], [rs, re])) .collect(); - self.inner.set_aligned_blocks(block_pairs, is_reverse); + let blocks = RustAlignedBlocks::new(block_pairs, self.inner.read_length) + .with_query_offset(query_offset); + self.inner.set_aligned_blocks_raw(blocks, is_reverse); + } + + /// Where SEQ starts in the annotation frame: 0 unless the annotations + /// describe the full read and this record is a hard-clipped part of it. + /// Query coordinates from get_coords / get_ref_coords / iter_type are in + /// the annotation frame; subtract this before indexing the sequence. + #[getter] + pub fn query_offset(&self) -> u32 { + self.inner.query_offset() } /// Check if aligned blocks are set. diff --git a/molecular-annotation/python/tests/test_molecular_annotation.py b/molecular-annotation/python/tests/test_molecular_annotation.py index 86a45c015..8db28c7d0 100644 --- a/molecular-annotation/python/tests/test_molecular_annotation.py +++ b/molecular-annotation/python/tests/test_molecular_annotation.py @@ -888,3 +888,86 @@ def test_add_annotation_then_write(self, tmp_path): assert len(annot2.iter_type("nuc")) == nuc_before # ...and the MM/ML base mods are still intact. assert len(annot2.iter_type("a")) == 1541 + + +class TestFullReadFrame: + """A hard-clipped record whose MA read length spans the hard-clipped + bases keeps its full-read annotations; only the lift carries the + leading hard clip (#136).""" + + def _rec(self, pysam, cigar, seq_len, ma, flag=0, mm=False): + r = TestAlignedBlocks()._record(pysam, cigar, seq_len) # ref_start 1000 + r.flag = flag + r.set_tag("Ma", ma) + if mm: + r.set_tag("MM", "A+a,0;") + r.set_tag("ML", [200]) + return r + + def test_hard_clips(self): + pysam = pytest.importorskip("pysam") + from molecular_annotation.pysam_utils import _hard_clips + + rec = TestAlignedBlocks()._record + assert _hard_clips(rec(pysam, "5H10S40M2D40M10S5H", 100)) == (5, 5) + assert _hard_clips(rec(pysam, "100M", 100)) == (0, 0) + assert _hard_clips(rec(pysam, "80M20H", 80)) == (0, 20) + + def test_forward_leading_hard_clip(self): + pysam = pytest.importorskip("pysam") + from molecular_annotation.pysam_utils import from_record + + # molecular [60,80) of a 150 bp read; SEQ is [50,150) + annot = from_record(self._rec(pysam, "50H100M", 100, "150;nuc.:61-20", mm=True)) + assert annot.read_length == 150 + assert annot.query_offset == 50 + assert annot.get_ref_coords("nuc") == [(60, 80, 1010, 1030)] + assert "a" not in annot.annotation_type_names(), "MM/ML must not be parsed" + + def test_reverse_trailing_hard_clip(self): + pysam = pytest.importorskip("pysam") + from molecular_annotation.pysam_utils import from_record + + annot = from_record(self._rec(pysam, "100M50H", 100, "150;nuc.:61-20", flag=16)) + assert annot.query_offset == 0 + assert annot.get_ref_coords("nuc") == [(70, 90, 1070, 1090)] + + def test_reverse_leading_hard_clip(self): + pysam = pytest.importorskip("pysam") + from molecular_annotation.pysam_utils import from_record + + annot = from_record(self._rec(pysam, "50H100M", 100, "150;nuc.:61-20", flag=16)) + assert annot.query_offset == 50 + assert annot.get_ref_coords("nuc") == [(70, 90, 1020, 1040)] + + def test_clipped_frame_has_no_offset(self): + pysam = pytest.importorskip("pysam") + from molecular_annotation.pysam_utils import from_record + + annot = from_record(self._rec(pysam, "50H100M", 100, "100;nuc.:11-20", mm=True)) + assert annot.query_offset == 0 + assert annot.get_ref_coords("nuc") == [(10, 30, 1010, 1030)] + assert "a" in annot.annotation_type_names(), "MA vouches for MM/ML" + + def test_before_seq_does_not_lift(self): + pysam = pytest.importorskip("pysam") + from molecular_annotation.pysam_utils import from_record + + annot = from_record(self._rec(pysam, "50H100M", 100, "150;nuc.:11-20")) + assert annot.get_ref_coords("nuc") == [(10, 30, None, None)] + + def test_set_aligned_blocks_default_offset(self): + annot = MolecularAnnotations(80) + annot.set_aligned_blocks([((0, 80), (1000, 1080))], is_reverse=False) + assert annot.query_offset == 0 + annot.set_aligned_blocks([((0, 80), (1000, 1080))], query_offset=10) + assert annot.query_offset == 10 + + def test_round_trip_keeps_full_read_frame(self): + pysam = pytest.importorskip("pysam") + from molecular_annotation.pysam_utils import from_record, to_record + + r = self._rec(pysam, "50H100M", 100, "150;nuc.:61-20") + annot = from_record(r) + to_record(annot, r) + assert from_record(r).to_ma_string().startswith("150;") diff --git a/molecular-annotation/src/coords.rs b/molecular-annotation/src/coords.rs index 709232dad..1e35b72bc 100644 --- a/molecular-annotation/src/coords.rs +++ b/molecular-annotation/src/coords.rs @@ -59,6 +59,26 @@ impl MolecularAnnotations { self.aligned_blocks.as_ref() } + /// Where SEQ starts in the annotation frame (see + /// [`AlignedBlocks::query_offset`]); 0 without aligned blocks or when + /// the annotation frame is SEQ. BAM-orientation query coordinates from + /// `get_coords`, `get_ref_coords` and `iter_type` are in the annotation + /// frame: subtract this (checked) before indexing SEQ with them, and + /// bound the result by SEQ's length. `project_query` already subtracts + /// it. + pub fn query_offset(&self) -> u32 { + self.aligned_blocks.as_ref().map_or(0, |b| b.query_offset()) + } + + /// Set where SEQ starts in the annotation frame on the aligned blocks + /// (see [`query_offset`](Self::query_offset)). Call after the blocks are + /// set; without blocks there is nothing to lift and this is a no-op. + pub fn set_query_offset(&mut self, offset: u32) { + if let Some(b) = self.aligned_blocks.take() { + self.aligned_blocks = Some(b.with_query_offset(offset)); + } + } + /// Check if the read is reverse-aligned. pub fn is_reverse_aligned(&self) -> bool { self.is_reverse_aligned @@ -89,8 +109,10 @@ impl MolecularAnnotations { /// Get coordinates for a specific annotation type in BAM orientation. /// /// Query coordinates are returned in **BAM orientation** (forward-oriented, matching - /// the sequence as stored in the BAM file). For reverse-aligned reads, this means - /// the coordinates are flipped from the original molecular orientation. + /// the sequence as stored in the BAM file, or, on a hard-clipped record whose + /// annotations describe the full read, that full read in BAM orientation: SEQ + /// starts at [`query_offset`](Self::query_offset)). For reverse-aligned reads, + /// this means the coordinates are flipped from the original molecular orientation. /// /// This is analogous to pysam's `modified_bases` which returns positions relative /// to the BAM sequence. For original molecular orientation, use diff --git a/molecular-annotation/src/decode.rs b/molecular-annotation/src/decode.rs index 888450645..c7808388d 100644 --- a/molecular-annotation/src/decode.rs +++ b/molecular-annotation/src/decode.rs @@ -60,6 +60,10 @@ impl MolecularAnnotations { /// from the record itself, so liftover-based getters (`ref_coords`, etc.) /// work without additional setup. /// + /// On a hard-clipped record whose MA read length spans the hard-clipped + /// bases, the aligned blocks carry the leading hard clip as query offset + /// and MM/ML are not parsed (see [`crate::full_read_query_offset`]). + /// /// **Idempotency:** this is the only public entry point that parses MM/ML. /// Each call constructs a fresh `MolecularAnnotations`; there is no API /// to re-parse into an existing object. @@ -95,8 +99,23 @@ impl MolecularAnnotations { None => Self::new(record.seq_len() as u32), }; - annot.aligned_blocks = Some(crate::AlignedBlocks::from_record(record)); + // Annotation frame vs SEQ. Aligners that hard-clip (minimap2 and + // dorado aligner without -Y) copy the full-length read's tags onto + // the clipped supplementary record. The MA family is then still + // correct: its coordinates are molecular coordinates of the full + // read, and only the lift needs to know that SEQ starts H_lead + // bases into that frame. MM/ML are SEQ-relative deltas that cannot + // be recovered without the clipped bases, so in that frame they are + // not parsed (writers strip them). + let full_read_offset = crate::liftover::full_read_query_offset(annot.read_length, record); + annot.aligned_blocks = Some( + crate::AlignedBlocks::from_record(record) + .with_query_offset(full_read_offset.unwrap_or(0)), + ); annot.is_reverse_aligned = record.is_reverse(); + if full_read_offset.is_some() { + return annot; + } // Extract MM/ML and the forward-oriented sequence off the record, then // hand the raw slices to the htslib-free parser. MM is always written diff --git a/molecular-annotation/src/iter.rs b/molecular-annotation/src/iter.rs index 94d90d030..45cd11b26 100644 --- a/molecular-annotation/src/iter.rs +++ b/molecular-annotation/src/iter.rs @@ -131,14 +131,16 @@ impl MolecularAnnotations { /// Project annotations into a query-coordinate system anchored at 0. /// - /// Each annotation's coords (in BAM orientation, matching `iter_full`) are - /// shifted by `-anchor`. If `flip` is true, each interval is reversed - /// around 0: `[a, b)` → `[-(b-1), -(a-1))`. + /// Each annotation's coords (in BAM orientation, matching `iter_full`, + /// then made SEQ-relative by subtracting [`query_offset`](Self::query_offset)) + /// are shifted by `-anchor`, so `anchor` is a SEQ position. If `flip` is + /// true, each interval is reversed around 0: `[a, b)` → `[-(b-1), -(a-1))`. pub fn project_query( &self, anchor: i64, flip: bool, ) -> impl Iterator> + '_ { + let offset = self.query_offset() as i64; self.annotation_types.iter().flat_map(move |t| { t.annotations.iter().map(move |a| { let (qs, qe) = if self.is_reverse_aligned { @@ -146,7 +148,8 @@ impl MolecularAnnotations { } else { (a.start, a.end()) }; - let (start, end) = project_interval(qs as i64, qe as i64, anchor, flip); + let (start, end) = + project_interval(qs as i64 - offset, qe as i64 - offset, anchor, flip); ProjectedAnnotation { type_name: &t.name, start, diff --git a/molecular-annotation/src/lib.rs b/molecular-annotation/src/lib.rs index d49d3010d..36575d462 100644 --- a/molecular-annotation/src/lib.rs +++ b/molecular-annotation/src/lib.rs @@ -26,6 +26,10 @@ //! | `get_coords()` | BAM | Coordinates as stored in BAM (flipped for reverse reads) | //! | `get_forward_coords()` | Molecular | Original read orientation (never flipped) | //! +//! Query coordinates are in the annotation frame. On a hard-clipped record +//! whose annotations describe the full read, SEQ starts `query_offset()` +//! bases into that frame; only the liftover applies the offset. +//! //! # Module layout //! //! The public surface lives behind re-exports from this crate root. Internally @@ -131,6 +135,8 @@ mod basemods; #[cfg(feature = "htslib")] pub use decode::ma_family_tags; +#[cfg(feature = "htslib")] +pub use liftover::{full_read_query_offset, hard_clips, query_span}; pub use liftover::{AlignedBlock, AlignedBlocks}; pub use types::{ Annotation, AnnotationInfo, AnnotationType, Encoding, LiftedCoords, MaParts, MmGroup, diff --git a/molecular-annotation/src/liftover.rs b/molecular-annotation/src/liftover.rs index f5fe8df01..d7fdf24fc 100644 --- a/molecular-annotation/src/liftover.rs +++ b/molecular-annotation/src/liftover.rs @@ -108,6 +108,14 @@ pub struct AlignedBlocks { blocks: Vec, /// Length of the query sequence pub query_len: u32, + /// Where SEQ starts in the annotation frame. Block query coordinates + /// are always SEQ-relative (from the CIGAR). Callers pass query + /// coordinates in the frame the annotations were made in; when that + /// frame is the full read and this record holds a hard-clipped part of + /// it, SEQ starts `query_offset` bases in (the leading hard clip, in + /// BAM/CIGAR orientation). `lift_to_reference` subtracts it and + /// `lift_to_query` adds it. 0 when the annotation frame is SEQ. + query_offset: u32, } impl AlignedBlocks { @@ -132,7 +140,11 @@ impl AlignedBlocks { .into_iter() .map(|([q_st, q_en], [r_st, r_en])| AlignedBlock::new(q_st, q_en, r_st, r_en)) .collect(); - Self { blocks, query_len } + Self { + blocks, + query_len, + query_offset: 0, + } } /// Create `AlignedBlocks` from an iterator of block pairs. @@ -146,7 +158,26 @@ impl AlignedBlocks { let blocks = iter .map(|([q_st, q_en], [r_st, r_en])| AlignedBlock::new(q_st, q_en, r_st, r_en)) .collect(); - Self { blocks, query_len } + Self { + blocks, + query_len, + query_offset: 0, + } + } + + /// Set where SEQ starts in the annotation frame (see `query_offset`). + /// Use the leading hard clip when the annotations describe the full + /// read and this record is a hard-clipped part of it. + pub fn with_query_offset(mut self, query_offset: u32) -> Self { + self.query_offset = query_offset; + self + } + + /// Where SEQ starts in the annotation frame; 0 unless the annotations + /// describe the full read and SEQ is a hard-clipped part of it. + #[inline] + pub fn query_offset(&self) -> u32 { + self.query_offset } /// Check if there are any aligned blocks. @@ -167,7 +198,8 @@ impl AlignedBlocks { /// Lift a range from query to reference coordinates. /// /// # Arguments - /// * `start` - 0-based query start position (inclusive) + /// * `start` - 0-based query start position (inclusive), in the + /// annotation frame (see `query_offset`) /// * `end` - 0-based query end position (exclusive) /// /// # Returns @@ -194,6 +226,13 @@ impl AlignedBlocks { /// assert_eq!(re, Some(1050)); /// ``` pub fn lift_to_reference(&self, start: u32, end: u32) -> (Option, Option) { + // Annotation frame -> SEQ-relative. Bases before SEQ (inside the + // leading hard clip) clamp to 0: a range that straddles the clip + // snaps forward into the first aligned base, exactly as it would + // across a soft clip; a range that ends before SEQ collapses to + // `start >= end` and lifts to nothing. + let start = start.saturating_sub(self.query_offset); + let end = end.saturating_sub(self.query_offset); if self.blocks.is_empty() || start >= end { return (None, None); } @@ -233,7 +272,8 @@ impl AlignedBlocks { /// * `end` - 0-based reference end position (exclusive) /// /// # Returns - /// Tuple of `(query_start, query_end)` as 0-based half-open interval. + /// Tuple of `(query_start, query_end)` as 0-based half-open interval, in + /// the annotation frame (see `query_offset`). /// Returns `(None, None)` if the range cannot be lifted. /// /// # Behavior @@ -250,6 +290,7 @@ impl AlignedBlocks { if range_len == 1 { // 1bp interval: require exact match if let Some(query_pos) = self.lift_exact_to_query(start) { + let query_pos = query_pos + self.query_offset; return (Some(query_pos), Some(query_pos + 1)); } return (None, None); @@ -267,7 +308,9 @@ impl AlignedBlocks { let query_end = self.lift_end_to_query(start, end); match (query_start, query_end) { - (Some(qs), Some(qe)) if qs < qe => (Some(qs), Some(qe)), + (Some(qs), Some(qe)) if qs < qe => { + (Some(qs + self.query_offset), Some(qe + self.query_offset)) + } _ => (None, None), } } @@ -438,6 +481,62 @@ impl AlignedBlocks { } } +/// Leading and trailing hard clips of a record, in BAM/CIGAR orientation. +/// The CIGAR is written in reference orientation, so "leading" is the same +/// end on both strands: the bases before the first base of SEQ. +#[cfg(feature = "htslib")] +pub fn hard_clips(record: &rust_htslib::bam::Record) -> (u32, u32) { + // A record with no CIGAR (unmapped, or a bare `Record::new()`) has no + // data to read; `cigar()` on it trips a debug assertion. + if record.cigar_len() == 0 { + return (0, 0); + } + let cigar = record.cigar(); + ( + cigar.leading_hardclips() as u32, + cigar.trailing_hardclips() as u32, + ) +} + +/// Query bases the record holds: the SEQ length, or, when SEQ is absent +/// (`*`), the CIGAR's query span (M/I/S/=/X; hard clips excluded). +#[cfg(feature = "htslib")] +pub fn query_span(record: &rust_htslib::bam::Record) -> u32 { + use rust_htslib::bam::record::Cigar; + let seq_len = record.seq_len() as u32; + if seq_len > 0 || record.cigar_len() == 0 { + return seq_len; + } + record + .cigar() + .iter() + .map(|c| match c { + Cigar::Match(l) + | Cigar::Ins(l) + | Cigar::SoftClip(l) + | Cigar::Equal(l) + | Cigar::Diff(l) => *l, + _ => 0, + }) + .sum() +} + +/// The query offset to lift with when `read_length` (the frame the +/// annotations were made in) is the full read and `record` is a +/// hard-clipped part of it: `Some(H_lead)` iff the record has a hard clip +/// and `read_length == query_span + H_lead + H_trail`. `None` when the +/// frame is SEQ (no hard clips, or tags computed after clipping) or when +/// it matches neither (the caller decides what to do with those). +#[cfg(feature = "htslib")] +pub fn full_read_query_offset(read_length: u32, record: &rust_htslib::bam::Record) -> Option { + let (lead, trail) = hard_clips(record); + if lead + trail > 0 && read_length == query_span(record) + lead + trail { + Some(lead) + } else { + None + } +} + #[cfg(test)] mod tests { use super::*; @@ -843,4 +942,35 @@ mod tests { let b = test_blocks(); assert_eq!(b.lift_to_query(103, 108), (Some(5), Some(7))); } + + // A full-read annotation frame on a hard-clipped record: SEQ starts + // `query_offset` bases in. Lifts subtract it, reverse lifts add it. + #[test] + fn test_query_offset_forward() { + let b = AlignedBlocks::new(vec![([0, 80], [1000, 1080])], 80).with_query_offset(10); + assert_eq!(b.query_offset(), 10); + assert_eq!(b.lift_to_reference(20, 50), (Some(1010), Some(1040))); + // straddling the clip snaps into the first aligned base + assert_eq!(b.lift_to_reference(5, 30), (Some(1000), Some(1020))); + // entirely inside the clip lifts to nothing, no underflow + assert_eq!(b.lift_to_reference(0, 10), (None, None)); + assert_eq!(b.lift_to_reference(7, 8), (None, None)); + assert_eq!(b.lift_to_query(1010, 1040), (Some(20), Some(50))); + assert_eq!(b.lift_to_query(1000, 1001), (Some(10), Some(11))); + assert_eq!( + AlignedBlocks::new(vec![([0, 80], [1000, 1080])], 80).query_offset(), + 0 + ); + assert_eq!(AlignedBlocks::default().query_offset(), 0); + } + + #[test] + fn test_query_offset_with_gaps() { + let b = test_blocks().with_query_offset(3); + // same as the unshifted [0,2) + assert_eq!(b.lift_to_reference(3, 5), (Some(100), Some(102))); + // unshifted [2,3): the insertion gap + assert_eq!(b.lift_to_reference(5, 6), (None, None)); + assert_eq!(b.lift_to_query(107, 111), (Some(9), Some(13))); + } } diff --git a/molecular-annotation/src/tests.rs b/molecular-annotation/src/tests.rs index 922c267d5..83e0a4c41 100644 --- a/molecular-annotation/src/tests.rs +++ b/molecular-annotation/src/tests.rs @@ -607,6 +607,91 @@ fn test_get_ref_coords_reverse_outside_aligned_region() { assert_eq!(coords[0].3, None); } +// A full-read annotation frame on a hard-clipped record (#136): query +// coordinates stay in that frame, only the lift and project_query subtract +// the leading hard clip. +#[test] +fn test_query_offset_forward_lifts_through_leading_hard_clip() { + // 150 bp read, CIGAR 50H100M: SEQ is molecular [50,150), ref [1000,1100) + let mut a = MolecularAnnotations::new(150); + a.add_annotation_type("nuc", QualitySpec::none(), Encoding::Ma) + .add(60, 20, Strand::Unknown, vec![], None) // [60,80): inside SEQ + .add(10, 20, Strand::Unknown, vec![], None) // [10,30): entirely in the clip + .add(40, 20, Strand::Unknown, vec![], None); // [40,60): straddles the clip + a.set_aligned_blocks(vec![([0, 100], [1000, 1100])], false); + assert_eq!(a.query_offset(), 0); + a.set_query_offset(50); + assert_eq!(a.query_offset(), 50); + let c = a.get_ref_coords("nuc").unwrap(); + assert_eq!( + c[0], + (60, 80, Some(1010), Some(1030)), + "query coords stay full-frame, ref uses the offset" + ); + assert_eq!( + c[1], + (10, 30, None, None), + "annotation before SEQ does not lift and does not panic" + ); + assert_eq!( + c[2], + (40, 60, Some(1000), Some(1010)), + "straddling annotation snaps like a soft clip" + ); + let infos: Vec<_> = a.iter_type("nuc").unwrap().collect(); + assert_eq!( + ( + infos[0].query_start, + infos[0].query_end, + infos[0].ref_start, + infos[0].ref_end + ), + (60, 80, Some(1010), Some(1030)) + ); + let pq: Vec<(i64, i64)> = a + .project_query(0, false) + .map(|p| (p.start, p.end)) + .collect(); + assert_eq!( + pq, + vec![(10, 30), (-40, -20), (-10, 10)], + "project_query is SEQ-relative" + ); + let pr: Vec<(i64, i64)> = a + .project_reference(1000, false) + .map(|p| (p.start, p.end)) + .collect(); + assert_eq!(pr, vec![(10, 30), (0, 10)]); + // the container's raw lift is in the annotation frame too + assert_eq!(a.lift_to_reference(60, 80), Some((Some(1010), Some(1030)))); + assert_eq!(a.lift_to_query(1010, 1030), Some((Some(60), Some(80)))); +} + +#[test] +fn test_query_offset_reverse_aligned() { + // molecular [60,80) of a 150 bp reverse read -> BAM [70,90) + let mut a = MolecularAnnotations::new(150); + a.add_annotation_type("nuc", QualitySpec::none(), Encoding::Ma) + .add(60, 20, Strand::Unknown, vec![], None); + a.set_aligned_blocks(vec![([0, 100], [1000, 1100])], true); + // CIGAR 100M50H (trailing clip): no offset, today's result + assert_eq!( + a.get_ref_coords("nuc").unwrap(), + vec![(70, 90, Some(1070), Some(1090))] + ); + // CIGAR 50H100M on a reverse read: BAM [70,90) is SEQ [20,40) + a.set_query_offset(50); + assert_eq!( + a.get_ref_coords("nuc").unwrap(), + vec![(70, 90, Some(1020), Some(1040))] + ); + let pq: Vec<(i64, i64)> = a + .project_query(0, false) + .map(|p| (p.start, p.end)) + .collect(); + assert_eq!(pq, vec![(20, 40)]); +} + #[test] fn test_flip_range() { let annotations = MolecularAnnotations::new(1000); @@ -2268,3 +2353,134 @@ fn ma_family_tags_accept_both_spellings() { } assert!(matches!(record.aux(b"Ma"), Ok(Aux::String(_)))); } + +/// A mapped record at reference position 1000 with the given SEQ, CIGAR and +/// flags, for the hard-clip frame tests below. +#[cfg(feature = "htslib")] +fn aligned_record(seq: &[u8], cigar: &str, flags: u16) -> rust_htslib::bam::Record { + use rust_htslib::bam::record::CigarString; + let mut r = rust_htslib::bam::Record::new(); + let cigar = CigarString::try_from(cigar).unwrap(); + r.set(b"read", Some(&cigar), seq, &vec![255u8; seq.len()]); + r.set_flags(flags); + r.set_tid(0); + r.set_pos(1000); + r +} + +#[cfg(feature = "htslib")] +fn with_ma_and_mm(seq: &[u8], cigar: &str, flags: u16, ma: &str) -> rust_htslib::bam::Record { + use rust_htslib::bam::record::Aux; + let mut r = aligned_record(seq, cigar, flags); + r.push_aux(b"Ma", Aux::String(ma)).unwrap(); + r.push_aux(b"MM", Aux::String("A+a,0;")).unwrap(); + r.push_aux(b"ML", Aux::ArrayU8((&[200u8][..]).into())) + .unwrap(); + r +} + +// MA read length == SEQ + hard clips: the tags describe the full read. The +// lift carries the leading hard clip; MM/ML are not parsed. +#[cfg(feature = "htslib")] +#[test] +fn from_record_full_read_frame_forward() { + let seq = vec![b'A'; 80]; + // molecular [20,50) of a 100 bp read; SEQ is [10,90) + let r = with_ma_and_mm(&seq, "10H80M10H", 0, "100;msp.:21-30"); + let annot = MolecularAnnotations::from_record(&r); + assert_eq!(annot.read_length, 100); + assert_eq!(annot.query_offset(), 10); + assert_eq!( + annot.get_ref_coords("msp"), + Some(vec![(20, 50, Some(1010), Some(1040))]) + ); + assert!(annot.get_type("a").is_none(), "MM/ML must not be parsed"); + assert!(annot.to_ma_string().starts_with("100;"), "frame kept as is"); +} + +#[cfg(feature = "htslib")] +#[test] +fn from_record_full_read_frame_reverse() { + let seq = vec![b'A'; 80]; + // molecular [20,50) -> flip with L=100 -> BAM [50,80) -> SEQ [45,75) + let r = with_ma_and_mm(&seq, "5H80M15H", 16, "100;msp.:21-30"); + let annot = MolecularAnnotations::from_record(&r); + assert_eq!(annot.query_offset(), 5); + assert_eq!(annot.get_coords("msp"), Some(vec![(50, 80)])); + assert_eq!( + annot.get_ref_coords("msp"), + Some(vec![(50, 80, Some(1045), Some(1075))]) + ); + // molecular [0,10) -> BAM [90,100) -> SEQ [85,95): past the 80 bp SEQ + let r = with_ma_and_mm(&seq, "5H80M15H", 16, "100;msp.:1-10"); + let annot = MolecularAnnotations::from_record(&r); + assert_eq!(annot.get_ref_coords("msp").unwrap()[0].2, None); + // trailing clip only: still the full-read frame, offset 0 + let r = with_ma_and_mm(&seq, "80M20H", 16, "100;msp.:31-10"); + let annot = MolecularAnnotations::from_record(&r); + assert_eq!(annot.query_offset(), 0); + assert_eq!( + annot.get_ref_coords("msp"), + Some(vec![(60, 70, Some(1060), Some(1070))]) + ); + assert!(annot.get_type("a").is_none()); +} + +// Every other shape takes the old path: no offset, MM/ML parsed. +#[cfg(feature = "htslib")] +#[test] +fn from_record_seq_frame_is_unchanged() { + use rust_htslib::bam::record::Aux; + let seq = vec![b'A'; 80]; + // tags computed after clipping + let r = with_ma_and_mm(&seq, "10H80M10H", 0, "80;msp.:21-30"); + let annot = MolecularAnnotations::from_record(&r); + assert_eq!(annot.query_offset(), 0); + assert_eq!( + annot.get_ref_coords("msp"), + Some(vec![(20, 50, Some(1020), Some(1050))]) + ); + assert!(annot.get_type("a").is_some()); + // soft clip + let r = with_ma_and_mm(&seq, "10S70M", 0, "80;msp.:21-30"); + let annot = MolecularAnnotations::from_record(&r); + assert_eq!(annot.query_offset(), 0); + assert_eq!( + annot.get_ref_coords("msp"), + Some(vec![(20, 50, Some(1010), Some(1040))]) + ); + // no MA tag with hard clips + let mut r = aligned_record(&seq, "10H80M10H", 0); + r.push_aux(b"MM", Aux::String("A+a,0;")).unwrap(); + r.push_aux(b"ML", Aux::ArrayU8((&[200u8][..]).into())) + .unwrap(); + let annot = MolecularAnnotations::from_record(&r); + assert_eq!((annot.read_length, annot.query_offset()), (80, 0)); + assert!(annot.get_type("a").is_some()); + // a read length matching neither: the library leaves it to the caller + let r = with_ma_and_mm(&seq, "10H80M10H", 0, "90;msp.:21-30"); + assert_eq!(MolecularAnnotations::from_record(&r).query_offset(), 0); +} + +#[cfg(feature = "htslib")] +#[test] +fn hard_clip_frame_helpers() { + use crate::{full_read_query_offset, hard_clips, query_span}; + use rust_htslib::bam::record::CigarString; + let seq = vec![b'A'; 80]; + let r = aligned_record(&seq, "10H80M10H", 0); + assert_eq!(hard_clips(&r), (10, 10)); + assert_eq!(query_span(&r), 80); + assert_eq!(full_read_query_offset(100, &r), Some(10)); + assert_eq!(full_read_query_offset(80, &r), None); + assert_eq!(full_read_query_offset(90, &r), None); + // SEQ-less: the CIGAR query span stands in for SEQ + let mut seqless = rust_htslib::bam::Record::new(); + let cigar = CigarString::try_from("10H70M10S10H").unwrap(); + seqless.set(b"read", Some(&cigar), b"", &[]); + assert_eq!(query_span(&seqless), 80); + assert_eq!(full_read_query_offset(100, &seqless), Some(10)); + let r = aligned_record(&seq, "80M", 0); + assert_eq!(hard_clips(&r), (0, 0)); + assert_eq!(full_read_query_offset(80, &r), None); +} diff --git a/src/fiber.rs b/src/fiber.rs index 71a59b3ba..e228f11ef 100644 --- a/src/fiber.rs +++ b/src/fiber.rs @@ -185,26 +185,41 @@ impl FiberseqData { pub fn callable_state(&self) -> (CallableState, i64, i64) { // A tag whose recorded read length no longer matches the record is // stale: something rewrote the read after tagging. Treat as Untagged. - // A SEQ-less record (seq_len 0, SEQ dropped to save space) is NOT - // stale: the MA read length is the source of truth there. - if crate::utils::ma_io::read_length_is_stale( - self.annotations.read_length, - self.record.seq_len(), - ) { + // The frame is SEQ, or SEQ plus the hard clips for a full-read-frame + // record. A SEQ-less record (seq_len 0, SEQ dropped to save space) + // is NOT stale: the MA read length is the source of truth there. + if crate::utils::ma_io::frame_is_stale(&self.annotations, &self.record) { return (CallableState::Untagged, 0, 0); } + // Decide the 1-0 NotCallable marker on the raw model, before any + // liftover: the marker carries no position, and on a full-read-frame + // record position 0 sits inside the leading hard clip. + let Some(t) = self.annotations.get_type(FIBERSEQ_CALLABLE_TYPE) else { + return (CallableState::Untagged, 0, 0); + }; + let Some(raw) = t.annotations.first() else { + return (CallableState::Untagged, 0, 0); + }; + if raw.length == 0 { + return (CallableState::NotCallable, 0, 0); + } + // A hard-clipped part of a tagged read has no m6A (MM/ML were + // dropped), so it is never callable, whatever its inherited tag says. + if self.is_full_read_frame() { + return (CallableState::NotCallable, 0, 0); + } let view = self.fiberseq_callable(); let infos = view.infos(); let Some(a) = infos.first() else { - return (CallableState::Untagged, 0, 0); + return (CallableState::NotCallable, 0, 0); }; let len = self.frame_length() as i64; let (cs, ce) = ( (a.query_start as i64).clamp(0, len), (a.query_end as i64).clamp(0, len), ); - // Empty after clamping covers both the 1-0 NotCallable marker and - // a corrupt foreign span lying outside the frame. + // Empty after clamping covers a corrupt foreign span lying outside + // the frame. if ce <= cs { return (CallableState::NotCallable, cs, cs); } @@ -224,6 +239,14 @@ impl FiberseqData { } } + /// True when the annotations describe the full-length read of which + /// this record's SEQ is a hard-clipped part (#136). Query coordinates + /// from the views are then in that frame, SEQ starts + /// `annotations.query_offset()` bases in, and the read has no m6A. + pub fn is_full_read_frame(&self) -> bool { + crate::utils::ma_io::model_is_full_read_frame(&self.annotations, &self.record) + } + /// True when the read carries a non-empty `fiberseq_callable` span. pub fn is_callable(&self) -> bool { self.callable_state().0 == CallableState::Callable @@ -386,7 +409,13 @@ impl FiberseqData { } else { ct = &name; start = 0; - end = self.record.seq_len() as i64; + // The blocks are in the annotation frame: SEQ, or the full read + // this record was hard-clipped from. + end = if self.is_full_read_frame() { + self.annotations.read_length as i64 + } else { + self.record.seq_len() as i64 + }; } let score = self.ec.round() as i64; let strand = if self.record.is_reverse() { '-' } else { '+' }; @@ -489,7 +518,13 @@ impl FiberseqData { // PB features let name = std::str::from_utf8(self.record.qname()).unwrap(); let score = self.ec.round() as i64; - let q_len = self.record.seq_len() as i64; + // The molecular columns are in the annotation frame, which on a + // hard-clipped full-read-frame record is the whole read, not SEQ. + let q_len = if self.is_full_read_frame() { + self.annotations.read_length as i64 + } else { + self.record.seq_len() as i64 + }; let rq = match self.get_rq() { Some(x) => format!("{x}"), None => ".".to_string(), diff --git a/src/main.rs b/src/main.rs index c62ffa0d4..50626d451 100644 --- a/src/main.rs +++ b/src/main.rs @@ -146,6 +146,8 @@ pub fn main() -> Result<(), Error> { } None => {} }; + utils::ma_io::report_stale_frames(); + utils::ma_io::report_full_frames(); let duration = pg_start.elapsed(); log::info!( "{} done! Time elapsed: {}", diff --git a/src/subcommands/ddda_to_m6a.rs b/src/subcommands/ddda_to_m6a.rs index 5b2e22d2c..6a39521f2 100644 --- a/src/subcommands/ddda_to_m6a.rs +++ b/src/subcommands/ddda_to_m6a.rs @@ -71,6 +71,17 @@ pub fn ddda_to_m6a_record(record: &mut Record, opts: &DddaToM6aOptions) { ); MolecularAnnotations::from_record(record) }); + if ma_io::model_is_full_read_frame(&annot, record) { + // The Y/R calls below are positions in SEQ; the inherited nuc/msp + // live in the full read's frame and cannot share a tag with them. + // Restart from SEQ. + annot.annotation_types.clear(); + annot.read_length = record.seq_len() as u32; + annot.set_aligned_blocks_raw( + molecular_annotation::AlignedBlocks::from_record(record), + record.is_reverse(), + ); + } // Verdict from the calling-time set, before the m6A rebuild. ma_io::sync_fiberseq_callable(&mut annot, record, &opts.input.filters); annot diff --git a/src/subcommands/fire.rs b/src/subcommands/fire.rs index 202b7295f..0f26d68e1 100644 --- a/src/subcommands/fire.rs +++ b/src/subcommands/fire.rs @@ -9,12 +9,31 @@ use itertools::Itertools; use rayon::prelude::*; use utils::fire::*; +/// True when the record carries at least one m6A call in memory. A record +/// with msp but no m6A (a hard-clipped full-read-frame record whose MM/ML +/// were dropped, or one stripped on purpose) cannot be scored: FireFeats +/// indexes SEQ by m6A and MSP coordinates and needs both. +fn has_m6a(rec: &FiberseqData) -> bool { + rec.annotations + .get_type(crate::utils::basemods::M6A_TYPE) + .is_some_and(|t| !t.annotations.is_empty()) +} + pub fn add_fire_to_rec( rec: &mut FiberseqData, fire_opts: &FireOptions, model: &GBDT, precision_table: &MapPrecisionValues, ) { + if !has_m6a(rec) { + // No m6A, so no FIRE features. Write the model back unchanged: + // nuc/msp (and any fire calls made while the read still had m6A) + // are kept; the callable state synced at read time (NotCallable + // for a full-read-frame record) is written with them. + log::debug!("FIRE: no m6A on {}; writing it unscored", rec.get_qname()); + rec.serialize_annotations(); + return; + } let fire_feats = FireFeats::new(rec, fire_opts); let mut precisions = fire_feats.predict_with_xgb(model, precision_table); // FIRE produces precisions in MSP-iteration (BAM) order. Convert to @@ -30,7 +49,11 @@ pub fn add_fire_to_rec( // and their paired precisions, keeping only entries with p > 0. let (fire_starts, fire_lens, fire_quals): (Vec, Vec, Vec) = { let Some(msp) = rec.annotations.get_type(ma_io::MSP_TYPE) else { - log::warn!("FIRE: no msp annotations on record; skipping"); + // Nothing to score. Still write the model back so a record whose + // stale tags were dropped by the reader leaves as an honest + // untagged read instead of passing those tags on. + log::debug!("FIRE: no msp annotations on record; writing it unscored"); + rec.serialize_annotations(); return; }; if msp.annotations.len() != precisions.len() { @@ -98,6 +121,7 @@ pub fn add_fire_to_bam(fire_opts: &mut FireOptions) -> Result<(), anyhow::Error> let chunk: Vec = chunk.collect(); let feats: Vec = chunk .par_iter() + .filter(|r| has_m6a(r)) .map(|r| FireFeats::new(r, fire_opts)) .collect(); feats.iter().for_each(|f| { diff --git a/src/subcommands/footprint.rs b/src/subcommands/footprint.rs index a277b6fda..170a6ae1b 100644 --- a/src/subcommands/footprint.rs +++ b/src/subcommands/footprint.rs @@ -127,8 +127,12 @@ impl<'a> Footprint<'a> { pub fn new(motif: &'a ReferenceMotif, in_fibers: &'a Vec) -> Self { let mut fibers = vec![]; for fiber in in_fibers { - // add if fiber spans the footprint - if motif.spans(fiber.record.reference_start(), fiber.record.reference_end()) { + // add if fiber spans the footprint; a full-read-frame record has + // no m6A, so it was never measured and must not count as a + // fully footprinted fiber + if !fiber.is_full_read_frame() + && motif.spans(fiber.record.reference_start(), fiber.record.reference_end()) + { fibers.push(fiber); } } diff --git a/src/subcommands/predict_m6a.rs b/src/subcommands/predict_m6a.rs index acfa276a5..40b37018e 100644 --- a/src/subcommands/predict_m6a.rs +++ b/src/subcommands/predict_m6a.rs @@ -250,6 +250,15 @@ where && t.name != ma_io::MSP_TYPE && t.name != ma_io::FIRE_TYPE }); + if ma_io::model_is_full_read_frame(&annot, record) { + // Full-read tags cannot be re-derived without kinetics: + // leave an honest NotCallable record in the SEQ frame. + annot.annotation_types.clear(); + annot.set_aligned_blocks_raw( + molecular_annotation::AlignedBlocks::from_record(record), + record.is_reverse(), + ); + } // Sync the frame or the marker reads back as Untagged. annot.read_length = record.seq_len() as u32; ma_io::set_fiberseq_callable( diff --git a/src/subcommands/qc.rs b/src/subcommands/qc.rs index 8b90914ee..237370d5c 100644 --- a/src/subcommands/qc.rs +++ b/src/subcommands/qc.rs @@ -121,6 +121,9 @@ pub struct QcStats<'a> { // reads whose tag was stale (read_length mismatch); folded into // Untagged, tracked for a distinct warning stale_tags: i64, + // hard-clipped reads whose tags are in the full-read frame: nuc/msp + // kept, m6A dropped, counted as NotCallable; tracked for a warning + full_frame_reads: i64, // phasing information phased_reads: HashMap, phased_bp: HashMap, @@ -156,6 +159,7 @@ impl<'a> QcStats<'a> { .map(|k| (k, Counts::default())) .collect(), stale_tags: 0, + full_frame_reads: 0, qc_opts, phased_reads: HashMap::new(), phased_bp: HashMap::new(), @@ -183,13 +187,13 @@ impl<'a> QcStats<'a> { }; bump(&mut self.fiberseq_callable, state_key, 1, f.inc()); if state == fiber::CallableState::Untagged - && crate::utils::ma_io::read_length_is_stale( - fiber.annotations.read_length, - fiber.record.seq_len(), - ) + && crate::utils::ma_io::record_frame_reason(&fiber.record).is_some() { self.stale_tags += 1; } + if fiber.is_full_read_frame() { + self.full_frame_reads += 1; + } // add auto-correlation of m6a self.add_m6a_starts_for_acf(fiber, passes); @@ -224,7 +228,9 @@ impl<'a> QcStats<'a> { // add the m6a to the working queue let mut m6a_vec: Vec = vec![0.0; fiber.record.seq_len()]; for m6a in fiber.m6a().starts().iter() { - m6a_vec[*m6a as usize] = 1.0; + if let Some(v) = m6a_vec.get_mut(*m6a as usize) { + *v = 1.0; + } } let elem = AcfRead { m6a: m6a_vec, @@ -267,7 +273,14 @@ impl<'a> QcStats<'a> { // have zero nucleosomes (the callable state needs m6A and MSPs, and the // nuc view is post-pruning), so the filtered side skips those // reads rather than admit an inf key. - let read_length = fiber.frame_length() as f32 / nuc.len() as f32; + // On a full-read-frame record the nucleosomes span the whole read, + // so divide that length, not SEQ. + let frame_len = if fiber.is_full_read_frame() { + fiber.annotations.read_length as f32 + } else { + fiber.frame_length() as f32 + }; + let read_length = frame_len / nuc.len() as f32; bump( &mut self.read_length_per_nuc, ordered_float_10k_round(read_length), @@ -540,18 +553,28 @@ pub fn run_qc(opts: &mut QcOpts) -> Result<(), anyhow::Error> { if untagged > 0 { log::warn!( "{untagged} reads have no fiberseq_callable state (no nuc/msp \ - calls, no SEQ to derive from, or a stale tag). These reads never enter \ - count_filtered. Run ft add-nucleosomes or ft predict-m6a to \ - call them." + calls, no SEQ to derive from, or tags from a hard-clipped alignment). \ + These reads never enter count_filtered. Reads that were never called: \ + run ft add-nucleosomes or ft predict-m6a." ); } if stats.stale_tags > 0 { log::warn!( - "{} reads carry a stale fiberseq_callable tag: the recorded read \ - length does not match the record. These reads count as Untagged.", + "{} of them carry Fiber-seq tags from a longer read than their SEQ \ + (hard-clipped supplementary alignments); their calls were dropped. \ + See the warning above for how to realign.", stats.stale_tags ); } + if stats.full_frame_reads > 0 { + log::warn!( + "{} reads are hard-clipped alignments whose nuc/msp tags are in the frame of the \ + full-length read. Their nuc/msp are counted, but their m6A (MM/ML) was dropped, so \ + they are NotCallable and never enter count_filtered. See the warning above \ + for how to realign.", + stats.full_frame_reads + ); + } let mut out = bio_io::writer(&opts.out)?; stats.write(&mut out)?; stats.write_m6a_acf(&mut out)?; diff --git a/src/utils/bio_io.rs b/src/utils/bio_io.rs index 875a36da2..c201da571 100644 --- a/src/utils/bio_io.rs +++ b/src/utils/bio_io.rs @@ -295,16 +295,9 @@ where for r in self.bam.by_ref().take(self.chunk_size) { pulled += 1; let r = r.unwrap(); - let has_mm_and_ml = r.aux(b"MM").is_ok() && r.aux(b"ML").is_ok(); - if has_mm_and_ml - && (r.cigar().leading_hardclips() > 0 || r.cigar().trailing_hardclips() > 0) - { - log::warn!( - "Skipping read ({}) because it has been hard clipped and has ML and MM tags. This read will be excluded from calculations and any output.", - String::from_utf8_lossy(r.qname()) - ); - continue; - } + // Hard-clipped records with stale tags are not skipped here: + // ma_io::read_record clears their annotations and the + // writers strip the tags, so every command sees one policy. // filter by bit flag if r.flags() & self.bit_flag_filter != 0 { continue; diff --git a/src/utils/input_bam.rs b/src/utils/input_bam.rs index 302a6996f..6f902f1d5 100644 --- a/src/utils/input_bam.rs +++ b/src/utils/input_bam.rs @@ -272,6 +272,8 @@ impl FiberFilters { #[derive(Debug, Args)] pub struct InputBam { /// Input BAM file. If no path is provided stdin is used. For m6A prediction, this should be a HiFi bam file with kinetics data. For other commands, this should be a bam file with m6A calls. + /// + /// Fiber-seq tags describe the full read, so aligned input must be soft-clipped: pbmm2 never hard-clips; for ONT use `dorado aligner --mm2-opts "-Y"` or `minimap2 -Y -y`. On hard-clipped records the nucleosome/MSP/FIRE tags are kept and lifted; the m6A cannot be recovered and is dropped, so such reads are not fiberseq-callable. #[clap(default_value = "-", value_hint = ValueHint::AnyPath)] pub bam: String, #[clap(flatten)] diff --git a/src/utils/ma_io.rs b/src/utils/ma_io.rs index 4f476abb7..308896f6e 100644 --- a/src/utils/ma_io.rs +++ b/src/utils/ma_io.rs @@ -20,6 +20,11 @@ //! (NotCallable); only never-processed reads have none. Non-default //! minimums name the annotation in the AN tag (e.g. "m20a10"). //! +//! Hard-clipped records whose tags describe the full-length read (an aligner +//! copied the primary's tags onto a supplementary; see [`record_frame`]) keep +//! their nuc/msp/fire in that frame, lifted with the leading hard clip, and +//! lose their MM/ML. They are always NotCallable: they have no m6A. +//! //! `m6a` and `cpg` types may appear *in memory* on a [`MolecularAnnotations`] //! populated by the library's MM/ML parser. Their on-disk source of truth //! is `MM`/`ML`; they carry `Encoding::MmMl` (set at construction, whether read @@ -27,7 +32,10 @@ //! library serializes them into MM/ML rather than the MA tag set. use anyhow::{bail, Result}; -use molecular_annotation::{ma_family_tags, Encoding, MolecularAnnotations, QualitySpec, Strand}; +use molecular_annotation::{ + hard_clips, ma_family_tags, query_span, AlignedBlocks, Encoding, MolecularAnnotations, + QualitySpec, Strand, +}; use rust_htslib::bam::{self, record::Aux}; /// Annotation type names used by fibertools-rs. @@ -49,68 +57,468 @@ type MspInput<'a> = (&'a [u32], &'a [u32], Option<&'a [u8]>); /// /// Delegates to the library's combined MA-spec + MM/ML parser. If the /// library returns an empty annotation set and the record carries legacy -/// `ns`/`nl`/`as`/`al`/`aq` tags, falls back to `read_legacy_nuc_msp`. +/// `ns`/`nl`/`as`/`al`/`aq` tags, falls back to the legacy reader. /// Tolerant of malformed MM/ML — the library handles those internally /// without panicking. +/// +/// The record's frame ([`record_frame`]) decides what is parsed: a stale +/// record comes back empty, a full-read-frame record keeps its MA/legacy +/// annotations under the full read length (lifted with the leading hard +/// clip) and never parses MM/ML, and everything else parses as before. pub fn read_record(record: &bam::Record) -> Result { - let mut annot = MolecularAnnotations::from_record(record); - // If MA tag is absent, also ingest legacy nuc/msp tags. The library - // already populates basemod types (m6a/cpg) from MM/ML, so we merge - // legacy-derived nuc/msp into whatever the library produced rather - // than gating on `annotation_types.is_empty()` (which would skip the - // legacy fallback whenever MM/ML is present). - let has_ma = ma_family_tags(record).is_some(); - if !has_ma { - // Provenance: the gates are type-checked (`has_legacy_nuc_msp`, - // `has_legacy_fibertig`), so a foreign tool reusing these two-letter - // names with a different aux type is never parsed as fibertools data. - // The write path (`strip_consumed_legacy_tags`) removes legacy tags - // under pair-level gates mirroring `read_legacy_nuc_msp`'s - // consumption — i.e. only what this reader ingested — so keep the - // two in sync. - if has_legacy_nuc_msp(record) { - merge_missing_types(&mut annot, read_legacy_nuc_msp(record)?); + let mut annot = match record_frame(record) { + // A record whose tags describe another read is untagged, and parsing + // its MM/ML would only produce truncation noise: decide before parsing. + Frame::Stale(why) => { + let mut annot = MolecularAnnotations::new(0); + drop_stale_frame(&mut annot, record, why); + return Ok(annot); + } + Frame::FullRead { + h_lead, + read_length, + } => { + let annot = if ma_family_tags(record).is_some() { + // The library sets the query offset from the MA read length + // and skips MM/ML itself. + MolecularAnnotations::from_record(record) + } else { + // Legacy arrays record no frame, and every producer wrote + // them on the full read: build that frame directly so MM/ML + // are never decoded against the wrong SEQ. + let mut annot = MolecularAnnotations::new(read_length); + annot.set_aligned_blocks_raw( + AlignedBlocks::from_record(record).with_query_offset(h_lead), + record.is_reverse(), + ); + read_legacy_nuc_msp_into(record, &mut annot)?; + read_legacy_fibertig_into(record, &mut annot)?; + annot + }; + annot } - // Legacy fibertig (`fs`/`fl`/`fa`): the pre-MA fibertig wire format, - // dropped from the writer in favour of the MA-spec `AN` tag. Kept - // readable so older fibertig BAMs stay consumable by `extract`. - if has_legacy_fibertig(record) { - merge_missing_types(&mut annot, read_legacy_fibertig(record)?); + Frame::Seq => { + let mut annot = MolecularAnnotations::from_record(record); + // The MA tag vouches for the frame, but an MN tag that disagrees + // with SEQ says the base mods were copied from a longer read: keep + // the calls, drop only the m6A/CpG (the writer strips MM/ML/MN). + if mn_disagrees(record) { + annot.annotation_types.retain(|t| !t.is_mm_ml()); + } + // If MA tag is absent, also ingest legacy nuc/msp tags. The library + // already populates basemod types (m6a/cpg) from MM/ML, so we add + // legacy-derived nuc/msp to whatever the library produced rather + // than gating on `annotation_types.is_empty()` (which would skip + // the legacy fallback whenever MM/ML is present). + if ma_family_tags(record).is_none() { + // Provenance: the gates are type-checked (`has_legacy_nuc_msp`, + // `has_legacy_fibertig`), so a foreign tool reusing these + // two-letter names with a different aux type is never parsed + // as fibertools data. The write path + // (`strip_consumed_legacy_tags`) removes legacy tags under + // pair-level gates mirroring `read_legacy_nuc_msp_into`'s + // consumption, i.e. only what this reader ingested, so keep + // the two in sync. + if has_legacy_nuc_msp(record) { + read_legacy_nuc_msp_into(record, &mut annot)?; + } + // Legacy fibertig (`fs`/`fl`/`fa`): the pre-MA fibertig wire + // format, dropped from the writer in favour of the MA-spec + // `AN` tag. Kept readable so older fibertig BAMs stay + // consumable by `extract`. + if has_legacy_fibertig(record) { + read_legacy_fibertig_into(record, &mut annot)?; + } + } + annot + } + }; + // Backstop on the parsed model: a read length that disagrees with the + // frame, or any annotation ending past it. + if let Some(why) = model_frame_reason(&annot, record) { + drop_stale_frame(&mut annot, record, why); + } else if model_is_full_read_frame(&annot, record) { + // Warn and count only when this pass really drops m6A. The writer + // strips MM/ML, so a second run over ft's own output stays quiet. + if matches!(record.aux(b"MM"), Ok(Aux::String(_))) { + note_full_read_frame(record); + } else { + log::debug!( + "full-read frame for {} (no m6A to drop)", + String::from_utf8_lossy(record.qname()) + ); } } Ok(annot) } +/// True when an MN tag is present and names a different length than SEQ: +/// the MM/ML next to it were written for another read. +fn mn_disagrees(record: &bam::Record) -> bool { + let seq_len = record.seq_len(); + seq_len > 0 + && matches!(record.aux(b"MM"), Ok(Aux::String(_))) + && record + .aux(b"MN") + .ok() + .and_then(aux_as_usize) + .is_some_and(|mn| mn != seq_len) +} + +/// How many stale-frame records get a WARN before the rest drop to DEBUG. +/// ONT BAMs can hold thousands of hard-clipped supplementary reads. A total +/// is printed at exit by [`report_stale_frames`]. +const STALE_FRAME_WARN_LIMIT: usize = 10; + +/// Records whose annotations were dropped because their tags did not fit SEQ. +static STALE_FRAMES: std::sync::atomic::AtomicUsize = std::sync::atomic::AtomicUsize::new(0); + +/// What to tell a user who hit a stale frame. Aligners that hard-clip +/// supplementary alignments (minimap2 and dorado aligner without -Y) copy the +/// full-length read's tags onto the clipped record; the tags cannot be +/// recovered, only avoided. +pub const HARD_CLIP_REMEDY: &str = "Realign with soft clipping: pbmm2 align (PacBio; it never \ + hard-clips), dorado aligner --mm2-opts \"-Y\", or minimap2 -Y -y. Or drop supplementary \ + alignments with -F 2048."; + +/// Treat a record whose tags come from another frame as untagged (#136, #31). +/// Hard-clipped supplementary alignments keep the full-length read's tags, +/// which would index past SEQ: consumers panic, or emit misplaced and +/// u32-wrapped coordinates. Lives here so every path that parses a record +/// (the fiber reader, convert-tags, strip-basemods, ddda-to-m6a, predict-m6a, +/// fibertig) sees the same thing, and [`write_record`] strips the stale tags +/// on the way out so nothing downstream trusts them. +/// +/// Not reframed on purpose: nuc/msp could be shifted by the hard clip, but +/// MM/ML cannot be recovered without the clipped bases, and half a record +/// (nucleosomes, no m6A) is worse than an honest untagged one. The spec's +/// answer is soft clipping (`HARD_CLIP_REMEDY`). +fn drop_stale_frame(annot: &mut MolecularAnnotations, record: &bam::Record, why: String) { + use std::sync::atomic::Ordering; + let n = STALE_FRAMES.fetch_add(1, Ordering::Relaxed); + let qname = String::from_utf8_lossy(record.qname()); + if n == 0 { + // One message, so the cause never prints after the symptom when + // several threads hit this at once. + log::warn!( + "some records carry Fiber-seq tags that describe the full-length read, not their SEQ \ + (hard-clipped supplementary alignments). Their annotations are dropped and they are \ + treated as untagged reads. {HARD_CLIP_REMEDY}\n\ + dropping annotations for {qname}: {why}" + ); + } else if n < STALE_FRAME_WARN_LIMIT { + log::warn!("dropping annotations for {qname}: {why}"); + if n + 1 == STALE_FRAME_WARN_LIMIT { + log::warn!( + "further such records are logged at debug level; a total is printed at exit" + ); + } + } else { + log::debug!("dropping annotations for {qname}: {why}"); + } + annot.annotation_types.clear(); + // An untagged read's frame is its SEQ, or its CIGAR query span without SEQ, + // so the writers persist a clean `Ma:Z:` instead of the stale length. + // The liftover follows the frame (no query offset). + annot.read_length = query_span(record); + annot.set_aligned_blocks_raw(AlignedBlocks::from_record(record), record.is_reverse()); +} + +/// Records whose Fiber-seq tags describe the full-length read while they +/// hold a hard-clipped part of it (`Frame::FullRead`). +static FULL_FRAMES: std::sync::atomic::AtomicUsize = std::sync::atomic::AtomicUsize::new(0); + +/// What happens to the base mods of a full-read-frame record. Part of the +/// once-per-run warning, so a test can look for it. +pub const FULL_FRAME_M6A_DROPPED: &str = + "their m6A (MM/ML) describes bases this record does not carry and was dropped"; + +/// Warn once per run that a record is in the full-read frame; later records +/// are logged at debug level and totalled at exit by [`report_full_frames`]. +fn note_full_read_frame(record: &bam::Record) { + use std::sync::atomic::Ordering; + let n = FULL_FRAMES.fetch_add(1, Ordering::Relaxed); + let qname = String::from_utf8_lossy(record.qname()); + if n == 0 { + log::warn!( + "some hard-clipped records carry Fiber-seq tags computed on the full-length read \ + (supplementary alignments; the aligner copied them from the primary). Their \ + nucleosome/MSP/FIRE calls are kept and lifted with the hard-clip offset; \ + {FULL_FRAME_M6A_DROPPED}, so these reads are not fiberseq-callable. \ + {HARD_CLIP_REMEDY}\n\ + full-read frame for {qname}" + ); + } else { + log::debug!("full-read frame for {qname}"); + } +} + +/// Log the number of full-read-frame records. Called once at exit by +/// `main`, next to [`report_stale_frames`]. +pub fn report_full_frames() { + let n = FULL_FRAMES.load(std::sync::atomic::Ordering::Relaxed); + if n > 0 { + let records = if n == 1 { "record" } else { "records" }; + log::warn!( + "kept nucleosome/MSP calls on {n} hard-clipped {records} whose Fiber-seq tags \ + describe the full-length read; they have no m6A and are not fiberseq-callable. \ + {HARD_CLIP_REMEDY}" + ); + } +} + +/// Log the number of records whose annotations were dropped for a stale +/// frame. Called once at exit by `main`. +pub fn report_stale_frames() { + let n = STALE_FRAMES.load(std::sync::atomic::Ordering::Relaxed); + if n > 0 { + let records = if n == 1 { "record" } else { "records" }; + log::warn!( + "dropped annotations on {n} {records} whose Fiber-seq tags did not match their SEQ \ + (hard-clipped supplementary alignments). {HARD_CLIP_REMEDY}" + ); + } +} + +/// True when the model is in the full-read frame of a hard-clipped record: +/// SEQ is present, the record has hard clips, and the model's read length is +/// SEQ plus both clips. Such a model has no m6A (the reader never parses +/// MM/ML in this frame) and its query coordinates run past SEQ; the lift +/// carries the leading hard clip as `annot.query_offset()`. +pub(crate) fn model_is_full_read_frame(annot: &MolecularAnnotations, record: &bam::Record) -> bool { + let span = query_span(record) as usize; + if span == 0 { + return false; + } + let (lead, trail) = hard_clips(record); + lead + trail > 0 && annot.read_length as usize == span + lead as usize + trail as usize +} + +/// True when a present SEQ disagrees with the model's frame: the recorded +/// read length must equal SEQ, or SEQ plus the hard clips for a full-read +/// frame (something else rewrote the read after tagging). SEQ-less records +/// are never stale: the MA read length is the frame. +pub(crate) fn frame_is_stale(annot: &MolecularAnnotations, record: &bam::Record) -> bool { + let seq_len = record.seq_len(); + seq_len > 0 && annot.read_length as usize != seq_len && !model_is_full_read_frame(annot, record) +} + +/// Which read a record's Fiber-seq tags describe; see [`record_frame`]. +#[derive(Debug, Clone, PartialEq, Eq)] +pub(crate) enum Frame { + /// The tags describe SEQ (unclipped, soft-clipped, SEQ-less, or computed + /// after clipping). + Seq, + /// The tags describe the full-length read of which SEQ is a hard-clipped + /// part. `h_lead` is the leading hard clip (BAM/CIGAR orientation, the + /// same end on both strands), `read_length` the full read's length. + FullRead { h_lead: u32, read_length: u32 }, + /// The tags describe some other read: drop them. + Stale(String), +} + +/// Which read a record's Fiber-seq tags describe, from signals read straight +/// off the record. Each tag family has its own frame signal: +/// - MA: the tag's own read length. Equal to SEQ (a hard-clipped record whose +/// tags were computed after clipping included): `Seq`. Equal to SEQ plus +/// the hard clips: `FullRead`, the aligner copied the full-length read's +/// tags onto this clipped part of it. Anything else: `Stale`. +/// - legacy `ns`/`nl`/`as`/`al` and fibertig `fs`/`fl`: no frame is recorded +/// and every producer wrote them on the full read, so any hard clip means +/// `FullRead` with `read_length = seq_len + H_lead + H_trail`; a missing +/// SEQ is `Stale`, since nothing then anchors them. +/// - MM/ML (with no MA or legacy tags deciding first): the SAM `MN` tag when +/// present (the spec's frame for exactly this case), else hard clips. +/// +/// In the full-read frame the MA-family coordinates are still right (they +/// are molecular coordinates of the full read; only the lift needs the +/// leading hard clip), but MM/ML are SEQ-relative and cannot be recovered: +/// the reader drops them and the writer strips them. +/// +/// SEQ-less MA records are never stale: the MA read length is the frame. +/// Both the reader ([`read_record`]) and the writer ([`write_record`]) use +/// this, so a stale record is cleared on the way in and cleaned on the way +/// out. +pub(crate) fn record_frame(record: &bam::Record) -> Frame { + let seq_len = record.seq_len(); + let (lead, trail) = hard_clips(record); + let hard_clipped = lead + trail > 0; + // The bases this record holds: SEQ, or the CIGAR query span when SEQ was + // dropped. A hard-clipped record's frame is judged against that span so + // a SEQ-less supplementary still gets the offset lift. + let span = query_span(record) as usize; + let full = span + lead as usize + trail as usize; + if let Some((ma, _, _)) = ma_family_tags(record) { + if let Some(read_length) = ma.split(';').next().and_then(|s| s.parse::().ok()) { + if hard_clipped && span > 0 && read_length == full { + return Frame::FullRead { + h_lead: lead, + read_length: read_length as u32, + }; + } + if seq_len > 0 && read_length != seq_len { + if hard_clipped && read_length == full { + return Frame::FullRead { + h_lead: lead, + read_length: read_length as u32, + }; + } + if hard_clipped { + return Frame::Stale(format!( + "MA read length {read_length} matches neither the {seq_len} bp sequence \ + nor the {full} bp read it was hard-clipped from" + )); + } + return Frame::Stale(format!( + "MA read length {read_length} does not match the {seq_len} bp sequence" + )); + } + // The MA tag was written for this SEQ, so it vouches for the + // MM/ML next to it too: fibertools' own writers emit MA and + // MM/ML together, and a hard-clipped record it produced is fine. + return Frame::Seq; + } + } else if has_legacy_nuc_msp(record) || has_legacy_fibertig(record) { + if seq_len == 0 { + return Frame::Stale("legacy nuc/msp tags on a record without SEQ".to_string()); + } + if hard_clipped { + // Any MM/ML next to them are SEQ-relative copies from the full + // read, dropped with the frame; they are not a stale signal. + return Frame::FullRead { + h_lead: lead, + read_length: full as u32, + }; + } + } + if matches!(record.aux(b"MM"), Ok(Aux::String(_))) { + if seq_len == 0 { + // Positions are implicit in SEQ; without it there is nothing to + // decode against (minimap2 -y writes SEQ-less secondaries). + return Frame::Stale("MM/ML on a record without SEQ".to_string()); + } + match record.aux(b"MN").ok().and_then(aux_as_usize) { + Some(mn) => { + if mn != seq_len { + return Frame::Stale(format!( + "MN {mn} does not match the {seq_len} bp sequence" + )); + } + } + None => { + if hard_clipped { + return Frame::Stale( + "MM/ML on a hard-clipped alignment without an MN tag".to_string(), + ); + } + } + } + } + Frame::Seq +} + +/// Why a record's tags describe a different read than its SEQ, or `None` +/// when they describe SEQ or the full read it was hard-clipped from. +pub(crate) fn record_frame_reason(record: &bam::Record) -> Option { + match record_frame(record) { + Frame::Stale(why) => Some(why), + _ => None, + } +} + +fn aux_as_usize(aux: Aux) -> Option { + match aux { + Aux::I8(v) => usize::try_from(v).ok(), + Aux::U8(v) => Some(v as usize), + Aux::I16(v) => usize::try_from(v).ok(), + Aux::U16(v) => Some(v as usize), + Aux::I32(v) => usize::try_from(v).ok(), + Aux::U32(v) => Some(v as usize), + _ => None, + } +} + +/// The parsed model disagrees with its frame: a read length that matches +/// neither SEQ nor the full read SEQ was hard-clipped from, or an annotation +/// ending past the frame. Backstop behind [`record_frame`] for tags that +/// carry no frame signal of their own. +fn model_frame_reason(annot: &MolecularAnnotations, record: &bam::Record) -> Option { + let seq_len = record.seq_len(); + if seq_len == 0 { + return None; + } + if frame_is_stale(annot, record) { + return Some(format!( + "MA read length {} does not match the {seq_len} bp sequence", + annot.read_length + )); + } + // SEQ, or the full read in the full-read frame. + let frame = annot.read_length as usize; + annot + .annotation_types + .iter() + .find(|t| { + t.annotations + .iter() + .any(|a| a.start as usize + a.length as usize > frame) + }) + .map(|t| format!("{} coordinates extend past the {frame} bp read", t.name)) +} + +/// Why a record's annotations do not fit its SEQ, or `None` when they do: +/// the record-level signals first, then the parsed model as a backstop. +#[cfg(test)] +pub(crate) fn stale_frame_reason( + annot: &MolecularAnnotations, + record: &bam::Record, +) -> Option { + record_frame_reason(record).or_else(|| model_frame_reason(annot, record)) +} + /// Make the fiberseq_callable annotation agree with the CLI minimums. /// Runs right after parsing, before any consumer-side pruning. Derives /// when the tag is absent (backfill) or the minimums are non-default /// (recalculation); otherwise the on-disk tag stands. +/// +/// A full-read-frame record (hard clips, tags from the full read) has no +/// m6A: the reader dropped its MM/ML. FIRE cannot score it, so it is +/// NotCallable whatever its on-disk tag or the minimums say. That is the +/// only "no m6A" this function knows about; m6A is never inspected, so a +/// read whose MM/ML were stripped on purpose (ft strip-basemods) keeps the +/// verdict from calling time. pub fn sync_fiberseq_callable( annot: &mut MolecularAnnotations, record: &bam::Record, filters: &crate::utils::input_bam::FiberFilters, ) { + let (min_msp, min_ave) = filters.callable_minimums(); + if model_is_full_read_frame(annot, record) { + // No m6A means not callable. Without nuc/msp there is nothing to + // judge either, so a never-called read gets no marker, but a copied + // Callable span must not survive the frame it no longer describes. + if has_calls(annot) || annot.get_type(FIBERSEQ_CALLABLE_TYPE).is_some() { + set_fiberseq_callable(annot, None, callable_minimums_name(min_msp, min_ave)); + } + return; + } let needed = filters.callable_minimums_are_custom() || annot.get_type(FIBERSEQ_CALLABLE_TYPE).is_none(); if needed && can_derive_callable(annot, record) { - let (min_msp, min_ave) = filters.callable_minimums(); derive_fiberseq_callable(annot, min_msp, min_ave); } } -/// True when a present SEQ disagrees with the recorded read length -/// (something rewrote the read after tagging). SEQ-less records are never -/// stale: the MA read length is the frame. -pub(crate) fn read_length_is_stale(read_length: u32, seq_len: usize) -> bool { - seq_len > 0 && read_length as usize != seq_len +/// True when calling ran: nuc or msp present. +fn has_calls(annot: &MolecularAnnotations) -> bool { + annot.get_type(NUC_TYPE).is_some() || annot.get_type(MSP_TYPE).is_some() } /// True when the callable state can be derived: calling ran (nuc or msp /// present) and the frame is not stale. Derivation is pure MA-tag /// arithmetic, so SEQ-less records derive fine. fn can_derive_callable(annot: &MolecularAnnotations, record: &bam::Record) -> bool { - !read_length_is_stale(annot.read_length, record.seq_len()) - && (annot.get_type(NUC_TYPE).is_some() || annot.get_type(MSP_TYPE).is_some()) + !frame_is_stale(annot, record) && has_calls(annot) } /// AN name recording non-default minimums, e.g. "m20a10". `None` for the @@ -146,10 +554,11 @@ pub fn set_fiberseq_callable( /// Derive the callable span and state from the nuc/MSP annotations; the /// single source of truth for every writer. The span is the extent of the -/// surviving calls; callable requires >= `min_msp` MSPs with mean length -/// >= `min_ave_msp_size`. m6A is never inspected: the caller emits no MSP -/// without m6A, and MA-only derivation is what lets SEQ-less records -/// derive. +/// surviving calls; callable requires at least `min_msp` MSPs with a mean +/// length of at least `min_ave_msp_size`. m6A is never inspected: the caller +/// emits no MSP without m6A, and MA-only derivation is what lets SEQ-less +/// records derive. A full-read-frame record (no m6A) is forced NotCallable +/// by [`sync_fiberseq_callable`], not here. pub fn derive_fiberseq_callable( annot: &mut MolecularAnnotations, min_msp: usize, @@ -184,23 +593,6 @@ pub fn derive_fiberseq_callable( ); } -/// Merge every annotation type from `src` into `dst`, skipping any type whose -/// name already exists on `dst`. Lets the legacy-tag fallbacks layer their -/// annotations on top of whatever the MA/MM/ML library parser already produced -/// rather than clobbering it (e.g. legacy nuc/msp alongside library-parsed -/// base mods). -fn merge_missing_types(dst: &mut MolecularAnnotations, src: MolecularAnnotations) { - for t in src.annotation_types.into_iter() { - if dst.get_type(&t.name).is_some() { - continue; - } - let new_t = dst.add_annotation_type(&t.name, t.quality_spec.clone(), t.encoding); - for a in t.annotations.into_iter() { - new_t.add_shared(a.start, a.length, a.strand, a.qualities, a.name); - } - } -} - /// Writes MA-family tags (MA/AQ/AN) to a BAM record, **preserving the /// record's existing MM/ML bytes**. /// @@ -238,16 +630,52 @@ fn merge_missing_types(dst: &mut MolecularAnnotations, src: MolecularAnnotations /// Producers that create or modify base mods must instead call /// [`write_record_with_basemods`], which canonically re-emits MM/ML. pub fn write_record(record: &mut bam::Record, annot: &MolecularAnnotations) { - strip_consumed_legacy_tags(record); + match record_frame(record) { + // The reader cleared this record's annotations (drop_stale_frame); + // leave no stale tag behind for another tool to trust. + Frame::Stale(_) => strip_all_fiber_tags(record), + Frame::FullRead { .. } => { + // The MA tag is rewritten from the model below, in the full + // read's frame. The SEQ-relative base mods describe bases this + // record does not carry: the reader dropped them, so leave none + // behind for another tool to trust. + for tag in [b"MM", b"ML", b"MN"] { + record.remove_aux(tag).ok(); + } + strip_consumed_legacy_tags(record); + } + Frame::Seq => { + if mn_disagrees(record) { + // The reader dropped these base mods (see read_record). + for tag in [b"MM", b"ML", b"MN"] { + record.remove_aux(tag).ok(); + } + } + strip_consumed_legacy_tags(record) + } + } annot.to_record(record); } +/// Every Fiber-seq tag fibertools knows how to read: legacy nuc/msp/fibertig +/// arrays, MM/ML/MN base mods, and both spellings of the MA family. Used only +/// for records whose frame is stale (`record_frame`), where none of them +/// describe this SEQ. +fn strip_all_fiber_tags(record: &mut bam::Record) { + for tag in [ + b"ns", b"nl", b"as", b"al", b"aq", b"fs", b"fl", b"fa", b"MM", b"ML", b"MN", b"Ma", b"Aq", + b"An", b"MA", b"AQ", b"AN", + ] { + record.remove_aux(tag).ok(); + } +} + /// Provenance rule shared by every MA write: strip exactly the legacy tags /// that [`read_record`] consumed as the source of the annotation model being /// written — they are superseded by the MA-family tags (v0.9 replace /// semantics; otherwise legacy readers silently see stale calls forever). /// -/// The gates mirror [`read_legacy_nuc_msp`]'s consumption at PAIR level +/// The gates mirror [`read_legacy_nuc_msp_into`]'s consumption at PAIR level /// (ns+nl, as+al, aq only inside the msp pair, fa only with fs+fl), which is /// what makes this provenance-safe: a tag is only removed when the reader /// ingested it. Records that already carry an MA-family main tag were not @@ -265,10 +693,13 @@ fn strip_consumed_legacy_tags(record: &mut bam::Record) { // writing an empty model — nothing was consumed, so nothing may be // removed. Gating on the same parse keeps strip and read consumption // identical by construction. - if read_legacy_nuc_msp(record).is_err() || read_legacy_fibertig(record).is_err() { + let mut scratch = MolecularAnnotations::new(0); + if read_legacy_nuc_msp_into(record, &mut scratch).is_err() + || read_legacy_fibertig_into(record, &mut scratch).is_err() + { return; } - // Pair-level gates mirroring read_legacy_nuc_msp: ns+nl only as a pair, + // Pair-level gates mirroring read_legacy_nuc_msp_into: ns+nl only as a pair, // as+al only as a pair, aq only inside a valid msp pair (its count is // validated by the parse above). A lone or orphan tag was never // consumed and so is never removed. @@ -284,7 +715,7 @@ fn strip_consumed_legacy_tags(record: &mut bam::Record) { } } // fibertig: fs+fl as a pair; fa is only consumed (and so only stripped) - // when the pair is non-empty, mirroring read_legacy_fibertig's early + // when the pair is non-empty, mirroring read_legacy_fibertig_into's early // return on empty fs. let fs = u32_array(record, b"fs"); if fs.is_some() && u32_array(record, b"fl").is_some() { @@ -328,6 +759,14 @@ fn has_legacy_fibertig(record: &bam::Record) -> bool { pub fn write_record_with_basemods(record: &mut bam::Record, annot: &MolecularAnnotations) { write_record(record, annot); annot.write_mm_ml(record); + // MN is the SAM spec's frame for MM/ML: the SEQ length they were written + // against. With it, any consumer (samtools, modkit, this reader) can tell + // when a later hard clip has made them stale. + record.remove_aux(b"MN").ok(); + let seq_len = record.seq_len(); + if seq_len > 0 && matches!(record.aux(b"MM"), Ok(Aux::String(_))) { + record.push_aux(b"MN", Aux::I32(seq_len as i32)).ok(); + } } /// Read annotations from a BAM record. @@ -348,8 +787,11 @@ fn read_ma_tags(record: &bam::Record) -> Result> { }; let mut annot = MolecularAnnotations::from_tags(&ma, aq.as_deref(), an.as_deref()) .map_err(|e| anyhow::anyhow!("MA tag parse error: {e}"))?; + // Same frame rule as the library's from_record: a full-read MA tag on a + // hard-clipped record lifts with the leading hard clip. + let offset = molecular_annotation::full_read_query_offset(annot.read_length, record); annot.set_aligned_blocks_raw( - molecular_annotation::AlignedBlocks::from_record(record), + AlignedBlocks::from_record(record).with_query_offset(offset.unwrap_or(0)), record.is_reverse(), ); Ok(Some(annot)) @@ -357,7 +799,13 @@ fn read_ma_tags(record: &bam::Record) -> Result> { fn read_legacy_nuc_msp(record: &bam::Record) -> Result { let mut annot = MolecularAnnotations::from_record(record); + read_legacy_nuc_msp_into(record, &mut annot)?; + Ok(annot) +} +/// Add the legacy `ns`/`nl`/`as`/`al`/`aq` tags to `annot` as nuc/msp/fire +/// types. The model's frame (read length, aligned blocks) is the caller's. +fn read_legacy_nuc_msp_into(record: &bam::Record, annot: &mut MolecularAnnotations) -> Result<()> { let ns = u32_array(record, b"ns"); let nl = u32_array(record, b"nl"); let a_starts = u32_array(record, b"as"); @@ -426,7 +874,7 @@ fn read_legacy_nuc_msp(record: &bam::Record) -> Result { } } - Ok(annot) + Ok(()) } /// Read the legacy fibertig `fs`/`fl`/`fa` tags into a [`FIBERTIG_TYPE`] @@ -439,13 +887,11 @@ fn read_legacy_nuc_msp(record: &bam::Record) -> Result { /// those older BAMs consumable. Mirrors the in-memory shape the MA path /// produces (`Strand::Forward`, no quality, `Encoding::Ma`) so downstream /// liftover to reference coordinates is identical either way. -fn read_legacy_fibertig(record: &bam::Record) -> Result { +fn read_legacy_fibertig_into(record: &bam::Record, annot: &mut MolecularAnnotations) -> Result<()> { use crate::utils::fibertig::FIBERTIG_TYPE; - let mut annot = MolecularAnnotations::from_record(record); - let (Some(fs), Some(fl)) = (u32_array(record, b"fs"), u32_array(record, b"fl")) else { - return Ok(annot); + return Ok(()); }; if fs.len() != fl.len() { bail!( @@ -455,7 +901,7 @@ fn read_legacy_fibertig(record: &bam::Record) -> Result { ); } if fs.is_empty() { - return Ok(annot); + return Ok(()); } // `fa` is optional; when present it must have one `|`-separated segment per @@ -483,7 +929,7 @@ fn read_legacy_fibertig(record: &bam::Record) -> Result { let name = names.as_ref().and_then(|v| v[i].clone()); t.add(*s, *l, Strand::Forward, vec![], name); } - Ok(annot) + Ok(()) } /// Convenience for callers that still operate on raw `i64` arrays of @@ -941,7 +1387,9 @@ mod tests { #[test] fn rewrite_replaces_ma_tag_instead_of_appending() { - let mut record = synth_record(b"ATCGATCGAT"); + // 300 bp so the msps below fit SEQ; read_record drops annotations + // that run past the sequence (#136). + let mut record = synth_record(&b"ATCGATCGAT".repeat(30)); let qspec_q = "Q".parse::().unwrap(); // First write: one msp at 100..150. @@ -980,4 +1428,483 @@ mod tests { assert_eq!(msp.annotations.len(), 1); assert_eq!(msp.annotations[0].start, 200, "read back stale MA tag"); } + + /// A mapped synthetic record with the given SEQ, CIGAR and flags. + fn synth_aligned(seq: &[u8], cigar: &str, flags: u16) -> bam::Record { + use rust_htslib::bam::record::CigarString; + let mut record = bam::Record::new(); + let qual = vec![60u8; seq.len()]; + let cigar = CigarString::try_from(cigar).expect("cigar parses"); + record.set(b"frame_test", Some(&cigar), seq, &qual); + record.set_flags(flags); + record.set_tid(0); + record.set_pos(0); + record + } + + fn legacy(record: &mut bam::Record, starts: &[u32], lens: &[u32]) { + record + .push_aux(b"ns", Aux::ArrayU32(starts.into())) + .unwrap(); + record.push_aux(b"nl", Aux::ArrayU32(lens.into())).unwrap(); + } + + fn nuc_starts(record: &bam::Record) -> Vec { + let annot = read_record(record).expect("read_record"); + annot + .get_type(NUC_TYPE) + .map(|t| t.annotations.iter().map(|a| a.start).collect()) + .unwrap_or_default() + } + + // MA carries its own frame: a mismatch is stale on either strand, a + // match is fine even with hard clips (tags computed after clipping). + #[test] + fn stale_frame_ma_read_length() { + let seq = b"ACGT".repeat(50); // 200 bp + for flags in [0u16, 16] { + let mut r = synth_aligned(&seq, "200M", flags); + r.push_aux(b"Ma", Aux::String("300;nuc.:10-40")).unwrap(); + assert!(record_frame_reason(&r).is_some(), "flags {flags}"); + assert!(nuc_starts(&r).is_empty()); + assert_eq!( + read_record(&r).unwrap().read_length, + 200, + "frame reset to SEQ" + ); + } + let mut r = synth_aligned(&seq, "50H200M", 2048); + r.push_aux(b"Ma", Aux::String("200;nuc.:10-40")).unwrap(); + assert!(record_frame_reason(&r).is_none()); + assert_eq!(nuc_starts(&r), vec![9], "MA text is 1-based"); + } + + // Legacy tags carry no frame and were always written on the full read: + // a hard clip puts them in the full-read frame, a soft clip does not, + // and no SEQ is stale. + #[test] + fn stale_frame_legacy_uses_hard_clips() { + let seq = b"ACGT".repeat(50); + let mut fits = synth_aligned(&seq, "30H200M", 2048); + legacy(&mut fits, &[10, 100], &[20, 20]); + assert_eq!( + record_frame(&fits), + Frame::FullRead { + h_lead: 30, + read_length: 230 + } + ); + assert!(record_frame_reason(&fits).is_none()); + assert_eq!( + nuc_starts(&fits), + vec![10, 100], + "legacy coords are full-read molecular coords" + ); + assert_eq!(read_record(&fits).unwrap().read_length, 230); + + let mut soft = synth_aligned(&seq, "30S170M", 2048); + legacy(&mut soft, &[10, 100], &[20, 20]); + assert!(record_frame_reason(&soft).is_none()); + assert_eq!(nuc_starts(&soft), vec![10, 100]); + + let mut exact = synth_aligned(&seq, "200M", 0); + legacy(&mut exact, &[180], &[20]); + assert!(stale_frame_reason(&read_record(&exact).unwrap(), &exact).is_none()); + assert_eq!(nuc_starts(&exact), vec![180]); + + let mut past = synth_aligned(&seq, "200M", 0); + legacy(&mut past, &[180], &[21]); + assert!(nuc_starts(&past).is_empty()); + + let mut seqless = synth_aligned(b"", "200M", 256); + legacy(&mut seqless, &[10], &[20]); + assert!(record_frame_reason(&seqless).is_some()); + let annot = read_record(&seqless).unwrap(); + assert!(annot.annotation_types.is_empty()); + assert_eq!( + annot.read_length, 200, + "frame from the CIGAR when SEQ is absent" + ); + } + + // MM/ML: MN is the frame when present, hard clips otherwise. + #[test] + fn stale_frame_mm_ml_uses_mn_then_hard_clips() { + let seq = b"ACGT".repeat(50); + let mm = |r: &mut bam::Record| { + r.push_aux(b"MM", Aux::String("A+a.,0;")).unwrap(); + r.push_aux(b"ML", Aux::ArrayU8((&[200u8][..]).into())) + .unwrap(); + }; + let mut mn_mismatch = synth_aligned(&seq, "200M", 0); + mm(&mut mn_mismatch); + mn_mismatch.push_aux(b"MN", Aux::I32(500)).unwrap(); + assert!(record_frame_reason(&mn_mismatch).is_some()); + assert!(read_record(&mn_mismatch) + .unwrap() + .annotation_types + .is_empty()); + + let mut mn_ok_clipped = synth_aligned(&seq, "50H200M", 2048); + mm(&mut mn_ok_clipped); + mn_ok_clipped.push_aux(b"MN", Aux::I32(200)).unwrap(); + assert!(record_frame_reason(&mn_ok_clipped).is_none()); + assert!(!read_record(&mn_ok_clipped) + .unwrap() + .annotation_types + .is_empty()); + + let mut no_mn_clipped = synth_aligned(&seq, "50H200M", 2048); + mm(&mut no_mn_clipped); + assert!(record_frame_reason(&no_mn_clipped).is_some()); + + let mut no_mn_soft = synth_aligned(&seq, "50S150M", 2048); + mm(&mut no_mn_soft); + assert!(record_frame_reason(&no_mn_soft).is_none()); + + // An MA tag written for this SEQ vouches for the MM/ML next to it: + // that is the shape fibertools' own writers produce. + let mut ma_vouches = synth_aligned(&seq, "50H200M", 2048); + mm(&mut ma_vouches); + ma_vouches + .push_aux(b"Ma", Aux::String("200;nuc.:10-40")) + .unwrap(); + assert!(record_frame_reason(&ma_vouches).is_none()); + let annot = read_record(&ma_vouches).unwrap(); + assert!(annot.get_type(NUC_TYPE).is_some()); + assert!(annot.get_type(crate::utils::basemods::M6A_TYPE).is_some()); + + let mut seqless = synth_aligned(b"", "200M", 256); + mm(&mut seqless); + assert!(record_frame_reason(&seqless).is_some()); + } + + // fibertools' own base-mod writer records the frame (MN) so its output + // survives a later hard clip check, including on a hard-clipped record. + #[test] + fn write_record_with_basemods_writes_mn_and_reads_back() { + use crate::utils::basemods::{canonical_header, M6A_TYPE}; + let seq = b"ACGT".repeat(50); + let mut r = synth_aligned(&seq, "50H200M", 2048); + let mut annot = read_record(&r).unwrap(); + let qspec = "Q".parse::().unwrap(); + let header = canonical_header(M6A_TYPE, b'A').unwrap().to_string(); + annot + .add_annotation_type(M6A_TYPE, qspec, Encoding::mm_ml()) + .add(0, 1, Strand::Forward, vec![200], Some(header)); + write_record_with_basemods(&mut r, &annot); + assert!( + matches!(r.aux(b"MN"), Ok(Aux::I32(200))), + "MN = {:?}", + r.aux(b"MN") + ); + assert!(record_frame_reason(&r).is_none()); + let back = read_record(&r).unwrap(); + assert!(back.get_type(M6A_TYPE).is_some(), "m6a lost on re-read"); + } + + // A stale record leaves the writer as an honest untagged read: no legacy + // arrays, no MM/ML/MN, and an MA tag whose frame is SEQ. + #[test] + fn write_record_strips_every_tag_of_a_stale_record() { + let seq = b"ACGT".repeat(50); + let mut r = synth_aligned(&seq, "200M", 0); + legacy(&mut r, &[10], &[20]); + r.push_aux(b"Ma", Aux::String("300;nuc.:11-20")).unwrap(); + r.push_aux(b"MM", Aux::String("A+a.,0;")).unwrap(); + r.push_aux(b"ML", Aux::ArrayU8((&[200u8][..]).into())) + .unwrap(); + let annot = read_record(&r).unwrap(); + assert!(annot.annotation_types.is_empty()); + write_record(&mut r, &annot); + for tag in [b"ns", b"nl", b"MM", b"ML", b"MN"] { + assert!( + r.aux(tag).is_err(), + "{} survived", + String::from_utf8_lossy(tag) + ); + } + let ma = r.aux(b"Ma"); + assert!( + matches!(ma, Ok(Aux::String(s)) if s.split(';').next() == Some("200")), + "Ma = {ma:?}" + ); + } + + // MA read length decides the frame three ways; the offset is the leading + // hard clip on both strands. + #[test] + fn frame_rule_ma_three_way() { + let seq = b"ACGT".repeat(50); // 200 bp + let mut full = synth_aligned(&seq, "50H200M50H", 2048); + // molecular [60,100), [200,240) + full.push_aux(b"Ma", Aux::String("300;nuc.:61-40,201-40")) + .unwrap(); + assert_eq!( + record_frame(&full), + Frame::FullRead { + h_lead: 50, + read_length: 300 + } + ); + assert!(record_frame_reason(&full).is_none()); + let annot = read_record(&full).unwrap(); + assert_eq!(annot.read_length, 300); + assert_eq!(annot.query_offset(), 50); + assert_eq!(nuc_starts(&full), vec![60, 200], "tag kept as is"); + assert_eq!( + annot.get_ref_coords(NUC_TYPE).unwrap(), + vec![ + (60, 100, Some(10), Some(50)), + (200, 240, Some(150), Some(190)) + ] + ); + + let mut clipped = synth_aligned(&seq, "50H200M50H", 2048); + clipped + .push_aux(b"Ma", Aux::String("200;nuc.:11-40")) + .unwrap(); + assert_eq!(record_frame(&clipped), Frame::Seq); + let annot = read_record(&clipped).unwrap(); + assert_eq!((annot.read_length, annot.query_offset()), (200, 0)); + assert_eq!( + annot.get_ref_coords(NUC_TYPE).unwrap(), + vec![(10, 50, Some(10), Some(50))], + "clipped frame lifts as today" + ); + + let mut stale = synth_aligned(&seq, "50H200M50H", 2048); + stale + .push_aux(b"Ma", Aux::String("250;nuc.:11-40")) + .unwrap(); + assert!(matches!(record_frame(&stale), Frame::Stale(_))); + assert!(nuc_starts(&stale).is_empty()); + assert_eq!(read_record(&stale).unwrap().read_length, 200); + + // reverse strand, trailing clip: offset 0; molecular [200,240) -> BAM [60,100) + let mut rev = synth_aligned(&seq, "200M100H", 2064); + rev.push_aux(b"Ma", Aux::String("300;nuc.:201-40")).unwrap(); + let annot = read_record(&rev).unwrap(); + assert_eq!(annot.query_offset(), 0); + assert_eq!( + annot.get_ref_coords(NUC_TYPE).unwrap(), + vec![(60, 100, Some(60), Some(100))] + ); + // reverse strand, leading clip: molecular [100,140) -> BAM [160,200) -> SEQ [60,100) + let mut rev_lead = synth_aligned(&seq, "100H200M", 2064); + rev_lead + .push_aux(b"Ma", Aux::String("300;nuc.:101-40")) + .unwrap(); + let annot = read_record(&rev_lead).unwrap(); + assert_eq!(annot.query_offset(), 100); + assert_eq!( + annot.get_ref_coords(NUC_TYPE).unwrap(), + vec![(160, 200, Some(60), Some(100))] + ); + // an annotation past the full read length is still stale + let mut past = synth_aligned(&seq, "50H200M50H", 2048); + past.push_aux(b"Ma", Aux::String("300;nuc.:281-40")) + .unwrap(); + assert!(nuc_starts(&past).is_empty()); + } + + // A full-read frame keeps nuc/msp, drops m6A, and the writer strips + // MM/ML/MN. + #[test] + fn full_frame_drops_mm_ml_and_strips_on_write() { + use crate::utils::basemods::M6A_TYPE; + let seq = b"ACGT".repeat(50); + // MN of the full read (dorado copy) or of SEQ: dropped either way + for mn in [250i32, 200] { + let mut r = synth_aligned(&seq, "50H200M", 2048); + r.push_aux(b"Ma", Aux::String("250;nuc.:61-40;msp.:101-20")) + .unwrap(); + r.push_aux(b"MM", Aux::String("A+a.,0;")).unwrap(); + r.push_aux(b"ML", Aux::ArrayU8((&[200u8][..]).into())) + .unwrap(); + r.push_aux(b"MN", Aux::I32(mn)).unwrap(); + let annot = read_record(&r).unwrap(); + assert!(annot.get_type(NUC_TYPE).is_some() && annot.get_type(MSP_TYPE).is_some()); + assert!( + annot.get_type(M6A_TYPE).is_none(), + "MN {mn}: m6A must be dropped on a full-read frame" + ); + write_record(&mut r, &annot); + for tag in [b"MM", b"ML", b"MN"] { + assert!( + r.aux(tag).is_err(), + "MN {mn}: {} survived", + String::from_utf8_lossy(tag) + ); + } + let ma = r.aux(b"Ma"); + assert!( + matches!(ma, Ok(Aux::String(s)) if s.starts_with("250;") && s.contains("nuc") && s.contains("msp")), + "Ma = {ma:?}" + ); + // write_record_with_basemods on the same model writes no MM/MN either + write_record_with_basemods(&mut r, &annot); + assert!(r.aux(b"MM").is_err() && r.aux(b"MN").is_err()); + // and the rewritten record reads back in the same frame + let back = read_record(&r).unwrap(); + assert_eq!((back.read_length, back.query_offset()), (250, 50)); + assert_eq!(back.get_forward_coords(NUC_TYPE), Some(vec![(60, 100)])); + } + } + + // Legacy ns/nl on a hard clip: read_length = seq + H_lead + H_trail, + // full-read frame. + #[test] + fn legacy_on_hard_clip_is_a_full_read_frame() { + let seq = b"ACGT".repeat(50); + let mut r = synth_aligned(&seq, "30H200M", 2048); + // [10,30) before SEQ, [20,60) straddles, [100,120) inside + legacy(&mut r, &[10, 20, 100], &[20, 40, 20]); + let annot = read_record(&r).unwrap(); + assert_eq!((annot.read_length, annot.query_offset()), (230, 30)); + assert_eq!( + annot.get_ref_coords(NUC_TYPE).unwrap(), + vec![ + (10, 30, None, None), + (20, 60, Some(0), Some(30)), + (100, 120, Some(70), Some(90)) + ] + ); + write_record(&mut r, &annot); + assert!( + r.aux(b"ns").is_err() && r.aux(b"nl").is_err(), + "consumed legacy tags stripped" + ); + assert!( + matches!(r.aux(b"Ma"), Ok(Aux::String(s)) if s.starts_with("230;") && s.contains("nuc.:11-20,21-40,101-20")), + "Ma = {:?}", + r.aux(b"Ma") + ); + let mut trail = synth_aligned(&seq, "200M30H", 0); + legacy(&mut trail, &[100], &[20]); + let annot = read_record(&trail).unwrap(); + assert_eq!((annot.read_length, annot.query_offset()), (230, 0)); + assert_eq!( + annot.get_ref_coords(NUC_TYPE).unwrap(), + vec![(100, 120, Some(100), Some(120))] + ); + // molecular [100,120) -> BAM [110,130) + let mut rev = synth_aligned(&seq, "200M30H", 2064); + legacy(&mut rev, &[100], &[20]); + assert_eq!( + read_record(&rev).unwrap().get_ref_coords(NUC_TYPE).unwrap(), + vec![(110, 130, Some(110), Some(130))] + ); + } + + // No m6A on a full-read frame means NotCallable, whatever the on-disk + // span says; an unprocessed full-frame record stays Untagged; a clipped + // frame keeps its span. + #[test] + fn full_frame_records_derive_not_callable() { + use crate::fiber::{CallableState, FiberseqData}; + let seq = b"ACGT".repeat(50); + let filters = crate::utils::input_bam::FiberFilters::default(); + // 10 MSPs, mean 15, end 150 + let msp = (0..10) + .map(|i| format!("{}-15", 1 + i * 15)) + .collect::>() + .join(","); + let state = |r: &bam::Record| { + FiberseqData::new(r.clone(), None, &filters) + .callable_state() + .0 + }; + let marker = |r: &bam::Record| { + let mut annot = read_record(r).unwrap(); + sync_fiberseq_callable(&mut annot, r, &filters); + annot + .get_type(FIBERSEQ_CALLABLE_TYPE) + .map(|t| (t.annotations[0].start, t.annotations[0].length)) + }; + let mut full = synth_aligned(&seq, "50H200M", 2048); + full.push_aux( + b"Ma", + Aux::String(&format!("250;msp.:{msp};fiberseq_callable.:1-150")), + ) + .unwrap(); + assert_eq!( + marker(&full), + Some((0, 0)), + "no m6A: NotCallable marker replaces the span" + ); + assert_eq!(state(&full), CallableState::NotCallable); + let mut clipped = synth_aligned(&seq, "50H200M", 2048); + clipped + .push_aux( + b"Ma", + Aux::String(&format!("200;msp.:{msp};fiberseq_callable.:1-150")), + ) + .unwrap(); + assert_eq!(marker(&clipped), Some((0, 150))); + assert_eq!(state(&clipped), CallableState::Callable); + let mut bare = synth_aligned(&seq, "50H200M", 2048); + bare.push_aux(b"Ma", Aux::String("250")).unwrap(); + assert_eq!(marker(&bare), None); + assert_eq!(state(&bare), CallableState::Untagged); + } + + // A SEQ-less hard-clipped record with a full-read MA tag is judged by + // its CIGAR span, so it gets the offset lift like a record with SEQ. + #[test] + fn seqless_hard_clipped_ma_is_full_read_frame() { + let mut r = synth_aligned(b"", "50H200M", 2048); + r.push_aux(b"Ma", Aux::String("250;nuc.:60-40")).unwrap(); + assert_eq!( + record_frame(&r), + Frame::FullRead { + h_lead: 50, + read_length: 250 + } + ); + let annot = read_record(&r).unwrap(); + assert!(annot.get_type(NUC_TYPE).is_some()); + assert!(model_is_full_read_frame(&annot, &r)); + } + + // An MN tag that disagrees with SEQ next to an MA tag that matches it: + // the calls stay, only the base mods go, on read and on write. + #[test] + fn mn_disagreement_drops_only_base_mods() { + let seq = b"ACGT".repeat(50); + let mut r = synth_aligned(&seq, "200M", 0); + r.push_aux(b"Ma", Aux::String("200;nuc.:10-40")).unwrap(); + r.push_aux(b"MM", Aux::String("A+a.,0;")).unwrap(); + r.push_aux(b"ML", Aux::ArrayU8((&[200u8][..]).into())) + .unwrap(); + r.push_aux(b"MN", Aux::I32(500)).unwrap(); + assert_eq!(record_frame(&r), Frame::Seq); + let annot = read_record(&r).unwrap(); + assert!(annot.get_type(NUC_TYPE).is_some()); + assert!(annot.annotation_types.iter().all(|t| !t.is_mm_ml())); + write_record(&mut r, &annot); + for tag in [b"MM", b"ML", b"MN"] { + assert!( + r.aux(tag).is_err(), + "{} survived", + String::from_utf8_lossy(tag) + ); + } + assert!(r.aux(b"Ma").is_ok()); + } + + // A Callable span copied onto a full-frame record with no calls of its + // own is replaced by the NotCallable marker rather than kept. + #[test] + fn full_frame_copied_callable_span_is_replaced() { + let seq = b"ACGT".repeat(50); + let filters = crate::utils::input_bam::FiberFilters::default(); + let mut r = synth_aligned(&seq, "50H200M", 2048); + r.push_aux(b"Ma", Aux::String("250;fiberseq_callable.:1-200")) + .unwrap(); + let mut annot = read_record(&r).unwrap(); + sync_fiberseq_callable(&mut annot, &r, &filters); + let t = annot.get_type(FIBERSEQ_CALLABLE_TYPE).expect("marker kept"); + assert_eq!(t.annotations[0].length, 0, "copied Callable span survived"); + } } diff --git a/src/utils/nucleosome.rs b/src/utils/nucleosome.rs index 47eb37a13..f42de93fd 100644 --- a/src/utils/nucleosome.rs +++ b/src/utils/nucleosome.rs @@ -201,6 +201,20 @@ pub fn add_nucleosomes_to_annotations( if record.seq_len() == 0 { return; } + // A full-read-frame record (hard clips, tags from the full read) has no + // m6A to call from; re-calling would erase its full-read nuc/msp and + // rewrite the read length to SEQ. Keep the record as it is; it is + // NotCallable (sync_fiberseq_callable) and cannot be scored by FIRE. + if ma_io::model_is_full_read_frame(annot, record) { + static FULL_FRAME_LOG: std::sync::Once = std::sync::Once::new(); + FULL_FRAME_LOG.call_once(|| { + log::warn!( + "add-nucleosomes skips hard-clipped reads whose tags are in the frame of the \ + full-length read (no m6A to call from); their nuc/msp are kept as is" + ); + }); + return; + } // The annotations must describe THIS record. An inherited MA field 0 // from a differently-sized input would make the tag we are about to // write read back as stale (Untagged) forever. diff --git a/tests/data/ont_hardclip_full_frame.bam b/tests/data/ont_hardclip_full_frame.bam new file mode 100644 index 000000000..329f0daa3 Binary files /dev/null and b/tests/data/ont_hardclip_full_frame.bam differ diff --git a/tests/data/ont_hardclip_full_frame.bam.bai b/tests/data/ont_hardclip_full_frame.bam.bai new file mode 100644 index 000000000..0b969b5bc Binary files /dev/null and b/tests/data/ont_hardclip_full_frame.bam.bai differ diff --git a/tests/data/ont_hardclip_full_frame.center.bed b/tests/data/ont_hardclip_full_frame.center.bed new file mode 100644 index 000000000..f73db0db1 --- /dev/null +++ b/tests/data/ont_hardclip_full_frame.center.bed @@ -0,0 +1,2 @@ +chr1_MATERNAL 2171 2172 nuc_a 0 + +chr1_MATERNAL 2171 2172 nuc_a_minus 0 - diff --git a/tests/data/ont_hardclip_mmml.bam b/tests/data/ont_hardclip_mmml.bam new file mode 100644 index 000000000..a9c3f1248 Binary files /dev/null and b/tests/data/ont_hardclip_mmml.bam differ diff --git a/tests/data/ont_hardclip_supplementary.bam b/tests/data/ont_hardclip_supplementary.bam new file mode 100644 index 000000000..edee745d0 Binary files /dev/null and b/tests/data/ont_hardclip_supplementary.bam differ diff --git a/tests/molecular_annotation.rs b/tests/molecular_annotation.rs index 7bd09e462..03b4b82b1 100644 --- a/tests/molecular_annotation.rs +++ b/tests/molecular_annotation.rs @@ -675,3 +675,63 @@ fn callable_does_not_require_a_decodable_m6a_type() { let fiber = FiberseqData::new(nomm, None, &explicit_defaults); assert_eq!(fiber.callable_state().0, CallableState::Callable); } + +/// A hard-clipped copy of a tagged record whose tags describe the full read +/// (#136). `set` keeps the aux tags, so the MA read length stays the full +/// length: the full-read frame. MM/ML stay too and the reader must drop them. +fn full_frame_copy(h: u32) -> (bam::Record, u32) { + use fibertools_rs::utils::input_bam::FiberFilters; + use fibertools_rs::utils::ma_io::sync_fiberseq_callable; + use rust_htslib::bam::record::{Cigar, CigarString}; + let record = read_records("msp_nuc.bam").into_iter().next().unwrap(); + let mut annot = read_record(&record).unwrap(); + sync_fiberseq_callable(&mut annot, &record, &FiberFilters::default()); + let mut tagged = record.clone(); + write_record(&mut tagged, &annot); + let len = tagged.seq_len() as u32; + let seq = tagged.seq().as_bytes()[h as usize..].to_vec(); + let qual = tagged.qual()[h as usize..].to_vec(); + let cigar = CigarString(vec![Cigar::HardClip(h), Cigar::Match(len - h)]); + let mut clipped = tagged.clone(); + clipped.set(tagged.qname(), Some(&cigar), &seq, &qual); + (clipped, len) +} + +#[test] +fn full_frame_record_is_not_callable() { + use fibertools_rs::fiber::{CallableState, FiberseqData}; + use fibertools_rs::utils::input_bam::FiberFilters; + let (clipped, len) = full_frame_copy(500); + for filters in [ + FiberFilters::default(), + FiberFilters { + min_msp: Some(1), + min_ave_msp_size: Some(1), + ..FiberFilters::default() + }, + ] { + let fiber = FiberseqData::new(clipped.clone(), None, &filters); + assert!(fiber.is_full_read_frame()); + assert_eq!(fiber.callable_state(), (CallableState::NotCallable, 0, 0)); + assert!(!fiber.is_callable()); + assert!(fiber.callable_reference_range().is_none()); + assert!(fiber.m6a().is_empty(), "m6A must be dropped"); + assert!(!fiber.nuc().is_empty() && !fiber.msp().is_empty()); + assert_eq!(fiber.frame_length(), (len - 500) as usize); + assert_eq!(fiber.annotations.query_offset(), 500); + } +} + +#[test] +fn full_frame_record_without_calls_stays_untagged() { + use fibertools_rs::fiber::{CallableState, FiberseqData}; + use fibertools_rs::utils::input_bam::FiberFilters; + let (mut clipped, len) = full_frame_copy(500); + clipped.remove_aux(b"Ma").unwrap(); + clipped + .push_aux(b"Ma", Aux::String(&len.to_string())) + .unwrap(); + let fiber = FiberseqData::new(clipped, None, &FiberFilters::default()); + assert!(fiber.is_full_read_frame()); + assert_eq!(fiber.callable_state().0, CallableState::Untagged); +} diff --git a/tests/nucleosome.rs b/tests/nucleosome.rs index 95192a114..3dc6cf423 100644 --- a/tests/nucleosome.rs +++ b/tests/nucleosome.rs @@ -189,3 +189,29 @@ fn seqless_record_passes_through_untouched() { assert_eq!(t.annotations[0].start, 100); assert_eq!(t.annotations[0].length, 400); } + +/// A hard-clipped record whose tags are in the full-read frame has no m6A +/// to call from: the producer leaves it alone instead of erasing its +/// full-read msp and rewriting the read length to SEQ (#136). +#[test] +fn add_nucleosomes_skips_full_frame_records() { + use fibertools_rs::utils::ma_io::{self, MSP_TYPE}; + use rust_htslib::bam::record::{Cigar, CigarString}; + let o = fibertools_rs::cli::NucleosomeParameters::default(); + let mut r = rec(1000); + let cigar = CigarString(vec![Cigar::HardClip(200), Cigar::Match(1000)]); + r.set(b"test", Some(&cigar), &vec![b'A'; 1000], &vec![255u8; 1000]); + r.set_tid(0); + r.set_pos(0); + let mut annot = MolecularAnnotations::new(1200); + ma_io::add_msp_annotations(&mut annot, &[300], &[50], None); + ma_io::write_record(&mut r, &annot); + let mut annot = ma_io::read_record(&r).unwrap(); + assert_eq!((annot.read_length, annot.query_offset()), (1200, 200)); + add_nucleosomes_to_annotations(&r, &mut annot, &[], &o, (10, 10)); + assert_eq!(annot.read_length, 1200); + let msp = annot.get_type(MSP_TYPE).expect("msp kept"); + assert_eq!(msp.annotations.len(), 1); + assert_eq!(msp.annotations[0].start, 300); + assert!(annot.get_type(FIBERSEQ_CALLABLE_TYPE).is_none()); +} diff --git a/tests/regression/center.rs b/tests/regression/center.rs index f37ff8d34..da691a5c5 100644 --- a/tests/regression/center.rs +++ b/tests/regression/center.rs @@ -26,3 +26,76 @@ fn center_default() { ] )); } + +// Centering on ref 2171 (the first retained nucleosome of the forward +// full-frame record). Rows of the two 2400 bp full-frame records +// (query_length 2400; the primary is 5376) must use the clip offset: +// molecular mode reproduces the primary's rows for the forward record and +// the reverse record's flipped ones; reference mode lists only the +// nucleosomes that lift. Expected values from pysam (#136). +#[test] +fn center_full_frame_records_use_the_clip_offset() { + let bam = fixture("ont_hardclip_full_frame.bam"); + let bed = fixture("ont_hardclip_full_frame.center.bed"); + let rows = |extra: &[&str], query_length: &str| -> Vec<(i64, i64)> { + let mut args = vec![ + "center", + bam.to_str().unwrap(), + "--bed", + bed.to_str().unwrap(), + "--dist", + "200", + ]; + args.extend_from_slice(extra); + let out = run(&args); + let mut lines = out.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + let col = |n: &str| header.iter().position(|h| *h == n).unwrap(); + let (strand, qlen, ty, st, en, cqs, cqe) = ( + col("strand"), + col("query_length"), + col("centered_position_type"), + col("centered_start"), + col("centered_end"), + col("centered_query_start"), + col("centered_query_end"), + ); + assert!( + lines + .clone() + .any(|l| l.split('\t').nth(strand) == Some("-")), + "minus-strand centering must produce rows" + ); + let mut v: Vec<(i64, i64)> = lines + .map(|l| l.split('\t').collect::>()) + .filter(|f| f[strand] == "+" && f[ty] == "nuc" && f[qlen] == query_length) + .inspect(|f| { + // molecular mode: leading columns are SEQ-relative (anchor at + // SEQ position 37, SEQ length 2400) + if query_length == "2400" && extra.is_empty() { + assert_eq!( + (f[cqs], f[cqe]), + ("-37", "2363"), + "SEQ frame leading columns" + ); + } + }) + .map(|f| (f[st].parse().unwrap(), f[en].parse().unwrap())) + .collect(); + v.sort(); + v + }; + // the primary, for reference (unchanged behaviour) + assert_eq!(rows(&[], "5376"), vec![(-163, -62), (0, 138)]); + assert_eq!(rows(&["--reference"], "5376"), vec![(-164, -62), (0, 137)]); + // forward full-frame record: same rows as the primary; reverse record: (36,116),(163,280) + assert_eq!( + rows(&[], "2400"), + vec![(-163, -62), (0, 138), (36, 116), (163, 280)] + ); + // reference mode: forward keeps only its lifted nucleosome, reverse its two + assert_eq!( + rows(&["--reference"], "2400"), + vec![(0, 137), (36, 116), (162, 278)] + ); +} diff --git a/tests/regression/common.rs b/tests/regression/common.rs index a41b35b2c..8b4cc7581 100644 --- a/tests/regression/common.rs +++ b/tests/regression/common.rs @@ -95,6 +95,24 @@ pub fn run(args: &[&str]) -> String { String::from_utf8(out.stdout).expect("non-UTF8 stdout") } +/// Run ft; return (stdout, stderr). Panics on non-zero exit. +pub fn run_capture(args: &[&str]) -> (String, String) { + let out = Command::new(ft()) + .args(args) + .output() + .expect("failed to spawn ft"); + assert!( + out.status.success(), + "ft exited {}\nstderr: {}", + out.status, + String::from_utf8_lossy(&out.stderr) + ); + ( + String::from_utf8(out.stdout).expect("non-UTF8 stdout"), + String::from_utf8_lossy(&out.stderr).into_owned(), + ) +} + /// Run `ft add-nucleosomes` on a fixture into a temp BAM, so tests get a /// BAM with the fiberseq_callable tag on disk. pub fn tagged_bam(name: &str) -> tempfile::NamedTempFile { diff --git a/tests/regression/convert_tags.rs b/tests/regression/convert_tags.rs index c120c4b59..80d96564b 100644 --- a/tests/regression/convert_tags.rs +++ b/tests/regression/convert_tags.rs @@ -172,3 +172,111 @@ fn convert_tags_migrates_uppercase_to_canonical() { assert_eq!(ml(b), ml(a), "ML changed"); } } + +// Hard-clipped supplementary reads keep the full-length read's tags, which do +// not match SEQ (#136). read_record drops those annotations and write_record +// strips every stale tag, so convert-tags writes an honest untagged record: an +// MA tag with no sections and no legacy arrays or MM/ML/MN. Primaries convert +// normally and keep their MM/ML. +fn assert_stale_records_cleaned(bam: &str, n_supp: usize) { + let out = NamedTempFile::with_suffix(".bam").unwrap(); + convert(&fixture(bam), out.path()); + let mut seen_supp = 0; + for rec in records(out.path()) { + let tag = ma(&rec).expect("every record gets an MA tag"); + // MA is ";": no sections means no annotations + let has_annotations = tag.trim_end_matches(';').contains(';'); + assert_eq!(!has_annotations, rec.is_supplementary(), "{bam} Ma {tag:?}"); + assert_eq!( + tag.split(';').next().unwrap(), + rec.seq_len().to_string(), + "{bam}: MA frame must be SEQ" + ); + for legacy in LEGACY_TAGS { + assert!(rec.aux(legacy).is_err(), "{bam}: legacy tag survived"); + } + let has_mm = rec.aux(b"MM").is_ok(); + assert_eq!( + has_mm, + !rec.is_supplementary(), + "{bam}: MM/ML on a stale record" + ); + if rec.is_supplementary() { + assert!(rec.aux(b"MN").is_err(), "{bam}: MN on a stale record"); + seen_supp += 1; + } + } + assert_eq!(seen_supp, n_supp, "{bam}"); +} + +/// Full-read-frame supplementaries (#136): convert-tags keeps nuc/msp under +/// the full read length, marks them NotCallable (no m6A), and strips the +/// consumed legacy tags and MM/ML/MN. +fn assert_full_frame_records_kept(bam: &str, expected: &[(&str, u16, u32)]) { + let out = NamedTempFile::with_suffix(".bam").unwrap(); + convert(&fixture(bam), out.path()); + let mut seen = 0; + for rec in records(out.path()) { + let tag = ma(&rec).expect("every record gets an MA tag"); + for legacy in LEGACY_TAGS { + assert!(rec.aux(legacy).is_err(), "{bam}: legacy tag survived"); + } + if !rec.is_supplementary() { + assert_eq!(tag.split(';').next().unwrap(), rec.seq_len().to_string()); + assert!(rec.aux(b"MM").is_ok(), "{bam}: primary lost MM"); + continue; + } + let qname = String::from_utf8_lossy(rec.qname()).to_string(); + let e = expected + .iter() + .find(|e| qname.starts_with(e.0) && rec.flags() == e.1) + .unwrap_or_else(|| panic!("{bam}: unexpected {qname} {}", rec.flags())); + assert_eq!( + tag.split(';').next().unwrap(), + e.2.to_string(), + "{bam} {qname}: MA frame must be the full read" + ); + assert!( + tag.contains(";nuc") && tag.contains(";msp"), + "{bam} {qname}: nuc/msp lost: {tag}" + ); + assert!( + tag.contains("fiberseq_callable.:1-0"), + "{bam} {qname}: no m6A means NotCallable: {tag}" + ); + for t in [b"MM", b"ML", b"MN"] { + assert!( + rec.aux(t).is_err(), + "{bam} {qname}: {} on a full-read frame", + String::from_utf8_lossy(t) + ); + } + seen += 1; + } + assert_eq!(seen, expected.len(), "{bam}"); +} + +#[test] +fn convert_tags_keeps_full_frame_legacy_records() { + assert_full_frame_records_kept( + "ont_hardclip_supplementary.bam", + &[("8ac3be13", 2048, 29940), ("4bd15181", 2064, 33088)], + ); +} + +#[test] +fn convert_tags_keeps_full_frame_ma_records() { + assert_full_frame_records_kept( + "ont_hardclip_full_frame.bam", + &[ + ("f2009f4d", 2048, 5376), + ("f2009f4d", 2064, 5376), + ("7b40cfd0", 2048, 9693), + ], + ); +} + +#[test] +fn convert_tags_cleans_hard_clipped_mm_ml_records() { + assert_stale_records_cleaned("ont_hardclip_mmml.bam", 1); +} diff --git a/tests/regression/extract.rs b/tests/regression/extract.rs index c09aa7f5a..859f014d0 100644 --- a/tests/regression/extract.rs +++ b/tests/regression/extract.rs @@ -1,4 +1,4 @@ -use super::common::{fixture, run, select_bed12_cols, select_tsv_cols}; +use super::common::{fixture, run, run_capture, select_bed12_cols, select_tsv_cols}; use tempfile::NamedTempFile; // bed12 columns worth snapshotting: locator + per-record feature data. @@ -101,3 +101,275 @@ fn extract_reads_ma_spelled_fixture() { let out = std::fs::read_to_string(tmp.path()).unwrap(); insta::assert_snapshot!(select_bed12_cols(&out, BED12_COLS)); } + +// The reader drops annotations that do not fit SEQ (hard-clipped supplementary +// reads keep the full-length read's tags, #136). Before this, extract printed +// misplaced coordinates for the forward read and u32-wrapped ones for the +// reverse read. Supplementaries must report `.` for nuc, msp and m6a; the +// primaries must not. +fn assert_supplementaries_untagged(bam: &str, n_primary: usize, n_supp: usize) { + let out = run(&["extract", "--all", "-", fixture(bam).to_str().unwrap()]); + let mut lines = out.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + let col = |name: &str| header.iter().position(|h| *h == name).unwrap(); + let (flag, nuc, msp, m6a) = ( + col("sam_flag"), + col("nuc_starts"), + col("msp_starts"), + col("m6a"), + ); + let (mut seen_primary, mut seen_supp) = (0, 0); + for line in lines { + let f: Vec<&str> = line.split('\t').collect(); + let supplementary = f[flag].parse::().unwrap() & 2048 != 0; + for c in [nuc, msp, m6a] { + assert_eq!(f[c] == ".", supplementary, "{bam} {}: {}", header[c], f[c]); + } + if supplementary { + seen_supp += 1; + } else { + seen_primary += 1; + } + } + assert_eq!((seen_primary, seen_supp), (n_primary, n_supp), "{bam}"); +} + +/// A hard-clipped supplementary whose tags describe the full-length read +/// (#136): nuc/msp are kept in the full read's frame and lifted through the +/// hard clip, m6A is dropped. Expected values were computed with pysam +/// get_aligned_pairs, independently of ft. +struct FullFrame { + /// qname prefix + qname: &'static str, + flag: u16, + /// SEQ length + /// The annotation frame: on a full-read-frame record the whole read, + /// not SEQ, so it matches the molecular nuc/msp columns and the + /// molecular-mode BED12 end. + fiber_length: i64, + /// every nucleosome of the full read stays in the tag + n_nuc: usize, + /// nuc_starts[0]: BAM orientation, full-read frame + first_nuc_start: i64, + /// ref_nuc_starts entries != -1, in output order (a prefix when shorter + /// than n_lifted) + lifted: &'static [i64], + n_lifted: usize, + same_nucs_as_primary: bool, +} + +fn assert_full_frame_supplementaries(bam: &str, n_primary: usize, expected: &[FullFrame]) { + let out = run(&["extract", "--all", "-", fixture(bam).to_str().unwrap()]); + let mut lines = out.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + let col = |name: &str| header.iter().position(|h| *h == name).unwrap(); + let (fiber, flag, len, nuc, ref_nuc, m6a) = ( + col("fiber"), + col("sam_flag"), + col("fiber_length"), + col("nuc_starts"), + col("ref_nuc_starts"), + col("m6a"), + ); + let ints = |s: &str| -> Vec { + s.trim_end_matches(',') + .split(',') + .map(|x| x.parse().unwrap()) + .collect() + }; + let mut primary_nucs = std::collections::HashMap::new(); + let mut supp_rows = Vec::new(); + for line in lines { + let f: Vec = line.split('\t').map(str::to_string).collect(); + if f[flag].parse::().unwrap() & 2048 == 0 { + assert_ne!(f[m6a], ".", "{bam} {}: primary lost m6A", f[fiber]); + primary_nucs.insert(f[fiber].clone(), f[nuc].clone()); + } else { + supp_rows.push(f); + } + } + assert_eq!( + (primary_nucs.len(), supp_rows.len()), + (n_primary, expected.len()), + "{bam}" + ); + for f in supp_rows { + let sam_flag: u16 = f[flag].parse().unwrap(); + let e = expected + .iter() + .find(|e| f[fiber].starts_with(e.qname) && e.flag == sam_flag) + .unwrap_or_else(|| { + panic!( + "{bam}: unexpected supplementary {} flag {sam_flag}", + f[fiber] + ) + }); + assert_eq!( + f[len].parse::().unwrap(), + e.fiber_length, + "{bam} {}", + e.qname + ); + assert_eq!( + f[m6a], ".", + "{bam} {} {}: m6A must be dropped on a full-read frame", + e.qname, e.flag + ); + let nucs = ints(&f[nuc]); + assert_eq!( + nucs.len(), + e.n_nuc, + "{bam} {} {}: nuc_starts", + e.qname, + e.flag + ); + assert_eq!( + nucs[0], e.first_nuc_start, + "{bam} {} {}: nuc_starts[0]", + e.qname, e.flag + ); + assert!( + nucs.windows(2).all(|w| w[0] < w[1]), + "{bam}: nuc_starts not ascending" + ); + if e.same_nucs_as_primary { + assert_eq!( + f[nuc], primary_nucs[&f[fiber]], + "{bam} {}: molecular nuc_starts must be unchanged", + e.qname + ); + } + let refs = ints(&f[ref_nuc]); + assert_eq!(refs.len(), e.n_nuc); + let lifted: Vec = refs.into_iter().filter(|r| *r != -1).collect(); + assert_eq!( + lifted.len(), + e.n_lifted, + "{bam} {} {}: lifted nucleosomes", + e.qname, + e.flag + ); + assert_eq!( + &lifted[..e.lifted.len()], + e.lifted, + "{bam} {} {}: ref_nuc_starts", + e.qname, + e.flag + ); + } +} + +#[test] +fn extract_lifts_full_frame_ma_records() { + assert_full_frame_supplementaries( + "ont_hardclip_full_frame.bam", + 2, + &[ + FullFrame { + qname: "f2009f4d", + flag: 2048, + fiber_length: 5376, + n_nuc: 31, + first_nuc_start: 64, + n_lifted: 13, + lifted: &[ + 2171, 2392, 2557, 2698, 2862, 3006, 3209, 3438, 3626, 3837, 3978, 4218, 4380, + ], + same_nucs_as_primary: true, + }, + FullFrame { + qname: "f2009f4d", + flag: 2064, + fiber_length: 5376, + n_nuc: 31, + first_nuc_start: 73, + n_lifted: 13, + lifted: &[ + 2207, 2333, 2575, 2721, 2925, 3093, 3231, 3498, 3666, 3806, 3982, 4166, 4359, + ], + same_nucs_as_primary: false, + }, + FullFrame { + qname: "7b40cfd0", + flag: 2048, + fiber_length: 9693, + n_nuc: 44, + first_nuc_start: 84, + n_lifted: 6, + lifted: &[102454, 102616, 102810, 102955, 103409, 103554], + same_nucs_as_primary: true, + }, + ], + ); +} + +// Legacy ns/nl/as/al on hard clips (dorado aligner shape): full-read frame +// with read_length = SEQ + H_lead + H_trail (29940 = 7440+22492+8; +// 33088 = 9910+0+23178). The 8ac3be13 nucleosome that straddles the leading +// clip snaps to the first aligned base (3834036), like across a soft clip. +#[test] +fn extract_lifts_full_frame_legacy_records() { + assert_full_frame_supplementaries( + "ont_hardclip_supplementary.bam", + 2, + &[ + FullFrame { + qname: "8ac3be13", + flag: 2048, + fiber_length: 29940, + n_nuc: 150, + first_nuc_start: 172, + n_lifted: 34, + lifted: &[ + 3834036, 3834112, 3834330, 3834492, 3834646, 3834860, 3834979, 3835452, + 3835701, 3835898, 3836264, 3836417, 3836689, 3836805, + ], + same_nucs_as_primary: false, + }, + FullFrame { + qname: "4bd15181", + flag: 2064, + fiber_length: 33088, + n_nuc: 162, + first_nuc_start: 217, + n_lifted: 48, + lifted: &[ + 10475654, 10475848, 10476012, 10476197, 10476540, 10476719, 10476912, 10477059, + 10477293, 10477492, 10477716, 10477893, 10478079, + ], + same_nucs_as_primary: false, + }, + ], + ); +} + +// One WARN per run that m6A was dropped, with the realignment remedy; no +// stale-frame noise. +#[test] +fn full_frame_m6a_dropped_warns_once() { + let (_, err) = run_capture(&[ + "extract", + "--all", + "-", + fixture("ont_hardclip_full_frame.bam").to_str().unwrap(), + ]); + assert_eq!( + err.matches("their m6A (MM/ML) describes bases this record does not carry and was dropped") + .count(), + 1, + "{err}" + ); + assert!(err.contains("minimap2 -Y -y"), "remedy missing: {err}"); + assert!( + !err.contains("dropping annotations for"), + "full-frame records are not stale: {err}" + ); +} + +// MM/ML/MN copied verbatim onto a 2376H hard-clipped supplementary (the +// plain minimap2 shape): caught by MN != SEQ length. Until 0.14 this record +// was deleted from every output with a per-record warning instead. +#[test] +fn extract_drops_mm_ml_on_hard_clipped_reads() { + assert_supplementaries_untagged("ont_hardclip_mmml.bam", 1, 1); +} diff --git a/tests/regression/fire.rs b/tests/regression/fire.rs index ab7e621e6..3a451552f 100644 --- a/tests/regression/fire.rs +++ b/tests/regression/fire.rs @@ -1,4 +1,4 @@ -use super::common::{fixture, run, select_tsv_cols, tagged_bam}; +use super::common::{fixture, run, run_capture, select_tsv_cols, tagged_bam}; use rust_htslib::bam::{self, Read}; use tempfile::NamedTempFile; @@ -27,6 +27,184 @@ fn fire_on_legacy_input_strips_consumed_legacy_tags() { } } +// Hard-clipped supplementary alignments keep the full-length read's tags +// (#136). The reader drops those annotations, so `ft fire` has nothing to +// score; it writes the record as an untagged read (MA tag with no sections, +// stale legacy arrays and MM/ML stripped) instead of panicking or passing the +// stale tags on. The fixtures hold scorable primary reads plus hard-clipped +// supplementaries from TEnCATS ONT data: one with legacy tags only, one with +// MM/ML/MN copied verbatim. +fn assert_fire_cleans_stale_records(bam: &str, n_scored: usize, n_cleaned: usize) { + let scored = NamedTempFile::with_suffix(".bam").unwrap(); + run(&[ + "fire", + "--ont", + fixture(bam).to_str().unwrap(), + scored.path().to_str().unwrap(), + ]); + let mut reader = bam::Reader::from_path(scored.path()).unwrap(); + let (mut seen_scored, mut seen_cleaned) = (0, 0); + for rec in reader.records() { + let rec = rec.unwrap(); + let ma = match rec.aux(b"Ma") { + Ok(bam::record::Aux::String(s)) => s.to_string(), + _ => panic!("{bam}: record without Ma tag"), + }; + if rec.is_supplementary() { + assert!( + !ma.trim_end_matches(';').contains(';'), + "{bam}: stale record kept annotations: {ma}" + ); + for tag in [b"as", b"ns", b"MM", b"ML", b"MN"] { + assert!( + rec.aux(tag).is_err(), + "{bam}: stale {} survived", + String::from_utf8_lossy(tag) + ); + } + seen_cleaned += 1; + } else { + assert!( + ma.contains("msp"), + "{bam}: scorable record missing msp in {ma}" + ); + seen_scored += 1; + } + } + assert_eq!((seen_scored, seen_cleaned), (n_scored, n_cleaned), "{bam}"); +} + +// Records with msp but no m6A cannot be scored: fire writes the model back +// with the full read length, no fire section, the NotCallable marker, and no +// MM/ML/MN (#136). +fn assert_fire_keeps_full_frame_records(bam: &str, n_scored: usize, expected: &[(&str, u16, u32)]) { + let scored = NamedTempFile::with_suffix(".bam").unwrap(); + run(&[ + "fire", + "--ont", + fixture(bam).to_str().unwrap(), + scored.path().to_str().unwrap(), + ]); + let mut reader = bam::Reader::from_path(scored.path()).unwrap(); + let (mut seen_scored, mut seen_full) = (0, 0); + for rec in reader.records() { + let rec = rec.unwrap(); + let ma = match rec.aux(b"Ma") { + Ok(bam::record::Aux::String(s)) => s.to_string(), + _ => panic!("{bam}: record without Ma tag"), + }; + if rec.is_supplementary() { + let qname = String::from_utf8_lossy(rec.qname()).to_string(); + let e = expected + .iter() + .find(|e| qname.starts_with(e.0) && rec.flags() == e.1) + .unwrap_or_else(|| panic!("{bam}: unexpected {qname}")); + assert_eq!( + ma.split(';').next().unwrap(), + e.2.to_string(), + "{bam} {qname}: Ma frame" + ); + assert!( + ma.contains(";nuc") && ma.contains(";msp"), + "{bam} {qname}: nuc/msp lost: {ma}" + ); + assert!( + !ma.contains(";fire"), + "{bam} {qname}: scored without m6A: {ma}" + ); + assert!( + ma.contains("fiberseq_callable.:1-0"), + "{bam} {qname}: not NotCallable: {ma}" + ); + for tag in [b"as", b"ns", b"MM", b"ML", b"MN"] { + assert!( + rec.aux(tag).is_err(), + "{bam} {qname}: {} survived", + String::from_utf8_lossy(tag) + ); + } + seen_full += 1; + } else { + assert!( + ma.contains("msp"), + "{bam}: scorable record missing msp in {ma}" + ); + seen_scored += 1; + } + } + assert_eq!( + (seen_scored, seen_full), + (n_scored, expected.len()), + "{bam}" + ); + // and the scored BAM counts them as NotCallable, never Untagged + let qc = run(&["qc", scored.path().to_str().unwrap()]); + let count = |state: &str| { + qc.lines() + .find(|l| l.starts_with(&format!("fiberseq_callable\t{state}\t"))) + .unwrap() + .split('\t') + .nth(2) + .unwrap() + .parse::() + .unwrap() + }; + assert_eq!( + (count("Callable"), count("NotCallable"), count("Untagged")), + (n_scored, expected.len(), 0), + "{bam}" + ); +} + +#[test] +fn fire_keeps_full_frame_legacy_records() { + assert_fire_keeps_full_frame_records( + "ont_hardclip_supplementary.bam", + 2, + &[("8ac3be13", 2048, 29940), ("4bd15181", 2064, 33088)], + ); +} + +#[test] +fn fire_keeps_full_frame_ma_records() { + assert_fire_keeps_full_frame_records( + "ont_hardclip_full_frame.bam", + 2, + &[ + ("f2009f4d", 2048, 5376), + ("f2009f4d", 2064, 5376), + ("7b40cfd0", 2048, 9693), + ], + ); +} + +// FireFeats::new slices SEQ by MSP coordinates: full-frame records must be +// skipped before it runs (feats-to-text) and contribute no feature rows. +#[test] +fn fire_feats_to_text_skips_full_frame_records() { + let bam = fixture("ont_hardclip_full_frame.bam"); + let all = run(&["fire", "--ont", "--feats-to-text", bam.to_str().unwrap()]); + let primaries = run(&[ + "fire", + "--ont", + "--feats-to-text", + "-F", + "2048", + bam.to_str().unwrap(), + ]); + assert_eq!( + all, primaries, + "full-frame records must add no feature rows" + ); + // must not panic + run(&["fire", "--ont", "--extract", bam.to_str().unwrap()]); +} + +#[test] +fn fire_cleans_hard_clipped_mm_ml_records() { + assert_fire_cleans_stale_records("ont_hardclip_mmml.bam", 1, 1); +} + fn extract_fdrs(out: &str) -> Vec { out.lines() .map(|l| l.split('\t').nth(9).unwrap().parse().unwrap()) @@ -215,3 +393,25 @@ fn fire_bam_mode_coverage_keeps_every_read_and_drop_removes() { "--drop must remove uncallable reads ({n_drop} vs {n_in})" ); } + +// ft fire strips the copied MM/ML from full-frame records, so a second pass +// over its own output has no m6A to drop and must stay quiet. +#[test] +fn full_frame_second_run_is_silent() { + let scored = NamedTempFile::with_suffix(".bam").unwrap(); + let (_, first) = run_capture(&[ + "fire", + "--ont", + fixture("ont_hardclip_full_frame.bam").to_str().unwrap(), + scored.path().to_str().unwrap(), + ]); + assert!( + first.contains("was dropped"), + "first run should warn: {first}" + ); + let (_, second) = run_capture(&["extract", "--all", "-", scored.path().to_str().unwrap()]); + assert!( + !second.contains("was dropped") && !second.contains("hard-clipped"), + "second run warned again: {second}" + ); +} diff --git a/tests/regression/pileup.rs b/tests/regression/pileup.rs index be69919b6..8f29e1f88 100644 --- a/tests/regression/pileup.rs +++ b/tests/regression/pileup.rs @@ -191,3 +191,61 @@ fn pileup_callable_fibers_shrinks_and_intersects() { assert_eq!(callable, fire, "--fire-coverage is the same flag"); assert_eq!(callable, fire_filter, "--fire-filter is the same flag"); } + +/// Weighted column sums of a pileup: (coverage, fire_coverage, nuc_coverage) bp. +fn pileup_sums(args: &[&str]) -> (i64, i64, i64) { + let tmp = NamedTempFile::new().unwrap(); + let mut a = vec!["pileup"]; + a.extend_from_slice(args); + a.extend_from_slice(&["-o", tmp.path().to_str().unwrap()]); + run(&a); + let out = std::fs::read_to_string(tmp.path()).unwrap(); + let mut lines = out.lines(); + let header: Vec<&str> = lines.next().unwrap().split('\t').collect(); + let col = |n: &str| header.iter().position(|h| *h == n).unwrap(); + let (s, e, cov, fire, nuc) = ( + col("start"), + col("end"), + col("coverage"), + col("fire_coverage"), + col("nuc_coverage"), + ); + lines.fold((0, 0, 0), |acc, l| { + let f: Vec = l.split('\t').map(|x| x.parse().unwrap_or(0)).collect(); + let w = f[e] - f[s]; + (acc.0 + f[cov] * w, acc.1 + f[fire] * w, acc.2 + f[nuc] * w) + }) +} + +// Full-read-frame supplementaries add their lifted nucleosomes to the +// nucleosome track (1737 + 1737 + 867 bp on top of the primaries' 8579) but +// never FIRE coverage, and --callable-fibers drops them from the denominator +// entirely (#136). +#[test] +fn pileup_full_frame_nucleosomes_count_fire_does_not() { + let bam = fixture("ont_hardclip_full_frame.bam"); + let bam = bam.to_str().unwrap(); + assert_eq!( + pileup_sums(&[bam, "-F", "2048"]), + (14257, 0, 8579), + "primaries alone" + ); + assert_eq!(pileup_sums(&[bam]), (20337, 0, 8579 + 1737 + 1737 + 867)); + assert_eq!( + pileup_sums(&[bam, "--callable-fibers"]), + (14069, 0, 8579), + "NotCallable reads leave the FIRE denominator" + ); + let scored = NamedTempFile::with_suffix(".bam").unwrap(); + run(&["fire", "--ont", bam, scored.path().to_str().unwrap()]); + index(scored.path()); + let s = scored.path().to_str().unwrap(); + let (_, fire_all, nuc_all) = pileup_sums(&[s]); + let (_, fire_prim, _) = pileup_sums(&[s, "-F", "2048"]); + assert!(fire_prim > 0, "primaries score some FIRE (615 bp today)"); + assert_eq!( + fire_all, fire_prim, + "fire coverage must not change when full-frame records are removed" + ); + assert_eq!(nuc_all, 12920); +} diff --git a/tests/regression/qc.rs b/tests/regression/qc.rs index 856375cdd..1b7ec28a7 100644 --- a/tests/regression/qc.rs +++ b/tests/regression/qc.rs @@ -326,3 +326,28 @@ fn qc_custom_minimums_reach_state_rows() { let filt: i64 = rows_for(&out, "phased_reads").iter().map(|r| r.2).sum(); assert_eq!(filt, 0, "filtered column agrees with the statet rows"); } + +fn callable_counts(bam: &str) -> (i64, i64, i64) { + let out = run(&["qc", fixture(bam).to_str().unwrap()]); + let get = |state: &str| -> i64 { + out.lines() + .find(|l| l.starts_with(&format!("fiberseq_callable\t{state}\t"))) + .unwrap_or_else(|| panic!("{bam}: no {state} row")) + .split('\t') + .nth(2) + .unwrap() + .parse() + .unwrap() + }; + (get("Callable"), get("NotCallable"), get("Untagged")) +} + +// Full-read-frame records have nuc/msp but no m6A: NotCallable, not Untagged +// (#136). MM/ML copied verbatim with no MA and no legacy tags stay stale +// (Untagged). +#[test] +fn qc_counts_full_frame_records_as_not_callable() { + assert_eq!(callable_counts("ont_hardclip_supplementary.bam"), (2, 2, 0)); + assert_eq!(callable_counts("ont_hardclip_full_frame.bam"), (2, 3, 0)); + assert_eq!(callable_counts("ont_hardclip_mmml.bam"), (1, 0, 1)); +}