Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
110 changes: 87 additions & 23 deletions src/lib/commands/dedup.rs
Original file line number Diff line number Diff line change
Expand Up @@ -288,7 +288,7 @@ pub(crate) struct CollectedDedupCounts {
/// accumulate this in `serialize_fn`: that closure also runs once per group,
/// but workers execute it in parallel completion order, not coordinate order.
#[derive(Default)]
struct DuplicationLadderRecorder {
pub(crate) struct DuplicationLadderRecorder {
/// Snapshot interval in cumulative templates (`--ladder-interval`).
interval: u64,
/// Per-library running totals and next snapshot threshold.
Expand Down Expand Up @@ -323,7 +323,7 @@ struct LadderLibraryState {
}

impl DuplicationLadderRecorder {
fn new(interval: u64) -> Self {
pub(crate) fn new(interval: u64) -> Self {
Self { interval, per_library: AHashMap::new(), rows: Vec::new() }
}

Expand All @@ -341,7 +341,7 @@ impl DuplicationLadderRecorder {
///
/// MUST be called in serial/coordinate order — see the ordering note on
/// [`DuplicationLadderRecorder`].
fn record(&mut self, library_idx: u16, group_counts: &DedupCounts) {
pub(crate) fn record(&mut self, library_idx: u16, group_counts: &DedupCounts) {
if group_counts.total_templates == 0 {
return;
}
Expand Down Expand Up @@ -381,7 +381,7 @@ impl DuplicationLadderRecorder {

/// Emits a final snapshot per library at its true total, unless the last
/// interval snapshot already landed exactly on that total.
fn finish(&mut self) {
pub(crate) fn finish(&mut self) {
for (&library_idx, state) in &self.per_library {
if state.templates_seen > state.last_emitted_at {
self.rows.push((
Expand Down Expand Up @@ -1074,7 +1074,7 @@ pub(crate) fn process_position_group(
//////////////////////////////////////////////////////////////////////////////

/// UMI-aware duplicate marking command.
#[derive(Debug, Parser)]
#[derive(Debug, Clone, Parser)]
#[command(
name = "dedup",
about = "\x1b[38;5;151m[DEDUP]\x1b[0m \x1b[36mMark or remove PCR duplicates using UMI information\x1b[0m",
Expand Down Expand Up @@ -1327,9 +1327,13 @@ impl Command for MarkDuplicates {
bail!("--no-umi cannot be used with --strategy paired");
}

// Handle --no-umi mode: force identity strategy
// Handle --no-umi mode: force identity strategy. The effective-strategy
// computation itself must run on both paths (validate_index_threshold
// below consumes it), but the override log is gated to the non-chain
// path: the --threads chain re-emits it via add_dedup, so logging here
// too would double-log.
let (effective_strategy, no_umi_edits_override) = if self.no_umi {
if !matches!(self.strategy, Strategy::Identity) {
if !matches!(self.strategy, Strategy::Identity) && self.threading.threads.is_none() {
info!("--no-umi mode: overriding strategy to identity");
}
(Strategy::Identity, true)
Expand Down Expand Up @@ -1357,6 +1361,19 @@ impl Command for MarkDuplicates {
// Validate the input exists (stdin paths are exempt).
self.io.validate()?;

// --threads N: run the dedup stage on the declarative chain builder.
// Dispatch here — after the reader-free pre-flight validations above
// (output collisions, strategy/min-umi combos, index-threshold,
// input existence), which must run for both paths — but BEFORE the
// timer/banner/reader below. `execute_chain` → `add_dedup` re-emits the
// timer, banner, and threading log lines and opens its own source, so
// running them here too would double-log and pre-consume stdin
// (breaking stdin + --threads). The no-`--threads` path keeps the
// hand-rolled unified pipeline below.
if self.threading.threads.is_some() {
return self.execute_chain(command_line);
}

let min_mapq: u8 = self.min_map_q.unwrap_or(0);

let timer = OperationTimer::new("Marking duplicates");
Expand Down Expand Up @@ -1669,28 +1686,52 @@ impl Command for MarkDuplicates {
);
}

// --metrics is optional, so it must not be the only channel reporting
// dropped templates: a run that filters everything otherwise logs a bare
// "0 templates" and is indistinguishable from empty input. Reasons are
// named by their column suffix so the log points at the metrics column.
let filtered = final_counts.filter_counts.total_rejected_templates();
if filtered > 0 {
let reasons: Vec<String> = TemplateFilterReason::ALL
.into_iter()
.filter_map(|reason| {
let count = final_counts.filter_counts.rejected_templates(reason);
(count > 0).then(|| format!("{count} {}", reason.column_suffix()))
})
.collect();
info!("Filtered out {filtered} templates before marking: {}", reasons.join(", "));
}
log_filtered_templates(&final_counts.filter_counts);

timer.log_completion(final_counts.total_reads);

Ok(())
}
}

impl MarkDuplicates {
/// Run the dedup stage on the declarative chain builder (the `--threads N`
/// path).
///
/// Replaces the hand-rolled unified-pipeline construction in `execute` for
/// the threaded case. The chain opens its own source, validates the
/// template-coordinate sort order, injects `@PG`, assigns `MoleculeId`s
/// deterministically, writes the output BAM, and writes the metrics /
/// family-size histogram / duplication ladder via its finalize hook — all
/// through the same shared helpers as the non-chain path, so the two
/// orchestrations stay in parity. The no-`--threads` path keeps its own
/// unified pipeline in `execute`, which is the in-process parity oracle for
/// this one (see `test_dedup_chain_matches_single_threaded`).
fn execute_chain(&self, command_line: &str) -> Result<()> {
use crate::pipeline::chains::{
ChainSpec, SingleStageContext, Stage, StageOptionsBag, build_for,
};

// `add_dedup` re-emits the timer/banner/threading log lines but not the
// CRC-verify status; emit it here so the --threads path reports it once,
// matching the non-chain path. (dedup has no memory-debug knobs, so
// unlike group there is no warn block to add here.)
self.io.log_effective_check_crc();

let stage_opts = StageOptionsBag { dedup: Some(self.clone()), ..Default::default() };
let ctx = SingleStageContext {
io: &self.io,
threading: &self.threading,
compression: &self.compression,
scheduler: &self.scheduler_opts,
queue_memory: &self.queue_memory,
command_line,
};
let spec = ChainSpec::single_stage(Stage::Dedup, stage_opts, &ctx);
build_for(spec)?.run()
}
}

//////////////////////////////////////////////////////////////////////////////
// Metrics writing
//////////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -1766,7 +1807,30 @@ pub(crate) fn write_family_size_histogram(

/// Writes the `--duplication-ladder` TSV: one row per (library, snapshot),
/// sorted deterministically by library name then ascending `templates_seen`.
fn write_duplication_ladder(
/// Log the "Filtered out N templates before marking" diagnostic, if any were
/// rejected. Called by both the non-chain `execute` tail and the chain's
/// `DedupFinalizeHook` so `--threads` and no-`--threads` report filtered
/// templates identically.
///
/// `--metrics` is optional, so this must not be the only channel reporting
/// dropped templates: a run that filters everything otherwise logs a bare
/// "0 templates" and is indistinguishable from empty input. Reasons are named
/// by their column suffix so the log points at the metrics column.
pub(crate) fn log_filtered_templates(filter_counts: &TemplateFilterCounts) {
let filtered = filter_counts.total_rejected_templates();
if filtered > 0 {
let reasons: Vec<String> = TemplateFilterReason::ALL
.into_iter()
.filter_map(|reason| {
let count = filter_counts.rejected_templates(reason);
(count > 0).then(|| format!("{count} {}", reason.column_suffix()))
})
.collect();
info!("Filtered out {filtered} templates before marking: {}", reasons.join(", "));
}
}

pub(crate) fn write_duplication_ladder(
recorder: &DuplicationLadderRecorder,
library_index: &LibraryIndex,
path: &PathBuf,
Expand Down
17 changes: 16 additions & 1 deletion src/lib/pipeline/chains/builder.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4295,7 +4295,18 @@ impl<'a> ChainBuilder<'a> {
dedup.include_unmapped,
accumulators_for_process,
);
let mi_assign_step = build_mi_assign_step(self.tuning.per_step_byte_limit);
// Duplication-saturation ladder recorder (`--duplication-ladder`).
// `Some` only when the flag is set; shared (via `Arc`) between the serial
// MI-assign step, which accumulates it in coordinate order, and the
// finalize hook, which writes it. `None` means the MI-assign step does
// zero extra work.
let ladder_recorder = dedup.duplication_ladder.as_ref().map(|_| {
Arc::new(parking_lot::Mutex::new(
crate::commands::dedup::DuplicationLadderRecorder::new(dedup.ladder_interval),
))
});
let mi_assign_step =
build_mi_assign_step(self.tuning.per_step_byte_limit, ladder_recorder.clone());
let serialize_step = build_serialize_step(
self.tuning.per_step_byte_limit,
dedup.remove_duplicates,
Expand All @@ -4321,6 +4332,10 @@ impl<'a> ChainBuilder<'a> {
accumulators,
metrics_path: dedup.metrics.clone(),
family_size_histogram_path: dedup.family_size_histogram.clone(),
// Pair the path with its recorder: both are derived from
// `dedup.duplication_ladder`, so `zip` yields `Some` exactly when the
// flag is set — the invariant is now in the type, not an assert.
duplication_ladder: dedup.duplication_ladder.clone().zip(ladder_recorder),
library_index: fgumi_bam_io::LibraryIndex::from_header(&self.header),
sample: crate::commands::dedup::resolve_sample(&self.header, dedup.sample.as_deref()),
timer,
Expand Down
45 changes: 43 additions & 2 deletions src/lib/pipeline/chains/commands/dedup.rs
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,9 @@ use ahash::AHashMap;
use anyhow::{Result, bail};

use crate::commands::dedup::{
BatchedProcessedDedupGroups, CollectedDedupMetrics, DedupCounts, ProcessedDedupGroup,
process_position_group, write_dedup_metrics, write_family_size_histogram,
BatchedProcessedDedupGroups, CollectedDedupMetrics, DedupCounts, DuplicationLadderRecorder,
ProcessedDedupGroup, process_position_group, write_dedup_metrics, write_duplication_ladder,
write_family_size_histogram,
};
use crate::logging::OperationTimer;
use crate::per_thread_accumulator::PerThreadAccumulator;
Expand Down Expand Up @@ -81,6 +82,16 @@ pub(crate) struct DedupFinalizeHook {
pub(crate) accumulators: Arc<PerThreadAccumulator<CollectedDedupMetrics>>,
pub(crate) metrics_path: Option<std::path::PathBuf>,
pub(crate) family_size_histogram_path: Option<std::path::PathBuf>,
/// `--duplication-ladder` output path paired with the recorder accumulated
/// by the serial MI-assign step; `Some` exactly when the flag is set. A
/// single `Option<(path, recorder)>` (rather than two coupled `Option`s)
/// makes the "both or neither" invariant unrepresentable-when-broken, so
/// finalize needs no runtime assert. The recorder is held behind an `Arc`
/// shared with the MI-assign step, so the ladder is read here **through the
/// lock** rather than by `Arc::try_unwrap` — the step's clone may still be
/// alive when finalize hooks run.
pub(crate) duplication_ladder:
Option<(std::path::PathBuf, Arc<parking_lot::Mutex<DuplicationLadderRecorder>>)>,
/// Library index resolved from the input header, needed for the per-library
/// metrics rows.
pub(crate) library_index: LibraryIndex,
Expand All @@ -95,6 +106,7 @@ impl FinalizeHook for DedupFinalizeHook {
accumulators,
metrics_path,
family_size_histogram_path,
duplication_ladder,
library_index,
sample,
timer,
Expand Down Expand Up @@ -136,6 +148,15 @@ impl FinalizeHook for DedupFinalizeHook {
write_family_size_histogram(&final_family_sizes, path)?;
}

// Write the duplication-saturation ladder if requested. Accessed through
// the lock (not `Arc::try_unwrap`): the serial MI-assign step holds a
// clone of this `Arc` that may still be alive here.
if let Some((path, recorder)) = duplication_ladder {
let mut guard = recorder.lock();
guard.finish();
write_duplication_ladder(&guard, &library_index, &path)?;
}

// Log summary banner (verbatim from the original execute).
log::info!(
"Deduplication complete: {} templates ({} unique, {} duplicates, {:.2}% duplicate rate)",
Expand All @@ -160,6 +181,10 @@ impl FinalizeHook for DedupFinalizeHook {
);
}

// Report filtered templates identically to the non-chain path (shared
// helper, so the two paths cannot drift).
crate::commands::dedup::log_filtered_templates(&final_counts.filter_counts);

timer.log_completion(final_counts.total_reads);

Ok(())
Expand Down Expand Up @@ -242,10 +267,20 @@ pub(crate) fn build_process_step(
/// Build the `MiAssignDedup` step: serial, `ByItemOrdinal`, assigns
/// monotonically increasing MI offsets to each batch.
///
/// When `ladder_recorder` is `Some`, this step also records the
/// `--duplication-ladder` saturation curve — per position group, in the same
/// serial/coordinate-order seam the non-chain path records it (`mi_assign_fn`),
/// **not** in the parallel serialize step. `MiAssign` is `Serial` +
/// `ByItemOrdinal`, so batches reach this closure in input-record order and the
/// groups within a batch are in coordinate order, exactly reproducing the
/// non-chain per-group stream — see the ordering note on
/// [`crate::commands::dedup::DuplicationLadderRecorder`].
///
/// `pub(crate)` — consumed only by `ChainBuilder::add_dedup`.
#[allow(clippy::type_complexity)]
pub(crate) fn build_mi_assign_step(
limit_bytes: u64,
ladder_recorder: Option<Arc<parking_lot::Mutex<DuplicationLadderRecorder>>>,
) -> MiAssign<
BatchedProcessedDedupGroups,
BatchedProcessedDedupGroups,
Expand All @@ -259,6 +294,12 @@ pub(crate) fn build_mi_assign_step(
limit_bytes,
move |next_mi: &mut u64, mut item: BatchedProcessedDedupGroups| {
assign_mi_offsets(next_mi, &mut item)?;
if let Some(recorder) = &ladder_recorder {
let mut guard = recorder.lock();
for processed in &item.groups {
guard.record(processed.library_idx, &processed.dedup_counts);
}
}
Ok(item)
},
)
Expand Down
Loading
Loading