fix: skip FIRE scoring when tag coordinates run past the sequence

.github/workflows/CI.yml

@@ -2,9 +2,9 @@ name: CI
 
 on:
   push:
-    branches: [main, master, dev, development]
+    branches: [main, master, dev, development, "release/**"]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, "release/**"]
 
 jobs:
   Formatting:

.github/workflows/release-plz.yml

@@ -7,9 +7,14 @@ name: release-plz
 # GitHub releases. The only cross-workflow hop — firing cargo-dist to build the
 # fibertools-rs binaries — uses workflow_dispatch, which the default
 # GITHUB_TOKEN is allowed to trigger, so no PAT is required.
+# This copy lives on the release/v0.13 maintenance branch: it runs on pushes
+# to that branch (push-triggered workflows use the pushed branch's file), so
+# merging a fix opens the patch release PR and merging that PR publishes it,
+# the same flow main has.
 on:
   push:
-    branches: [main]
+    branches: ["release/**"]
+  workflow_dispatch:
 
 permissions:
   contents: write
@@ -71,7 +76,9 @@ jobs:
           tag=$(echo "$RELEASES" | jq -r '.[] | select(.package_name=="fibertools-rs") | .tag')
           if [ -n "$tag" ] && [ "$tag" != "null" ]; then
             echo "Dispatching cargo-dist for $tag"
-            gh workflow run release.yml --ref main -f tag="$tag"
+            # Dispatch on the tag itself (not main) so the binaries are
+            # built from the patched 0.13.x source, not main's 0.14 code.
+            gh workflow run release.yml --ref "$tag" -f tag="$tag"
           else
             echo "fibertools-rs not in this release; skipping cargo-dist"
           fi

src/subcommands/fire.rs

@@ -15,6 +15,11 @@ pub fn add_fire_to_rec(
     model: &GBDT,
     precision_table: &MapPrecisionValues,
 ) {
+    // Skip (and pass through unchanged) records whose annotations cannot be
+    // scored; fire_coords_fit_seq logs the reason (#136).
+    if !fire_coords_fit_seq(rec) {
+        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
@@ -82,6 +87,7 @@ pub fn add_fire_to_bam(fire_opts: &mut FireOptions) -> Result<(), anyhow::Error>
             let chunk: Vec<FiberseqData> = chunk.collect();
             let feats: Vec<FireFeats> = chunk
                 .par_iter()
+                .filter(|r| fire_coords_fit_seq(r))
                 .map(|r| FireFeats::new(r, fire_opts))
                 .collect();
             feats.iter().for_each(|f| {

src/utils/fire.rs

@@ -61,6 +61,38 @@ fn get_mid_point(start: i64, end: i64) -> i64 {
     (start + end) / 2
 }
 
+/// Check that the annotations FIRE reads (msp, m6a, cpg) all fit inside the
+/// stored sequence, warning and returning false if they do not.
+///
+/// Hard-clipped supplementary alignments can keep nuc/msp tag coordinates
+/// from the full-length read, so positions run past the clipped SEQ and
+/// scoring the record indexes out of bounds (#136). The check uses the raw
+/// molecular-orientation coordinates: building a BAM-oriented view flips
+/// coordinates through `read_length - end`, which itself overflows on these
+/// records.
+pub fn fire_coords_fit_seq(rec: &FiberseqData) -> bool {
+    use crate::utils::basemods::{CPG_TYPE, M6A_TYPE};
+    use crate::utils::ma_io::MSP_TYPE;
+    let seq_len = rec.record.seq_len() as u64;
+    for type_name in [MSP_TYPE, M6A_TYPE, CPG_TYPE] {
+        let Some(t) = rec.annotations.get_type(type_name) else {
+            continue;
+        };
+        for a in &t.annotations {
+            if a.start as u64 + a.length as u64 > seq_len {
+                log::warn!(
+                    "skipping FIRE for {}: {} coordinates extend past the {} bp sequence (hard-clipped supplementary alignment?)",
+                    String::from_utf8_lossy(rec.record.qname()),
+                    type_name,
+                    seq_len
+                );
+                return false;
+            }
+        }
+    }
+    true
+}
+
 /// ```
 /// use fibertools_rs::utils::fire::get_bins;
 /// let bins = get_bins(50, 5, 20, 200);

tests/regression/fire.rs

@@ -27,6 +27,43 @@ fn fire_on_legacy_input_strips_consumed_legacy_tags() {
     }
 }
 
+// Hard-clipped supplementary alignments can keep nuc/msp tag coordinates
+// from the full-length read, so positions run past the clipped SEQ (and wrap
+// below zero when flipped on reverse-strand records). `ft fire` must skip
+// scoring these records instead of panicking, and still write them to the
+// output unchanged (#136). The fixture holds two scorable primary reads plus
+// a forward and a reverse hard-clipped supplementary read from TEnCATS ONT
+// data.
+#[test]
+fn fire_skips_records_whose_coords_exceed_the_sequence() {
+    let scored = NamedTempFile::with_suffix(".bam").unwrap();
+    run(&[
+        "fire",
+        "--ont",
+        fixture("ont_hardclip_supplementary.bam").to_str().unwrap(),
+        scored.path().to_str().unwrap(),
+    ]);
+    let mut reader = bam::Reader::from_path(scored.path()).unwrap();
+    let mut n_scored = 0;
+    let mut n_skipped = 0;
+    for rec in reader.records() {
+        let rec = rec.unwrap();
+        if rec.is_supplementary() {
+            assert!(rec.aux(b"Ma").is_err(), "unscorable record got a Ma tag");
+            assert!(
+                rec.aux(b"as").is_ok(),
+                "skipped record lost its original tags"
+            );
+            n_skipped += 1;
+        } else {
+            assert!(rec.aux(b"Ma").is_ok(), "scorable record missing Ma tag");
+            n_scored += 1;
+        }
+    }
+    assert_eq!(n_scored, 2);
+    assert_eq!(n_skipped, 2);
+}
+
 fn extract_fdrs(out: &str) -> Vec<f64> {
     out.lines()
         .map(|l| l.split('\t').nth(9).unwrap().parse().unwrap())