From b35e92ca809b953fbfb2c118e9e4330e7e1893fe Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Wed, 5 Aug 2026 00:03:37 -0700 Subject: [PATCH 1/3] Apply indexed edge filters consistently --- modkit-core/src/mod_bam.rs | 2 +- modkit-core/src/sample_probs/mod.rs | 280 +++++++++++++++++++++++----- 2 files changed, 238 insertions(+), 44 deletions(-) diff --git a/modkit-core/src/mod_bam.rs b/modkit-core/src/mod_bam.rs index ef6a80e7..460d7968 100644 --- a/modkit-core/src/mod_bam.rs +++ b/modkit-core/src/mod_bam.rs @@ -1665,7 +1665,7 @@ pub fn base_mod_probs_from_record( extract_mod_probs(record, &forward_seq, &mm, &ml, &converter) } -#[derive(new, Debug)] +#[derive(new, Debug, Clone, Copy)] pub struct EdgeFilter { pub edge_filter_start: usize, pub edge_filter_end: usize, diff --git a/modkit-core/src/sample_probs/mod.rs b/modkit-core/src/sample_probs/mod.rs index a0080af9..8ce7e145 100644 --- a/modkit-core/src/sample_probs/mod.rs +++ b/modkit-core/src/sample_probs/mod.rs @@ -841,8 +841,26 @@ pub(crate) struct ProbsExtractor { sample: bool, sample_frac: f64, motif_bases: [DnaBase; 4], - edge_filter_start: usize, - edge_filter_end: usize, + edge_filter: Option, +} + +#[inline] +fn indexed_query_position_is_kept( + edge_filter: Option<&EdgeFilter>, + query_position: usize, + read_length: usize, + reverse: bool, +) -> bool { + if query_position >= read_length { + return false; + } + let forward_position = + if reverse { read_length - 1 - query_position } else { query_position }; + edge_filter + .map(|filter| { + filter.keep_position(forward_position, read_length).unwrap_or(false) + }) + .unwrap_or(true) } impl ProbsExtractor { @@ -854,8 +872,7 @@ impl ProbsExtractor { should_count_record: bool, start_pos: u32, end_pos: u32, - edge_filter_start: usize, - edge_filter_end: usize, + edge_filter: Option, ) -> anyhow::Result> + use<'a>> { let reverse = record.is_reverse(); @@ -872,8 +889,12 @@ impl ProbsExtractor { *rpos >= start_pos && *rpos < end_pos }) .filter_ok(move |(qpos, _rpos, _mod_state)| { - (*qpos as usize) >= edge_filter_start - && (*qpos as usize) < (read_length - edge_filter_end) + indexed_query_position_is_kept( + edge_filter.as_ref(), + *qpos as usize, + read_length, + reverse, + ) }) .filter_map_ok(move |(qpos, rpos, mod_state)| { let rpos = rpos @@ -1059,7 +1080,7 @@ impl ExtractsMleProbs for ProbsExtractor { sample: bool, sample_frac: f64, motif_bases: [DnaBase; 4], - _edge_filter: Option<&EdgeFilter>, + edge_filter: Option<&EdgeFilter>, ) -> Self { let rng = SmallRng::seed_from_u64(seed); Self { @@ -1067,8 +1088,7 @@ impl ExtractsMleProbs for ProbsExtractor { sample, sample_frac, motif_bases, - edge_filter_start: 0usize, - edge_filter_end: 0usize, + edge_filter: edge_filter.copied(), } } @@ -1099,16 +1119,25 @@ impl ExtractsMleProbs for ProbsExtractor { let mut modbase_iter = BaseModsAdapter::<16>::new(record)?; let bases = modbase_iter.primary_bases_in_record(); byte_to_bool_positions(bases, records_with_base_mods); + let read_length = record.seq_len(); + let reverse = record.is_reverse(); loop { match modbase_iter.next_modified_position_no_thresh() { Ok(Some(mod_state)) => { - increment_base_counts( - mod_state, - explicit_canonical_probs, - base_hist, - base_totals, - &mut self.rng, - ); + if indexed_query_position_is_kept( + self.edge_filter.as_ref(), + mod_state.mod_position, + read_length, + reverse, + ) { + increment_base_counts( + mod_state, + explicit_canonical_probs, + base_hist, + base_totals, + &mut self.rng, + ); + } } Ok(None) => break, Err(e) => { @@ -1126,7 +1155,7 @@ impl ExtractsMleProbs for ProbsExtractor { sample: bool, sample_frac: f64, motif_bases: [DnaBase; 4], - _edge_filter: Option<&EdgeFilter>, + edge_filter: Option<&EdgeFilter>, ) -> Self { let rng = SmallRng::seed_from_u64(seed); Self { @@ -1134,8 +1163,7 @@ impl ExtractsMleProbs for ProbsExtractor { sample, sample_frac, motif_bases, - edge_filter_start: 0usize, - edge_filter_end: 0usize, + edge_filter: edge_filter.copied(), } } @@ -1165,17 +1193,26 @@ impl ExtractsMleProbs for ProbsExtractor { let mut modbase_iter = BaseModsAdapter::<16>::new(record)?; let bases = modbase_iter.primary_bases_in_record(); byte_to_bool_positions(bases, records_with_base_mods); + let read_length = record.seq_len(); + let reverse = record.is_reverse(); loop { match modbase_iter.next_modified_position_no_thresh() { Ok(Some(mod_state)) => { - increment_mods_counts( - mod_state, - &mut self.rng, - explicit_canonical_probs, - base_hist, - base_totals, - mods_hists, - ); + if indexed_query_position_is_kept( + self.edge_filter.as_ref(), + mod_state.mod_position, + read_length, + reverse, + ) { + increment_mods_counts( + mod_state, + &mut self.rng, + explicit_canonical_probs, + base_hist, + base_totals, + mods_hists, + ); + } } Ok(None) => break, Err(e) => { @@ -1196,17 +1233,12 @@ impl ExtractsMleProbs for ProbsExtractor { edge_filter: Option<&EdgeFilter>, ) -> Self { let rng = SmallRng::seed_from_u64(seed); - let edge_filter_start = - edge_filter.map(|ef| ef.edge_filter_start).unwrap_or(0usize); - let edge_filter_end = - edge_filter.map(|ef| ef.edge_filter_end).unwrap_or(0usize); Self { rng, sample, sample_frac, motif_bases, - edge_filter_start, - edge_filter_end, + edge_filter: edge_filter.copied(), } } @@ -1247,8 +1279,7 @@ impl ExtractsMleProbs for ProbsExtractor { should_count_record, start_pos, end_pos, - self.edge_filter_start, - self.edge_filter_end, + self.edge_filter, )?; for res in modstate_iter { @@ -1281,17 +1312,12 @@ impl ExtractsMleProbs for ProbsExtractor { edge_filter: Option<&EdgeFilter>, ) -> Self { let rng = SmallRng::seed_from_u64(seed); - let edge_filter_start = - edge_filter.map(|ef| ef.edge_filter_start).unwrap_or(0usize); - let edge_filter_end = - edge_filter.map(|ef| ef.edge_filter_end).unwrap_or(0usize); Self { rng, sample, sample_frac, motif_bases, - edge_filter_start, - edge_filter_end, + edge_filter: edge_filter.copied(), } } @@ -1332,8 +1358,7 @@ impl ExtractsMleProbs for ProbsExtractor { should_count_record, start_pos, end_pos, - self.edge_filter_start, - self.edge_filter_end, + self.edge_filter, )?; for res in modstate_iter { match res { @@ -1759,3 +1784,172 @@ fn byte_to_bool_positions(b: u8, agg: &mut [u32; SIZE]) { agg[i] = count.saturating_add(1u32); } } + +#[cfg(test)] +mod tests { + use bitvec::bitvec; + use rust_htslib::bam::record::{Aux, Cigar, CigarString}; + + use super::*; + + fn make_five_base_record(reverse: bool) -> bam::Record { + let mut record = bam::Record::new(); + let cigar = CigarString(vec![Cigar::Match(5)]); + let seq = if reverse { b"TTTTT" } else { b"AAAAA" }; + record.set(b"read", Some(&cigar), seq, &[255; 5]); + record.set_tid(0); + record.set_pos(0); + record.push_aux(b"MM", Aux::String("A+a?,0,0,0,0,0;")).unwrap(); + let ml = vec![201u8, 202, 203, 204, 205]; + record.push_aux(b"ML", Aux::ArrayU8((&ml).into())).unwrap(); + if reverse { + record.set_reverse(); + } + record + } + + fn retained_qualities( + record: &bam::Record, + edge_filter: Option<&EdgeFilter>, + ) -> Vec + where + ProbsExtractor: ExtractsMleProbs, + { + let mut extractor = >::new( + 42, + false, + 0.0, + [DnaBase::A; 4], + edge_filter, + ); + let chrom_coords = ChromCoordinates::new( + 0, + 0, + 5, + FocusPositions2::MaskedPositions { mask: bitvec![1; 10] }, + true, + ); + let mut hist = QualHist::default(); + let processed = + >::process_record( + &mut extractor, + record, + &chrom_coords, + &mut hist.explicit_canonical_probs, + &mut hist.hist, + &mut hist.base_totals, + &mut hist.mods_hists, + &mut hist.num_records_with_base_mods, + ) + .unwrap(); + assert!(processed); + + let mut qualities = hist + .hist + .iter() + .chain(hist.mods_hists.iter().map(|mod_hist| &mod_hist.hist)) + .flat_map(|counts| { + counts.iter().enumerate().flat_map(|(qual, count)| { + std::iter::repeat(qual as u8).take(*count as usize) + }) + }) + .collect::>(); + qualities.sort_unstable(); + qualities + } + + macro_rules! handler_results { + ($marker:ty, $forward:expr, $reverse:expr, $filter:expr) => { + ( + retained_qualities::<$marker>($forward, None), + retained_qualities::<$marker>($forward, Some($filter)), + retained_qualities::<$marker>($reverse, None), + retained_qualities::<$marker>($reverse, Some($filter)), + ) + }; + } + + #[test] + fn indexed_handlers_apply_asymmetric_edge_filter() { + let forward = make_five_base_record(false); + let reverse = make_five_base_record(true); + let edge_filter = EdgeFilter::new(1, 2, false); + let actual = [ + handler_results!( + AlignedBaseArgmaxProbs, + &forward, + &reverse, + &edge_filter + ), + handler_results!( + AlignedBaseAndModArgmaxProbs, + &forward, + &reverse, + &edge_filter + ), + handler_results!(BaseArgmaxProbs, &forward, &reverse, &edge_filter), + handler_results!( + BaseAndModArgmaxProbs, + &forward, + &reverse, + &edge_filter + ), + ]; + let all_qualities = vec![201, 202, 203, 204, 205]; + let retained_qualities = vec![202, 203]; + + for (forward_all, forward_filtered, reverse_all, reverse_filtered) in + actual + { + assert_eq!(forward_all, all_qualities); + assert_eq!(forward_filtered, retained_qualities); + assert_eq!(reverse_all, all_qualities); + assert_eq!(reverse_filtered, retained_qualities); + } + } + + #[test] + fn indexed_handlers_treat_overlong_edge_filter_as_empty() { + let record = make_five_base_record(false); + let edge_filter = EdgeFilter::new(1, 6, false); + let actual = [ + retained_qualities::(&record, Some(&edge_filter)), + retained_qualities::( + &record, + Some(&edge_filter), + ), + retained_qualities::( + &record, + Some(&edge_filter), + ), + retained_qualities::( + &record, + Some(&edge_filter), + ), + ]; + + for qualities in actual { + assert!(qualities.is_empty()); + } + } + + #[test] + fn indexed_query_filter_preserves_inversion() { + let edge_filter = EdgeFilter::new(1, 2, true); + let retained = |reverse| { + (0..5) + .filter(|qpos| { + indexed_query_position_is_kept( + Some(&edge_filter), + *qpos, + 5, + reverse, + ) + }) + .collect::>() + }; + + assert_eq!(retained(false), vec![0, 3, 4]); + assert_eq!(retained(true), vec![0, 1, 4]); + } +} From ced01a5533f260fecd268132788508ce47e7c5e6 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Thu, 6 Aug 2026 01:44:11 -0700 Subject: [PATCH 2/3] Add edge-filter threshold consistency regressions --- modkit-core/src/sample_probs/mod.rs | 203 ++++++++++++++++++++++++---- modkit/tests/test_extract.rs | 139 +++++++++++++++++++ 2 files changed, 312 insertions(+), 30 deletions(-) diff --git a/modkit-core/src/sample_probs/mod.rs b/modkit-core/src/sample_probs/mod.rs index 8ce7e145..d5cb84f1 100644 --- a/modkit-core/src/sample_probs/mod.rs +++ b/modkit-core/src/sample_probs/mod.rs @@ -1788,9 +1788,15 @@ fn byte_to_bool_positions(b: u8, agg: &mut [u32; SIZE]) { #[cfg(test)] mod tests { use bitvec::bitvec; - use rust_htslib::bam::record::{Aux, Cigar, CigarString}; + use rust_htslib::bam::{ + header::HeaderRecord, + record::{Aux, Cigar, CigarString}, + Format, Header, Read, Writer, + }; + use tempfile::tempdir; use super::*; + use crate::position_filter::{GenomeIntervals, Iv}; fn make_five_base_record(reverse: bool) -> bam::Record { let mut record = bam::Record::new(); @@ -1844,6 +1850,10 @@ mod tests { .unwrap(); assert!(processed); + qualities_from_hist(&hist) + } + + fn qualities_from_hist(hist: &QualHist) -> Vec { let mut qualities = hist .hist .iter() @@ -1858,6 +1868,55 @@ mod tests { qualities } + fn write_bam(bam_path: &std::path::Path, record: &bam::Record) { + let mut header = Header::new(); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", &"chr1"); + sq.push_tag(b"LN", &5); + header.push_record(&sq); + + let mut writer = + Writer::from_path(bam_path, &header, Format::Bam).unwrap(); + writer.write(record).unwrap(); + } + + fn all_positions_filter() -> StrandedPositionFilter<()> { + let intervals = + || GenomeIntervals::new(vec![Iv { start: 0, stop: 5, val: () }]); + let mut pos_positions = FxHashMap::default(); + pos_positions.insert(0, intervals()); + let mut neg_positions = FxHashMap::default(); + neg_positions.insert(0, intervals()); + StrandedPositionFilter { pos_positions, neg_positions } + } + + fn serial_retained_qualities( + record: &bam::Record, + edge_filter: Option<&EdgeFilter>, + with_position_filter: bool, + ) -> Vec { + let temp_dir = tempdir().unwrap(); + let bam_path = temp_dir.path().join("record.bam"); + write_bam(&bam_path, record); + let mut reader = bam::Reader::from_path(&bam_path).unwrap(); + let position_filter = with_position_filter.then(all_positions_filter); + let histogram = QualHist::from_records( + reader.records(), + position_filter, + None, + None, + None, + edge_filter, + false, + false, + false, + &MultiProgress::new(), + ) + .unwrap(); + + qualities_from_hist(&histogram) + } + macro_rules! handler_results { ($marker:ty, $forward:expr, $reverse:expr, $filter:expr) => { ( @@ -1910,46 +1969,130 @@ mod tests { #[test] fn indexed_handlers_treat_overlong_edge_filter_as_empty() { - let record = make_five_base_record(false); - let edge_filter = EdgeFilter::new(1, 6, false); + for reverse in [false, true] { + let record = make_five_base_record(reverse); + for inverted in [false, true] { + let edge_filter = EdgeFilter::new(1, 6, inverted); + let actual = [ + retained_qualities::( + &record, + Some(&edge_filter), + ), + retained_qualities::( + &record, + Some(&edge_filter), + ), + retained_qualities::( + &record, + Some(&edge_filter), + ), + retained_qualities::( + &record, + Some(&edge_filter), + ), + ]; + + for qualities in actual { + assert!(qualities.is_empty()); + } + } + } + } + + #[test] + fn indexed_handlers_preserve_asymmetric_filter_inversion() { + let forward = make_five_base_record(false); + let reverse = make_five_base_record(true); + let edge_filter = EdgeFilter::new(1, 2, true); let actual = [ - retained_qualities::(&record, Some(&edge_filter)), - retained_qualities::( - &record, - Some(&edge_filter), + handler_results!( + AlignedBaseArgmaxProbs, + &forward, + &reverse, + &edge_filter ), - retained_qualities::( - &record, - Some(&edge_filter), + handler_results!( + AlignedBaseAndModArgmaxProbs, + &forward, + &reverse, + &edge_filter ), - retained_qualities::( - &record, - Some(&edge_filter), + handler_results!(BaseArgmaxProbs, &forward, &reverse, &edge_filter), + handler_results!( + BaseAndModArgmaxProbs, + &forward, + &reverse, + &edge_filter ), ]; + let all_qualities = vec![201, 202, 203, 204, 205]; + let retained_qualities = vec![201, 204, 205]; - for qualities in actual { - assert!(qualities.is_empty()); + for (forward_all, forward_filtered, reverse_all, reverse_filtered) in + actual + { + assert_eq!(forward_all, all_qualities); + assert_eq!(forward_filtered, retained_qualities); + assert_eq!(reverse_all, all_qualities); + assert_eq!(reverse_filtered, retained_qualities); } } #[test] - fn indexed_query_filter_preserves_inversion() { - let edge_filter = EdgeFilter::new(1, 2, true); - let retained = |reverse| { - (0..5) - .filter(|qpos| { - indexed_query_position_is_kept( + fn serial_sampling_applies_edge_filters_in_molecular_orientation() { + let all_qualities = vec![201, 202, 203, 204, 205]; + let interior_qualities = vec![202, 203]; + let edge_qualities = vec![201, 204, 205]; + let normal = EdgeFilter::new(1, 2, false); + let inverted = EdgeFilter::new(1, 2, true); + + for reverse in [false, true] { + let record = make_five_base_record(reverse); + for with_position_filter in [false, true] { + assert_eq!( + serial_retained_qualities( + &record, + None, + with_position_filter, + ), + all_qualities, + ); + assert_eq!( + serial_retained_qualities( + &record, + Some(&normal), + with_position_filter, + ), + interior_qualities, + ); + assert_eq!( + serial_retained_qualities( + &record, + Some(&inverted), + with_position_filter, + ), + edge_qualities, + ); + } + } + } + + #[test] + fn serial_sampling_treats_overlong_edge_filters_as_empty() { + for reverse in [false, true] { + let record = make_five_base_record(reverse); + for with_position_filter in [false, true] { + for inverted in [false, true] { + let edge_filter = EdgeFilter::new(1, 6, inverted); + let qualities = serial_retained_qualities( + &record, Some(&edge_filter), - *qpos, - 5, - reverse, - ) - }) - .collect::>() - }; + with_position_filter, + ); - assert_eq!(retained(false), vec![0, 3, 4]); - assert_eq!(retained(true), vec![0, 1, 4]); + assert!(qualities.is_empty()); + } + } + } } } diff --git a/modkit/tests/test_extract.rs b/modkit/tests/test_extract.rs index 729153bb..24331fc6 100644 --- a/modkit/tests/test_extract.rs +++ b/modkit/tests/test_extract.rs @@ -3,11 +3,19 @@ use crate::common::{ }; use anyhow::{anyhow, Context}; use common::{check_legal_csv, run_modkit, ExtractFullRecord}; +use rust_htslib::bam::{ + self, + header::HeaderRecord, + record::{Aux, Cigar, CigarString}, + Format, Header, Writer, +}; +use serde::Deserialize; use std::collections::{HashMap, HashSet}; use std::fs::File; use std::io::BufRead; use std::io::BufReader; use std::path::{Path, PathBuf}; +use tempfile::tempdir; mod common; @@ -27,6 +35,137 @@ fn test_extract_help() { .unwrap(); } +fn make_probability_edge_filter_record( + read_name: &[u8], + mapped: bool, +) -> bam::Record { + let mut record = bam::Record::new(); + let cigar = mapped.then(|| CigarString(vec![Cigar::Match(9)])); + record.set(read_name, cigar.as_ref(), b"AAAAAAAAA", &[30; 9]); + record.push_aux(b"MM", Aux::String("A+a?,0,0,0,0,0,0,0,0,0;")).unwrap(); + let ml = vec![149u8, 150, 200, 200, 200, 200, 200, 150, 1]; + record.push_aux(b"ML", Aux::ArrayU8((&ml).into())).unwrap(); + + if mapped { + record.set_tid(0); + record.set_pos(0); + record.set_mapq(60); + } else { + record.set_tid(-1); + record.set_pos(-1); + record.set_unmapped(); + } + record +} + +fn write_probability_edge_filter_bam(bam_path: &Path) { + let mut header = Header::new(); + let mut hd = HeaderRecord::new(b"HD"); + hd.push_tag(b"VN", &"1.6"); + hd.push_tag(b"SO", &"coordinate"); + header.push_record(&hd); + let mut sq = HeaderRecord::new(b"SQ"); + sq.push_tag(b"SN", &"chr1"); + sq.push_tag(b"LN", &9); + header.push_record(&sq); + + let mut writer = Writer::from_path(bam_path, &header, Format::Bam).unwrap(); + writer + .write(&make_probability_edge_filter_record(b"mapped", true)) + .unwrap(); + writer + .write(&make_probability_edge_filter_record(b"unmapped", false)) + .unwrap(); +} + +#[derive(Debug, Deserialize, Eq, PartialEq, Ord, PartialOrd)] +struct ProbabilityEdgeFilterCall { + read_id: String, + forward_read_position: usize, + call_prob: String, + call_code: String, + fail: bool, +} + +fn assert_probability_edge_filter_calls(input_bam: &Path, output_tsv: &Path) { + run_modkit(&[ + "extract", + "calls", + input_bam.to_str().unwrap(), + output_tsv.to_str().unwrap(), + "--edge-filter", + "2", + "--invert-edge-filter", + "--filter-percentile", + "0.5", + "--threads", + "1", + "--io-threads", + "1", + "--suppress-progress", + "--force", + ]) + .unwrap(); + + let mut reader = csv::ReaderBuilder::new() + .delimiter(b'\t') + .from_path(output_tsv) + .unwrap(); + let mut calls = reader + .deserialize::() + .map(|result| result.unwrap()) + .collect::>(); + calls.sort_unstable(); + + let mut expected = ["mapped", "unmapped"] + .into_iter() + .flat_map(|read_id| { + [ + (0, "0.5839844", "a", true), + (1, "0.5878906", "a", false), + (7, "0.5878906", "a", false), + (8, "0.9941406", "-", false), + ] + .into_iter() + .map(move |(position, call_prob, call_code, fail)| { + ProbabilityEdgeFilterCall { + read_id: read_id.to_owned(), + forward_read_position: position, + call_prob: call_prob.to_owned(), + call_code: call_code.to_owned(), + fail, + } + }) + }) + .collect::>(); + expected.sort_unstable(); + + assert_eq!(calls, expected); +} + +#[test] +fn test_extract_inverted_edge_filter_stream_threshold_matches_output() { + let temp_dir = tempdir().unwrap(); + let unindexed_bam = temp_dir.path().join("unindexed.bam"); + write_probability_edge_filter_bam(&unindexed_bam); + assert_probability_edge_filter_calls( + &unindexed_bam, + &temp_dir.path().join("unindexed.tsv"), + ); +} + +#[test] +fn test_extract_inverted_edge_filter_indexed_unmapped_fallback_threshold() { + let temp_dir = tempdir().unwrap(); + let indexed_bam = temp_dir.path().join("indexed.bam"); + write_probability_edge_filter_bam(&indexed_bam); + bam::index::build(&indexed_bam, None, bam::index::Type::Bai, 1).unwrap(); + assert_probability_edge_filter_calls( + &indexed_bam, + &temp_dir.path().join("indexed.tsv"), + ); +} + fn parse_bed_file(fp: &PathBuf) -> HashMap> { let reader = BufReader::new(File::open(fp).unwrap()); reader From e043ada73a77c4a5bc156b8cc211a6184223ca30 Mon Sep 17 00:00:00 2001 From: SuhasSrinivasan <32346517+SuhasSrinivasan@users.noreply.github.com> Date: Thu, 6 Aug 2026 01:45:29 -0700 Subject: [PATCH 3/3] Apply inverted edge filters during probability sampling --- modkit-core/src/sample_probs/mod.rs | 35 ++++++++++++++--------------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/modkit-core/src/sample_probs/mod.rs b/modkit-core/src/sample_probs/mod.rs index d5cb84f1..d2d4a25f 100644 --- a/modkit-core/src/sample_probs/mod.rs +++ b/modkit-core/src/sample_probs/mod.rs @@ -403,10 +403,6 @@ impl QualHist { mapped_only: bool, multi_progress: &MultiProgress, ) -> anyhow::Result { - let edge_filter_start = - edge_filter.map(|ef| ef.edge_filter_start).unwrap_or(0usize); - let edge_filter_end = - edge_filter.map(|ef| ef.edge_filter_end).unwrap_or(0usize); let record_counter = multi_progress.add(get_ticker()); let mut rng = SmallRng::seed_from_u64(seed.unwrap_or(42u64)); record_counter.set_message("records processed"); @@ -438,11 +434,7 @@ impl QualHist { continue 'records; } let read_length = record.seq_len(); - let (left_edge_filter, right_edge_filter) = if record.is_reverse() { - (edge_filter_end, read_length.saturating_sub(edge_filter_start)) - } else { - (edge_filter_start, read_length.saturating_sub(edge_filter_end)) - }; + let reverse = record.is_reverse(); if let Some(sampling_frac) = sampling_frac { if !rng.gen_bool(sampling_frac) { continue 'records; @@ -475,8 +467,12 @@ impl QualHist { let aligned_pairs_iter = aligned_pairs_iter .into_iter() .filter_ok(move |(qpos, _rpos, _mod_state)| { - (*qpos as usize) >= left_edge_filter - && (*qpos as usize) < right_edge_filter + edge_filter_keeps_query_position( + edge_filter, + *qpos as usize, + read_length, + reverse, + ) }) .filter_ok(|(_, rpos, _)| { spf.contains(chrom_id, *rpos as u64, strand) @@ -525,9 +521,12 @@ impl QualHist { 'mods: loop { match modbase_iter.next_modified_position_no_thresh() { Ok(Some(mod_state)) => { - if mod_state.mod_position >= left_edge_filter - && mod_state.mod_position < right_edge_filter - { + if edge_filter_keeps_query_position( + edge_filter, + mod_state.mod_position, + read_length, + reverse, + ) { if collect_mod_probs { increment_mods_counts( mod_state, @@ -845,7 +844,7 @@ pub(crate) struct ProbsExtractor { } #[inline] -fn indexed_query_position_is_kept( +fn edge_filter_keeps_query_position( edge_filter: Option<&EdgeFilter>, query_position: usize, read_length: usize, @@ -889,7 +888,7 @@ impl ProbsExtractor { *rpos >= start_pos && *rpos < end_pos }) .filter_ok(move |(qpos, _rpos, _mod_state)| { - indexed_query_position_is_kept( + edge_filter_keeps_query_position( edge_filter.as_ref(), *qpos as usize, read_length, @@ -1124,7 +1123,7 @@ impl ExtractsMleProbs for ProbsExtractor { loop { match modbase_iter.next_modified_position_no_thresh() { Ok(Some(mod_state)) => { - if indexed_query_position_is_kept( + if edge_filter_keeps_query_position( self.edge_filter.as_ref(), mod_state.mod_position, read_length, @@ -1198,7 +1197,7 @@ impl ExtractsMleProbs for ProbsExtractor { loop { match modbase_iter.next_modified_position_no_thresh() { Ok(Some(mod_state)) => { - if indexed_query_position_is_kept( + if edge_filter_keeps_query_position( self.edge_filter.as_ref(), mod_state.mod_position, read_length,