From 90de2e0cdedd9229db74bdb1d1d61d0af81d2f0b Mon Sep 17 00:00:00 2001 From: ebanks Date: Thu, 24 Sep 2009 15:29:50 +0000 Subject: [PATCH] Added ability to specify whether you want to use a point estimate or fair coin test calculation; for now you can use either but fair coin test is still experimental as it needs to be parametrized correctly. This job will hopefully be done by the future Bioinformatic Analyst... git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1716 348d0f76-0448-11de-a6fe-93d51630548a --- .../gatk/walkers/filters/RatioFilter.java | 29 +++++++++++++------ .../walkers/filters/VECAlleleBalance.java | 8 +++++ .../filters/VECOnOffGenotypeRatio.java | 8 +++++ 3 files changed, 36 insertions(+), 9 deletions(-) diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java index 7aa8a84b4..3f923401a 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java @@ -1,8 +1,7 @@ package org.broadinstitute.sting.gatk.walkers.filters; import org.broadinstitute.sting.gatk.contexts.VariantContext; -import org.broadinstitute.sting.gatk.refdata.RodGeliText; -import org.broadinstitute.sting.gatk.refdata.AllelicVariant; +import org.broadinstitute.sting.gatk.refdata.*; import org.broadinstitute.sting.utils.*; import org.apache.log4j.Logger; @@ -21,6 +20,11 @@ public abstract class RatioFilter implements VariantExclusionCriterion { protected double lowThreshold = -1; protected double highThreshold = -1; + protected enum AnalysisType { POINT_ESTIMATE, FAIR_COIN_TEST }; + protected AnalysisType analysis = AnalysisType.POINT_ESTIMATE; + protected double integralPValueThreshold = 0.05; + private final static double SEARCH_INCREMENT = 0.01; + protected boolean exclude = false; public RatioFilter(final String name, Class myClass, Tail tail ) { @@ -37,6 +41,10 @@ public abstract class RatioFilter implements VariantExclusionCriterion { highThreshold = threshold; } + protected void setIntegralPvalue(double pvalue) { + integralPValueThreshold = pvalue; + } + protected abstract Pair getRatioCounts(char ref, ReadBackedPileup pileup, RodGeliText variant); protected abstract boolean excludeHetsOnly(); @@ -51,8 +59,13 @@ public abstract class RatioFilter implements VariantExclusionCriterion { Pair counts = getRatioCounts(ref, pileup, variant); boolean highGenotypeConfidence = variant.getConsensusConfidence() > minGenotypeConfidenceToTest; - boolean excludable = !excludeHetsOnly() || GenotypeUtils.isHet((AllelicVariant)variant); - exclude = excludable && highGenotypeConfidence && pointEstimateExclude(counts); + boolean excludable = !excludeHetsOnly() || GenotypeUtils.isHet(variant); + + if ( analysis == AnalysisType.POINT_ESTIMATE ) + exclude = excludable && highGenotypeConfidence && pointEstimateExclude(counts); + else if ( analysis == AnalysisType.FAIR_COIN_TEST ) + exclude = excludable && highGenotypeConfidence && integralExclude(counts); + // // for printing only // @@ -63,10 +76,7 @@ public abstract class RatioFilter implements VariantExclusionCriterion { value, exclude, variant.getLocation(), pileup.getBases())); } - private final static double SEARCH_INCREMENT = 0.01; - private final static double integralPValueThreshold = 0.05; - - // TODO - this whole calculation needs to be redone correctly + // TODO - this calculation needs to be parameterized correctly private boolean pointEstimateExclude(Pair counts) { int n = counts.first + counts.second; if ( n < minDepthOfCoverage ) @@ -90,6 +100,7 @@ public abstract class RatioFilter implements VariantExclusionCriterion { } double percentExcluded = sumExclude / sumP; + return 1 - percentExcluded <= integralPValueThreshold ; } @@ -105,4 +116,4 @@ public abstract class RatioFilter implements VariantExclusionCriterion { return true; } } -} \ No newline at end of file +} diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java index 40a3a3230..fd1889713 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java @@ -17,6 +17,14 @@ public class VECAlleleBalance extends RatioFilter { } public void initialize(HashMap args) { + if ( args.get("analysis") != null ) { + if ( args.get("analysis").equalsIgnoreCase("fair_coin_test") ) + analysis = AnalysisType.FAIR_COIN_TEST; + else if ( args.get("analysis").equalsIgnoreCase("point_estimate") ) + analysis = AnalysisType.POINT_ESTIMATE; + } + if ( args.get("pvalue") != null ) + setIntegralPvalue(Double.valueOf(args.get("pvalue"))); if ( args.get("low") != null ) lowThreshold = Double.valueOf(args.get("low")); if ( args.get("high") != null ) diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java index 4b8b1669d..ec314be0e 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java @@ -14,6 +14,14 @@ public class VECOnOffGenotypeRatio extends RatioFilter { } public void initialize(HashMap args) { + if ( args.get("analysis") != null ) { + if ( args.get("analysis").equalsIgnoreCase("fair_coin_test") ) + analysis = AnalysisType.FAIR_COIN_TEST; + else if ( args.get("analysis").equalsIgnoreCase("point_estimate") ) + analysis = AnalysisType.POINT_ESTIMATE; + } + if ( args.get("pvalue") != null ) + setIntegralPvalue(Double.valueOf(args.get("pvalue"))); if ( args.get("threshold") != null ) threshold = Double.valueOf(args.get("threshold")); if ( args.get("confidence") != null )