Fix clipping issues when reads overlap

Author: clintval

src/main/scala/com/fulcrumgenomics/bam/ClipBam.scala

@@ -126,6 +126,8 @@ class ClipBam
         case (Some(r1), Some(r2)) =>
           clipPair(r1=r1, r2=r2, r1Metric=metricsMap.get(ReadType.ReadOne), r2Metric=metricsMap.get(ReadType.ReadTwo))
           SamPairUtil.setMateInfo(r1.asSam, r2.asSam, true)
+          // SamPairUtil.setMateInfo does not clear the proper-pair (0x2) flag when one read becomes unmapped
+          if (r1.unmapped || r2.unmapped) { r1.properlyPaired = false; r2.properlyPaired = false }
           template.r1Supplementals.foreach(s => SamPairUtil.setMateInformationOnSupplementalAlignment(s.asSam, r2.asSam, true))
           template.r2Supplementals.foreach(s => SamPairUtil.setMateInformationOnSupplementalAlignment(s.asSam, r1.asSam, true))
         case (Some(frag), None) =>

src/main/scala/com/fulcrumgenomics/bam/SamRecordClipper.scala

@@ -314,9 +314,14 @@ class SamRecordClipper(val mode: ClippingMode, val autoClipAttributes: Boolean)
         else retval
       }
       val readEnd   = rec.readPosAtRefPos(pos=midPoint, returnLastBaseIfDeleted=true)
-      val mateStart = { // NB: need to be careful if the midpoint falls in a deletion
-        val retval = mate.readPosAtRefPos(pos=midPoint + 1, returnLastBaseIfDeleted=false)
-        if (retval != 0) retval else mate.readPosAtRefPos(pos=midPoint + 1, returnLastBaseIfDeleted=true) + 1
+      val mateStart = { // NB: need to be careful if the midpoint falls in a deletion or past the end of the mate
+        if (midPoint >= mate.end) {
+          // The midpoint is at or past the mate's last aligned base: clip all of the mate.
+          mate.length + 1
+        } else {
+          val retval = mate.readPosAtRefPos(pos=midPoint + 1, returnLastBaseIfDeleted=false)
+          if (retval != 0) retval else mate.readPosAtRefPos(pos=midPoint + 1, returnLastBaseIfDeleted=true) + 1
+        }
       }
       val numOverlappingBasesRead = this.clip3PrimeEndOfRead(rec, rec.cigar.trailingHardClippedBases + rec.length - readEnd)
       val numOverlappingBasesMate = this.clip3PrimeEndOfRead(mate, mate.cigar.leadingHardClippedBases + mateStart - 1)
@@ -354,7 +359,7 @@ class SamRecordClipper(val mode: ClippingMode, val autoClipAttributes: Boolean)
     * @return the additional number of bases clipped (3' end in sequencing order) for the read and mate respectively
     */
   def clipExtendingPastMateEnds(rec: SamRecord, mate: SamRecord): (Int, Int) = {
-    if (rec.isFrPair) {
+    if (rec.isFrPair && rec.mapped && mate.mapped) {
       val basesClipped1 = clipExtendingPastMateEnd(
         rec                    = rec,
         mateUnSoftClippedStart = mate.unSoftClippedStart,

src/test/scala/com/fulcrumgenomics/bam/ClipBamTest.scala

@@ -618,6 +618,240 @@ class ClipBamTest extends UnitSpec with ErrorLogLevel with OptionValues {
     r2.end shouldBe 90 + 200 - 1 // 289: the original end of R2
   }
 
+  /** Helper that asserts two reads have the same start, end, and CIGAR. */
+  private def assertSameClipping(got: SamRecord, expected: SamRecord, label: String): Unit = {
+    withClue(s"$label start:") { got.start shouldBe expected.start }
+    withClue(s"$label end:")   { got.end shouldBe expected.end }
+    withClue(s"$label cigar:") { got.cigar.toString shouldBe expected.cigar.toString }
+    ()
+  }
+
+  // Tests for the ordering issue raised in https://github.com/fulcrumgenomics/fgbio/issues/1141#issuecomment-3998369518
+  // The order of --clip-overlapping-reads and --clip-bases-past-mate should not affect the final result.
+  // Each test runs the two operations in both orderings and asserts they produce identical results.
+
+  private val samClipper = new SamRecordClipper(mode = ClippingMode.Hard, autoClipAttributes = false)
+
+  it should "produce the same result for clip-overlapping-reads then clip-bases-past-mate as the reverse order for symmetrically dovetailed reads" in {
+    // R1 and R2 dovetail symmetrically: each extends 10 bp past the other's 5' start
+    val Seq(r1a, r2a) = new SamBuilder(readLength = 100).addPair(start1 = 100, start2 = 90)
+    val Seq(r1b, r2b) = new SamBuilder(readLength = 100).addPair(start1 = 100, start2 = 90)
+
+    // Order A (current production order): overlapping → past-mate
+    samClipper.clipOverlappingReads(r1a, r2a)
+    samClipper.clipExtendingPastMateEnds(r1a, r2a)
+
+    // Order B (reversed): past-mate → overlapping
+    samClipper.clipExtendingPastMateEnds(r1b, r2b)
+    samClipper.clipOverlappingReads(r1b, r2b)
+
+    assertSameClipping(r1a, r1b, "R1")
+    assertSameClipping(r2a, r2b, "R2")
+  }
+
+  it should "produce the same result for clip-overlapping-reads then clip-bases-past-mate as the reverse order for heavily dovetailed reads" in {
+    // R1 and R2 heavily dovetail: R2 starts 50 bp before R1's 5' start
+    val Seq(r1a, r2a) = new SamBuilder(readLength = 100).addPair(start1 = 100, start2 = 50)
+    val Seq(r1b, r2b) = new SamBuilder(readLength = 100).addPair(start1 = 100, start2 = 50)
+
+    // Order A: overlapping → past-mate
+    samClipper.clipOverlappingReads(r1a, r2a)
+    samClipper.clipExtendingPastMateEnds(r1a, r2a)
+
+    // Order B: past-mate → overlapping
+    samClipper.clipExtendingPastMateEnds(r1b, r2b)
+    samClipper.clipOverlappingReads(r1b, r2b)
+
+    assertSameClipping(r1a, r1b, "R1")
+    assertSameClipping(r2a, r2b, "R2")
+  }
+
+  it should "produce the same result for clip-overlapping-reads then clip-bases-past-mate as the reverse order when R2 is much longer and extends past both ends of R1" in {
+    // R2 is much longer than R1 and dovetails, extending past both R1.start and R1.end.
+    // This is the scenario eboyden described where the shorter read could be clipped to zero length.
+    // R1 is 20 bases long; R2 is 100 bases long.
+    val bases20  = "A" * 20
+    val quals20  = "F" * 20
+    val bases100 = "A" * 100
+    val quals100 = "F" * 100
+    val Seq(r1a, r2a) = new SamBuilder(readLength = 100).addPair(
+      start1 = 100, start2 = 80, cigar1 = "20M", cigar2 = "100M",
+      bases1 = bases20, quals1 = quals20, bases2 = bases100, quals2 = quals100,
+    )
+    val Seq(r1b, r2b) = new SamBuilder(readLength = 100).addPair(
+      start1 = 100, start2 = 80, cigar1 = "20M", cigar2 = "100M",
+      bases1 = bases20, quals1 = quals20, bases2 = bases100, quals2 = quals100,
+    )
+
+    // Order A: overlapping → past-mate
+    samClipper.clipOverlappingReads(r1a, r2a)
+    samClipper.clipExtendingPastMateEnds(r1a, r2a)
+
+    // Order B: past-mate → overlapping
+    samClipper.clipExtendingPastMateEnds(r1b, r2b)
+    samClipper.clipOverlappingReads(r1b, r2b)
+
+    assertSameClipping(r1a, r1b, "R1")
+    assertSameClipping(r2a, r2b, "R2")
+  }
+
+  it should "produce the same result for clip-overlapping-reads then clip-bases-past-mate as the reverse order when R1 is much longer and contains a short R2" in {
+    // R1 is much longer; R2 is short and starts late within R1's range.
+    // The concern is that clipOverlappingReads clips the midpoint to well before R2.start,
+    // dramatically shortening R1 but leaving R2 nearly intact.
+    val bases100 = "A" * 100
+    val quals100 = "F" * 100
+    val bases20  = "A" * 20
+    val quals20  = "F" * 20
+    val Seq(r1a, r2a) = new SamBuilder(readLength = 100).addPair(
+      start1 = 100, start2 = 180, cigar1 = "100M", cigar2 = "20M",
+      bases1 = bases100, quals1 = quals100, bases2 = bases20, quals2 = quals20,
+    )
+    val Seq(r1b, r2b) = new SamBuilder(readLength = 100).addPair(
+      start1 = 100, start2 = 180, cigar1 = "100M", cigar2 = "20M",
+      bases1 = bases100, quals1 = quals100, bases2 = bases20, quals2 = quals20,
+    )
+
+    // Order A: overlapping → past-mate
+    samClipper.clipOverlappingReads(r1a, r2a)
+    samClipper.clipExtendingPastMateEnds(r1a, r2a)
+
+    // Order B: past-mate → overlapping
+    samClipper.clipExtendingPastMateEnds(r1b, r2b)
+    samClipper.clipOverlappingReads(r1b, r2b)
+
+    assertSameClipping(r1a, r1b, "R1")
+    assertSameClipping(r2a, r2b, "R2")
+  }
+
+  it should "produce the same result for clip-overlapping-reads then clip-bases-past-mate as the reverse order with 5-prime hard-clipped reads (issue 1141 scenario)" in {
+    // Reproduces the asymmetric hard-clipped scenario from issue #1141:
+    //   - R1 has 25H at its 5' (left) end
+    //   - R2 has 25H at its 5' (right/trailing-CIGAR) end
+    //   - Both reads dovetail: R2 starts before R1's reference start
+    // In this scenario the mate's hard-clipped 5' region must be ignored when determining
+    // clipping boundaries (fixed in PR #1065). The order of the two operations should not matter.
+    val readLen = 89 // bases after removing 25 hard-clipped
+    val Seq(r1a, r2a) = new SamBuilder(readLength = readLen).addPair(
+      start1  = 100,
+      start2  = 75,
+      cigar1  = s"25H${readLen}M",
+      cigar2  = s"${readLen}M25H",
+    )
+    val Seq(r1b, r2b) = new SamBuilder(readLength = readLen).addPair(
+      start1  = 100,
+      start2  = 75,
+      cigar1  = s"25H${readLen}M",
+      cigar2  = s"${readLen}M25H",
+    )
+
+    // Order A: overlapping → past-mate
+    samClipper.clipOverlappingReads(r1a, r2a)
+    samClipper.clipExtendingPastMateEnds(r1a, r2a)
+
+    // Order B: past-mate → overlapping
+    samClipper.clipExtendingPastMateEnds(r1b, r2b)
+    samClipper.clipOverlappingReads(r1b, r2b)
+
+    assertSameClipping(r1a, r1b, "R1")
+    assertSameClipping(r2a, r2b, "R2")
+  }
+
+  it should "produce the same result for clip-overlapping-reads then clip-bases-past-mate as the reverse order for soft-clipped dovetailed reads" in {
+    // R2 has 10 soft-clipped bases at its 3' (left reference) end extending the footprint to
+    // the left. Both reads dovetail: R2's soft-clipped bases reach R1.start.
+    // cigar2 = "10S90M" has lengthOnQuery=100, matching readLength=100.
+    val Seq(r1a, r2a) = new SamBuilder(readLength = 100).addPair(
+      start1  = 100,
+      start2  = 110, // R2 aligned start is 110; the leading 10S extend to ref position 100
+      cigar1  = "100M",
+      cigar2  = "10S90M",
+    )
+    val Seq(r1b, r2b) = new SamBuilder(readLength = 100).addPair(
+      start1  = 100,
+      start2  = 110,
+      cigar1  = "100M",
+      cigar2  = "10S90M",
+    )
+
+    // Order A: overlapping → past-mate
+    samClipper.clipOverlappingReads(r1a, r2a)
+    samClipper.clipExtendingPastMateEnds(r1a, r2a)
+
+    // Order B: past-mate → overlapping
+    samClipper.clipExtendingPastMateEnds(r1b, r2b)
+    samClipper.clipOverlappingReads(r1b, r2b)
+
+    assertSameClipping(r1a, r1b, "R1")
+    assertSameClipping(r2a, r2b, "R2")
+  }
+
+  it should "produce a 1bp R1 and an unmapped R2 when both 5-prime ends are at the same reference position and R2 has leading soft-clips" in {
+    // R1 is 1bp at ref 100. R2 has 5 soft-clipped bases followed by 1 aligned base, also starting at ref 100.
+    // The midpoint falls at mate.end (ref 100), so clipOverlappingReads clips all of R2, unmapping it.
+    // clipExtendingPastMateEnds is then a no-op because the mate is already unmapped.
+    val builder = new SamBuilder(readLength = 1)
+    val Seq(r1, r2) = builder.addPair(
+      start1 = 100,
+      start2 = 100,
+      cigar1 = "1M",
+      cigar2 = "5S1M",
+      bases1 = "A",
+      quals1 = "F",
+      bases2 = "AAAAAA",
+      quals2 = "FFFFFF",
+    )
+
+    samClipper.clipOverlappingReads(r1, r2)
+    samClipper.clipExtendingPastMateEnds(r1, r2)
+
+    r1.start    shouldBe 100
+    r1.end      shouldBe 100
+    r2.unmapped shouldBe true
+  }
+
+  it should "produce a 1bp R1 and an unmapped R2 when both reads are 1bp and fully co-located" in {
+    // When two 1bp reads are at the same reference position, the midpoint equals mate.end,
+    // so clipOverlappingReads clips all of R2, resulting in R2 being unmapped.
+    val Seq(r1a, r2a) = new SamBuilder(readLength = 1).addPair(start1 = 100, start2 = 100)
+    val Seq(r1b, r2b) = new SamBuilder(readLength = 1).addPair(start1 = 100, start2 = 100)
+
+    // Order A: overlapping → past-mate
+    samClipper.clipOverlappingReads(r1a, r2a)
+    samClipper.clipExtendingPastMateEnds(r1a, r2a)
+
+    // Order B: past-mate → overlapping
+    samClipper.clipExtendingPastMateEnds(r1b, r2b)
+    samClipper.clipOverlappingReads(r1b, r2b)
+
+    r1a.start    shouldBe 100
+    r1a.end      shouldBe 100
+    r2a.unmapped shouldBe true
+
+    r1b.start    shouldBe 100
+    r1b.end      shouldBe 100
+    r2b.unmapped shouldBe true
+  }
+
+  it should "clear the proper-pair flag on the surviving read when its mate is clipped to unmapped" in {
+    val builder = new SamBuilder(readLength = 1, sort = Some(Queryname))
+    builder.addPair(start1 = 100, start2 = 100)
+
+    val out = makeTempFile("out.", ".bam")
+    new ClipBam(input = builder.toTempFile(), output = out, ref = ref, clipOverlappingReads = true).execute()
+
+    val clipped = SamSource(out).toSeq
+    clipped should have length 2
+
+    val mapped   = clipped.find(_.mapped).value
+    val unmapped = clipped.find(_.unmapped).value
+
+    unmapped.unmapped       shouldBe true
+    mapped.mateUnmapped     shouldBe true
+    mapped.properlyPaired   shouldBe false
+    unmapped.properlyPaired shouldBe false
+  }
+
   it should "unmap reads when the hard clipping length requested is greater than the length of the reads" in {
     val builder = new SamBuilder(readLength = 100, sort = Some(Queryname))
     val clipper = new ClipBam(