[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 <[email protected]>
fulcrumgenomics · Jul 28, 2022 · da9ecbc · da9ecbc
1 parent 74d920e
commit da9ecbc
Show file tree

Hide file tree

Showing 2 changed files with 34 additions and 3 deletions.
diff --git 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
@@ -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
+    }
+  }
 }