feat: auto-detect sequencing platform per read in FIRE

src/cli/fire_opts.rs

@@ -9,8 +9,11 @@ pub struct FireOptions {
     /// Output file (BAM by default, table of MSP features if `--feats-to-text` is used, and bed9 + if `--extract`` is used)
     #[clap(default_value = "-")]
     pub out: String,
-    /// Use a ONT heuristic adjustment for FIRE calling.
+    /// Force the ONT heuristic adjustment for every read in FIRE calling.
     /// This adjusts the observed number of m6A counts by adding pseudo counts to account for the single stranded nature of ONT data.
+    /// ONT reads are auto-detected per read (via aux tags / read name), so this
+    /// flag is only needed to override detection; reads matching no known
+    /// platform convention default to PacBio.
     #[clap(long, env)]
     pub ont: bool,
     /// Use a human model for FIRE calling

src/fiber.rs

@@ -183,6 +183,12 @@ impl FiberseqData {
         }
     }
 
+    /// Detect the sequencing platform (PacBio vs ONT) for this read from its
+    /// aux tags, falling back to the read name if those were stripped.
+    pub fn platform(&self) -> crate::utils::platform::SeqPlatform {
+        crate::utils::platform::platform_from_record(&self.record)
+    }
+
     pub fn get_hp(&self) -> String {
         match self.record.aux(b"HP") {
             Ok(Aux::U8(v)) => format!("H{v}"),

src/utils.rs

@@ -8,5 +8,6 @@ pub mod input_bam;
 pub mod ma_io;
 pub mod nucleosome;
 pub mod panspec;
+pub mod platform;
 // test modules for expressions
 pub mod ftexpression;

src/utils/fire.rs

@@ -1,6 +1,7 @@
 use crate::cli::FireOptions;
 use crate::fiber::FiberseqData;
 use crate::utils::bamannotations::AnnotationTypeView;
+use crate::utils::platform::SeqPlatform;
 use crate::*;
 use anyhow;
 use derive_builder::Builder;
@@ -172,6 +173,9 @@ pub struct FireFeats<'a> {
     frac_m6a: f32,
     //frac_m6a_in_msps: f32,
     fire_opts: &'a FireOptions,
+    /// Whether the ONT single-strand heuristic applies to this read. Set from
+    /// `--ont` (forces it for all reads) or per-read platform detection.
+    is_ont: bool,
     seq: Vec<u8>,
     fire_feats: Vec<(i64, i64, Vec<f32>)>,
 }
@@ -181,6 +185,11 @@ impl<'a> FireFeats<'a> {
         let seq_len = rec.record.seq_len();
         let seq = rec.record.seq().as_bytes();
 
+        // Apply the ONT single-strand m6A heuristic when the user forces it
+        // (`--ont`) or when this read is detected as ONT. Reads matching no
+        // known platform convention default to PacBio (no heuristic).
+        let is_ont = fire_opts.ont || matches!(rec.platform(), SeqPlatform::Ont);
+
         let mut rtn = Self {
             rec,
             m6a_view: rec.m6a(),
@@ -190,6 +199,7 @@ impl<'a> FireFeats<'a> {
             frac_m6a: 0.0,
             //frac_m6a_in_msps,
             fire_opts,
+            is_ont,
             seq,
             fire_feats: vec![],
         };
@@ -204,7 +214,7 @@ impl<'a> FireFeats<'a> {
         };
 
         rtn.get_fire_features();
-        if rtn.fire_opts.ont {
+        if rtn.is_ont {
             rtn.validate_that_ont_is_single_strand();
         }
         rtn
@@ -246,7 +256,7 @@ impl<'a> FireFeats<'a> {
         let mut m6a_count = self.m6a_view.count_query_in(start, end);
 
         // estimate what the count would be if we sequenced the other strand
-        if self.fire_opts.ont {
+        if self.is_ont {
             let mut sequenced_bp = self.get_bp_count(start, end, b'A');
             let mut un_sequenced_bp = self.get_bp_count(start, end, b'T');
             if self.rec.record.is_reverse() {
@@ -318,7 +328,7 @@ impl<'a> FireFeats<'a> {
         if msp_len < self.fire_opts.min_msp_length_for_positive_fire_call {
             return vec![];
         }
-        let ccs_passes = if self.fire_opts.ont { 4.0 } else { self.rec.ec };
+        let ccs_passes = if self.is_ont { 4.0 } else { self.rec.ec };
 
         // find the 100bp window within the range with the most m6a
         let mut max_m6a_count = 0;

src/utils/platform.rs

@@ -0,0 +1,174 @@
+//! Detect the sequencing platform (PacBio vs ONT) per read, so downstream
+//! commands can auto-select models/parameters. Aux tags are checked first;
+//! when they've been stripped, the read name is used as a fallback.
+
+use lazy_static::lazy_static;
+use regex::Regex;
+use rust_htslib::bam::Record;
+
+/// PacBio-only per-read aux tags: `zm` ZMW, `np` passes, `rq` read quality,
+/// `sn` SNR. Never emitted by ONT or re-used by fibertools.
+/// See <https://pacbiofileformats.readthedocs.io/en/13.1/BAM.html#use-of-read-tags-for-per-read-information>.
+const PACBIO_READ_TAGS: [&[u8; 2]; 4] = [b"zm", b"np", b"rq", b"sn"];
+
+/// ONT-only per-read aux tags: `mx` mux, `ch` channel, `st` start time,
+/// `du` duration, `ts` trimmed samples. Excludes `ns`/`nl`/`fn`/`rn`, which
+/// fibertools re-uses for its own annotation arrays.
+/// See <https://nanoporetech.github.io/ont-output-specifications/latest/read_formats/bam/#read-tags>.
+const ONT_READ_TAGS: [&[u8; 2]; 5] = [b"mx", b"ch", b"st", b"du", b"ts"];
+
+lazy_static! {
+    /// PacBio read names are `<movie>/<zmw>/<type>`, e.g.
+    /// `m84081_231207_200206_s1/240456095/ccs`.
+    /// See <https://pacbiofileformats.readthedocs.io/en/13.1/BAM.html#qname-convention>.
+    static ref PACBIO_QNAME: Regex = Regex::new(r"^m\w+/\d+/").unwrap();
+    /// ONT read names are the per-read UUID (`read_id`), e.g.
+    /// `d78c9b31-cec7-4381-bc1c-5f1513d070b9`.
+    /// See <https://nanoporetech.github.io/ont-output-specifications/latest/read_formats/bam/>.
+    static ref ONT_QNAME: Regex = Regex::new(
+        r"^[0-9a-fA-F]{8}-[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{12}$"
+    ).unwrap();
+}
+
+/// The sequencing platform a single read originated from.
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum SeqPlatform {
+    PacBio,
+    Ont,
+    /// No usable signal on this read.
+    Unknown,
+}
+
+/// Detect the platform of a single read: aux tags first, then the read name as
+/// a fallback for reads whose basecaller tags have been stripped.
+pub fn platform_from_record(rec: &Record) -> SeqPlatform {
+    match platform_from_aux(rec) {
+        SeqPlatform::Unknown => platform_from_qname(rec.qname()),
+        known => known,
+    }
+}
+
+/// Detect from per-read aux tag signatures. Presence is a positive signal only;
+/// absence means nothing (many tags are optional/flag-gated).
+pub fn platform_from_aux(rec: &Record) -> SeqPlatform {
+    let has_any = |tags: &[&[u8; 2]]| tags.iter().any(|t| rec.aux(&t[..]).is_ok());
+    let pacbio = has_any(&PACBIO_READ_TAGS);
+    let ont = has_any(&ONT_READ_TAGS);
+    match (pacbio, ont) {
+        (true, false) => SeqPlatform::PacBio,
+        (false, true) => SeqPlatform::Ont,
+        // Neither, or contradictory: don't guess.
+        _ => SeqPlatform::Unknown,
+    }
+}
+
+/// Detect from the read-name convention: PacBio `<movie>/<zmw>/...` vs an ONT
+/// UUID. Survives aux-tag stripping.
+pub fn platform_from_qname(qname: &[u8]) -> SeqPlatform {
+    let name = String::from_utf8_lossy(qname);
+    if PACBIO_QNAME.is_match(&name) {
+        SeqPlatform::PacBio
+    } else if ONT_QNAME.is_match(&name) {
+        SeqPlatform::Ont
+    } else {
+        SeqPlatform::Unknown
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use rust_htslib::bam::{self, Read};
+
+    fn first_record(path: &str) -> Record {
+        let mut reader = bam::Reader::from_path(path).unwrap();
+        let mut rec = Record::new();
+        reader
+            .read(&mut rec)
+            .expect("fixture has at least one record")
+            .expect("read ok");
+        rec
+    }
+
+    const PACBIO_FIXTURES: [&str; 5] = [
+        "tests/data/revio.bam",
+        "tests/data/two_two.bam",
+        "tests/data/three_two.bam",
+        "tests/data/all.bam",
+        "tests/data/ctcf.bam",
+    ];
+    const ONT_FIXTURE: &str = "tests/data/ONT.NAPA.bam";
+
+    #[test]
+    fn aux_detects_pacbio() {
+        for f in PACBIO_FIXTURES {
+            assert_eq!(
+                platform_from_aux(&first_record(f)),
+                SeqPlatform::PacBio,
+                "{f}"
+            );
+        }
+    }
+
+    #[test]
+    fn aux_detects_ont() {
+        assert_eq!(
+            platform_from_aux(&first_record(ONT_FIXTURE)),
+            SeqPlatform::Ont
+        );
+    }
+
+    #[test]
+    fn qname_detects_pacbio() {
+        for f in PACBIO_FIXTURES {
+            assert_eq!(
+                platform_from_qname(first_record(f).qname()),
+                SeqPlatform::PacBio,
+                "{f}"
+            );
+        }
+    }
+
+    #[test]
+    fn qname_detects_ont() {
+        assert_eq!(
+            platform_from_qname(first_record(ONT_FIXTURE).qname()),
+            SeqPlatform::Ont
+        );
+    }
+
+    #[test]
+    fn qname_recovers_platform_after_aux_stripping() {
+        // Strip the discriminator aux tags; the read name must still classify.
+        let mut rec = first_record(ONT_FIXTURE);
+        for t in ONT_READ_TAGS {
+            rec.remove_aux(t).ok();
+        }
+        assert_eq!(platform_from_aux(&rec), SeqPlatform::Unknown);
+        assert_eq!(platform_from_record(&rec), SeqPlatform::Ont);
+
+        let mut rec = first_record("tests/data/revio.bam");
+        for t in PACBIO_READ_TAGS {
+            rec.remove_aux(t).ok();
+        }
+        assert_eq!(platform_from_aux(&rec), SeqPlatform::Unknown);
+        assert_eq!(platform_from_record(&rec), SeqPlatform::PacBio);
+    }
+
+    #[test]
+    fn qname_patterns_and_non_matches() {
+        assert_eq!(
+            platform_from_qname(b"m84081_231207_200206_s1/240456095/ccs"),
+            SeqPlatform::PacBio
+        );
+        assert_eq!(
+            platform_from_qname(b"m54329U_210810_004956/120523046/0_5000"),
+            SeqPlatform::PacBio
+        );
+        assert_eq!(
+            platform_from_qname(b"d78c9b31-cec7-4381-bc1c-5f1513d070b9"),
+            SeqPlatform::Ont
+        );
+        assert_eq!(platform_from_qname(b"read_12345"), SeqPlatform::Unknown);
+    }
+}