diff --git a/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java b/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java index ac1c751ca..21971e189 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java +++ b/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java @@ -102,11 +102,14 @@ public class TraverseActiveRegions extends TraversalEngine extends Walker> activeRegionBindings = null; + @Advanced + @Argument(fullName="activeRegionExtension", shortName="activeRegionExtension", doc="The active region extension; if not provided defaults to Walker annotated default", required = false) + public Integer activeRegionExtension = null; + + @Advanced + @Argument(fullName="activeRegionMaxSize", shortName="activeRegionMaxSize", doc="The active region maximum size; if not provided defaults to Walker annotated default", required = false) + public Integer activeRegionMaxSize = null; + + @Advanced + @Argument(fullName="bandPassFilterSize", shortName="bandPassFilterSize", doc="The filter size of band pass filter; if not provided defaults to Walker annotated default", required = false) + public Integer bandPassFilterSize = null; + + @Advanced + @Argument(fullName="bandPassSigma", shortName="bandPassSigma", doc="The sigma of the band pass filter Gaussian kernel; if not provided defaults to Walker annotated default", required = false) + public Double bandPassSigma = null; + private GenomeLocSortedSet presetActiveRegions = null; @Override diff --git a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java index f1f0ab9b1..0c3ed87c0 100644 --- a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java @@ -651,7 +651,7 @@ public class MathUtils { final double sum = sum(array); final double[] normalized = new double[array.length]; - if ( sum < 0.0 || sum > 1.0 ) throw new IllegalArgumentException("Values in probability array sum to a negative number " + sum); + if ( sum < 0.0 ) throw new IllegalArgumentException("Values in probability array sum to a negative number " + sum); for ( int i = 0; i < array.length; i++ ) { normalized[i] = array[i] / sum; } diff --git a/public/java/src/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfile.java b/public/java/src/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfile.java index 1a8bac086..5c6389c26 100644 --- a/public/java/src/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfile.java +++ b/public/java/src/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfile.java @@ -44,8 +44,10 @@ import java.util.LinkedList; */ public class BandPassActivityProfile extends ActivityProfile { public static final int DEFAULT_FILTER_SIZE = 80; + public static final double DEFAULT_SIGMA = 55.0; private final int filterSize; + private final double sigma; private final double[] GaussianKernel; /** @@ -53,7 +55,7 @@ public class BandPassActivityProfile extends ActivityProfile { * @param parser our genome loc parser */ public BandPassActivityProfile(final GenomeLocParser parser) { - this(parser, DEFAULT_FILTER_SIZE); + this(parser, DEFAULT_FILTER_SIZE, DEFAULT_SIGMA); } /** @@ -63,16 +65,18 @@ public class BandPassActivityProfile extends ActivityProfile { * side that are included in the band. So a filter size of 1 implies that the actual band * is 3 bp, 1 for the center site and 1 on each size. 2 => 5, etc. */ - public BandPassActivityProfile(final GenomeLocParser parser, final int filterSize) { + public BandPassActivityProfile(final GenomeLocParser parser, final int filterSize, final double sigma) { super(parser); if ( filterSize < 0 ) throw new IllegalArgumentException("Filter size must be greater than or equal to 0 but got " + filterSize); + if ( sigma < 0 ) throw new IllegalArgumentException("Sigma must be greater than or equal to 0 but got " + sigma); // setup the Gaussian kernel for the band pass filter this.filterSize = filterSize; + this.sigma = sigma; final double[] kernel = new double[getBandSize()]; for( int iii = 0; iii < 2* filterSize + 1; iii++ ) { - kernel[iii] = MathUtils.NormalDistribution(filterSize, 55.0, iii); + kernel[iii] = MathUtils.NormalDistribution(filterSize, sigma, iii); } this.GaussianKernel = MathUtils.normalizeFromRealSpace(kernel); } @@ -108,6 +112,15 @@ public class BandPassActivityProfile extends ActivityProfile { return filterSize; } + /** + * Get the Gaussian kernel sigma value + * @return a positive double + */ + @Ensures("result >= 0") + public double getSigma() { + return sigma; + } + /** * Get the kernel of this band pass filter. Do not modify returned result * @return the kernel used in this band pass filter diff --git a/public/java/test/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfileUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfileUnitTest.java index a2a85f1d0..0a71bad14 100644 --- a/public/java/test/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfileUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/activeregion/BandPassActivityProfileUnitTest.java @@ -64,11 +64,13 @@ public class BandPassActivityProfileUnitTest extends BaseTest { for ( boolean precedingIsActive : Arrays.asList(true, false) ) { for ( int precedingSites: Arrays.asList(0, 1, 10, 100) ) { for ( int bandPassSize : Arrays.asList(0, 1, 10, 100) ) { + for ( double sigma : Arrays.asList(1.0, 2.0, BandPassActivityProfile.DEFAULT_SIGMA) ) { // for ( int start : Arrays.asList(10) ) { // for ( boolean precedingIsActive : Arrays.asList(false) ) { // for ( int precedingSites: Arrays.asList(0) ) { // for ( int bandPassSize : Arrays.asList(1) ) { - tests.add(new Object[]{ start, precedingIsActive, precedingSites, bandPassSize }); + tests.add(new Object[]{ start, precedingIsActive, precedingSites, bandPassSize, sigma }); + } } } } @@ -78,10 +80,12 @@ public class BandPassActivityProfileUnitTest extends BaseTest { } @Test(dataProvider = "BandPassBasicTest") - public void testBandPass(final int start, final boolean precedingIsActive, final int nPrecedingSites, final int bandPassSize) { - final BandPassActivityProfile profile = new BandPassActivityProfile(genomeLocParser, bandPassSize); + public void testBandPass(final int start, final boolean precedingIsActive, final int nPrecedingSites, final int bandPassSize, final double sigma) { + final BandPassActivityProfile profile = new BandPassActivityProfile(genomeLocParser, bandPassSize, sigma); final int expectedBandSize = bandPassSize * 2 + 1; + Assert.assertEquals(profile.getFilteredSize(), bandPassSize, "Wrong filter size"); + Assert.assertEquals(profile.getSigma(), sigma, "Wrong sigma"); Assert.assertEquals(profile.getBandSize(), expectedBandSize, "Wrong expected band size"); final String contig = genomeLocParser.getContigs().getSequences().get(0).getSequenceName(); @@ -132,7 +136,7 @@ public class BandPassActivityProfileUnitTest extends BaseTest { @Test( dataProvider = "BandPassComposition") public void testBandPassComposition(final int bandPassSize, final int integrationLength) { final int start = 1; - final BandPassActivityProfile profile = new BandPassActivityProfile(genomeLocParser, bandPassSize); + final BandPassActivityProfile profile = new BandPassActivityProfile(genomeLocParser, bandPassSize, BandPassActivityProfile.DEFAULT_SIGMA); final double[] rawActiveProbs = new double[integrationLength + bandPassSize * 2]; // add a buffer so that we can get all of the band pass values