feat!: keep nuc/MSP/FIRE calls on hard-clipped alignments, drop only their m6A, never panic

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`

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

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
 

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.

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;")

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

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

molecular-annotation/src/iter.rs

@@ -131,22 +131,25 @@ 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<Item = ProjectedAnnotation<'_>> + '_ {
+        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 {
                     self.flip_range(a.start, a.end())
                 } 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,

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,