Change how per-base errors are counted when calling consensuses. (#166)

* Change how per-base errors are counted when calling consensuses.

When the final consensus call is a no-call (N), we instead count the
number of errors relative to the most likely consensus call, whereas we
previously used the total depth.

Author: nh13

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

@@ -50,7 +50,8 @@ object ConsensusTags {
   object PerBase {
     /** The per-base number of raw-reads contributing to the consensus (stored as a short[]). */
     val RawReadCount    = "cd" // consensus depth
-    /** The number of bases at each position that disagreed with the final consensus call (stored as a short[]). */
+    /** The number of bases at each position that disagreed with the final consensus call (stored as a short[]). If the
+      * final consensus call is a no call (N), then we use the most likely consensus base instead of the final call.  */
     val RawReadErrors   = "ce" // consensus errors
 
     // Duplex versions of the above tags for the two single strand consensus reads

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

@@ -187,23 +187,19 @@ class VanillaUmiConsensusCaller(override val readNamePrefix: String,
         }
 
         // Call the consensus and do any additional filtering
+        val (rawBase, rawQual) = builder.call()
         val (base, qual) = {
-          if (builder.contributions < this.options.minReads) {
-            (NoCall, NotEnoughReadsQual)
-          }
-          else {
-            val (b,q) = builder.call()
-            if (q < this.options.minConsensusBaseQuality) (NoCall, TooLowQualityQual) else (b,q)
-
-          }
+          if (builder.contributions < this.options.minReads)       (NoCall, NotEnoughReadsQual)
+          else if (rawQual < this.options.minConsensusBaseQuality) (NoCall, TooLowQualityQual)
+          else (rawBase, rawQual)
         }
 
         consensusBases(positionInRead) = base
         consensusQuals(positionInRead) = qual
 
         // Generate the values for depth and count of errors
         val depth  = builder.contributions
-        val errors = if (base == NoCall) depth else depth - builder.observations(base)
+        val errors = if (rawBase == NoCall) depth else depth - builder.observations(rawBase)
         consensusDepths(positionInRead) = if (depth  > Short.MaxValue) Short.MaxValue else depth.toShort
         consensusErrors(positionInRead) = if (errors > Short.MaxValue) Short.MaxValue else errors.toShort
 

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

@@ -33,11 +33,12 @@ import com.fulcrumgenomics.util.NumericTypes._
 import htsjdk.samtools.util.CloserUtil
 import htsjdk.samtools.{SAMRecordSetBuilder, SAMUtils}
 import net.jafama.FastMath._
+import org.scalatest.OptionValues
 
 /**
   * Tests for ConsensusCaller.
   */
-class VanillaUmiConsensusCallerTest extends UnitSpec {
+class VanillaUmiConsensusCallerTest extends UnitSpec with OptionValues {
   /** Helper function to make a set of consensus caller options. */
   def cco = VanillaUmiConsensusCallerOptions
 
@@ -423,4 +424,18 @@ class VanillaUmiConsensusCallerTest extends UnitSpec {
     recs should have size 11
     recs.map(_.getCigarString).distinct.sorted shouldBe Seq("25M1D25M", "5S20M1D25M", "5S20M1D20M5H", "25M1D20M5S").sorted
   }
+
+  it should "calculate the # of errors relative to the most likely consensus call, even when the final call is an N" in {
+    // NB: missing last base on the first read, which causes an N no-call, but errors should still be 1/3
+    val call = cc(cco(minReads=4)).consensusCall(Seq(
+      src("GATTACAN", Array(20,20,20,20,20,20,20,20)),
+      src("GATTACAG", Array(20,20,20,20,20,20,20,20)),
+      src("GATTACAG", Array(20,20,20,20,20,20,20,20)),
+      src("GATTACAT", Array(20,20,20,20,20,20,20,20))
+    )).value
+
+    call.baseString shouldBe "GATTACAN"
+    call.depths should contain theSameElementsInOrderAs Seq(4,4,4,4,4,4,4,3)
+    call.errors should contain theSameElementsInOrderAs Seq(0,0,0,0,0,0,0,1)
+  }
 }