Skip to content

Commit cf9ecfd

Browse files
committed
fix(pileup, generic): actually use combine-strands with generic workers.
1 parent a0ad253 commit cf9ecfd

File tree

2 files changed

+41
-8
lines changed

2 files changed

+41
-8
lines changed

modkit-core/src/pileup/pileup_processor.rs

Lines changed: 40 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -456,6 +456,7 @@ pub struct GenericPileupWorker {
456456
motif_bases: Vec<MotifInfo>,
457457
caller: MultipleThresholdModCaller,
458458
pileup_numeric_options: PileupNumericOptions,
459+
combine_strands: bool,
459460
}
460461

461462
impl GenericPileupWorker {
@@ -465,6 +466,7 @@ impl GenericPileupWorker {
465466
motif_bases: Vec<MotifInfo>,
466467
caller: MultipleThresholdModCaller,
467468
pileup_numeric_options: PileupNumericOptions,
469+
combine_strands: bool,
468470
) -> anyhow::Result<Self> {
469471
let mut reader = bam::IndexedReader::from_path(in_bam_fp)?;
470472
if reader_is_cram(&reader) {
@@ -485,7 +487,13 @@ impl GenericPileupWorker {
485487
}
486488
};
487489

488-
Ok(Self { reader, motif_bases, caller, pileup_numeric_options })
490+
Ok(Self {
491+
reader,
492+
motif_bases,
493+
caller,
494+
pileup_numeric_options,
495+
combine_strands,
496+
})
489497
}
490498
}
491499
impl PileupWorker for GenericPileupWorker {
@@ -525,9 +533,9 @@ impl PileupWorker for GenericPileupWorker {
525533
FxHashMap::default();
526534
let mut add_to_tally = |position_strand: PositionStrand,
527535
call: Call,
528-
ref_base: Option<DnaBase>,
536+
motif_info: Option<&MotifInfo>,
529537
motif_idxs: u8| {
530-
let call = match ref_base {
538+
let call = match motif_info.map(|x| x.primary_base) {
531539
Some(x) if call.matches_dna_base(&x) => call,
532540
Some(_x) => match call {
533541
Call::Delete => call,
@@ -540,6 +548,25 @@ impl PileupWorker for GenericPileupWorker {
540548
},
541549
None => call,
542550
};
551+
let position_strand = match motif_info {
552+
Some(motif_info) if self.combine_strands => {
553+
let (rpos, strand) = position_strand;
554+
if strand == Strand::Negative {
555+
(
556+
rpos.saturating_sub(
557+
(motif_info.reverse_offset
558+
- motif_info.forward_offset)
559+
as u32,
560+
),
561+
Strand::Positive,
562+
)
563+
} else {
564+
position_strand
565+
}
566+
}
567+
_ => position_strand,
568+
};
569+
543570
chrom_features
544571
.entry(position_strand)
545572
.or_insert_with(Tally2::default)
@@ -625,13 +652,13 @@ impl PileupWorker for GenericPileupWorker {
625652
(
626653
qpos,
627654
rpos,
628-
Some(self.motif_bases[idx].primary_base),
655+
Some(&self.motif_bases[idx]),
629656
indices_to_byte(bs),
630657
)
631658
})
632659
}
633660
FocusPositions2::AllPositions => {
634-
Some((qpos, rpos, Option::<DnaBase>::None, 0u8))
661+
Some((qpos, rpos, Option::<&MotifInfo>::None, 0u8))
635662
}
636663
}
637664
});
@@ -875,7 +902,12 @@ impl PileupWorker for GenericPileupWorker {
875902
o @ _ => o,
876903
})
877904
.flat_map(|((ref_pos, strand), tally)| {
878-
tally.into_counts(start_pos, ref_pos, combine_mods, strand)
905+
tally.into_counts(
906+
start_pos,
907+
ref_pos,
908+
combine_mods,
909+
if self.combine_strands { '.' } else { strand.to_char() },
910+
)
879911
})
880912
.filter(|x| x.is_valid())
881913
.collect::<Vec<PileupFeatureCounts2>>();
@@ -2114,7 +2146,7 @@ impl Tally2 {
21142146
start_pos: u32,
21152147
relative_position: u32,
21162148
combine_mods: bool,
2117-
strand: Strand,
2149+
strand: char,
21182150
) -> Vec<PileupFeatureCounts2> {
21192151
self.modcall_counts
21202152
.into_iter()
@@ -2175,7 +2207,7 @@ impl Tally2 {
21752207
ModifiedBase::Modified(mod_code) => {
21762208
Some(PileupFeatureCounts2::new(
21772209
relative_position.saturating_add(start_pos),
2178-
strand.to_char(),
2210+
strand,
21792211
filtered_coverage,
21802212
mod_code,
21812213
n_canonical,

modkit-core/src/pileup/subcommand.rs

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1229,6 +1229,7 @@ impl ModBamPileup {
12291229
motif_infos.clone(),
12301230
threshold_caller.clone(),
12311231
pileup_options.clone(),
1232+
self.combine_strands,
12321233
)
12331234
})
12341235
.collect::<anyhow::Result<Vec<_>>>()?;

0 commit comments

Comments
 (0)