diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java index f5f707690..a7c2ff23c 100755 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java @@ -27,6 +27,7 @@ package org.broadinstitute.sting.gatk.walkers.haplotypecaller; import com.google.java.contract.Ensures; import net.sf.picard.reference.IndexedFastaSequenceFile; +import org.broadinstitute.sting.utils.activeregion.ActivityProfileResult; import org.broadinstitute.sting.utils.help.DocumentedGATKFeature; import org.broadinstitute.sting.commandline.*; import org.broadinstitute.sting.gatk.CommandLineGATK; @@ -303,8 +304,8 @@ public class HaplotypeCaller extends ActiveRegionWalker implem public boolean wantsNonPrimaryReads() { return true; } @Override - @Ensures({"result >= 0.0", "result <= 1.0"}) - public double isActive( final RefMetaDataTracker tracker, final ReferenceContext ref, final AlignmentContext context ) { + @Ensures({"result.isActiveProb >= 0.0", "result.isActiveProb <= 1.0"}) + public ActivityProfileResult isActive( final RefMetaDataTracker tracker, final ReferenceContext ref, final AlignmentContext context ) { if( UG_engine.getUAC().GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ) { for( final VariantContext vc : tracker.getValues(UG_engine.getUAC().alleles, ref.getLocus()) ) { @@ -313,15 +314,15 @@ public class HaplotypeCaller extends ActiveRegionWalker implem } } if( tracker.getValues(UG_engine.getUAC().alleles, ref.getLocus()).size() > 0 ) { - return 1.0; + return new ActivityProfileResult(1.0); } } if( USE_ALLELES_TRIGGER ) { - return ( tracker.getValues(UG_engine.getUAC().alleles, ref.getLocus()).size() > 0 ? 1.0 : 0.0 ); + return new ActivityProfileResult( tracker.getValues(UG_engine.getUAC().alleles, ref.getLocus()).size() > 0 ? 1.0 : 0.0 ); } - if( context == null ) { return 0.0; } + if( context == null ) { return new ActivityProfileResult(0.0); } final List noCall = new ArrayList(); // used to noCall all genotypes until the exact model is applied noCall.add(Allele.NO_CALL); @@ -362,7 +363,7 @@ public class HaplotypeCaller extends ActiveRegionWalker implem alleles.add( FAKE_REF_ALLELE ); alleles.add( FAKE_ALT_ALLELE ); final VariantCallContext vcOut = UG_engine_simple_genotyper.calculateGenotypes(new VariantContextBuilder("HCisActive!", context.getContig(), context.getLocation().getStart(), context.getLocation().getStop(), alleles).genotypes(genotypes).make(), GenotypeLikelihoodsCalculationModel.Model.INDEL); - return ( vcOut == null ? 0.0 : QualityUtils.qualToProb( vcOut.getPhredScaledQual() ) ); + return new ActivityProfileResult( vcOut == null ? 0.0 : QualityUtils.qualToProb( vcOut.getPhredScaledQual() ) ); } //--------------------------------------------------------------------------------------------------------------- diff --git a/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java b/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java index 6ff9f3bd5..7d035d208 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java +++ b/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java @@ -279,7 +279,6 @@ public class LocusIteratorByState extends LocusIterator { */ private void lazyLoadNextAlignmentContext() { while (nextAlignmentContext == null && readStates.hasNext()) { - // this call will set hasExtendedEvents to true if it picks up a read with indel right before the current position on the ref: readStates.collectPendingReads(); final GenomeLoc location = getLocation(); @@ -378,7 +377,7 @@ public class LocusIteratorByState extends LocusIterator { CigarOperator op = state.stepForwardOnGenome(); if (op == null) { // we discard the read only when we are past its end AND indel at the end of the read (if any) was - // already processed. Keeping the read state that retunred null upon stepForwardOnGenome() is safe + // already processed. Keeping the read state that returned null upon stepForwardOnGenome() is safe // as the next call to stepForwardOnGenome() will return null again AND will clear hadIndel() flag. it.remove(); // we've stepped off the end of the object } 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 1b9c12fb0..845c4eacf 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java +++ b/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java @@ -13,6 +13,7 @@ import org.broadinstitute.sting.gatk.walkers.Walker; import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.GenomeLocSortedSet; import org.broadinstitute.sting.utils.activeregion.ActivityProfile; +import org.broadinstitute.sting.utils.activeregion.ActivityProfileResult; import org.broadinstitute.sting.utils.pileup.PileupElement; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; @@ -69,8 +70,7 @@ public class TraverseActiveRegions extends TraversalEngine extends TraversalEngine extends TraversalEngine walker, + private final ActivityProfileResult walkerActiveProb(final ActiveRegionWalker walker, final RefMetaDataTracker tracker, final ReferenceContext refContext, final AlignmentContext locus, final GenomeLoc location) { if ( walker.hasPresetActiveRegions() ) { - return walker.presetActiveRegions.overlaps(location) ? 1.0 : 0.0; + return new ActivityProfileResult(walker.presetActiveRegions.overlaps(location) ? 1.0 : 0.0); } else { return walker.isActive( tracker, refContext, locus ); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java index e38e166ea..b2975cbbf 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java @@ -12,6 +12,7 @@ import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.GenomeLocParser; import org.broadinstitute.sting.utils.GenomeLocSortedSet; +import org.broadinstitute.sting.utils.activeregion.ActivityProfileResult; import org.broadinstitute.sting.utils.interval.IntervalMergingRule; import org.broadinstitute.sting.utils.interval.IntervalSetRule; import org.broadinstitute.sting.utils.interval.IntervalUtils; @@ -72,7 +73,7 @@ public abstract class ActiveRegionWalker extends Walker isActiveList; + final List isActiveList; private GenomeLoc lastLoc = null; private static final int FILTER_SIZE = 65; - private static final Double[] GaussianKernel; + private static final double[] GaussianKernel; static { - GaussianKernel = new Double[2*FILTER_SIZE + 1]; + GaussianKernel = new double[2*FILTER_SIZE + 1]; for( int iii = 0; iii < 2*FILTER_SIZE + 1; iii++ ) { GaussianKernel[iii] = MathUtils.NormalDistribution(FILTER_SIZE, 40.0, iii); } @@ -63,22 +62,22 @@ public class ActivityProfile { // todo -- add unit tests // TODO -- own preset regions public ActivityProfile(final GenomeLocParser parser, final boolean presetRegions) { - this(parser, presetRegions, new ArrayList(), null); + this(parser, presetRegions, new ArrayList(), null); } - protected ActivityProfile(final GenomeLocParser parser, final boolean presetRegions, final List isActiveList, final GenomeLoc regionStartLoc) { + protected ActivityProfile(final GenomeLocParser parser, final boolean presetRegions, final List isActiveList, final GenomeLoc regionStartLoc) { this.parser = parser; this.presetRegions = presetRegions; this.isActiveList = isActiveList; this.regionStartLoc = regionStartLoc; } - public void add(final GenomeLoc loc, final double score) { + public void add(final GenomeLoc loc, final ActivityProfileResult result) { if ( loc.size() != 1 ) throw new ReviewedStingException("Bad add call to ActivityProfile: loc " + loc + " size != 1" ); if ( lastLoc != null && loc.getStart() != lastLoc.getStop() + 1 ) throw new ReviewedStingException("Bad add call to ActivityProfile: lastLoc added " + lastLoc + " and next is " + loc); - isActiveList.add(score); + isActiveList.add(result); if( regionStartLoc == null ) { regionStartLoc = loc; } @@ -93,22 +92,33 @@ public class ActivityProfile { * @return a new ActivityProfile that's the band-pass filtered version of this profile */ public ActivityProfile bandPassFilter() { - final Double[] activeProbArray = isActiveList.toArray(new Double[isActiveList.size()]); - final Double[] filteredProbArray = new Double[activeProbArray.length]; + final double[] activeProbArray = new double[isActiveList.size()]; + int iii = 0; + for( final ActivityProfileResult result : isActiveList ) { + activeProbArray[iii++] = result.isActiveProb; + } + final double[] filteredProbArray = new double[activeProbArray.length]; if( !presetRegions ) { - for( int iii = 0; iii < activeProbArray.length; iii++ ) { - final Double[] kernel = (Double[]) ArrayUtils.subarray(GaussianKernel, Math.max(FILTER_SIZE-iii, 0), Math.min(GaussianKernel.length,FILTER_SIZE + activeProbArray.length - iii)); - final Double[] activeProbSubArray = (Double[]) ArrayUtils.subarray(activeProbArray, Math.max(0,iii - FILTER_SIZE), Math.min(activeProbArray.length,iii + FILTER_SIZE + 1)); + for( iii = 0; iii < activeProbArray.length; iii++ ) { + final double[] kernel = ArrayUtils.subarray(GaussianKernel, Math.max(FILTER_SIZE-iii, 0), Math.min(GaussianKernel.length,FILTER_SIZE + activeProbArray.length - iii)); + final double[] activeProbSubArray = ArrayUtils.subarray(activeProbArray, Math.max(0,iii - FILTER_SIZE), Math.min(activeProbArray.length,iii + FILTER_SIZE + 1)); filteredProbArray[iii] = MathUtils.dotProduct(activeProbSubArray, kernel); } } - return new ActivityProfile(parser, presetRegions, Arrays.asList(filteredProbArray), regionStartLoc); + iii = 0; + for( final double prob : filteredProbArray ) { + final ActivityProfileResult result = isActiveList.get(iii++); + result.isActiveProb = prob; + result.resultState = ActivityProfileResult.ActivityProfileResultState.NONE; + result.resultValue = null; + } + return new ActivityProfile(parser, presetRegions, isActiveList, regionStartLoc); } /** * Partition this profile into active regions - * @param activeRegionExtension - * @return + * @param activeRegionExtension the amount of margin overlap in the active region + * @return the list of active regions */ public List createActiveRegions( final int activeRegionExtension, final int maxRegionSize ) { final double ACTIVE_PROB_THRESHOLD = 0.002; // TODO: needs to be set-able by the walker author @@ -119,14 +129,14 @@ public class ActivityProfile { return Collections.emptyList(); } else if( isActiveList.size() == 1 ) { // there's a single element, it's either active or inactive - boolean isActive = isActiveList.get(0) > ACTIVE_PROB_THRESHOLD; + boolean isActive = isActiveList.get(0).isActiveProb > ACTIVE_PROB_THRESHOLD; returnList.addAll(createActiveRegion(isActive, 0, 0, activeRegionExtension, maxRegionSize)); } else { // there are 2+ elements, divide these up into regions - boolean isActive = isActiveList.get(0) > ACTIVE_PROB_THRESHOLD; + boolean isActive = isActiveList.get(0).isActiveProb > ACTIVE_PROB_THRESHOLD; int curStart = 0; for(int iii = 1; iii < isActiveList.size(); iii++ ) { - final boolean thisStatus = isActiveList.get(iii) > ACTIVE_PROB_THRESHOLD; + final boolean thisStatus = isActiveList.get(iii).isActiveProb > ACTIVE_PROB_THRESHOLD; if( isActive != thisStatus ) { returnList.addAll(createActiveRegion(isActive, curStart, iii - 1, activeRegionExtension, maxRegionSize)); isActive = thisStatus; @@ -143,7 +153,7 @@ public class ActivityProfile { * @param isActive should the region be active? * @param curStart offset (0-based) from the start of this region * @param curEnd offset (0-based) from the start of this region - * @param activeRegionExtension + * @param activeRegionExtension the amount of margin overlap in the active region * @return a fully initialized ActiveRegion with the above properties */ private final List createActiveRegion(final boolean isActive, final int curStart, final int curEnd, final int activeRegionExtension, final int maxRegionSize) { @@ -161,7 +171,7 @@ public class ActivityProfile { final int size = curEnd - curStart + 1; for( int iii = curStart + (int)(size*0.25); iii < curEnd - (int)(size*0.25); iii++ ) { - if( isActiveList.get(iii) < minProb ) { minProb = isActiveList.get(iii); cutPoint = iii; } + if( isActiveList.get(iii).isActiveProb < minProb ) { minProb = isActiveList.get(iii).isActiveProb; cutPoint = iii; } } final List leftList = createActiveRegion(isActive, curStart, cutPoint, activeRegionExtension, maxRegionSize, new ArrayList()); final List rightList = createActiveRegion(isActive, cutPoint+1, curEnd, activeRegionExtension, maxRegionSize, new ArrayList()); diff --git a/public/java/src/org/broadinstitute/sting/utils/activeregion/ActivityProfileResult.java b/public/java/src/org/broadinstitute/sting/utils/activeregion/ActivityProfileResult.java new file mode 100644 index 000000000..8dc29aa3c --- /dev/null +++ b/public/java/src/org/broadinstitute/sting/utils/activeregion/ActivityProfileResult.java @@ -0,0 +1,31 @@ +package org.broadinstitute.sting.utils.activeregion; + +/** + * Created with IntelliJ IDEA. + * User: rpoplin + * Date: 7/27/12 + */ + +public class ActivityProfileResult { + public double isActiveProb; + public ActivityProfileResultState resultState; + public Number resultValue; + + public enum ActivityProfileResultState { + NONE, + HIGH_QUALITY_SOFT_CLIPS + } + + public ActivityProfileResult( final double isActiveProb ) { + this.isActiveProb = isActiveProb; + this.resultState = ActivityProfileResultState.NONE; + this.resultValue = null; + } + + public ActivityProfileResult( final double isActiveProb, final ActivityProfileResultState resultState, final Number resultValue ) { + this.isActiveProb = isActiveProb; + this.resultState = resultState; + this.resultValue = resultValue; + } + +} diff --git a/public/java/test/org/broadinstitute/sting/utils/activeregion/ActivityProfileUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/activeregion/ActivityProfileUnitTest.java index 282f19d8a..f7c564c74 100644 --- a/public/java/test/org/broadinstitute/sting/utils/activeregion/ActivityProfileUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/activeregion/ActivityProfileUnitTest.java @@ -123,12 +123,12 @@ public class ActivityProfileUnitTest extends BaseTest { for ( int i = 0; i < cfg.probs.size(); i++ ) { double p = cfg.probs.get(i); GenomeLoc loc = genomeLocParser.createGenomeLoc(cfg.regionStart.getContig(), cfg.regionStart.getStart() + i, cfg.regionStart.getStart() + i); - profile.add(loc, p); + profile.add(loc, new ActivityProfileResult(p)); } Assert.assertEquals(profile.regionStartLoc, genomeLocParser.createGenomeLoc(cfg.regionStart.getContig(), cfg.regionStart.getStart(), cfg.regionStart.getStart() )); Assert.assertEquals(profile.size(), cfg.probs.size()); - Assert.assertEquals(profile.isActiveList, cfg.probs); + assertProbsAreEqual(profile.isActiveList, cfg.probs); assertRegionsAreEqual(profile.createActiveRegions(0, 100), cfg.expectedRegions); } @@ -140,5 +140,12 @@ public class ActivityProfileUnitTest extends BaseTest { } } + private void assertProbsAreEqual(List actual, List expected) { + Assert.assertEquals(actual.size(), expected.size()); + for ( int i = 0; i < actual.size(); i++ ) { + Assert.assertEquals(actual.get(i).isActiveProb, expected.get(i)); + } + } + // todo -- test extensions } \ No newline at end of file