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
This commit is contained in:
ebanks 2009-09-24 15:29:50 +00:00
parent d262cbd41c
commit 90de2e0cde
3 changed files with 36 additions and 9 deletions

View File

@ -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<Integer, Integer> getRatioCounts(char ref, ReadBackedPileup pileup, RodGeliText variant);
protected abstract boolean excludeHetsOnly();
@ -51,8 +59,13 @@ public abstract class RatioFilter implements VariantExclusionCriterion {
Pair<Integer, Integer> 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<Integer, Integer> 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;
}
}
}
}

View File

@ -17,6 +17,14 @@ public class VECAlleleBalance extends RatioFilter {
}
public void initialize(HashMap<String,String> 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 )

View File

@ -14,6 +14,14 @@ public class VECOnOffGenotypeRatio extends RatioFilter {
}
public void initialize(HashMap<String,String> 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 )