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

Added a new constructor in which priors for hom-ref, het, and hom-var can be specified. Otherwise, it uses the default values of 0.999, 1e-3, and 1e-5 respectively. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@991 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kiran 2009-06-11 20:33:45 +00:00
parent 909fefa40a
commit a12009e9e7
1 changed files with 20 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java
View File

@ -37,6 +37,11 @@ public class GenotypeLikelihoods {
public double[] likelihoods; public double[] likelihoods;
public String[] genotypes; public String[] genotypes;
// The genotype priors;
private double priorHomRef;
private double priorHet;
private double priorHomVar;
// Store the 2nd-best base priors for on-genotype primary bases // Store the 2nd-best base priors for on-genotype primary bases
HashMap<String, Double> onNextBestBasePriors = new HashMap<String, Double>(); HashMap<String, Double> onNextBestBasePriors = new HashMap<String, Double>();
@ -44,6 +49,18 @@ public class GenotypeLikelihoods {
HashMap<String, Double> offNextBestBasePriors = new HashMap<String, Double>(); HashMap<String, Double> offNextBestBasePriors = new HashMap<String, Double>();
public GenotypeLikelihoods() { public GenotypeLikelihoods() {
initialize(1.0 - 1e-3, 1e-3, 1e-5);
}
public GenotypeLikelihoods(double priorHomRef, double priorHet, double priorHomVar) {
initialize(priorHomRef, priorHet, priorHomVar);
}
private void initialize(double priorHomRef, double priorHet, double priorHomVar) {
this.priorHomRef = priorHomRef;
this.priorHet = priorHet;
this.priorHomVar = priorHomVar;
likelihoods = new double[10]; likelihoods = new double[10];
genotypes = new String[10]; genotypes = new String[10];
@ -153,17 +170,17 @@ public class GenotypeLikelihoods {
if ((genotypes[i].charAt(0) == ref) && (genotypes[i].charAt(1) == ref)) if ((genotypes[i].charAt(0) == ref) && (genotypes[i].charAt(1) == ref))
{ {
// hom-ref // hom-ref
likelihoods[i] += Math.log10(1.0 - 1e-3); likelihoods[i] += Math.log10(priorHomRef);
} }
else if ((genotypes[i].charAt(0) != ref) && (genotypes[i].charAt(1) != ref)) else if ((genotypes[i].charAt(0) != ref) && (genotypes[i].charAt(1) != ref))
{ {
// hom-nonref // hom-nonref
likelihoods[i] += Math.log10(1e-5); likelihoods[i] += Math.log10(priorHomVar);
} }
else else
{ {
// het // het
likelihoods[i] += Math.log10(1e-3); likelihoods[i] += Math.log10(priorHet);
} }
if (Double.isInfinite(likelihoods[i])) { likelihoods[i] = -1000; } if (Double.isInfinite(likelihoods[i])) { likelihoods[i] = -1000; }
} }
@ -231,15 +248,7 @@ public class GenotypeLikelihoods {
double offPrior = MathUtils.binomialProbability(offIsGenotypic, offTotal, offNextBestBasePriors.get(genotypes[genotypeIndex])); double offPrior = MathUtils.binomialProbability(offIsGenotypic, offTotal, offNextBestBasePriors.get(genotypes[genotypeIndex]));
double onPrior = MathUtils.binomialProbability(onIsGenotypic, onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex])); double onPrior = MathUtils.binomialProbability(onIsGenotypic, onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex]));
//System.out.printf("%c%c [%d %d %f %e] [%d %d %f %e]\n",
// firstAllele, secondAllele,
// offIsGenotypic, offTotal, offNextBestBasePriors.get(genotypes[genotypeIndex]), offPrior,
// onIsGenotypic, onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex]), onPrior);
//System.out.printf("%f\n", (new cern.jet.random.Binomial(onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex]), cern.jet.random.engine.RandomEngine.makeDefault())).pdf(onIsGenotypic));
likelihoods[genotypeIndex] += Math.log10(offPrior) + Math.log10(onPrior); likelihoods[genotypeIndex] += Math.log10(offPrior) + Math.log10(onPrior);
//System.out.println(genotypes[genotypeIndex] + " " + offNextBestBasePriors.get(genotypes[genotypeIndex]) + " " + offIsGenotypic + " " + offTotal + " " + (((double) offIsGenotypic)/((double) offTotal)) + " " + offPrior);
} }
this.sort(); this.sort();
} }