diff --git a/src/cli/fire_opts.rs b/src/cli/fire_opts.rs index 63160359..47fd8cf8 100644 --- a/src/cli/fire_opts.rs +++ b/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 diff --git a/src/fiber.rs b/src/fiber.rs index f584be7b..7b16bfb7 100644 --- a/src/fiber.rs +++ b/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}"), diff --git a/src/utils.rs b/src/utils.rs index d4d33260..b7b38f59 100644 --- a/src/utils.rs +++ b/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; diff --git a/src/utils/fire.rs b/src/utils/fire.rs index 7ca70d79..ba9b1340 100644 --- a/src/utils/fire.rs +++ b/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, fire_feats: Vec<(i64, i64, Vec)>, } @@ -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; diff --git a/src/utils/platform.rs b/src/utils/platform.rs new file mode 100644 index 00000000..10d94216 --- /dev/null +++ b/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 . +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 . +const ONT_READ_TAGS: [&[u8; 2]; 5] = [b"mx", b"ch", b"st", b"du", b"ts"]; + +lazy_static! { + /// PacBio read names are `//`, e.g. + /// `m84081_231207_200206_s1/240456095/ccs`. + /// See . + 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 . + 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 `//...` 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); + } +}