Make PileupBuilder.includeMapPositionsOutsideFrInsert intuitively cor…

…rect
fulcrumgenomics · Apr 24, 2024 · 7f76630 · 7f76630
1 parent 3a74fd2
commit 7f76630
Show file tree

Hide file tree

Showing 3 changed files with 50 additions and 11 deletions.
diff --git a/src/main/scala/com/fulcrumgenomics/bam/pileup/PileupBuilder.scala b/src/main/scala/com/fulcrumgenomics/bam/pileup/PileupBuilder.scala
@@ -116,14 +116,11 @@ object PileupBuilder {
     )
   }
 
-  /** Returns true if <rec> is in a mapped FR pair but the position <pos> is outside the insert coordinates of <rec>.
-    * Returns false if <rec> is in a mapped FR pair and the position <pos> is inside the insert coordinates of <rec> or
-    * <rec> is not in a mapped FR pair.
-    */
-  private def positionIsOutsideFrInsert(rec: SamRecord, refIndex: Int, pos: Int): Boolean = {
-    rec.isFrPair && {
+  /** Returns true if <rec> is in a mapped FR pair and the position <pos> is inside the insert coordinates of <rec>. */
+  private def positionIsInsideFrInsert(rec: SamRecord, refIndex: Int, pos: Int): Boolean = {
+    refIndex == rec.refIndex && rec.isFrPair && {
       val (start, end) = Bams.insertCoordinates(rec)
-      rec.refIndex == refIndex && pos >= start && pos <= end
+      pos >= start && pos <= end
     }
   }
 
@@ -200,8 +197,8 @@ trait PileupBuilder extends PileupParameters {
         if (compare) compare = this.acceptRecord(rec)
         if (compare) compare = rec.end >= pos
         if (compare) compare = rec.start <= pos || PileupBuilder.startsWithInsertion(rec)
-        if (compare) compare = if (!includeMapPositionsOutsideFrInsert && rec.isFrPair) {
-          PileupBuilder.positionIsOutsideFrInsert(rec, refIndex = refIndex, pos = pos)
+        if (compare) compare = if (rec.paired && !includeMapPositionsOutsideFrInsert) {
+          PileupBuilder.positionIsInsideFrInsert(rec, refIndex = refIndex, pos = pos)
         } else { compare }
         compare
       }

diff --git a/src/main/scala/com/fulcrumgenomics/testing/SamBuilder.scala b/src/main/scala/com/fulcrumgenomics/testing/SamBuilder.scala
@@ -27,15 +27,17 @@ package com.fulcrumgenomics.testing
 import java.nio.file.Files
 import java.text.DecimalFormat
 import java.util.concurrent.atomic.AtomicLong
-
 import com.fulcrumgenomics.FgBioDef._
 import com.fulcrumgenomics.alignment.Cigar
+import com.fulcrumgenomics.commons.CommonsDef.IteratorToJavaCollectionsAdapter
 import com.fulcrumgenomics.bam.api.SamRecord.{MissingBases, MissingQuals}
 import com.fulcrumgenomics.bam.api.{SamOrder, SamRecord, SamSource, SamWriter}
 import com.fulcrumgenomics.fasta.{SequenceDictionary, SequenceMetadata}
 import com.fulcrumgenomics.testing.SamBuilder._
+import htsjdk.samtools.SamPairUtil.PairOrientation
 import htsjdk.samtools._
 
+import java.util
 import scala.collection.mutable.ArrayBuffer
 import scala.util.Random
 
@@ -173,7 +175,8 @@ class SamBuilder(val readLength: Int=100,
     r2("RG") = this.rg.getId
     attrs.foreach { case (key, value) => r2(key) = value }
 
-    SamPairUtil.setMateInfo(r1.asSam, r2.asSam, true)
+    val orientations = util.Arrays.asList(PairOrientation.FR, PairOrientation.RF, PairOrientation.TANDEM)
+    SamPairUtil.setProperPairAndMateInfo(r1.asSam, r2.asSam, orientations, true)
     val recs = Seq(r1, r2)
     this.records ++= recs
     recs

diff --git a/src/test/scala/com/fulcrumgenomics/bam/pileup/PileupBuilderTest.scala b/src/test/scala/com/fulcrumgenomics/bam/pileup/PileupBuilderTest.scala
@@ -289,6 +289,26 @@ class PileupBuilderTest extends UnitSpec {
       piler.safelyClose()
     }
 
+    it should "not filter out single-end records when we are removing records outside the insert of FR pairs" in {
+      val builder = new SamBuilder(readLength = ReadLength, sd = Some(TestSequenceDictionary), sort = Some(Coordinate))
+
+      builder.addFrag(name = "q1", contig = Chr1Index, start = 100)
+
+      val source = builder.toSource
+      val piler  = PileupBuilder(
+        source,
+        accessPattern                      = accessPattern,
+        mappedPairsOnly                    = false,
+        includeMapPositionsOutsideFrInsert = true
+      )
+
+      piler.pileup(Chr1, 100).depth shouldBe 1
+      piler.pileup(Chr1, 100 + ReadLength - 1).depth shouldBe 1
+
+      source.safelyClose()
+      piler.safelyClose()
+    }
+
     it should "filter out records where a position is outside the insert for an FR pair" in {
       val builder = new SamBuilder(readLength = ReadLength, sd = Some(TestSequenceDictionary), sort = Some(Coordinate))
 
@@ -323,6 +343,23 @@ class PileupBuilderTest extends UnitSpec {
       piler.safelyClose()
     }
 
+    it should "exclude records that appear to be in an FR pair but are on different chromosomes" in {
+      val builder = new SamBuilder(readLength = ReadLength, sd = Some(TestSequenceDictionary), sort = Some(Coordinate))
+
+      builder.addPair(name = "q1", contig = Chr1Index, contig2 = Some(Chr2Index), start1 = 100, start2 = 125)
+
+      val source = builder.toSource
+      val piler  = PileupBuilder(source, accessPattern = accessPattern, includeMapPositionsOutsideFrInsert = false)
+
+      piler.pileup(Chr1, 100).depth shouldBe 0
+      piler.pileup(Chr1, 101).depth shouldBe 0
+      piler.pileup(Chr1, 100 + ReadLength - 1).depth shouldBe 0
+      piler.pileup(Chr1, 101 + ReadLength - 1).depth shouldBe 0
+
+      source.safelyClose()
+      piler.safelyClose()
+    }
+
     it should "not filter out records where a position is outside what might look like an 'insert' for a non-FR pair" in {
       val builder = new SamBuilder(readLength = ReadLength, sd = Some(TestSequenceDictionary), sort = Some(Coordinate))
 
@@ -340,6 +377,8 @@ class PileupBuilderTest extends UnitSpec {
       piler.safelyClose()
     }
 
+    //////
+
     it should "filter out records and entries with custom filters" in {
       val builder = new SamBuilder(readLength = ReadLength, sd = Some(TestSequenceDictionary), sort = Some(Coordinate))