feat: raise exception if CollectDuplexSeqMetrics run on consensus BAM

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

@@ -307,6 +307,18 @@ class CollectDuplexSeqMetrics
     // Build the iterator we'll use based on whether or not we're restricting to a set of intervals
     val in = SamSource(input)
     val _filteredIterator = in.iterator.filter(r => r.paired && r.mapped && r.mateMapped && r.firstOfPair && !r.secondary && !r.supplementary)
+      .bufferBetter
+    // Ensure the records are not consensus records
+    _filteredIterator.headOption.foreach {
+      r =>
+        def exceptionString = s"Input BAM file to CollectDuplexSeqMetrics ($input) appears to contain consensus reads." +
+          "CollectDuplexSeqMetrics cannot run on consensus BAMs, and instead requires the UMI-grouped BAM generated " +
+          "prior to consensus calling. The UMI-grouped BAM is the output of running GroupReadsByUmi." +
+          s"\nFirst record in $input has consensus tags present:\n${_filteredIterator.head}"
+
+        if (Umis.isConsensusRead(r)) throw new IllegalArgumentException(exceptionString)
+        else r
+    }
     val iterator = intervals match {
       case None       => _filteredIterator
       case Some(path) =>

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

@@ -26,6 +26,7 @@
 package com.fulcrumgenomics.umi
 
 import com.fulcrumgenomics.bam.api.SamRecord
+import com.fulcrumgenomics.umi.ConsensusTags
 import com.fulcrumgenomics.util.Sequences
 
 object Umis {
@@ -127,4 +128,14 @@ object Umis {
   @inline private def isValidUmiCharacter(ch: Char): Boolean = {
     ch == 'A' || ch == 'C' || ch == 'G' || ch == 'T' || ch == 'N' || ch == '-'
   }
+
+  /** Tests if a record is a consensus or not
+   *
+   * @param rec the record to test
+   * @return boolean indicating if the record is a consensus or not
+   */
+  def isConsensusRead(rec: SamRecord): Boolean = {
+    rec.contains(ConsensusTags.PerRead.RawReadCount) ||
+      (rec.contains(ConsensusTags.PerRead.AbRawReadCount) && rec.contains(ConsensusTags.PerRead.BaRawReadCount))
+  }
 }

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

@@ -26,12 +26,12 @@ package com.fulcrumgenomics.umi
 
 import java.nio.file.{Path, Paths}
 import java.util.Random
-
 import com.fulcrumgenomics.FgBioDef._
 import com.fulcrumgenomics.bam.api.SamOrder
 import com.fulcrumgenomics.commons.util.SimpleCounter
 import com.fulcrumgenomics.testing.SamBuilder.{Minus, Plus}
 import com.fulcrumgenomics.testing.{SamBuilder, UnitSpec}
+import com.fulcrumgenomics.umi.ConsensusTags
 import com.fulcrumgenomics.util.{Io, Metric, Rscript}
 import htsjdk.samtools.util.{Interval, IntervalList}
 import org.apache.commons.math3.distribution.NormalDistribution
@@ -41,6 +41,7 @@ import scala.math.{max, min}
 class CollectDuplexSeqMetricsTest extends UnitSpec {
   private val MI = "MI"
   private val RX = "RX"
+  private val aD = "aD"
 
   // Case class to hold pointers to all the outputs of CollectDuplexSeqMetrics
   private case class Outputs(families: Path, duplexFamilies: Path, umis: Path, duplexUmis: Path, yields: Path, plots: Path) {
@@ -312,6 +313,22 @@ class CollectDuplexSeqMetricsTest extends UnitSpec {
     duplexFamilies.find(f => f.ab_size == 6 && f.ba_size == 0).get.count shouldBe 1
   }
 
+  it should "raise an exception if duplex consensus records are provided" in {
+    val builder = new SamBuilder(readLength=100, sort=Some(SamOrder.TemplateCoordinate))
+    builder.addPair(
+      contig=1,
+      start1=1000,
+      start2=1100,
+      attrs=Map(
+        RX -> "AAA-GGG",
+        MI -> "1/A",
+        ConsensusTags.PerRead.AbRawReadCount -> 10,
+        ConsensusTags.PerRead.BaRawReadCount -> 10
+      )
+    )
+    an[IllegalArgumentException] shouldBe thrownBy { exec(builder) }
+  }
+
   "CollectDuplexSeqMetrics.updateUmiMetrics" should "not count duplex umis" in collector(duplex=false).foreach { c =>
     val builder = new SamBuilder(readLength=10)
     builder.addPair(start1=100, start2=200, attrs=Map(RX -> "AAA-CCC", MI -> "1/A"))

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

@@ -25,8 +25,9 @@
 
 package com.fulcrumgenomics.umi
 
-import com.fulcrumgenomics.bam.api.SamRecord
+import com.fulcrumgenomics.bam.api.{SamOrder, SamRecord}
 import com.fulcrumgenomics.testing.{SamBuilder, UnitSpec}
+import com.fulcrumgenomics.umi.ConsensusTags
 import org.scalatest.OptionValues
 
 class UmisTest extends UnitSpec with OptionValues {
@@ -111,4 +112,20 @@ class UmisTest extends UnitSpec with OptionValues {
     an[Exception] should be thrownBy copyUmiFromReadName(rec=rec("NAME:ACGT-CCKC"))
     an[Exception] should be thrownBy copyUmiFromReadName(rec=rec("NAME:CCKC-ACGT"))
   }
+
+  "Umis.isConsensusRead" should "return false for reads without consensus tags" in {
+    val builder = new SamBuilder(sort=Some(SamOrder.Coordinate), readLength=10, baseQuality=20)
+    builder.addFrag(start=100).exists(Umis.isConsensusRead) shouldBe false
+    builder.addPair(start1=100, start2=100, unmapped2=true).exists(Umis.isConsensusRead) shouldBe false
+  }
+
+  it should "return true for reads with consensus tags" in {
+    val builder = new SamBuilder(sort=Some(SamOrder.Coordinate), readLength=10, baseQuality=20)
+    builder.addFrag(start=10, attrs=Map(ConsensusTags.PerRead.RawReadCount -> 10))
+      .exists(Umis.isConsensusRead) shouldBe true
+    builder.addFrag(
+        start=10,
+        attrs=Map(ConsensusTags.PerRead.AbRawReadCount -> 10, ConsensusTags.PerRead.BaRawReadCount -> 10)
+    ).exists(Umis.isConsensusRead) shouldBe true
+  }
 }