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
30 changes: 28 additions & 2 deletions crates/fgumi-consensus/src/codec_caller.rs
Original file line number Diff line number Diff line change
Expand Up @@ -1454,8 +1454,26 @@ impl CodecConsensusCaller {
&mut self,
records: Vec<RawRecord>,
) -> std::result::Result<ConsensusOutput, CodecConsensusError> {
let result = self.consensus_reads_raw(&records)?;
// When a group fails to produce consensus, all its records are rejected
let result = match self.consensus_reads_raw(&records) {
Ok(r) => r,
// A recoverable high-duplex-disagreement reject. fgbio routes the
// source records to `--rejects` (RejectionReason.HighDuplexDisagreement)
// and bumps `consensusReadsFilteredHighDisagreement`. Mirror both here
// before propagating the typed Err: previously the `?` unwound before
// the reject-capture block below, so disagreeing reads were silently
// dropped from BOTH --output and --rejects, and the HDD counter never
// incremented (S9b-006).
Err(e) if e.is_duplex_disagreement() => {
self.stats.consensus_reads_rejected_hdd += 1;
if self.track_rejects {
self.rejected_reads.extend(records.into_iter().map(RawRecord::into_inner));
}
return Err(e);
}
Err(e) => return Err(e),
};
// When a group fails to produce consensus (without erroring), all its
// records are rejected.
if self.track_rejects && result.count == 0 && !records.is_empty() {
self.rejected_reads.extend(records.into_iter().map(RawRecord::into_inner));
}
Expand Down Expand Up @@ -2718,6 +2736,14 @@ mod tests {
"expected DuplexDisagreementCount, got: {err:?}"
);
assert!(err.is_duplex_disagreement());
// The disagreement arm bumps the high-duplex-disagreement counter
// (fgbio's `consensusReadsFilteredHighDisagreement`) exactly once for
// the single disagreeing molecule (S9b-006).
assert_eq!(
caller.statistics().consensus_reads_rejected_hdd,
1,
"one HDD rejection should be counted for the disagreeing molecule"
);
}

/// Verifies that `consensus_reads_typed` returns the typed
Expand Down
16 changes: 13 additions & 3 deletions src/lib/pipeline/chains/commands/codec.rs
Original file line number Diff line number Diff line change
Expand Up @@ -282,9 +282,19 @@ pub(crate) fn build_codec_consensus_step_with_rejects(
-> io::Result<Process2Output<DecompressedBlock, DecompressedBlock>> {
let (consensus, rejects) =
run_codec_consensus_batch(state, item, track_rejects, &accumulators, &progress)?;
match rejects {
Some(r) => Ok(Process2Output::both(consensus, r)),
None => Ok(Process2Output::only_a(consensus)),
if let Some(r) = rejects {
Ok(Process2Output::both(consensus, r))
} else {
// X5-001: an all-clean batch must still emit a (zero-byte) rejects
// block so the rejects branch's `ByItemOrdinal` reorder stage sees a
// dense serial sequence; pushing nothing (`only_a`) leaves a gap at
// this serial that wedges the rejects sink. The empty `Vec` does not
// allocate and produces no physical BGZF block.
let batch_serial = consensus.batch_serial;
Ok(Process2Output::both(
consensus,
DecompressedBlock { batch_serial, bytes: Vec::new() },
))
}
};

Expand Down
16 changes: 13 additions & 3 deletions src/lib/pipeline/chains/commands/duplex.rs
Original file line number Diff line number Diff line change
Expand Up @@ -370,9 +370,19 @@ pub(crate) fn build_duplex_consensus_step_with_rejects(
&accumulators,
&progress,
)?;
match rejects {
Some(r) => Ok(Process2Output::both(consensus, r)),
None => Ok(Process2Output::only_a(consensus)),
if let Some(r) = rejects {
Ok(Process2Output::both(consensus, r))
} else {
// X5-001: an all-clean batch must still emit a (zero-byte) rejects
// block so the rejects branch's `ByItemOrdinal` reorder stage sees a
// dense serial sequence; pushing nothing (`only_a`) leaves a gap at
// this serial that wedges the rejects sink. The empty `Vec` does not
// allocate and produces no physical BGZF block.
let batch_serial = consensus.batch_serial;
Ok(Process2Output::both(
consensus,
DecompressedBlock { batch_serial, bytes: Vec::new() },
))
}
};

Expand Down
18 changes: 15 additions & 3 deletions src/lib/pipeline/chains/commands/simplex.rs
Original file line number Diff line number Diff line change
Expand Up @@ -346,9 +346,21 @@ pub(crate) fn build_simplex_consensus_step_with_rejects(
&accumulators,
&progress,
)?;
match rejects {
Some(r) => Ok(Process2Output::both(consensus, r)),
None => Ok(Process2Output::only_a(consensus)),
if let Some(r) = rejects {
Ok(Process2Output::both(consensus, r))
} else {
// X5-001: an all-clean batch must still emit a (zero-byte) rejects
// block so the rejects branch's `ByItemOrdinal` reorder stage sees a
// dense serial sequence. Pushing nothing here (`only_a`) leaves a
// permanent gap at this serial that wedges `try_pop_in_order`, so
// the rejects sink never drains. The empty `Vec` does not allocate
// and produces no physical BGZF block (compress/write of `&[]` are
// no-ops; the single EOF is appended once at drain).
let batch_serial = consensus.batch_serial;
Ok(Process2Output::both(
consensus,
DecompressedBlock { batch_serial, bytes: Vec::new() },
))
}
};

Expand Down
13 changes: 10 additions & 3 deletions src/lib/pipeline/steps/correct/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -106,9 +106,16 @@ pub(crate) fn correct_step_with_rejects(
batch: BamTemplateBatch|
-> io::Result<Process2Output<BamTemplateBatch, DecompressedBlock>> {
let (kept, rejects_opt) = run_batch_with_rejects(state, &cfg_run, batch)?;
match rejects_opt {
Some(rejects) => Ok(Process2Output::both(kept, rejects)),
None => Ok(Process2Output::only_a(kept)),
if let Some(rejects) = rejects_opt {
Ok(Process2Output::both(kept, rejects))
} else {
// X5-001: an all-clean batch must still emit a (zero-byte) rejects
// block so the rejects branch's `ByItemOrdinal` reorder stage sees a
// dense serial sequence; pushing nothing (`only_a`) leaves a gap at
// this serial that wedges the rejects sink. The empty `Vec` does not
// allocate and produces no physical BGZF block.
let batch_serial = kept.batch_serial;
Ok(Process2Output::both(kept, DecompressedBlock { batch_serial, bytes: Vec::new() }))
}
};

Expand Down
17 changes: 17 additions & 0 deletions tests/integration/helpers/assertions.rs
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,23 @@ pub fn assert_rejects_header_matches_input(
);
}

/// Reads every record from a BAM file into owned `RecordBuf`s.
///
/// Useful for asserting full record identity (sequence, flags, tags) — e.g.
/// that a `--rejects` stream passed its source reads through byte-identically.
///
/// # Panics
///
/// Panics if the file cannot be opened, the header cannot be read, or any
/// record cannot be decoded.
pub fn read_record_bufs(path: &std::path::Path) -> Vec<RecordBuf> {
let mut reader = noodles::bam::io::Reader::new(
std::fs::File::open(path).expect("Failed to open BAM for record-buf read"),
);
let header = reader.read_header().expect("Failed to read BAM header");
reader.record_bufs(&header).map(|r| r.expect("Failed to read BAM record")).collect()
}

/// Reads a BAM file and returns the multiset of record (read) names it holds.
///
/// Used to assert *family routing* — e.g. that a `--min-reads`-rejected family
Expand Down
51 changes: 38 additions & 13 deletions tests/integration/test_codec_command.rs
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ use fgumi_lib::commands::command::Command;
use fgumi_lib::sam::SamTag;
use fgumi_raw_bam::{RawRecord, SamBuilder, flags as raw_flags};
use noodles::bam;
use noodles::sam::alignment::record_buf::RecordBuf;
use std::fs;
use std::path::PathBuf;
use tempfile::TempDir;
Expand Down Expand Up @@ -535,6 +536,7 @@ fn test_codec_command_recovers_from_duplex_disagreement() {
let temp_dir = TempDir::new().expect("Failed to create temp dir");
let input_bam = temp_dir.path().join("input.bam");
let output_bam = temp_dir.path().join("output.bam");
let rejects_bam = temp_dir.path().join("rejects.bam");

// Single offset pair → 20 single-strand positions = 20 disagreements,
// which trips the count threshold below.
Expand All @@ -547,6 +549,8 @@ fn test_codec_command_recovers_from_duplex_disagreement() {
input_bam.to_str().unwrap(),
"--output",
output_bam.to_str().unwrap(),
"--rejects",
rejects_bam.to_str().unwrap(),
"--min-reads",
"1",
"--min-duplex-length",
Expand All @@ -567,6 +571,25 @@ fn test_codec_command_recovers_from_duplex_disagreement() {
let _header = reader.read_header().unwrap();
let consensus_count = reader.records().count();
assert_eq!(consensus_count, 0, "Disagreeing molecule should produce no consensus");

// The disagreement-failed pair (R1 + R2) is routed to --rejects, matching
// fgbio (`RejectionReason.HighDuplexDisagreement`) and the threaded sibling
// test (S9b-006). Asserts the single-threaded path populates rejects too.
assert!(rejects_bam.exists(), "Rejects BAM file should be created");
// The disagreement-failed pair (R1 + R2) must reach --rejects byte-identical
// to its input records (passed through unchanged). Comparing full
// `RecordBuf`s catches sequence/quality/CIGAR/tag mutation AND distinguishes
// R1 from R2 via the segment flags — stronger than a (name, flags) check.
let expected: Vec<RecordBuf> = crate::helpers::assertions::read_record_bufs(&input_bam)
.into_iter()
.filter(|r| r.name().is_some_and(|n| n == "disagree_read"))
.collect();
assert_eq!(expected.len(), 2, "fixture sanity: the disagreeing pair is two records");
let observed = crate::helpers::assertions::read_record_bufs(&rejects_bam);
assert_eq!(
observed, expected,
"the disagreeing pair (R1+R2) must reach --rejects byte-identical to its input records"
);
Comment thread
nh13 marked this conversation as resolved.
}

/// Verifies that the parallel `Codec::execute_threads_mode` path recovers
Expand Down Expand Up @@ -611,20 +634,22 @@ fn test_codec_command_recovers_from_duplex_disagreement_threaded() {
let consensus_count = reader.records().count();
assert_eq!(consensus_count, 0, "Disagreeing molecule should produce no consensus");

// Exercising --rejects forces the parallel-mode `flush_byte_records`
// call inside the typed-disagreement arm (`is_duplex_disagreement()`
// branch in `execute_threads_mode`). The current behavior is that the
// disagreement short-circuit in `consensus_reads_typed` returns the
// typed error before populating `rejected_reads`, so the file is
// created but stays empty — match that here without baking the
// open question (whether disagreement-failed reads *should* land in
// --rejects) into the test.
// High-duplex-disagreement molecules are routed to --rejects, matching
// fgbio's `RejectionReason.HighDuplexDisagreement` (S9b-006). The single
// offset pair contributes its two source reads (R1 + R2), so the rejects
// BAM holds exactly 2 records.
assert!(rejects_bam.exists(), "Rejects BAM file should be created");
let mut rejects_reader = bam::io::Reader::new(fs::File::open(&rejects_bam).unwrap());
let _ = rejects_reader.read_header().unwrap();
let reject_count = rejects_reader.records().count();
// Full `RecordBuf` identity against the input records (passed through
// unchanged) — catches sequence/quality/CIGAR/tag mutation and distinguishes
// R1 from R2 via the segment flags, like the single-threaded sibling.
let expected: Vec<RecordBuf> = crate::helpers::assertions::read_record_bufs(&input_bam)
.into_iter()
.filter(|r| r.name().is_some_and(|n| n == "disagree_read_mt"))
.collect();
assert_eq!(expected.len(), 2, "fixture sanity: the disagreeing pair is two records");
let observed = crate::helpers::assertions::read_record_bufs(&rejects_bam);
assert_eq!(
reject_count, 0,
"Disagreement short-circuits before reject tracking; rejects file is empty"
observed, expected,
"the disagreeing pair (R1+R2) must reach --rejects byte-identical to its input records"
);
}
106 changes: 106 additions & 0 deletions tests/integration/test_correct_command.rs
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ use fgumi_lib::commands::correct::CorrectUmis;
use fgumi_raw_bam::RawRecord;
use noodles::bam;
use noodles::sam::alignment::io::Write as AlignmentWrite;
use noodles::sam::alignment::record_buf::RecordBuf;
use std::fs;
use std::path::PathBuf;
use tempfile::TempDir;
Expand Down Expand Up @@ -247,6 +248,111 @@ fn test_correct_command_rejects_streaming_threaded_integrity() {
assert_eq!(observed_names, expected_names, "unexpected reject-name set");
}

/// Regression test for X5-001: a fully-clean batch must not wedge the rejects
/// branch. The rejects fan-out's branch B uses a `ByItemOrdinal` reorder stage,
/// so every batch serial must produce a branch-B item. A batch with zero
/// rejects previously pushed nothing (`Process2Output::only_a`), leaving a
/// permanent gap in the serial sequence that stalls `try_pop_in_order` forever
/// — any `--rejects` run with at least one all-clean batch hangs (the deadlock
/// monitor eventually fails it). The fix emits an empty (zero-byte) rejects
/// block for clean batches so serials stay dense.
///
/// Construction: the correct step's input is batched by a byte budget
/// (`per_step_byte_limit`, 4 MiB), not a fixed template count, so the clean
/// prefix must exceed that budget to guarantee the first emitted batch(es) are
/// entirely correctable (zero rejects → the gap-triggering batches). 8000
/// correctable records with 600 bp sequences clear 4 MiB several times over,
/// followed by uncorrectable records whose rejects land in a later batch. The
/// existing streaming-integrity test above fits all its records in a single
/// byte-bounded batch, so it never exercises the all-clean-batch gap.
#[test]
fn test_correct_command_rejects_all_clean_batch_does_not_wedge() {
let temp_dir = TempDir::new().unwrap();
let input_bam = temp_dir.path().join("input.bam");
let output_bam = temp_dir.path().join("output.bam");
let rejects_bam = temp_dir.path().join("rejects.bam");
let whitelist = temp_dir.path().join("whitelist.txt");

// The pipeline batches the correct step's input by a byte budget
// (`per_step_byte_limit`, 4 MiB), not a fixed template count, so the clean
// prefix must exceed that budget to guarantee the first emitted batch is
// entirely correctable (zero rejects). Use long sequences so a few thousand
// records clear 4 MiB. Each record is its own template (unique QNAME via
// create_umi_family's "{base}_{i}").
let long_seq: String = "ACGT".repeat(150); // 600 bp
let clean_prefix_size: usize = 8000;
let reject_size: u32 = 200;
let clean = create_umi_family("ACGTACGT", clean_prefix_size, "clean", &long_seq, 30);
// UMI far from the whitelist → uncorrectable → rejected, landing in a batch
// after the all-clean leading batch.
let rej = create_umi_family("TTTTTTTT", reject_size as usize, "rej", &long_seq, 30);

let expected_order: Vec<String> = (0..reject_size).map(|i| format!("rej_{i}")).collect();
// Keep the generated reject inputs (converted the same way `create_umi_bam`
// writes them) so we can assert the rejects BAM preserves each record's full
// identity — sequence, flags, tags — not just its QNAME. A names-only check
// would miss mutation of a rejected record.
let expected_records: Vec<RecordBuf> = rej.iter().map(to_record_buf).collect();

create_umi_bam(&input_bam, vec![clean, rej]);
create_whitelist(&whitelist, &["ACGTACGT"]);

let cmd = CorrectUmis::try_parse_from([
"correct",
"--input",
input_bam.to_str().unwrap(),
"--output",
output_bam.to_str().unwrap(),
"--umi-files",
whitelist.to_str().unwrap(),
"--max-mismatches",
"1",
"--min-distance",
"1",
"--rejects",
rejects_bam.to_str().unwrap(),
"--threads",
"4",
// Shorten the deadlock-monitor window so that, if the fix regresses,
// the wedge is reported as a failure in ~18s (3s × the 6× fatal
// multiplier) rather than the ~60s default. The fixed path completes in
// well under a second, so this never fires in the passing case.
"--deadlock-timeout",
"3",
"--compression-level",
"1",
])
.expect("failed to parse correct args");
// Pre-fix this call wedges on the rejects sink (the run never completes;
// the deadlock monitor fails it after its timeout). Post-fix it returns.
cmd.execute("fgumi correct")
.expect("correct --rejects must complete even with an all-clean leading batch");

assert!(rejects_bam.exists(), "Rejects BAM not created");
crate::helpers::assertions::assert_has_bgzf_eof(&rejects_bam);

let mut reader = bam::io::Reader::new(fs::File::open(&rejects_bam).unwrap());
let header = reader.read_header().unwrap();
let observed_records: Vec<RecordBuf> =
reader.record_bufs(&header).map(|r| r.expect("read reject record")).collect();

// Names in input order — a clear message on drop/dupe/reorder.
let observed_order: Vec<String> = observed_records
.iter()
.map(|r| r.name().expect("reject record missing read name").to_string())
.collect();
assert_eq!(
observed_order, expected_order,
"exactly the uncorrectable records, in input order, must reach --rejects"
);
Comment thread
nh13 marked this conversation as resolved.
// Full record identity — catches mutation of a rejected read's sequence,
// flags, or tags that a names-only check would let through.
assert_eq!(
observed_records, expected_records,
"rejected records must be byte-identical to the uncorrectable inputs"
);
}

/// Smoke check — `fgumi correct` runs to completion via the typed-step
/// framework entry point (`execute_new_pipeline`) and emits its
/// distinctive log line. The legacy iterator path was deleted in Phase 1
Expand Down
Loading
Loading