diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 28eddc52d..f33309d02 100644 --- a/.github/workflows/CI.yml +++ b/.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: diff --git a/.github/workflows/release-plz.yml b/.github/workflows/release-plz.yml index 36775e20a..66d6b8999 100644 --- a/.github/workflows/release-plz.yml +++ b/.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 diff --git a/src/subcommands/fire.rs b/src/subcommands/fire.rs index d33cc87d8..5c95e5059 100644 --- a/src/subcommands/fire.rs +++ b/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 = chunk.collect(); let feats: Vec = chunk .par_iter() + .filter(|r| fire_coords_fit_seq(r)) .map(|r| FireFeats::new(r, fire_opts)) .collect(); feats.iter().for_each(|f| { diff --git a/src/utils/fire.rs b/src/utils/fire.rs index ba9b13401..cd37ce9bd 100644 --- a/src/utils/fire.rs +++ b/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); 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/regression/fire.rs b/tests/regression/fire.rs index f2c490f75..8ff6fa712 100644 --- a/tests/regression/fire.rs +++ b/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 { out.lines() .map(|l| l.split('\t').nth(9).unwrap().parse().unwrap())