Skip to content

Commit de8e5d4

Browse files
committed
fix(codec): avoid MatchError on dovetail FR pairs with asymmetric isFrPair
CodecConsensusCaller filtered primaries with a per-record `isFrPair` check then pattern-matched each read-name bucket as `Seq(rec, mate)`. For dovetail FR pairs whose aligned ends coincide (e.g. soft-clipping producing TLEN=+/-1), htsjdk's `SamPairUtil.getPairOrientation` returns FR for R1 but RF for R2, which kept one mate and dropped the other and produced a singleton bucket that blew up with `scala.MatchError`. Group primaries by read name first and keep only buckets that form a single primary FR pair, where the FR check is computed at the pair level from the unclipped 5' positions of both records and is therefore symmetric in its arguments. Templates that do not form a clean primary FR pair (singletons, tandems, RF, cross-contig) are now rejected cleanly with a new `NotPrimaryFrPair` reason instead of crashing. Preserve BAM-iteration order of templates when grouping by read name. The previous `groupBy(_.name).toSeq` produced a template order that depended on the JVM's `HashMap` layout, which propagated through `filterToMostCommonAlignment` and caused `maxBy(_.cigar.lengthOnTarget)` to break ties in hash order when selecting `longestR1Alignment` / `longestR2Alignment`. That could pair records from different templates, producing spurious `IndelErrorBetweenStrands` rejections when the resulting overlap geometry did not agree across the strands. A new `orderedPrimaryPairs` helper groups by read name but orders templates by first-occurrence in the input. Includes unit tests for the new helpers and a regression test that reproduces the exact geometry (68S53M8S / 28S48M53S, aligned ends coinciding) observed in real CODEC data.
1 parent baf6ad3 commit de8e5d4

3 files changed

Lines changed: 164 additions & 4 deletions

File tree

src/main/scala/com/fulcrumgenomics/umi/CodecConsensusCaller.scala

Lines changed: 51 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -170,10 +170,24 @@ class CodecConsensusCaller(readNamePrefix: String,
170170
Nil
171171
}
172172
else {
173-
// Extract just the primary alignments and then clip where they extend past the ends of their mates
174-
// TODO: Remove the isFrPair here and handle chimeric reads or reads with one unmapped etc.
175-
val primaries = pairs.filterNot(r => r.secondary || r.supplementary).filter(_.isFrPair)
176-
primaries.groupBy(_.name).foreach { case (_, Seq(rec, mate)) => clipper.clipExtendingPastMateEnds(rec, mate) }
173+
// Group the primary alignments by read name and retain only templates that form a single
174+
// primary FR pair. The FR check is performed at the pair level (with both records in hand)
175+
// because htsjdk's per-record `SamPairUtil.getPairOrientation` can disagree between mates
176+
// on dovetail pairs whose aligned ends coincide (TLEN=+/-1), which previously caused a
177+
// `scala.MatchError` on a singleton read-name bucket here.
178+
// Template order is preserved from BAM-iteration order (not hash order) so that downstream
179+
// `maxBy` tie-breaks on `cigar.lengthOnTarget` are stable across JVMs.
180+
// TODO: handle chimeric reads or reads with one unmapped etc.
181+
val primaryPairs = CodecConsensusCaller.orderedPrimaryPairs(
182+
pairs.filterNot(r => r.secondary || r.supplementary)
183+
)
184+
val (frPairs, nonFrPairs) = primaryPairs.partition {
185+
case (_, Seq(a, b)) => CodecConsensusCaller.isPrimaryFrPair(a, b)
186+
case _ => false
187+
}
188+
rejectRecords(nonFrPairs.flatMap(_._2), RejectionReason.NotPrimaryFrPair)
189+
frPairs.foreach { case (_, Seq(rec, mate)) => clipper.clipExtendingPastMateEnds(rec, mate) }
190+
val primaries = frPairs.flatMap(_._2)
177191

178192
val r1s = filterToMostCommonAlignment(primaries.filter(_.firstOfPair).map(toSourceReadForCodec))
179193
val r2s = filterToMostCommonAlignment(primaries.filter(_.secondOfPair).map(toSourceReadForCodec))
@@ -375,3 +389,36 @@ class CodecConsensusCaller(readNamePrefix: String,
375389
)
376390
}
377391
}
392+
393+
object CodecConsensusCaller {
394+
/** Groups primary records by read name while preserving BAM-iteration order of templates.
395+
*
396+
* Scala's `Seq.groupBy` returns a hash-keyed `Map`, so calling `.toSeq` on it produces a
397+
* template order that depends on the JVM's `String#hashCode` and the `HashMap` layout.
398+
* Downstream code selects the longest R1 and R2 via `maxBy(_.cigar.lengthOnTarget)`, which
399+
* returns the first element in iteration order on ties; relying on hash order there yields
400+
* non-obvious, JVM-dependent winners. This helper instead orders templates by the first
401+
* occurrence of each read name in the input (i.e. BAM-iteration order), so tie-breaks are
402+
* stable and match the order the records were observed.
403+
*/
404+
private[umi] def orderedPrimaryPairs(records: Seq[SamRecord]): Seq[(String, Seq[SamRecord])] = {
405+
val byName = records.groupBy(_.name)
406+
records.iterator.map(_.name).distinct.map(name => (name, byName(name))).toSeq
407+
}
408+
409+
/** Returns true if the two primary records form an FR pair. The check is symmetric in its
410+
* arguments so that both records of a pair always produce the same classification; htsjdk's
411+
* per-record `SamPairUtil.getPairOrientation` can disagree between mates on dovetail pairs
412+
* whose aligned ends coincide, which is why this is computed from unclipped 5' positions
413+
* rather than delegated to `SamRecord.isFrPair`.
414+
*/
415+
private[umi] def isPrimaryFrPair(a: SamRecord, b: SamRecord): Boolean = {
416+
if (!a.mapped || !b.mapped) false
417+
else if (a.refIndex != b.refIndex) false
418+
else if (a.positiveStrand == b.positiveStrand) false
419+
else {
420+
val (fwd, rev) = if (a.positiveStrand) (a, b) else (b, a)
421+
fwd.unclippedStart <= rev.unclippedEnd
422+
}
423+
}
424+
}

src/main/scala/com/fulcrumgenomics/umi/UmiConsensusCaller.scala

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,7 @@ object UmiConsensusCaller {
8585
case object IndelErrorBetweenStrands extends _RejectionReason("indel_error_between_strands", "Indel error between top/bottom strands", false, false, true)
8686
case object HighDuplexDisagreement extends _RejectionReason("high_duplex_disagreement", "Too many errors between top/bottoms strands", false, false, true)
8787
case object ClipOverlapFailed extends _RejectionReason("clip_overlap_failed", "See https://github.com/fulcrumgenomics/fgbio/issues/1090", false, false, true)
88+
case object NotPrimaryFrPair extends _RejectionReason("not_primary_fr_pair", "Template did not have a single primary FR pair of reads", false, false, true)
8889
}
8990

9091
/** Metric class for outputting consensus calling statistics. */

src/test/scala/com/fulcrumgenomics/umi/CodecConsensusCallerTest.scala

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -207,6 +207,38 @@ class CodecConsensusCallerTest extends UnitSpec with OptionValues {
207207
caller.consensusReadsFromSamRecords(rfPair) shouldBe Seq()
208208
}
209209

210+
it should "not throw on an FR dovetail pair whose aligned ends coincide (htsjdk asymmetric isFrPair)" in {
211+
// Reproduces a real template from a HEK293T CODEC dataset in which the aligned end of the
212+
// reverse-strand read coincides with the aligned start of the forward-strand read and
213+
// soft-clipping produces TLEN=+/-1. For this geometry htsjdk's
214+
// `SamPairUtil.getPairOrientation` returns FR for R1 and RF for R2, which previously caused
215+
// `CodecConsensusCaller` to drop R2 via a per-record `isFrPair` filter and then throw
216+
// `scala.MatchError` on the resulting singleton read-name bucket.
217+
val builder = new SamBuilder(readLength=129, baseQuality=35)
218+
val raw = builder.addPair(
219+
contig = 0, start1 = 96, start2 = 49,
220+
cigar1 = "68S53M8S", cigar2 = "28S48M53S",
221+
attrs = Map(("RX", "ACC-TGA"), ("MI", "hi"))
222+
).tapEach(setReadSequence)
223+
224+
val caller = new CodecConsensusCaller(readNamePrefix="codec", minReadsPerStrand=1, minDuplexLength=1)
225+
noException should be thrownBy caller.consensusReadsFromSamRecords(raw)
226+
}
227+
228+
it should "not throw on a template whose primary read-name bucket has only one record" in {
229+
// Exercise the defensive path for a degenerate template that reaches the CODEC caller with
230+
// only one primary record per read name (e.g. an orphan introduced upstream). Previously
231+
// the `Seq(rec, mate)` pattern match at line 176 crashed with MatchError.
232+
val builder = new SamBuilder(readLength=30, baseQuality=35)
233+
val pair = builder.addPair(
234+
contig=0, start1=1, start2=11, attrs=Map(("RX", "ACC-TGA"), ("MI", "hi"))
235+
).tapEach(setReadSequence)
236+
val onlyR1 = pair.filter(_.firstOfPair)
237+
238+
val caller = new CodecConsensusCaller(readNamePrefix="codec", minReadsPerStrand=1, minDuplexLength=1)
239+
caller.consensusReadsFromSamRecords(onlyR1) shouldBe Seq()
240+
}
241+
210242
it should "not emit a consensus when the reads are a cross-chromosomal chimeric pair" in {
211243
val builder = new SamBuilder(readLength=30, baseQuality=35)
212244
val caller = new CodecConsensusCaller(readNamePrefix="codec", minReadsPerStrand=1, minDuplexLength=1)
@@ -290,6 +322,86 @@ class CodecConsensusCallerTest extends UnitSpec with OptionValues {
290322
.consensusReadsFromSamRecords(raw) should have length 0
291323
}
292324

325+
//////////////////////////////////////////////////////////////////////////////
326+
// Tests for orderedPrimaryPairs
327+
//////////////////////////////////////////////////////////////////////////////
328+
329+
"CodecConsensusCaller.orderedPrimaryPairs" should "order templates by first-occurrence BAM-iteration order" in {
330+
// Names chosen so that Scala's `HashMap` iteration order on `groupBy(_.name).toSeq` is
331+
// extremely unlikely to coincide with insertion order; the helper must return them in the
332+
// order they were observed in the input regardless of hash layout.
333+
val builder = new SamBuilder(readLength=30, baseQuality=35)
334+
val names = Seq("zulu", "alpha", "mike", "bravo", "yankee")
335+
names.foreach { n =>
336+
builder.addPair(name=n, contig=0, start1=1, start2=11).tapEach(setReadSequence)
337+
}
338+
val ordered = CodecConsensusCaller.orderedPrimaryPairs(builder.toSeq)
339+
ordered.map(_._1) shouldBe names
340+
ordered.foreach { case (_, recs) => recs should have size 2 }
341+
}
342+
343+
it should "ignore secondary/supplementary records only when the caller pre-filters them" in {
344+
// The helper itself does not filter; callers pass pre-filtered primaries. This test documents
345+
// that contract by passing records that include a supplementary alignment and verifying it is
346+
// grouped alongside the primary under the same read name.
347+
val builder = new SamBuilder(readLength=30, baseQuality=35)
348+
val pair = builder.addPair(name="t1", contig=0, start1=1, start2=11).tapEach(setReadSequence)
349+
pair.head.supplementary = true
350+
val ordered = CodecConsensusCaller.orderedPrimaryPairs(builder.toSeq)
351+
ordered.map(_._1) shouldBe Seq("t1")
352+
ordered.head._2 should have size 2
353+
}
354+
355+
//////////////////////////////////////////////////////////////////////////////
356+
// Tests for isPrimaryFrPair
357+
//////////////////////////////////////////////////////////////////////////////
358+
359+
"CodecConsensusCaller.isPrimaryFrPair" should "return the same answer regardless of argument order" in {
360+
val builder = new SamBuilder(readLength=129, baseQuality=35)
361+
val pair = builder.addPair(
362+
contig=0, start1=96, start2=49, cigar1="68S53M8S", cigar2="28S48M53S"
363+
).tapEach(setReadSequence)
364+
val r1 = pair.head
365+
val r2 = pair.last
366+
CodecConsensusCaller.isPrimaryFrPair(r1, r2) shouldBe CodecConsensusCaller.isPrimaryFrPair(r2, r1)
367+
}
368+
369+
it should "classify a typical FR pair as FR" in {
370+
val builder = new SamBuilder(readLength=30, baseQuality=35)
371+
val Seq(r1, r2) = builder.addPair(contig=0, start1=10, start2=50).toSeq
372+
CodecConsensusCaller.isPrimaryFrPair(r1, r2) shouldBe true
373+
}
374+
375+
it should "classify an RF pair as not FR" in {
376+
val builder = new SamBuilder(readLength=30, baseQuality=35)
377+
val Seq(r1, r2) = builder.addPair(
378+
contig=0, start1=50, start2=100, strand1=Minus, strand2=Plus
379+
).toSeq
380+
CodecConsensusCaller.isPrimaryFrPair(r1, r2) shouldBe false
381+
}
382+
383+
it should "reject cross-chromosomal pairs" in {
384+
val builder = new SamBuilder(readLength=30, baseQuality=35)
385+
val Seq(r1, r2) = builder.addPair(contig=0, start1=10, start2=50).toSeq
386+
r1.refIndex = 0
387+
r2.refIndex = 1
388+
CodecConsensusCaller.isPrimaryFrPair(r1, r2) shouldBe false
389+
}
390+
391+
it should "reject a pair with one mate unmapped" in {
392+
val builder = new SamBuilder(readLength=30, baseQuality=35)
393+
val Seq(r1, r2) = builder.addPair(contig=0, start1=10, start2=10, unmapped2=true).toSeq
394+
CodecConsensusCaller.isPrimaryFrPair(r1, r2) shouldBe false
395+
}
396+
397+
it should "reject a tandem pair" in {
398+
val builder = new SamBuilder(readLength=30, baseQuality=35)
399+
val Seq(r1, r2) = builder.addPair(
400+
contig=0, start1=10, start2=50, strand1=Plus, strand2=Plus
401+
).toSeq
402+
CodecConsensusCaller.isPrimaryFrPair(r1, r2) shouldBe false
403+
}
404+
293405
//////////////////////////////////////////////////////////////////////////////
294406
// Tests for quality masking
295407
//////////////////////////////////////////////////////////////////////////////

0 commit comments

Comments
 (0)