fix(clip): keep MM/ML in step with hard-clipped and unmapped reads

crates/fgumi-consensus/src/filter.rs

@@ -15,6 +15,7 @@ use fgumi_raw_bam::{AsTagBytes, RawRecord, RawRecordView, SamTag};
 
 pub use crate::modifications::{
     drop_masked_modifications_raw, drop_modifications_at_raw, has_modification_tags,
+    trim_clipped_modifications_raw,
 };
 
 /// Expands a 1-3 element slice to a 3-element array, filling missing values from the last.
@@ -2764,6 +2765,132 @@ mod tests {
         }
     }
 
+    // -- trim_clipped_modifications_raw tests --
+
+    /// One `trim_clipped_modifications_raw` case: the record before and after clipping, and the
+    /// tags expected afterwards.
+    struct TrimCase {
+        pre_seq: &'static [u8],
+        pre_reverse: bool,
+        post_seq: &'static [u8],
+        post_flags: u16,
+        removed_start: Option<usize>,
+        mm: &'static str,
+        ml: &'static [u8],
+        mn: Option<i32>,
+        /// Expected `MM` (and `am`, written with the same value), `ML`, `MN`, and whether any tag
+        /// was removed.
+        expected: (Option<&'static str>, Option<&'static [u8]>, Option<i64>, bool),
+    }
+
+    /// A forward read `CACGTCAACG` (C at 0, 2, 5, 8) clipped to `post_seq`.
+    fn forward_trim(
+        post_seq: &'static [u8],
+        removed_start: Option<usize>,
+        mm: &'static str,
+        ml: &'static [u8],
+        mn: Option<i32>,
+        expected: (Option<&'static str>, Option<&'static [u8]>, Option<i64>, bool),
+    ) -> TrimCase {
+        TrimCase {
+            pre_seq: b"CACGTCAACG",
+            pre_reverse: false,
+            post_seq,
+            post_flags: R1_FWD,
+            removed_start,
+            mm,
+            ml,
+            mn,
+            expected,
+        }
+    }
+
+    /// Clipping removes bases from the ends of SEQ (or, unmapping a reverse read,
+    /// reverse-complements it): the calls outside the kept window are dropped, the skips are
+    /// recomputed over the window, and `MN` becomes the new length. Tags that cannot be placed
+    /// are removed.
+    #[rstest]
+    #[case::trimmed_start(forward_trim(b"CGTCAACG", Some(2), "C+m?,0,0,0;", &[10, 20, 30], Some(10), (Some("C+m?,0,0;"), Some(&[20, 30][..]), Some(8), false)))]
+    #[case::trimmed_both_ends(forward_trim(b"CGTCAA", Some(2), "C+m?,0,0,0;", &[10, 20, 30], Some(10), (Some("C+m?,0,0;"), Some(&[20, 30][..]), Some(6), false)))]
+    #[case::call_in_trimmed_tail(forward_trim(b"CGTCAA", Some(2), "C+m?,2,0;", &[10, 20], Some(10), (Some("C+m?,1;"), Some(&[10][..]), Some(6), false)))]
+    #[case::without_mn(forward_trim(b"CGTCAACG", Some(2), "C+m?,0,0,0;", &[10, 20, 30], None, (Some("C+m?,0,0;"), Some(&[20, 30][..]), None, false)))]
+    #[case::group_on_n(forward_trim(b"CGTCAACG", Some(2), "N+n?,0,0,0;", &[10, 20, 30], Some(10), (Some("N+n?,0;"), Some(&[30][..]), Some(8), false)))]
+    #[case::masked_not_trimmed(forward_trim(b"NNCGTCAACG", Some(0), "C+m?,0,0,0;", &[10, 20, 30], Some(10), (Some("C+m?,0,0;"), Some(&[20, 30][..]), Some(10), false)))]
+    #[case::unknown_offset(forward_trim(b"CGTCAACG", None, "C+m?,0,0,0;", &[10, 20, 30], Some(10), (None, None, None, true)))]
+    #[case::stale_before_clipping(forward_trim(b"CGTCAACG", Some(2), "C+m?,0,0,0;", &[10, 20, 30], Some(12), (None, None, None, true)))]
+    #[case::not_a_window_of_the_old_seq(forward_trim(b"CGTCAACC", Some(2), "C+m?,0,0,0;", &[10, 20, 30], Some(10), (None, None, None, true)))]
+    // Unmapping a reverse read reverse-complements SEQ into read orientation, where MM already
+    // pointed: nothing changes. `CACGTCAACG` reverse reads `CGTTGACGTG` (C at 0 and 6).
+    #[case::reverse_read_unmapped(TrimCase {
+        pre_seq: b"CACGTCAACG",
+        pre_reverse: true,
+        post_seq: b"CGTTGACGTG",
+        post_flags: flags::UNMAPPED,
+        removed_start: Some(0),
+        mm: "C+m?,0,0;",
+        ml: &[10, 20],
+        mn: Some(10),
+        expected: (Some("C+m?,0,0;"), Some(&[10, 20][..]), Some(10), false),
+    })]
+    // Hard clipping 4 bases from the start of a reverse read's stored SEQ removes the last 4
+    // bases of the read: `CGTTGA` keeps the C at 0 and drops the call on the C at 6.
+    #[case::reverse_read_trimmed(TrimCase {
+        pre_seq: b"CACGTCAACG",
+        pre_reverse: true,
+        post_seq: b"TCAACG",
+        post_flags: R1_REV,
+        removed_start: Some(4),
+        mm: "C+m?,0,0;",
+        ml: &[10, 20],
+        mn: Some(10),
+        expected: (Some("C+m?,0;"), Some(&[10][..]), Some(6), false),
+    })]
+    // A record without SEQ (`*`) is left alone, as htslib accepts `MM` there.
+    #[case::no_seq(TrimCase {
+        pre_seq: b"",
+        pre_reverse: false,
+        post_seq: b"",
+        post_flags: R1_FWD | flags::SECONDARY,
+        removed_start: Some(3),
+        mm: "C+m?,0;",
+        ml: &[10],
+        mn: None,
+        expected: (Some("C+m?,0;"), Some(&[10][..]), None, false),
+    })]
+    fn test_trim_clipped_modifications_raw(#[case] case: TrimCase) {
+        let mut b = RawSamBuilder::new();
+        b.flags(case.post_flags).ref_id(0).pos(0).mapq(60).sequence(case.post_seq);
+        if !case.post_seq.is_empty() && case.post_flags & flags::UNMAPPED == 0 {
+            b.cigar_ops(&[u32::try_from(case.post_seq.len()).expect("short fixture") << 4])
+                .qualities(&vec![30; case.post_seq.len()]);
+        } else if !case.post_seq.is_empty() {
+            b.qualities(&vec![30; case.post_seq.len()]);
+        }
+        b.add_string_tag(SamTag::MM, case.mm.as_bytes())
+            .add_array_u8(SamTag::ML, case.ml)
+            .add_string_tag(SamTag::AM_BASES, case.mm.as_bytes());
+        if let Some(mn) = case.mn {
+            b.add_int_tag(SamTag::MN, mn);
+        }
+        let mut record = b.build().as_ref().to_vec();
+
+        let removed = trim_clipped_modifications_raw(
+            &mut record,
+            case.pre_seq,
+            case.pre_reverse,
+            case.removed_start,
+        );
+
+        let aux = bam_fields::aux_data_slice(&record);
+        let ml = bam_fields::find_array_tag(aux, SamTag::ML).map(|a| a.data.to_vec());
+        let (want_mm, want_ml, want_mn, want_removed) = case.expected;
+        assert_eq!(string_tag(&record, SamTag::MM).as_deref(), want_mm, "MM");
+        assert_eq!(ml.as_deref(), want_ml, "ML");
+        assert_eq!(string_tag(&record, SamTag::AM_BASES).as_deref(), want_mm, "am");
+        assert_eq!(bam_fields::find_int_tag(aux, SamTag::MN), want_mn, "MN");
+        assert_eq!(removed, want_removed, "removed");
+    }
+
     // -- check_conversion_fraction_raw tests --
 
     /// OB-derived reads (R1 reverse, R2 forward) carry their evidence at reference G, so the

crates/fgumi-consensus/src/modifications.rs

@@ -176,6 +176,96 @@ pub fn drop_masked_modifications_raw(record: &mut Vec<u8>, pre_mask_seq: &[u8])
     })
 }
 
+/// Keeps `MM`/`ML`, `MN` and the per-strand `am`/`bm` consistent with SEQ after clipping.
+///
+/// Clipping can mask bases to `N` (`soft-with-mask`), remove bases from the ends of SEQ (hard
+/// clipping), and, when it unmaps a reverse-mapped read, reverse-complement SEQ. The tags index
+/// SEQ in original read orientation, so both sequences are first put in that orientation:
+/// `pre_clip_seq` (SEQ as stored before clipping) by `pre_clip_reverse` (the record's reverse
+/// flag then), and SEQ now by the record's reverse flag now. `removed_start` is the number of
+/// bases removed from the start of the stored pre-clip SEQ, or `None` when it is not known, in
+/// which case a shorter SEQ cannot be placed.
+///
+/// SEQ must be a window of the pre-clip SEQ, apart from bases masked to `N`. The calls outside
+/// the window or on a masked base are dropped and the skips recomputed over the bases that
+/// remain (a group on base `N` counts every base of the window), and `MN`, when present, is set
+/// to the new length. All the tags are removed instead when they did not fit SEQ before
+/// clipping (`MN` differs from its length) or SEQ is not such a window; a tag that cannot be
+/// edited (see [`drop_masked_modifications_raw`]) is removed alone. A record whose SEQ is `*`
+/// before and after is left alone. Returns whether any tag was removed.
+pub fn trim_clipped_modifications_raw(
+    record: &mut Vec<u8>,
+    pre_clip_seq: &[u8],
+    pre_clip_reverse: bool,
+    removed_start: Option<usize>,
+) -> bool {
+    if !has_modification_tags(record) {
+        return false;
+    }
+    let view = RawRecordView::new(record);
+    let post_reverse = view.is_reverse();
+    let post_clip_seq = view.sequence_vec();
+    if pre_clip_seq.is_empty() && post_clip_seq.is_empty() {
+        return false;
+    }
+    let mn = bam_fields::find_int_tag(bam_fields::aux_data_slice(record), SamTag::MN);
+    if mn.is_some_and(|mn| usize::try_from(mn).ok() != Some(pre_clip_seq.len())) {
+        return remove_modification_tags(record);
+    }
+    if pre_clip_reverse == post_reverse && post_clip_seq == pre_clip_seq {
+        return false;
+    }
+    let in_read_orientation = |seq: &[u8], reverse: bool| {
+        if reverse { fgumi_dna::dna::reverse_complement(seq) } else { seq.to_vec() }
+    };
+    let before = in_read_orientation(pre_clip_seq, pre_clip_reverse);
+    let after = in_read_orientation(&post_clip_seq, post_reverse);
+    // Where the window starts in `before`. Bases removed from the start of a reverse-mapped
+    // read's stored SEQ are removed from the end of the read.
+    let lo = if after.len() == before.len() {
+        Some(0)
+    } else {
+        removed_start.filter(|_| pre_clip_reverse == post_reverse).and_then(|start| {
+            if pre_clip_reverse {
+                before.len().checked_sub(start)?.checked_sub(after.len())
+            } else {
+                Some(start)
+            }
+        })
+    };
+    let window = lo.and_then(|lo| Some(lo..lo.checked_add(after.len())?));
+    let fits = window.as_ref().and_then(|w| before.get(w.clone())).is_some_and(|kept| {
+        kept.iter().zip(&after).all(|(&b, &a)| a == b || a.eq_ignore_ascii_case(&b'N'))
+    });
+    let Some(window) = window.filter(|_| fits) else {
+        return remove_modification_tags(record);
+    };
+    if after == before {
+        return false;
+    }
+    let removed = edit_modification_tags_raw(record, |tag, ml| {
+        rewrite_modifications(
+            tag,
+            ml,
+            &before,
+            |base, pos| {
+                window.contains(&pos) && (base == b'N' || after[pos - window.start] == base)
+            },
+            |_| true,
+        )
+    });
+    if after.len() != before.len()
+        && bam_fields::find_int_tag(bam_fields::aux_data_slice(record), SamTag::MN).is_some()
+    {
+        // `MN` is a 32-bit signed tag; a SEQ too long for it cannot be described.
+        let Ok(l_seq) = i32::try_from(after.len()) else {
+            return remove_modification_tags(record);
+        };
+        bam_fields::RawTagsEditor::from_vec(record).update_int(SamTag::MN, l_seq);
+    }
+    removed
+}
+
 /// Drops the methylation calls at `positions` (genomic orientation) from `MM`/`ML` and the
 /// per-strand `am`/`bm`, leaving SEQ unchanged.
 ///

crates/fgumi-sam/src/clipper.rs

@@ -12,6 +12,21 @@ const TAGS_TO_REVERSE_COMPLEMENT: [fgumi_raw_bam::SamTag; 2] =
 
 /// Quality-oriented aux tags that are reversed (with `QUAL`) when a read is
 /// reverse-complemented, mirroring htsjdk `SAMRecord.TAGS_TO_REVERSE`.
+/// Base modification tags. `MM`/`am`/`bm` strings and the `ML` array are not per-base, so
+/// `--auto-clip-attributes` must not slice them even when their length happens to equal the
+/// read's; the `clip` command keeps them in step with the clipped SEQ itself.
+const MODIFICATION_TAGS: [fgumi_raw_bam::SamTag; 4] = [
+    fgumi_raw_bam::SamTag::MM,
+    fgumi_raw_bam::SamTag::ML,
+    fgumi_raw_bam::SamTag::AM_BASES,
+    fgumi_raw_bam::SamTag::BM_BASES,
+];
+
+/// Whether `tag` is one of [`MODIFICATION_TAGS`].
+fn is_modification_tag(tag: [u8; 2]) -> bool {
+    MODIFICATION_TAGS.iter().any(|t| **t == tag)
+}
+
 const TAGS_TO_REVERSE: [fgumi_raw_bam::SamTag; 2] =
     [fgumi_raw_bam::SamTag::OQ, fgumi_raw_bam::SamTag::U2];
 
@@ -94,6 +109,9 @@ impl RawRecordClipper {
 
         for entry in view.iter_typed() {
             let (tag, value) = entry;
+            if is_modification_tag(tag) {
+                continue;
+            }
             match value {
                 TagValue::String(s) => {
                     if s.len() == old_length {
@@ -1044,6 +1062,9 @@ impl RawRecordClipper {
 
             for (tag, value) in view.iter_typed() {
                 use fgumi_raw_bam::TagValue;
+                if is_modification_tag(tag) {
+                    continue;
+                }
                 match value {
                     TagValue::String(s) => {
                         if s.len() == old_seq_len {

docs/src/guide/methylation.md

@@ -315,7 +315,7 @@ The methylation options of `filter` read the `cu`/`ct` counts written by methyla
 - **Simplex** (reference-anchored): an aligned reference cytosine of the read's own strand: reference C or G, depending on the read type and the direction it aligned. A site where the read shows a third base (a C>A or C>G variant) carries no call and is never masked. Simplex methylation lives in the sequence, so a failing site is masked to N.
 - **Duplex** (molecule-based): a C or G in the read's own sequence, which is the molecule's sequence. A site without counts carries no call (the strands did not confirm it, or they conflicted). A CpG is a C followed in the read by a G: an insertion between them splits it, and a deletion between them joins one. An N is never a site. Duplex methylation lives in `MM`/`ML`, so a failing site's calls are dropped from `MM`/`ML` (and `am`/`bm`) and the base is left alone: a base that `--min-reads` accepted is never lost for lack of a methylation call. No reference is needed, so unaligned duplex reads are filtered too.
 
-At every site a methylation filter rejects, filter also zeroes the counts (`cu`/`ct`, and `au`/`at`/`bu`/`bt` on duplex), so the counts never report a call that the sequence or `MM`/`ML` no longer carry. A dropped call becomes an unlisted base in its `?` group. Filter does not edit a group without `?`, where an unlisted base would read as unmethylated, and removes `MM`/`ML` instead. It also removes `MM`/`ML`/`MN`/`am`/`bm` when they no longer fit the sequence (for example after hard clipping changed its length). Records whose `cu`/`ct` are missing or no longer match the sequence are left unchecked.
+At every site a methylation filter rejects, filter also zeroes the counts (`cu`/`ct`, and `au`/`at`/`bu`/`bt` on duplex), so the counts never report a call that the sequence or `MM`/`ML` no longer carry. A dropped call becomes an unlisted base in its `?` group. Filter does not edit a group without `?`, where an unlisted base would read as unmethylated, and removes `MM`/`ML` instead. It also removes `MM`/`ML`/`MN`/`am`/`bm` when they no longer fit the sequence (for example after a tool other than `fgumi clip` hard-clipped it). Records whose `cu`/`ct` are missing or no longer match the sequence are left unchecked.
 
 At the end of the run, filter reports each of these counts: unmapped simplex reads skipped, records without `cu`/`ct`, records whose `cu`/`ct` do not match the sequence, records whose modification tags were removed, duplex calls dropped from `MM`/`ML`, and single-strand records that `--require-strand-methylation-agreement` could not be applied to. It says so plainly when no record was checked at all. Without `--reverse-per-base-tags`, it also reports how many reverse-mapped reads the methylation filters checked, whose counts were read at the wrong positions unless they had already been reversed.
 
@@ -368,7 +368,7 @@ With `--ref`, `filter` (and `clip`) recompute `NM`, `UQ` and `MD` from the read'
 
 `MD` always lists every difference, so the sequence, CIGAR and `MD` reconstruct the reference, as with `samtools calmd`; it can therefore differ from the `MD` a bisulfite-aware aligner wrote. A duplex read, whose sequence is the molecule's, is scored literally. The rule follows the read, not `--methylation-mode`.
 
-`clip --clipping-mode soft-with-mask` masks clipped bases to N and drops their calls from `MM`/`ML`; unlike filter, it leaves `cu`/`ct` unchanged.
+`clip` keeps `MM`/`ML` in step with the clipped read: `--clipping-mode soft-with-mask` masks the clipped bases to N and hard clipping removes them, and either way their calls are dropped from `MM`/`ML` (and `am`/`bm`); after hard clipping, `MN` (when present) is the new read length. `--auto-clip-attributes` leaves these tags to that step. Tags that cannot be kept in step are removed instead, as filter does, and clip reports how many records lost them at the end of the run: for example a group without the `?` flag that would lose a call, or a read that clipping shortened and then unmapped, which leaves no record of where the removed bases were. Unlike filter, clip does not zero `cu`/`ct`.
 
 ### Recommended Parameters