Skip to content

Preserve pileup calls around ambiguous query bases #714

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
SuhasSrinivasan wants to merge 3 commits into from
Closed

Preserve pileup calls around ambiguous query bases #714

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
aa15346
Handle ambiguous bases in optimized mod scanning
SuhasSrinivasan Aug 5, 2026
85d62ea
Preserve pileup reads around ambiguous query bases
SuhasSrinivasan Aug 5, 2026
bb20328
Reject reference-equal pileup records
SuhasSrinivasan Aug 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 41 additions & 2 deletions modkit-core/src/pileup/base_mods_adapter.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ impl<'a, const SIZE: usize> BaseModsAdapter<'a, SIZE> {
b'C' | b'c' => agg[2] += 1u32,
b'G' | b'g' => agg[1] += 1u32,
b'T' | b't' => agg[0] += 1u32,
_ => unreachable!(),
_ => {}
}
agg
},
Expand Down Expand Up @@ -465,11 +465,23 @@ fn parse_int<const DELIM: u8, const END: u8>(bs: &[u8]) -> (u32, usize) {

fn base_complement(base: u8) -> u8 {
match base {
b'=' => b'=',
b'A' => b'T',
b'C' => b'G',
b'M' => b'K',
b'G' => b'C',
b'R' => b'Y',
b'S' => b'S',
b'V' => b'B',
b'T' => b'A',
_ => panic!("not allowed base"),
b'W' => b'W',
b'Y' => b'R',
b'H' => b'D',
b'K' => b'M',
b'D' => b'H',
b'B' => b'V',
b'N' => b'N',
_ => base,
}
}

Expand Down Expand Up @@ -622,6 +634,33 @@ mod base_mods_adapter_tests {
assert!(mod_state.is_none());
}

#[test]
fn test_ambiguous_query_bases_preserve_mm_scanning() {
for (seq, reverse, expected) in [
("CARNCC", false, vec![(0, 235), (4, 225)]),
("GNYGG", true, vec![(3, 225), (4, 235)]),
] {
let record =
make_record("C+m?,0,0;", &[20, 30], seq, None, reverse);
let mut scanner = BaseModsAdapter::<1>::new(&record).unwrap();
let states = std::iter::from_fn(|| {
scanner.next_modified_position([0f32; 4], &[]).unwrap()
})
.map(|state| {
assert_eq!(
state.primary_base,
crate::mod_base_code::DnaBase::C
);
assert!(!state.modified);
assert!(!state.inferred);
(state.mod_position, state.mod_qual)
})
.collect::<Vec<_>>();

assert_eq!(states, expected, "reverse={reverse}");
}
}

#[test]
fn test_forward() {
let seq = "ATCATCATTCCTACCGCTATAGCCT";
Expand Down
2 changes: 2 additions & 0 deletions modkit-core/src/pileup/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,7 @@ pub struct ModBasePileup2 {
pub(crate) interval_width: usize,
pub(crate) stride: usize,
pub(crate) failed_records: usize,
pub(crate) reference_equal_records: usize,
pub(crate) phased_feature_counts: [Vec<PileupFeatureCounts2>; 2],
}

Expand All @@ -182,6 +183,7 @@ impl ModBasePileup2 {
interval_width: 0,
stride: 0,
failed_records: 0,
reference_equal_records: 0,
phased_feature_counts: [Vec::new(), Vec::new()],
}
}
Expand Down
162 changes: 90 additions & 72 deletions modkit-core/src/pileup/pileup_processor.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,48 @@ where
) -> anyhow::Result<ModBasePileup2>;
}

#[inline]
fn get_query_base(record: &bam::Record, qpos: usize) -> Option<DnaBase> {
DnaBase::try_from(record.seq()[qpos]).ok().map(|base| {
if record.is_reverse() {
base.complement()
} else {
base
}
})
}

#[inline]
fn record_uses_reference_equal_bases(record: &bam::Record) -> bool {
let seq = record.seq();
// Each complete byte stores two bases. For odd-length sequences, only
// inspect the final high nibble because the low nibble is zero padding.
let complete_bytes = seq.len() / 2;
seq.encoded[..complete_bytes]
.iter()
.any(|packed| packed & 0x0f == 0 || packed >> 4 == 0)
|| (seq.len() % 2 == 1 && seq.encoded[complete_bytes] >> 4 == 0)
}

#[cfg(test)]
mod reference_equal_base_tests {
use super::*;

#[test]
fn odd_length_padding_is_not_a_reference_equal_base() {
let mut canonical = bam::Record::new();
canonical.set(b"canonical", None, b"ACGTA", &[30; 5]);
let seq = canonical.seq();
assert_eq!(seq.len(), 5);
assert_eq!(seq.encoded.last().unwrap() & 0x0f, 0);
assert!(!record_uses_reference_equal_bases(&canonical));

let mut reference_equal = bam::Record::new();
reference_equal.set(b"reference-equal", None, b"AC=TA", &[30; 5]);
assert!(record_uses_reference_equal_bases(&reference_equal));
}
}

pub(super) struct DnaPileupWorker<
T,
M,
Expand Down Expand Up @@ -146,6 +188,7 @@ impl<
mut pileup_space: ModBasePileup2,
) -> anyhow::Result<ModBasePileup2> {
let mut erred_records = 0usize;
let mut reference_equal_records = 0usize;
let chrom_tid = item.chrom_tid;
let start_pos = item.start_pos;
let end_pos = item.end_pos;
Expand Down Expand Up @@ -192,6 +235,16 @@ impl<
});

'records: for (record, hp) in records {
// BAM code zero (`=` in SAM) means "same as the reference". Its
// identity cannot be recovered from SEQ alone, so reject the
// record before adding any of its observations to the matrix.
if record_uses_reference_equal_bases(&record) {
erred_records = erred_records.saturating_add(1);
reference_equal_records =
reference_equal_records.saturating_add(1);
continue 'records;
}

if self.allow_non_primary && record_is_not_primary(&record) {
if validate_mn_tag_on_record(&record).is_err() {
erred_records = erred_records.saturating_add(1);
Expand Down Expand Up @@ -319,17 +372,8 @@ impl<
// at this aligned position, but there is no modification
// call
(Some(q), Some(mp)) if q < mp.mod_position => {
let base = {
let Ok(tmp) = DnaBase::try_from(record.seq()[q])
else {
erred_records = erred_records.saturating_add(1);
continue 'records;
};
if record.is_reverse() {
tmp.complement()
} else {
tmp
}
let Some(base) = get_query_base(&record, q) else {
continue 'pileup;
};
self.matrix
.incr_diff_call(rpos, base, ref_base, reverse, hp);
Expand Down Expand Up @@ -362,22 +406,16 @@ impl<
break 'overran;
} else {
assert!(pos > q);
let base = {
let tmp = DnaBase::try_from(
record.seq()[q],
)
.unwrap();
if record.is_reverse() {
tmp.complement()
} else {
tmp
}
mod_state = Some(ms);
let Some(base) =
get_query_base(&record, q)
else {
break 'overran;
};
self.matrix.incr_diff_call(
rpos, base, ref_base, reverse,
hp,
);
mod_state = Some(ms);
break 'overran;
}
}
Expand All @@ -388,20 +426,14 @@ impl<
continue 'records;
}
Ok(None) => {
let base = {
let tmp =
DnaBase::try_from(record.seq()[q])
.unwrap();
if record.is_reverse() {
tmp.complement()
} else {
tmp
}
mod_state = None;
let Some(base) = get_query_base(&record, q)
else {
break 'overran;
};
self.matrix.incr_diff_call(
rpos, base, ref_base, reverse, hp,
);
mod_state = None;
break 'overran;
}
}
Expand All @@ -412,14 +444,8 @@ impl<
self.matrix.incr_delete(rpos, reverse, hp);
}
(Some(q), None) => {
let base = {
let tmp =
DnaBase::try_from(record.seq()[q]).unwrap();
if record.is_reverse() {
tmp.complement()
} else {
tmp
}
let Some(base) = get_query_base(&record, q) else {
continue 'pileup;
};
self.matrix
.incr_diff_call(rpos, base, ref_base, reverse, hp);
Expand Down Expand Up @@ -461,6 +487,7 @@ impl<
pileup_space.position_feature_counts = combined_counts;
pileup_space.phased_feature_counts = [hp1, hp2];
pileup_space.failed_records = erred_records;
pileup_space.reference_equal_records = reference_equal_records;
Ok(pileup_space)
} else {
pileup_space.chrom_name = chrom_name;
Expand All @@ -471,6 +498,7 @@ impl<
.filter(|x| x.is_valid())
.collect();
pileup_space.failed_records = erred_records;
pileup_space.reference_equal_records = reference_equal_records;
Ok(pileup_space)
}
}
Expand Down Expand Up @@ -555,6 +583,7 @@ impl PileupWorker for GenericPileupWorker {
mut pileup_space: ModBasePileup2,
) -> anyhow::Result<ModBasePileup2> {
let mut erred_records = 0usize;
let mut reference_equal_records = 0usize;
let chrom_tid = chrom_coordinates.chrom_tid;
let start_pos = chrom_coordinates.start_pos;
let end_pos = chrom_coordinates.end_pos;
Expand Down Expand Up @@ -630,6 +659,15 @@ impl PileupWorker for GenericPileupWorker {
let mut pos_base_mod_call = Option::<BaseModCall>::None;
let mut mod_strand = Option::<Strand>::None;
'records: for record in records {
// Keep this preflight before modification parsing and tallying so
// a reference-equal query cannot contribute a partial record.
if record_uses_reference_equal_bases(&record) {
erred_records = erred_records.saturating_add(1);
reference_equal_records =
reference_equal_records.saturating_add(1);
continue 'records;
}

let reverse = record.is_reverse();
let record_strand = if record.is_reverse() {
Strand::Negative
Expand Down Expand Up @@ -755,17 +793,8 @@ impl PileupWorker for GenericPileupWorker {
continue 'pileup;
}
(Some(q), Some(mp)) if q < mp => {
let base = {
let Ok(tmp) = DnaBase::try_from(record.seq()[q])
else {
erred_records = erred_records.saturating_add(1);
continue 'records;
};
if record.is_reverse() {
tmp.complement()
} else {
tmp
}
let Some(base) = get_query_base(&record, q) else {
continue 'pileup;
};
if implicit_bases.contains(&base) {
add_to_tally(
Expand Down Expand Up @@ -860,15 +889,13 @@ impl PileupWorker for GenericPileupWorker {
break 'overran;
} else {
assert!(pos > q);
let base = {
let tmp =
DnaBase::try_from(record.seq()[q])
.unwrap();
if record.is_reverse() {
tmp.complement()
} else {
tmp
}
mod_pos = Some(pos);
canonical_base = Some(can_base);
pos_base_mod_call = Some(call);
mod_strand = Some(pos_mod_strand);
let Some(base) = get_query_base(&record, q)
else {
break 'overran;
};
if implicit_bases.contains(&base) {
add_to_tally(
Expand All @@ -891,10 +918,6 @@ impl PileupWorker for GenericPileupWorker {
motif_idxs,
);
}
mod_pos = Some(pos);
canonical_base = Some(can_base);
pos_base_mod_call = Some(call);
mod_strand = Some(pos_mod_strand);
break 'overran;
}
}
Expand All @@ -910,14 +933,8 @@ impl PileupWorker for GenericPileupWorker {
);
}
(Some(q), None) => {
let base = {
let tmp =
DnaBase::try_from(record.seq()[q]).unwrap();
if record.is_reverse() {
tmp.complement()
} else {
tmp
}
let Some(base) = get_query_base(&record, q) else {
continue 'pileup;
};
if implicit_bases.contains(&base) {
add_to_tally(
Expand Down Expand Up @@ -971,6 +988,7 @@ impl PileupWorker for GenericPileupWorker {
pileup_space.interval_width = width;
pileup_space.position_feature_counts = position_feature_counts;
pileup_space.failed_records = erred_records;
pileup_space.reference_equal_records = reference_equal_records;
Ok(pileup_space)
}
}
Expand Down
Loading