From da9ecbcc55d746b752049a1fecd6f4eeacbad3cb Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Thu, 28 Jul 2022 14:09:12 -0700 Subject: [PATCH] [bugfix] fix issue #858 in CallDuplexConsensusReads (#864) * [bugfix] fix issue #858 in CallDuplexConsensusReads Fixes #858 where the AB strand single strand consensus read could have zero depth and all Ns, such that the consensus read produced has an aD value of zero but bD value > 0, causing issues in FilterConsensus Reads. Co-authored-by: Tim Fennell --- .../umi/DuplexConsensusCaller.scala | 9 ++++-- .../umi/DuplexConsensusCallerTest.scala | 28 ++++++++++++++++++- 2 files changed, 34 insertions(+), 3 deletions(-) diff --git a/src/main/scala/com/fulcrumgenomics/umi/DuplexConsensusCaller.scala b/src/main/scala/com/fulcrumgenomics/umi/DuplexConsensusCaller.scala index d402cfe40..00434b71c 100644 --- a/src/main/scala/com/fulcrumgenomics/umi/DuplexConsensusCaller.scala +++ b/src/main/scala/com/fulcrumgenomics/umi/DuplexConsensusCaller.scala @@ -340,11 +340,16 @@ class DuplexConsensusCaller(override val readNamePrefix: String, private[umi] def duplexConsensus(ab: Option[VanillaConsensusRead], ba: Option[VanillaConsensusRead], sourceReads: Seq[SourceRead]): Option[DuplexConsensusRead] = { - (ab, ba) match { + // Calculate the length that we'll create a duplex read as "the shorter of the available reads", + // and then filter each read to make sure it has some coverage in that section of the read + val len = min(ab.map(_.length).getOrElse(Int.MaxValue), ba.map(_.length).getOrElse(Int.MaxValue)) + val abX = ab.filter(_.depths.iterator.take(len).exists(_ > 0)) + val baX = ba.filter(_.depths.iterator.take(len).exists(_ > 0)) + + (abX, baX) match { case (Some(a), None) => Some(DuplexConsensusRead(id=a.id, a.bases, a.quals, a.errors, a, None)) case (None, Some(b)) => Some(DuplexConsensusRead(id=b.id, b.bases, b.quals, b.errors, b, None)) case (Some(a), Some(b)) => - val len = min(a.length, b.length) val id = a.id val bases = new Array[Byte](len) val quals = new Array[Byte](len) diff --git a/src/test/scala/com/fulcrumgenomics/umi/DuplexConsensusCallerTest.scala b/src/test/scala/com/fulcrumgenomics/umi/DuplexConsensusCallerTest.scala index 25047dcbc..e51fad2bd 100644 --- a/src/test/scala/com/fulcrumgenomics/umi/DuplexConsensusCallerTest.scala +++ b/src/test/scala/com/fulcrumgenomics/umi/DuplexConsensusCallerTest.scala @@ -28,8 +28,9 @@ import com.fulcrumgenomics.testing.SamBuilder.{Minus, Plus} import com.fulcrumgenomics.testing.{SamBuilder, UnitSpec} import com.fulcrumgenomics.umi.ConsensusTags.PerRead.{AbRawReadErrorRate, BaRawReadErrorRate, RawReadErrorRate} import com.fulcrumgenomics.util.NumericTypes.PhredScore +import org.scalatest.OptionValues -class DuplexConsensusCallerTest extends UnitSpec { +class DuplexConsensusCallerTest extends UnitSpec with OptionValues { // Function to create a caller without so much writing def caller(q: Int = 10, pre: Int = DuplexConsensusCaller.ErrorRatePreUmi, post: Int = DuplexConsensusCaller.ErrorRatePostUmi, minReads: Seq[Int] = Seq(1)) = new DuplexConsensusCaller(readNamePrefix="test", minInputBaseQuality=q.toByte, errorRatePreUmi=pre.toByte, errorRatePostUmi=post.toByte, @@ -468,4 +469,29 @@ class DuplexConsensusCallerTest extends UnitSpec { } } } + "DuplexConsensusCaller.duplexConsensus" should "swap AB and BA reads when AB reads have zero depth across the minimum length" in { + val cc = caller(minReads=Seq(1, 1, 0)) + val ab = VanillaConsensusRead( + id = "ab", + bases = "NNN".getBytes, + quals = "III".getBytes, + depths = Array(0, 0, 0), + errors = Array(0, 0, 0) + ) + val ba = ab.copy(depths=Array(1, 1, 1)) + + // ab has zero depth, so it should be swapped with BA + { + val read = cc.duplexConsensus(ab=Some(ab), ba=Some(ba), sourceReads=Seq.empty).value + read.abConsensus shouldBe ba + read.baConsensus.isEmpty shouldBe true + } + + // ba has zero depth, so it should be undefined. + { + val read = cc.duplexConsensus(ab=Some(ba), ba=Some(ab), sourceReads=Seq.empty).value + read.abConsensus shouldBe ba + read.baConsensus.isEmpty shouldBe true + } + } }