Skip to content

fix(umi): exclude unmapped reads from consensus calling - #1168

Draft
nh13 wants to merge 2 commits into
mainfrom
nh_test-consensus-cap-half-mapped-pair
Draft

nh13 wants to merge 2 commits into
mainfrom
nh_test-consensus-cap-half-mapped-pair

Conversation

@nh13

@nh13 nh13 commented Aug 10, 2026

Copy link
Copy Markdown
Member

Problem

Unmapped reads reach the consensus callers and contribute bases and depth to consensus reads, with no alignment support.

GroupReadsByUmi keeps a read whose mate is mapped, and ConsensusCallingIterator applies the same predicate, so the unmapped end of a half-mapped pair survives into consensus calling. An unmapped read has an empty cigar, and Cigar.isPrefixOf returns true for an empty left hand side — the loop never executes:

def isPrefixOf(that: Cigar): Boolean = {
  if (that.elems.size < this.elems.size) false
  else {
    val lastIndex = this.elems.size - 1   // -1 for an empty cigar
    ...                                   // while (isPrefix && i <= lastIndex) never runs

So in filterToMostCommonAlignment an empty cigar matches every alignment group instead of forming one of its own, and the read is retained rather than rejected as a minority alignment. Whether it is retained at all depends on where it lands in that method's length-descending sort: if it sorts first it seeds a group nothing can join and is rejected, and if it sorts later it joins the winning group. An unmapped read gets no mate-overlap trimming, so it is frequently the longest read in its group — the outcome flips on a read-length tiebreak.

Observed on main, for a tag family of three read pairs where one pair has an unmapped R2:

Tool Symptom
CallMolecularConsensusReads --max-reads 3 reports cD=3 for R2 when only two reads are eligible. --max-reads 2 retains the unmapped read in preference to a mapped one, and every base of the R2 consensus no-calls.
CallDuplexConsensusReads Reports aD=3 for the AB single-strand consensus of R2 when only two AB-R2s are mapped.
CallCodecConsensusReads Unaffected — isPrimaryFrPair rejects a template with an unmapped end as NotPrimaryFrPair before the R1/R2 split.

Fix

Require r.mapped in ConsensusCallingIterator instead of r.mapped || (r.paired && r.mateMapped).

All three tools construct this iterator (CallMolecularConsensusReads.scala:188, CallDuplexConsensusReads.scala:164, CallCodecConsensusReads.scala:151), so the change covers every one of them structurally rather than relying on each caller to filter.

This restores the intent of #948. That PR widened the predicate to stop discarding the mapped end of a half-mapped pair, after #940 had dropped both ends. That end is still kept — it satisfies r.mapped on its own. Only the unmapped end, admitted as a side effect of how the predicate was widened, is now dropped.

Why not subGroupRecords

Returning a fourth "unmapped" group was considered and rejected: CodecConsensusCaller extends DuplexConsensusCaller but never calls subGroupRecords, doing its own split at CodecConsensusCaller.scala:196-197, so a fourth element would miss the CODEC path entirely. It would also widen a protected 3-tuple across three call sites, two of which already discard the fragments slot as (_, abR1s, abR2s).

Tests

Written first, watched failing, then fixed — the two commits are test: then fix: in that order, so the first is red by design. Happy to squash if you would rather not carry a red commit.

Each of the three tools gets a regression test over the same family shape, so none can regress independently:

  • CallMolecularConsensusReadsTest — two cases, --max-reads 2 and 3, asserting R1 is capped at max-reads, that R2's depth is 2, and that R2's bases are those of the mapped reads. The bases assertion is what catches the --max-reads 2 case, where cD alone looks correct.
  • CallDuplexConsensusReadsTest — asserts aD/bD on both ends, with R1's aD=3 as a positive control that mapped reads are unaffected.
  • CallCodecConsensusReadsTest — passes before and after; a guard against losing the NotPrimaryFrPair protection.

Before: 3 failures (2 molecular, 1 duplex), codec green. After: full suite 1608 tests, 0 failures.

Deliberately not addressed

The root cause in filterToMostCommonAlignment is untouched. An empty cigar is still a prefix of every cigar. This PR makes that unreachable from the three consensus tools; any future path producing an empty cigar would hit it again. Fixing the grouping loop itself has a wider blast radius and wants its own change.

Accounting. Dropped reads are filtered before the caller sees them, so they are not counted in raw_reads_considered and do not reach --rejects — the same treatment secondary and supplementary reads already get here. Routing them through rejectRecords with a dedicated RejectionReason means touching the usedByVanilla/usedByDuplex/usedByCodec tables, which is the same shape of work as the raw_reads_used gap tracked in #1167. The codec and molecular tests assert only that every rejected record belongs to the half-mapped template, so they hold under either choice.

Duplex has a second, separate exposure not covered here: areAllSameStrand (DuplexConsensusCaller.scala:260) reads negativeStrand off an unmapped record, which can falsely reject an entire AB+BA molecule as PotentialCollision. This PR removes unmapped reads before that check, so it is no longer reachable from the tool, but the method remains unguarded.

nh13 added 2 commits August 10, 2026 09:36
Add a regression test to each of the three consensus calling tools covering a
tag family that contains a half-mapped read pair. Two of the three currently
fail.

GroupReadsByUmi keeps a read whose mate is mapped, and ConsensusCallingIterator
does the same, so the unmapped end of a half-mapped pair reaches the consensus
callers. An unmapped read has an empty cigar, and Cigar.isPrefixOf returns true
for an empty left hand side, so in filterToMostCommonAlignment it matches every
alignment group rather than forming one of its own. It is therefore retained and
contributes bases and depth to the consensus with no alignment support.

CallMolecularConsensusReads: a family of three pairs, one half-mapped, at
--max-reads 3 reports cD=3 for R2 when only two reads are eligible. At
--max-reads 2 the cap retains the unmapped read in preference to a mapped one
and every base of the R2 consensus no-calls.

CallDuplexConsensusReads: the same family shape reports aD=3 for the AB single
strand consensus of R2 when only two AB-R2s are mapped.

CallCodecConsensusReads: already unaffected, since isPrimaryFrPair rejects a
template with an unmapped end as NotPrimaryFrPair before the R1/R2 split. The
test is a guard against that protection being lost, and passes as written.
Require reads to be mapped in ConsensusCallingIterator rather than accepting any
read that is either mapped or has a mapped mate. All three consensus calling
tools construct this iterator, so the change reaches every one of them by
construction rather than by each caller remembering to filter.

This restores the intent of #948, which widened the predicate to stop discarding
the mapped end of a half-mapped pair after #940 had dropped both ends. That end
is still kept: it satisfies `r.mapped` on its own. Only the unmapped end, which
was admitted as a side effect of how the predicate was written, is now dropped.

A template with both ends unmapped was already discarded here, and GroupReadsByUmi
filters such templates too, so no supported input loses coverage. There is no
facility for consensus calling from wholly unaligned data.

Note this does not fix the underlying behaviour in filterToMostCommonAlignment,
where an empty cigar is a prefix of every cigar and so matches every alignment
group. It makes that unreachable from these three tools; any future path that
produces an empty cigar would hit it again.

The dropped reads are filtered before the caller sees them, so they are not
counted in raw_reads_considered and do not reach --rejects, matching how
secondary and supplementary reads are already handled here. Routing them through
rejectRecords with a dedicated RejectionReason would require touching the
usedByVanilla/usedByDuplex/usedByCodec tables and is left for the metrics work
tracked in #1167.

Full suite: 1608 tests, 0 failures.
@codecov

codecov Bot commented Aug 10, 2026 •

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.95%. Comparing base (798aca4) to head (3bc4526).
⚠️ Report is 5 commits behind head on main.

Additional details and impacted files
@@           Coverage Diff           @@
##             main    #1168   +/-   ##
=======================================
  Coverage   95.95%   95.95%           
=======================================
  Files         132      132           
  Lines        8347     8347           
  Branches      984      952   -32     
=======================================
  Hits         8009     8009           
  Misses        338      338           
Flag Coverage Δ
unittests 95.95% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@github-actions

Copy link
Copy Markdown
PR Preview Action v1.6.1

🚀 View preview at
https://fulcrumgenomics.github.io/fgbio/pr-preview/pr-1168/

Built to branch gh-pages at 2026-08-10 16:49 UTC.
Preview will be ready when the GitHub Pages deployment is complete.

nh13 added a commit to fulcrumgenomics/fgumi that referenced this pull request Aug 10, 2026
…mapped

Drop unmapped records in the consensus pre-group filter rather than keeping any
record that is mapped *or* has a mapped mate. This is fgumi's copy of fgbio's
`ConsensusCallingIterator` filter, shared by the simplex, duplex and codec
commands, so all three are covered by construction.

An unmapped record has no cigar, and `is_cigar_prefix` returns true for an empty
left hand side, so in `select_most_common_alignment_group` an empty cigar matched
every alignment group rather than forming one of its own. The read was retained and
contributed bases and depth to the consensus with no alignment supporting it.
Whether it was retained at all depended on where it landed in the length-descending
sort, so the outcome turned on a read-length tiebreak. The mapped end of a
half-mapped pair is unaffected: it is mapped in its own right.

This is a behaviour change for input containing half-mapped pairs. It restores
parity with fgbio, where the equivalent family produces no consensus at all rather
than an R2 built from unaligned reads (fulcrumgenomics/fgbio#1168).

`--allow-unmapped` continues to relax the rule so wholly unmapped input can be
consensus called, as `fgumi group --allow-unmapped` intends. Under that flag,
unmapped reads are admitted to grouping but still never displace mapped ones:
`drop_unmapped_if_any_mapped` prefers the mapped reads within each set being
consensus called, and leaves a wholly unmapped set intact. Reads it drops are
attributed to `RejectionReason::Unmapped`, which had no production call site until
now, and are added to the rejected-index set so they reach the rejects output.

The `unmapped_paired_mate_mapped` case in `test_consensus_pregroup_keep_flags` is
inverted to pin the new behaviour, and a case is added for the mapped end of a
half-mapped pair, which is still kept.

Full suite: 7247 tests, 0 failures. Formatting and clippy clean.
nh13 added a commit to fulcrumgenomics/fgumi that referenced this pull request Aug 12, 2026
…mapped

Drop unmapped records in the consensus pre-group filter rather than keeping any
record that is mapped *or* has a mapped mate. This is fgumi's copy of fgbio's
`ConsensusCallingIterator` filter, shared by the simplex, duplex and codec
commands, so all three are covered by construction.

An unmapped record has no cigar, and `is_cigar_prefix` returns true for an empty
left hand side, so in `select_most_common_alignment_group` an empty cigar matched
every alignment group rather than forming one of its own. The read was retained and
contributed bases and depth to the consensus with no alignment supporting it.
Whether it was retained at all depended on where it landed in the length-descending
sort, so the outcome turned on a read-length tiebreak. The mapped end of a
half-mapped pair is unaffected: it is mapped in its own right.

This is a behaviour change for input containing half-mapped pairs. It restores
parity with fgbio, where the equivalent family produces no consensus at all rather
than an R2 built from unaligned reads (fulcrumgenomics/fgbio#1168).

`--allow-unmapped` continues to relax the rule so wholly unmapped input can be
consensus called, as `fgumi group --allow-unmapped` intends. Under that flag,
unmapped reads are admitted to grouping but still never displace mapped ones:
`drop_unmapped_if_any_mapped` prefers the mapped reads within each set being
consensus called, and leaves a wholly unmapped set intact. Reads it drops are
attributed to `RejectionReason::Unmapped`, which had no production call site until
now, and are added to the rejected-index set so they reach the rejects output.

The `unmapped_paired_mate_mapped` case in `test_consensus_pregroup_keep_flags` is
inverted to pin the new behaviour, and a case is added for the mapped end of a
half-mapped pair, which is still kept.

Full suite: 7247 tests, 0 failures. Formatting and clippy clean.
nh13 added a commit to fulcrumgenomics/fgumi that referenced this pull request Aug 13, 2026
…mapped

Drop unmapped records in the consensus pre-group filter rather than keeping any
record that is mapped *or* has a mapped mate. This is fgumi's copy of fgbio's
`ConsensusCallingIterator` filter, shared by the simplex, duplex and codec
commands, so all three are covered by construction.

An unmapped record has no cigar, and `is_cigar_prefix` returns true for an empty
left hand side, so in `select_most_common_alignment_group` an empty cigar matched
every alignment group rather than forming one of its own. The read was retained and
contributed bases and depth to the consensus with no alignment supporting it.
Whether it was retained at all depended on where it landed in the length-descending
sort, so the outcome turned on a read-length tiebreak. The mapped end of a
half-mapped pair is unaffected: it is mapped in its own right.

This is a behaviour change for input containing half-mapped pairs. It restores
parity with fgbio, where the equivalent family produces no consensus at all rather
than an R2 built from unaligned reads (fulcrumgenomics/fgbio#1168).

`--allow-unmapped` continues to relax the rule so wholly unmapped input can be
consensus called, as `fgumi group --allow-unmapped` intends. Under that flag,
unmapped reads are admitted to grouping but still never displace mapped ones:
`drop_unmapped_if_any_mapped` prefers the mapped reads within each set being
consensus called, and leaves a wholly unmapped set intact. Reads it drops are
attributed to `RejectionReason::Unmapped`, which had no production call site until
now, and are added to the rejected-index set so they reach the rejects output.

The `unmapped_paired_mate_mapped` case in `test_consensus_pregroup_keep_flags` is
inverted to pin the new behaviour, and a case is added for the mapped end of a
half-mapped pair, which is still kept.

Full suite: 7247 tests, 0 failures. Formatting and clippy clean.
nh13 added a commit to fulcrumgenomics/fgumi that referenced this pull request Aug 14, 2026
…mapped (#737)

* test(consensus): pin that mapped reads are preferred over unmapped

Add regression tests to the three consensus commands covering a molecule that
contains a half-mapped read pair. The simplex and duplex cases currently fail.

`fgumi group --allow-unmapped` deliberately supports wholly unmapped templates,
but its gate is on templates with *both* reads unmapped, so a half-mapped pair is
kept regardless. The consensus pre-group filter then keeps any record that is
mapped *or* has a mapped mate, which admits the unmapped end as well. That end has
no cigar, and `is_cigar_prefix` returns true for an empty left hand side, so in
`select_most_common_alignment_group` it matches every alignment group rather than
forming one of its own. It is retained and contributes bases and depth to the
consensus with no alignment supporting it.

simplex: a family of three pairs, one half-mapped, reports cD=3 for R2 where only
two reads are eligible, and the unmapped read's bases reach the consensus.

duplex: the same molecule shape reports aD=3 for the AB single-strand consensus of
R2 where only two AB-R2s are mapped.

codec: already unaffected, since its alignment grouping is reached only for
templates forming a primary FR pair. The test guards that rather than reproducing
a defect, and passes as written.

Two further simplex cases cover `--allow-unmapped`, which must keep working: within
a set being consensus called mapped reads must still win, while a wholly unmapped
family must still produce a consensus.

* fix(consensus)!: require mapped reads for consensus unless --allow-unmapped

Drop unmapped records in the consensus pre-group filter rather than keeping any
record that is mapped *or* has a mapped mate. This is fgumi's copy of fgbio's
`ConsensusCallingIterator` filter, shared by the simplex, duplex and codec
commands, so all three are covered by construction.

An unmapped record has no cigar, and `is_cigar_prefix` returns true for an empty
left hand side, so in `select_most_common_alignment_group` an empty cigar matched
every alignment group rather than forming one of its own. The read was retained and
contributed bases and depth to the consensus with no alignment supporting it.
Whether it was retained at all depended on where it landed in the length-descending
sort, so the outcome turned on a read-length tiebreak. The mapped end of a
half-mapped pair is unaffected: it is mapped in its own right.

This is a behaviour change for input containing half-mapped pairs. It restores
parity with fgbio, where the equivalent family produces no consensus at all rather
than an R2 built from unaligned reads (fulcrumgenomics/fgbio#1168).

`--allow-unmapped` continues to relax the rule so wholly unmapped input can be
consensus called, as `fgumi group --allow-unmapped` intends. Under that flag,
unmapped reads are admitted to grouping but still never displace mapped ones:
`drop_unmapped_if_any_mapped` prefers the mapped reads within each set being
consensus called, and leaves a wholly unmapped set intact. Reads it drops are
attributed to `RejectionReason::Unmapped`, which had no production call site until
now, and are added to the rejected-index set so they reach the rejects output.

The `unmapped_paired_mate_mapped` case in `test_consensus_pregroup_keep_flags` is
inverted to pin the new behaviour, and a case is added for the mapped end of a
half-mapped pair, which is still kept.

Full suite: 7247 tests, 0 failures. Formatting and clippy clean.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant