fix: MA-refactor fallout — qc, track-decorators, filter expressions, mock-fire, stale legacy tags

molecular-annotation/README.md

@@ -16,10 +16,9 @@ molecular-annotation = { path = ".", features = ["htslib"] }
 
 ## Encoding Formats
 
-Two encoding formats are available via feature flags:
-
-- **`inline-lengths`** (default): Lengths in MA string: `MA:Z:1000;nuc+:100-50,200-60`
-- **`separate-lengths`**: Lengths in AL array: `MA:Z:1000;nuc+:100,200` + `AL:B:I,50,60`
+Lengths are encoded inline in the MA string (`MA:Z:1000;nuc+:100-50,200-60`).
+The retired separate-length encoding (`AL` array) is stripped on write and
+never emitted.
 
 ## Documentation
 

molecular-annotation/docs/mm-ml-per-group-passthrough.md

@@ -12,7 +12,7 @@ but non-canonical encodings (grouped multi-code `C+mh`, `N+a` wildcards).
 Today we work around this with two write functions, and the caller picks based
 on intent:
 
-- `write_record` / `to_record` — writes only MA-family tags (`MA`/`AL`/`AQ`/`AN`),
+- `write_record` / `to_record` — writes only MA-family tags (`MA`/`AQ`/`AN`),
   leaves the record's MM/ML bytes untouched. Used by structural editors
   (fire, add-nucs, footprint, pileup, extract, center, convert-tags) via
   `FiberseqData::serialize_annotations`.

molecular-annotation/examples/fiberseq_to_ma.rs

@@ -1,10 +1,7 @@
 //! Convert fiber-seq BAM tags (ns/nl, as/al/aq) to MolecularAnnotations format
 //!
 //! This example reads a CRAM/BAM file with fiber-seq annotations and converts
-//! them to the MA/AQ tag format (inline lengths, the default).
-//!
-//! It also writes M2/AL tags with the alternative separate format for compression
-//! comparison testing.
+//! them to the MA/AQ tag format (lengths are inline in the MA string).
 //!
 //! Run with:
 //!   cargo run --example fiberseq_to_ma

molecular-annotation/src/lib.rs

@@ -1,7 +1,7 @@
 //! Molecular Annotation Library
 //!
 //! This library provides types and functions for working with molecular annotations
-//! according to the MA/AL/AQ/AN tag specification for SAM/BAM files.
+//! according to the MA/AQ/AN tag specification for SAM/BAM files.
 //!
 //! # Coordinate Conventions
 //!

src/fiber.rs

@@ -157,6 +157,25 @@ impl FiberseqData {
         AnnotationTypeView::new(&self.annotations, FIRE_TYPE)
     }
 
+    /// Per-MSP FIRE quals, BAM-orient ascending (aligned with the `msp()`
+    /// view's accessors). FIRE quals live on the `fire` annotation type (the
+    /// subset of MSPs called as FIREs), not on the MSPs themselves — overlay
+    /// them onto the matching MSPs, keyed by BAM-orient query start (fire
+    /// entries are exact interval copies of their source MSPs). MSPs without
+    /// a fire entry fall back to their own qual, which is 0 for MA-era
+    /// records (FDR 100) and the legacy `aq` value for pre-MA BAMs.
+    pub fn msp_fire_quals(&self) -> Vec<u8> {
+        let fire = self.fire();
+        let fire_quals: std::collections::HashMap<i64, u8> =
+            fire.starts().into_iter().zip(fire.qual()).collect();
+        let msp = self.msp();
+        msp.qual()
+            .into_iter()
+            .zip(msp.starts())
+            .map(|(q, s)| fire_quals.get(&s).copied().unwrap_or(q))
+            .collect()
+    }
+
     /// Flush `self.annotations` onto the record's MA-family aux tags. The
     /// single write path for subcommands that edit nuc/msp/fire annotations;
     /// call this, then hand the record to the BAM writer.

src/subcommands/decorator.rs

@@ -166,12 +166,13 @@ pub fn fire_decorators(fiber: &FiberseqData) -> Vec<Decorator<'_>> {
         map.insert(color, vec![]);
     }
 
-    // FIRE is its own annotation type — color decorators by FIRE precision
-    // directly. Non-FIRE MSPs are not decorated here
-    let fire = fiber.fire();
-    let ref_starts = fire.reference_starts();
-    let ref_lengths = fire.reference_lengths();
-    let quals = fire.qual();
+    // FIRE quals live on the `fire` annotation type — overlay them onto the
+    // MSPs so every MSP is decorated: called FIREs by their precision, and
+    // non-FIRE MSPs as LINKER (qual 0 -> FDR 100), matching pre-MA output.
+    let msp = fiber.msp();
+    let ref_starts = msp.reference_starts();
+    let ref_lengths = msp.reference_lengths();
+    let quals = fiber.msp_fire_quals();
 
     for ((pos, length), qual) in ref_starts.iter().zip(ref_lengths.iter()).zip(quals.iter()) {
         if let (Some(p), Some(l)) = (pos, length) {
@@ -182,8 +183,11 @@ pub fn fire_decorators(fiber: &FiberseqData) -> Vec<Decorator<'_>> {
             map.get_mut(&fire_color).unwrap().push((p, l));
         }
     }
+    // iterate in FIRE_COLORS order so output is deterministic (HashMap
+    // iteration order is randomized per process)
     let mut rtn = vec![];
-    for (color, values) in map.into_iter() {
+    for color in FIRE_COLORS.iter().map(|(_, color)| color) {
+        let values = map.remove(&color).unwrap();
         let (starts, lengths): (Vec<Option<i64>>, Vec<Option<i64>>) = values.into_iter().unzip();
         let el_type = if *color == LINKER_COLOR {
             "LINKER"

src/subcommands/fire.rs

@@ -128,23 +128,7 @@ pub fn fire_to_bed9(fire_opts: &FireOptions, bam: &mut bam::Reader) -> Result<()
         let msp_ends = msp.reference_ends();
         let nuc_ends = nuc.reference_ends();
         let end_iter = msp_ends.iter().chain(nuc_ends.iter());
-        // FIRE quals live on the `fire` annotation type (the subset of MSPs
-        // called as FIREs), not on the MSPs themselves. Overlay them onto the
-        // matching MSPs (keyed by molecular start) so extraction reports the
-        // called FDR; MSPs without a fire entry fall back to their own qual
-        // (nonzero only for legacy-tag BAMs) and so report an FDR of 100.
-        let fire = rec.fire();
-        let fire_quals: std::collections::HashMap<i64, u8> = fire
-            .starts()
-            .into_iter()
-            .zip(fire.qual().into_iter())
-            .collect();
-        let msp_qual: Vec<u8> = msp
-            .qual()
-            .into_iter()
-            .zip(msp.starts().into_iter())
-            .map(|(q, s)| fire_quals.get(&s).copied().unwrap_or(q))
-            .collect();
+        let msp_qual = rec.msp_fire_quals();
         let nuc_qual = nuc.qual();
         let qual_iter = msp_qual.iter().chain(nuc_qual.iter());
         let n_msps = msp_starts.len();

src/subcommands/mock_fire.rs

@@ -129,13 +129,17 @@ fn create_mock_fire_record(
         quals.push(quality);
     }
 
-    // Build MA-spec annotations and emit. Mock FIRE only produces MSPs with
-    // the FIRE precision (Q-scaled) quality; no nucleosomes are populated.
+    // Build MA-spec annotations and emit. No nucleosomes are populated.
+    // Real fire-scored records carry every interval as an `msp` annotation
+    // (no quals) with the called subset duplicated into the `fire` type, so
+    // emit both — fire-only records are invisible to msp-driven consumers
+    // (pileup msp_coverage, ft extract, qc, filter expressions).
     // Using `ma_io::write_record_with_basemods` directly (rather than the
     // `FiberseqData::serialize_annotations` path used elsewhere) is intentional
     // here — we're synthesizing a record from BED, so there's no `FiberseqData`
     // to round-trip through.
     let mut annot = MolecularAnnotations::from_record(&record);
+    ma_io::add_msp_annotations(&mut annot, &starts, &lengths, None);
     ma_io::add_fire_annotations(&mut annot, &starts, &lengths, &quals);
     ma_io::write_record_with_basemods(&mut record, &annot);
 

src/subcommands/qc.rs

@@ -243,13 +243,12 @@ impl<'a> QcStats<'a> {
     fn m6a_per_msp(&mut self, fiber: &fiber::FiberseqData) {
         let msp = fiber.msp();
         let m6a_starts = fiber.m6a().starts();
-        for annotation in &msp {
+        // FIRE quals live on the `fire` annotation type, not the MSPs —
+        // msp_fire_quals() overlays them (both are BAM-orient ascending).
+        let quals = fiber.msp_fire_quals();
+        for (annotation, qual) in msp.infos().iter().zip(quals.into_iter()) {
             let st = annotation.query_start as i64;
             let en = annotation.query_end as i64;
-            let qual = crate::utils::bamannotations::primary_qual(
-                annotation.qualities,
-                annotation.type_name,
-            );
             let is_fire = qual >= 230;
             let msp_size = en - st;
             let m6a_count = m6a_starts

src/utils/ftexpression.rs

@@ -150,55 +150,62 @@ pub fn apply_filter_fsd(fsd: &mut FiberseqData, filt: &FiberFilters) -> Result<(
                     other => anyhow::bail!("Unknown feature name: {}", other),
                 };
                 match parser.fn_name.as_str() {
-                    "len" => fsd.annotations.retain(type_name, |a| {
-                        len(a.length as i64, &parser.op, &parser.threshold)
-                    }),
+                    "len" => {
+                        fsd.annotations.retain(type_name, |a| {
+                            len(a.length as i64, &parser.op, &parser.threshold)
+                        });
+                        // fire annotations are interval copies of their source
+                        // MSPs — apply the same predicate so they are not
+                        // orphaned when their parent MSP is removed.
+                        if type_name == "msp" {
+                            fsd.annotations.retain("fire", |a| {
+                                len(a.length as i64, &parser.op, &parser.threshold)
+                            });
+                        }
+                    }
                     "qual" if type_name == "msp" => {
                         // qual(msp) historically meant "filter MSPs by FIRE
                         // precision" (legacy fibertools wrote FIRE precisions
                         // onto the MSP `aq` tag). Post-MA, that quality lives
-                        // on the `fire` annotation type instead. Collect
-                        // per-MSP precisions in *molecular* order to align
-                        // with retain's iteration, then drop msp and fire
-                        // in lockstep.
+                        // on the `fire` annotation type, which holds the p>0
+                        // SUBSET of MSPs as exact interval copies. Overlay the
+                        // fire quals onto the MSPs keyed by molecular start;
+                        // MSPs without a fire entry fall back to their own
+                        // qual (legacy read path) and so evaluate as 0. Drop
+                        // msp and fire entries by the same kept-start set so
+                        // the two types stay paired.
                         let primary = crate::utils::bamannotations::primary_qual;
-                        let mol_quals: Vec<u8> = if let Some(f) =
-                            fsd.annotations.get_type("fire").filter(|f| {
-                                fsd.annotations
-                                    .get_type("msp")
-                                    .is_some_and(|m| m.annotations.len() == f.annotations.len())
-                            }) {
-                            f.annotations
-                                .iter()
-                                .map(|a| primary(&a.qualities, "fire"))
-                                .collect()
-                        } else if let Some(m) = fsd.annotations.get_type("msp") {
-                            m.annotations
-                                .iter()
-                                .map(|a| primary(&a.qualities, "msp"))
-                                .collect()
-                        } else {
-                            Vec::new()
-                        };
-                        let keep: Vec<bool> = mol_quals
-                            .iter()
-                            .map(|q| qual(*q, &parser.op, &parser.threshold))
-                            .collect();
-                        let has_fire = fsd.annotations.get_type("fire").is_some();
-                        let mut i = 0;
-                        fsd.annotations.retain("msp", |_| {
-                            let k = keep.get(i).copied().unwrap_or(false);
-                            i += 1;
-                            k
-                        });
-                        if has_fire {
-                            let mut i = 0;
-                            fsd.annotations.retain("fire", |_| {
-                                let k = keep.get(i).copied().unwrap_or(false);
-                                i += 1;
-                                k
-                            });
-                        }
+                        let fire_quals: std::collections::HashMap<u32, u8> = fsd
+                            .annotations
+                            .get_type("fire")
+                            .map(|f| {
+                                f.annotations
+                                    .iter()
+                                    .map(|a| (a.start, primary(&a.qualities, "fire")))
+                                    .collect()
+                            })
+                            .unwrap_or_default();
+                        let keep_starts: std::collections::HashSet<u32> = fsd
+                            .annotations
+                            .get_type("msp")
+                            .map(|m| {
+                                m.annotations
+                                    .iter()
+                                    .filter(|a| {
+                                        let q = fire_quals
+                                            .get(&a.start)
+                                            .copied()
+                                            .unwrap_or_else(|| primary(&a.qualities, "msp"));
+                                        qual(q, &parser.op, &parser.threshold)
+                                    })
+                                    .map(|a| a.start)
+                                    .collect()
+                            })
+                            .unwrap_or_default();
+                        fsd.annotations
+                            .retain("msp", |a| keep_starts.contains(&a.start));
+                        fsd.annotations
+                            .retain("fire", |a| keep_starts.contains(&a.start));
                     }
                     "qual" => fsd.annotations.retain(type_name, |a| {
                         qual(

src/utils/ma_io.rs

@@ -96,7 +96,7 @@ fn merge_missing_types(dst: &mut MolecularAnnotations, src: MolecularAnnotations
     }
 }
 
-/// Writes MA-family tags (MA/AL/AQ/AN) to a BAM record, **preserving the
+/// Writes MA-family tags (MA/AQ/AN) to a BAM record, **preserving the
 /// record's existing MM/ML bytes**.
 ///
 /// # Which write function do I call?

tests/regression.rs

@@ -4,12 +4,16 @@ mod center;
 mod common;
 #[path = "regression/convert_tags.rs"]
 mod convert_tags;
+#[path = "regression/decorator.rs"]
+mod decorator;
 #[path = "regression/extract.rs"]
 mod extract;
 #[path = "regression/fire.rs"]
 mod fire;
 #[path = "regression/footprint.rs"]
 mod footprint;
+#[path = "regression/mock_fire.rs"]
+mod mock_fire;
 #[path = "regression/pileup.rs"]
 mod pileup;
 #[path = "regression/predict_m6a.rs"]

tests/regression/decorator.rs

@@ -0,0 +1,35 @@
+use super::common::{fixture, run};
+use tempfile::NamedTempFile;
+
+// Every MSP must produce a decoration: called FIREs by their precision and
+// non-FIRE (precision-0) MSPs as LINKER, matching pre-MA output. FIRE quals
+// live on the `fire` annotation type, so the decorator must overlay them
+// onto the MSPs rather than iterate the fire type alone.
+#[test]
+fn track_decorators_emit_fire_and_linker() {
+    let scored = NamedTempFile::with_suffix(".bam").unwrap();
+    run(&[
+        "fire",
+        fixture("all.bam").to_str().unwrap(),
+        scored.path().to_str().unwrap(),
+    ]);
+    let bed12 = NamedTempFile::with_suffix(".bed").unwrap();
+    let out = run(&[
+        "track-decorators",
+        "--bed12",
+        bed12.path().to_str().unwrap(),
+        scored.path().to_str().unwrap(),
+    ]);
+    let count = |el: &str| {
+        out.lines()
+            .filter(|l| l.split('\t').any(|f| f == el))
+            .count()
+    };
+    let (fire, linker) = (count("FIRE"), count("LINKER"));
+    assert!(fire > 0, "no FIRE decorations emitted");
+    assert!(
+        linker > fire,
+        "expected precision-0 MSPs (the majority) to decorate as LINKER; got {linker} LINKER vs {fire} FIRE"
+    );
+    insta::assert_snapshot!(out);
+}