From 8bc925a2164dca6a81c2af3644c7077d7b10bfb7 Mon Sep 17 00:00:00 2001 From: kiran Date: Fri, 31 Jul 2009 21:18:30 +0000 Subject: [PATCH] Commit on the behalf of Mark: cleaning up some old and busted code in GenotypeLikelihood and associated objects. git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1361 348d0f76-0448-11de-a6fe-93d51630548a --- .../gatk/walkers/HLAcaller/CallHLAWalker.java | 8 +- .../gatk/walkers/MultiSampleCaller.java | 6 +- .../gatk/walkers/SingleSampleGenotyper.java | 37 ++--- .../playground/utils/GenotypeLikelihoods.java | 133 ++++-------------- 4 files changed, 49 insertions(+), 135 deletions(-) diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/HLAcaller/CallHLAWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/HLAcaller/CallHLAWalker.java index eb5c2bb0a..0e02475de 100644 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/HLAcaller/CallHLAWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/HLAcaller/CallHLAWalker.java @@ -35,12 +35,6 @@ public class CallHLAWalker extends LocusWalker>{ int[] HLAstoppos; int numHLAlleles = 0; double[][] Prob; String[][] Alleles; - boolean THREE_BASE_ERRORS = false; - String PRIORS_2ND_ON = "0.000,0.302,0.366,0.142,0.000,0.548,0.370,0.000,0.319,0.000"; - String PRIORS_2ND_OFF = "0.480,0.769,0.744,0.538,0.575,0.727,0.768,0.589,0.762,0.505"; - double[] p2ndon = priorsArray(PRIORS_2ND_ON); - double[] p2ndoff = priorsArray(PRIORS_2ND_OFF); - boolean keepQ0Bases = false; //HLAreads.add("Italian Riviera"); @@ -108,7 +102,7 @@ public class CallHLAWalker extends LocusWalker>{ } out.printf("%sAs\t%sCs\t%sTs\t%sGs\t",numAs,numCs,numTs,numGs); - GenotypeLikelihoods G = new GenotypeLikelihoods(THREE_BASE_ERRORS,0.999,0.000333,0.000667, p2ndon, p2ndoff, keepQ0Bases); + GenotypeLikelihoods G = new GenotypeLikelihoods(GenotypeLikelihoods.HUMAN_HETEROZYGOSITY); SSGGenotypeCall geno = (SSGGenotypeCall)G.callGenotypes(tracker, ref, pileup); diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java index 4e63cd29b..415c3c344 100644 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java @@ -150,7 +150,7 @@ public class MultiSampleCaller extends LocusWalker offNextBestBasePriors = new HashMap(); - public GenotypeLikelihoods() { - double[] p2ndon = {0.000, 0.302, 0.366, 0.142, 0.000, 0.548, 0.370, 0.000, 0.319, 0.000}; - double[] p2ndoff = {0.480, 0.769, 0.744, 0.538, 0.575, 0.727, 0.768, 0.589, 0.762, 0.505}; - keepQ0Bases = true; - initialize(false, 1.0 - 1e-3, 1e-3, 1e-5, p2ndon, p2ndoff); + public GenotypeLikelihoods(double heterozygosity) { + double[] vals = computePriors(heterozygosity); + initialize(vals[0], vals[1], vals[2]); } - public GenotypeLikelihoods(boolean threeBaseErrors , double priorHomRef, double priorHet, double priorHomVar) { - double[] p2ndon = {0.000, 0.302, 0.366, 0.142, 0.000, 0.548, 0.370, 0.000, 0.319, 0.000}; - double[] p2ndoff = {0.480, 0.769, 0.744, 0.538, 0.575, 0.727, 0.768, 0.589, 0.762, 0.505}; - keepQ0Bases = true; - initialize(threeBaseErrors, priorHomRef, priorHet, priorHomVar, p2ndon, p2ndoff); - } - - public GenotypeLikelihoods(boolean threeBaseErrors , double priorHomRef, double priorHet, double priorHomVar, double[] p2ndon, double[] p2ndoff, boolean keepQ0Bases) { - this.keepQ0Bases = keepQ0Bases; - initialize(threeBaseErrors, priorHomRef, priorHet, priorHomVar, p2ndon, p2ndoff); + public GenotypeLikelihoods(double priorHomRef, double priorHet, double priorHomVar) { + initialize(priorHomRef, priorHet, priorHomVar); } /** @@ -131,9 +140,8 @@ public class GenotypeLikelihoods implements GenotypeGenerator { this.filterQ0Bases = filterQ0Bases; } - private void initialize(boolean threeBaseErrors , double priorHomRef, double priorHet, double priorHomVar, double[] p2ndon, double[] p2ndoff) { - this.threeBaseErrors = threeBaseErrors ; - this.oneHalfMinusData = threeBaseErrors ? oneHalfMinusData3Base : oneHalfMinusDataArachne; + private void initialize(double priorHomRef, double priorHet, double priorHomVar) { + this.oneHalfMinusData = threeStateErrors ? oneHalfMinusData3Base : oneHalfMinusDataArachne; this.priorHomRef = priorHomRef; this.priorHet = priorHet; @@ -144,11 +152,6 @@ public class GenotypeLikelihoods implements GenotypeGenerator { coverage = 0; for (int i = 0; i < likelihoods.length; i++) { likelihoods[i] = Math.log10(0.1); } - - for (int genotypeIndex = 0; genotypeIndex < 10; genotypeIndex++) { - onNextBestBasePriors.put(genotypes[genotypeIndex], p2ndon[genotypeIndex]); - offNextBestBasePriors.put(genotypes[genotypeIndex], p2ndoff[genotypeIndex]); - } } public double getHomRefPrior() { @@ -175,38 +178,6 @@ public class GenotypeLikelihoods implements GenotypeGenerator { this.priorHomVar = priorHomVar; } - public double[] getOnGenotypeSecondaryPriors() { - double[] p2ndon = new double[10]; - - for (int genotypeIndex = 0; genotypeIndex < 10; genotypeIndex++) { - p2ndon[genotypeIndex] = onNextBestBasePriors.get(genotypes[genotypeIndex]); - } - - return p2ndon; - } - - public void setOnGenotypeSecondaryPriors(double[] p2ndon) { - for (int genotypeIndex = 0; genotypeIndex < 10; genotypeIndex++) { - onNextBestBasePriors.put(genotypes[genotypeIndex], p2ndon[genotypeIndex]); - } - } - - public double[] getOffGenotypeSecondaryPriors() { - double[] p2ndoff = new double[10]; - - for (int genotypeIndex = 0; genotypeIndex < 10; genotypeIndex++) { - p2ndoff[genotypeIndex] = offNextBestBasePriors.get(genotypes[genotypeIndex]); - } - - return p2ndoff; - } - - public void setOffGenotypeSecondaryPriors(double[] p2ndoff) { - for (int genotypeIndex = 0; genotypeIndex < 10; genotypeIndex++) { - offNextBestBasePriors.put(genotypes[genotypeIndex], p2ndoff[genotypeIndex]); - } - } - public int add(char ref, char read, byte qual) { if (qual <= 0) { @@ -253,7 +224,7 @@ public class GenotypeLikelihoods implements GenotypeGenerator { } else if ( (h1 != h2) && ((h1 == read) || (h2 == read)) ) { // het p_base = getOneHalfMinusQual(qual); - } else if ( this.threeBaseErrors ) { + } else if ( threeStateErrors ) { // error //System.out.printf("%s %b %f %f%n", genotype, h1 != h2, log10Of2_3, log10Of1_3 ); p_base = qual / -10.0 + ( h1 != h2 ? log10Of1_3 : log10Of1_3 ); @@ -329,48 +300,6 @@ public class GenotypeLikelihoods implements GenotypeGenerator { this.sort(); } - public void applySecondBaseDistributionPrior(String primaryBases, String secondaryBases) { - for (int genotypeIndex = 0; genotypeIndex < genotypes.length; genotypeIndex++) { - char firstAllele = genotypes[genotypeIndex].charAt(0); - char secondAllele = genotypes[genotypeIndex].charAt(1); - - int offIsGenotypic = 0; - int offTotal = 0; - - int onIsGenotypic = 0; - int onTotal = 0; - - for (int pileupIndex = 0; pileupIndex < primaryBases.length(); pileupIndex++) { - char primaryBase = primaryBases.charAt(pileupIndex); - - if (secondaryBases != null) { - char secondaryBase = secondaryBases.charAt(pileupIndex); - - if (primaryBase != firstAllele && primaryBase != secondAllele) { - if (secondaryBase == firstAllele || secondaryBase == secondAllele) { - offIsGenotypic++; - } - offTotal++; - } else { - if (secondaryBase == firstAllele || secondaryBase == secondAllele) { - onIsGenotypic++; - } - onTotal++; - } - } - } - - double offPrior = MathUtils.binomialProbability(offIsGenotypic, offTotal, offNextBestBasePriors.get(genotypes[genotypeIndex])); - double onPrior = MathUtils.binomialProbability(onIsGenotypic, onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex])); - - double logOffPrior = MathUtils.compareDoubles(offPrior, 0.0, 1e-10) == 0 ? Math.log10(Double.MIN_VALUE) : Math.log10(offPrior); - double logOnPrior = MathUtils.compareDoubles(onPrior, 0.0, 1e-10) == 0 ? Math.log10(Double.MIN_VALUE) : Math.log10(onPrior); - - likelihoods[genotypeIndex] += logOffPrior + logOnPrior; - } - this.sort(); - } - public double LodVsNextBest() { this.sort(); return sorted_likelihoods[0] - sorted_likelihoods[1]; @@ -438,14 +367,13 @@ public class GenotypeLikelihoods implements GenotypeGenerator { int offset = pileup.getOffsets().get(i); char base = read.getReadString().charAt(offset); byte qual = read.getBaseQualities()[offset]; - add(ref, base, qual); + add(ref, base, qual); } // save off the likelihoods if (likelihoods == null || likelihoods.length == 0) return null; // Apply the two calculations applyPrior(ref); - applySecondBaseDistributionPrior(pileup.getBases(), pileup.getSecondaryBasePileup()); // lets setup the locus List lst = new ArrayList(); @@ -454,5 +382,4 @@ public class GenotypeLikelihoods implements GenotypeGenerator { } return new SSGGenotypeCall(discoveryMode,ref,2,lst,likelihoods,pileup); } - }