diff --git a/src/main/scala/com/fulcrumgenomics/bam/pileup/PileupBuilder.scala b/src/main/scala/com/fulcrumgenomics/bam/pileup/PileupBuilder.scala index 1f0bba4b8..0db04bb5e 100644 --- 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 is in a mapped FR pair but the position is outside the insert coordinates of . - * Returns false if is in a mapped FR pair and the position is inside the insert coordinates of or - * is not in a mapped FR pair. - */ - private def positionIsOutsideFrInsert(rec: SamRecord, refIndex: Int, pos: Int): Boolean = { - rec.isFrPair && { + /** Returns true if is in a mapped FR pair and the position is inside the insert coordinates of . */ + 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 index 64508991e..f6afda284 100644 --- a/src/main/scala/com/fulcrumgenomics/testing/SamBuilder.scala +++ b/src/main/scala/com/fulcrumgenomics/testing/SamBuilder.scala @@ -24,18 +24,19 @@ 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.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.nio.file.Files +import java.text.DecimalFormat +import java.util +import java.util.concurrent.atomic.AtomicLong import scala.collection.mutable.ArrayBuffer import scala.util.Random @@ -173,7 +174,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.values(): _*) + 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 index 13ed46b2d..3dc82cd6a 100644 --- 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,24 @@ 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, 100 + ReadLength - 1).depth shouldBe 0 + + piler.pileup(Chr2, 125).depth shouldBe 0 + piler.pileup(Chr2, 125 + 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))