fix(umi): exclude unmapped reads from consensus calling

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

@@ -53,9 +53,13 @@ class ConsensusCallingIterator[ConsensusRead <: SimpleRead](sourceIterator: Iter
   private var collectedStats: Boolean = false
 
   protected val iter: Iterator[SamRecord] = {
+    // Only mapped reads may contribute to a consensus.  An unmapped read has no cigar, and an empty cigar is a prefix
+    // of every cigar, so `filterToMostCommonAlignment` would match it against whichever alignment group it was tested
+    // against first rather than rejecting it as a minority alignment.  Note that the mapped end of a half-mapped pair
+    // is still used; only the unmapped end is dropped.
     val filteredIterator = sourceIterator
       .filterNot (r => r.secondary || r.supplementary)
-      .filter(r => r.mapped || (r.paired && r.mateMapped))
+      .filter(r => r.mapped)
 
     // Wrap our input iterator in a progress logging iterator if we have a progress logger
     val progressIterator = progress match {

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

@@ -126,4 +126,27 @@ class CallCodecConsensusReadsTest extends UnitSpec with OptionValues {
       metric.value.toString.toInt shouldBe 2
     }
   }
+
+  it should "not use either end of a half-mapped pair when calling a consensus" in {
+    val builder = new SamBuilder(readLength=30, sort=Some(SamOrder.TemplateCoordinate))
+    val attrs   = Map(("RX", "ACC-TGA"), ("MI", "hi"))
+    builder.addPair(name="mapped1", start1=100, start2=100, bases1="AC" * 15, bases2="AC" * 15, attrs=attrs)
+    builder.addPair(name="mapped2", start1=100, start2=100, bases1="AC" * 15, bases2="AC" * 15, attrs=attrs)
+    // A template whose R2 is unmapped is not a primary FR pair, so neither end may contribute to the consensus.
+    builder.addPair(name="halfmapped", start1=100, start2=100, bases1="AC" * 15, bases2="GT" * 15, unmapped2=true, attrs=attrs)
+
+    val out = makeTempFile("codec.", ".bam")
+    val rej = makeTempFile("rejects.", ".bam")
+    new CallCodecConsensusReads(input=builder.toTempFile(), output=out, readGroupId="ZZ", rejects=Some(rej)).execute()
+
+    val recs = readBamRecs(out)
+    recs should have size 1
+
+    // Only the two fully mapped templates contribute to each single-strand consensus.
+    recs.head[Int](ConsensusTags.PerRead.AbRawReadCount) shouldBe 2
+    recs.head[Int](ConsensusTags.PerRead.BaRawReadCount) shouldBe 2
+
+    // Every rejected record belongs to the half-mapped template.
+    readBamRecs(rej).foreach { rec => rec.name shouldBe "halfmapped" }
+  }
 }

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

@@ -132,4 +132,31 @@ class CallDuplexConsensusReadsTest extends UnitSpec {
       }
     }
   }
+
+  it should "not use the unmapped end of a half-mapped pair when calling a single-strand consensus" in {
+    val builder = new SamBuilder(readLength=10, sort=Some(SamOrder.TemplateCoordinate))
+    // Three AB templates, the third of which has an unmapped R2 carrying different bases to the mapped R2s.
+    builder.addPair(name="ab1", start1=100, start2=100, attrs=Map(MI -> "1/A"), bases1="AAAAAAAAAA", bases2="AAAAAAAAAA")
+    builder.addPair(name="ab2", start1=100, start2=100, attrs=Map(MI -> "1/A"), bases1="AAAAAAAAAA", bases2="AAAAAAAAAA")
+    builder.addPair(name="ab3", start1=100, start2=100, attrs=Map(MI -> "1/A"), bases1="AAAAAAAAAA", bases2="CCCCCCCCCC", unmapped2=true)
+    // Two fully mapped BA templates.
+    builder.addPair(name="ba1", start1=100, start2=100, strand1=Minus, strand2=Plus, attrs=Map(MI -> "1/B"), bases1="AAAAAAAAAA", bases2="AAAAAAAAAA")
+    builder.addPair(name="ba2", start1=100, start2=100, strand1=Minus, strand2=Plus, attrs=Map(MI -> "1/B"), bases1="AAAAAAAAAA", bases2="AAAAAAAAAA")
+
+    val out = makeTempFile("duplex.", ".bam")
+    new CallDuplexConsensusReads(input=builder.toTempFile(), output=out, readGroupId="ZZ").execute()
+
+    val recs = readBamRecs(out)
+    recs should have size 2
+    val r1 = recs.find(_.firstOfPair).value
+    val r2 = recs.find(_.secondOfPair).value
+
+    // The AB-R1s are all mapped, so all three contribute to the AB single-strand consensus for R1.
+    r1[Int](ConsensusTags.PerRead.AbRawReadCount) shouldBe 3
+    r1[Int](ConsensusTags.PerRead.BaRawReadCount) shouldBe 2
+
+    // Only the two mapped AB-R2s may contribute to the AB single-strand consensus for R2; the unmapped one must not.
+    r2[Int](ConsensusTags.PerRead.AbRawReadCount) shouldBe 2
+    r2[Int](ConsensusTags.PerRead.BaRawReadCount) shouldBe 2
+  }
 }

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

@@ -110,4 +110,54 @@ class CallMolecularConsensusReadsTest extends UnitSpec {
       rec[String](specialCellTag) shouldBe "AB"
     }
   }
+
+  /** Builds a tag family of three read pairs at the same coordinates, two with both ends mapped and one whose R2 is
+    * unmapped.
+    *
+    * The two mapped R2s carry `G`s while the unmapped R2 carries `T`s, so the consensus R2 shows which reads
+    * contributed to it.  R2s are placed on the negative strand by [[SamBuilder]], so `toSourceRead` reverse
+    * complements them and a consensus built from the mapped R2s alone reads as `C`s. */
+  private def halfMappedFamily(rlen: Int): SamBuilder = {
+    val builder = new SamBuilder(baseQuality=30, readLength=rlen, readGroupId=Some("ABC"), sort=Some(SamOrder.TemplateCoordinate))
+    val attrs   = Map(DefaultTag -> "GATTACA:1")
+    builder.addPair(name="mapped:1",     start1=100, start2=300, bases1="A"*rlen, bases2="G"*rlen, attrs=attrs)
+    builder.addPair(name="mapped:2",     start1=100, start2=300, bases1="A"*rlen, bases2="G"*rlen, attrs=attrs)
+    builder.addPair(name="halfmapped:3", start1=100, start2=300, bases1="A"*rlen, bases2="T"*rlen, unmapped2=true, attrs=attrs)
+    builder
+  }
+
+  Seq(2, 3).foreach { maxReads =>
+    it should f"cap R1 at max-reads and build R2 from the mapped reads only with --max-reads=$maxReads" in {
+      val rlen    = 100
+      val output  = newBam
+      val rejects = newBam
+
+      new CallMolecularConsensusReads(
+        input       = halfMappedFamily(rlen).toTempFile(),
+        output      = output,
+        minReads    = 1,
+        maxReads    = Some(maxReads),
+        rejects     = Some(rejects),
+        readGroupId = "ABC"
+      ).execute()
+
+      // Exactly one consensus read pair is produced.
+      val records = readBamRecs(output)
+      records.size shouldBe 2
+      val r1 = records.find(_.firstOfPair).value
+      val r2 = records.find(_.secondOfPair).value
+
+      // All three R1s are mapped and eligible, so R1 is capped at max-reads.
+      r1[Int](ConsensusTags.PerRead.RawReadCount) shouldBe math.min(3, maxReads)
+      r1.basesString shouldBe "A" * rlen
+
+      // Only the two mapped R2s may contribute, whether or not the cap binds.  The unmapped R2's `C`s must not reach
+      // the consensus, and it must not be counted towards the consensus depth.
+      r2[Int](ConsensusTags.PerRead.RawReadCount) shouldBe 2
+      r2.basesString shouldBe "C" * rlen
+
+      // The unmapped R2 must not have been used, whether it was rejected by the caller or filtered before reaching it.
+      readBamRecs(rejects).foreach { rec => rec.name shouldBe "halfmapped:3" }
+    }
+  }
 }