zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix to priors to allow lower het values for mouse guys; no intergration test changes 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@2734 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
depristo 2010-01-29 15:36:12 +00:00
parent d57a86ad41
commit 9decd20f46
3 changed files with 15 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypePriors.java
View File

@ -17,9 +17,11 @@ public class DiploidGenotypePriors {
public static final double YRI_HETEROZYGOSITY = 1.0 / 850; public static final double YRI_HETEROZYGOSITY = 1.0 / 850;
/** /**
* Default value of the prob of seeing a true B/C het when the reference is A * Default value of the prob of seeing a reference error. Used to calculation the
* chance of seeing a true B/C het when the reference is A, which we assume is the product
* of the ref error rate and the het. Default value is Q60
*/ */
public static final double PROB_OF_TRISTATE_GENOTYPE = 1e-5; public static final double PROB_OF_REFERENCE_ERROR = 1e-6; // the reference is
private final static double[] flatPriors = new double[DiploidGenotype.values().length]; private final static double[] flatPriors = new double[DiploidGenotype.values().length];
@ -52,7 +54,7 @@ public class DiploidGenotypePriors {
double[] vals = heterozygosity2DiploidProbabilities(heterozygosity); double[] vals = heterozygosity2DiploidProbabilities(heterozygosity);
priors = getReferenceIndependentPriors(ref, vals[0], vals[1], vals[2]); priors = getReferenceIndependentPriors(ref, vals[0], vals[1], vals[2]);
} else { } else {
priors = getReferencePolarizedPriors(ref, heterozygosity, PROB_OF_TRISTATE_GENOTYPE); priors = getReferencePolarizedPriors(ref, heterozygosity, PROB_OF_REFERENCE_ERROR);
} }
} }
@ -198,16 +200,17 @@ public class DiploidGenotypePriors {
* *
* @param ref * @param ref
* @param heterozyosity * @param heterozyosity
* @param pTriStateGenotype * @param pRefError
*/ */
public static double[] getReferencePolarizedPriors(char ref, double heterozyosity, double pTriStateGenotype ) { public static double[] getReferencePolarizedPriors(char ref, double heterozyosity, double pRefError ) {
if ( ! MathUtils.isBounded(pTriStateGenotype, 0.0, 1.0) ) { if ( ! MathUtils.isBounded(pRefError, 0.0, 0.01) ) {
throw new RuntimeException(String.format("p Tristate genotype is out of bounds %f", pTriStateGenotype)); throw new RuntimeException(String.format("BUG: p Reference error is out of bounds (0.0 - 0.01) is allow range %f", pRefError));
} }
if ( pTriStateGenotype >= heterozyosity ) { double pTriStateGenotype = heterozyosity * PROB_OF_REFERENCE_ERROR;
throw new RuntimeException(String.format("p Tristate genotype %f is greater than the heterozygosity %f", pTriStateGenotype, heterozyosity)); // if ( pTriStateGenotype >= heterozyosity ) {
} // throw new RuntimeException(String.format("p Tristate genotype %f is greater than the heterozygosity %f", pTriStateGenotype, heterozyosity));
// }
double pHomRef = heterozygosity2HomRefProbability(heterozyosity); double pHomRef = heterozygosity2HomRefProbability(heterozyosity);
double pHet = heterozygosity2HetProbability(heterozyosity); double pHet = heterozygosity2HetProbability(heterozyosity);

2
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
View File

@ -297,7 +297,7 @@ public class UnifiedGenotyper extends LocusWalker<VariantCallContext, UnifiedGen
if ( stratifiedContexts == null ) if ( stratifiedContexts == null )
return null; return null;
DiploidGenotypePriors priors = new DiploidGenotypePriors(ref, UAC.heterozygosity, DiploidGenotypePriors.PROB_OF_TRISTATE_GENOTYPE); DiploidGenotypePriors priors = new DiploidGenotypePriors(ref, UAC.heterozygosity, DiploidGenotypePriors.PROB_OF_REFERENCE_ERROR);
VariantCallContext call = gcm.get().callLocus(tracker, ref, rawContext.getLocation(), stratifiedContexts, priors); VariantCallContext call = gcm.get().callLocus(tracker, ref, rawContext.getLocation(), stratifiedContexts, priors);
// annotate the call, if possible // annotate the call, if possible

2
java/src/org/broadinstitute/sting/playground/gatk/walkers/papergenotyper/GATKPaperGenotyper.java
View File

@ -45,7 +45,7 @@ public class GATKPaperGenotyper extends LocusWalker<SimpleCall, Integer> impleme
ReadBackedPileup pileup = context.getPileup(); ReadBackedPileup pileup = context.getPileup();
double likelihoods[] = DiploidGenotypePriors.getReferencePolarizedPriors(ref.getBase(), double likelihoods[] = DiploidGenotypePriors.getReferencePolarizedPriors(ref.getBase(),
DiploidGenotypePriors.HUMAN_HETEROZYGOSITY, DiploidGenotypePriors.HUMAN_HETEROZYGOSITY,
DiploidGenotypePriors.PROB_OF_TRISTATE_GENOTYPE); 0.01);
// get the bases and qualities from the pileup // get the bases and qualities from the pileup
byte bases[] = pileup.getBases(); byte bases[] = pileup.getBases();
byte quals[] = pileup.getQuals(); byte quals[] = pileup.getQuals();