From deb45865599df06e2fb96c177c377f6c982c479c Mon Sep 17 00:00:00 2001 From: Guillermo del Angel Date: Sat, 24 Mar 2012 21:49:43 -0400 Subject: [PATCH] Next intermediate commit for new pool caller structure: a) Bug fixes in pool GL computation. Now, correct GL's are returned per each pool to the UG engine. Work still needs to be done in redoing interface with exact model. b) Added unit tests for new MathUtils dot product and logDotProduct functions. c) Refactorings of UnifiedGentotyperEngine since N (size of prior/posterior arrays) is no longer necessarily nSamples+1 but, in general, nSamplesPerPool*nPools+1 --- build.xml | 4 ++-- .../gatk/walkers/genotyper/UnifiedGenotyper.java | 2 +- .../walkers/genotyper/UnifiedGenotyperEngine.java | 11 ++++++----- .../org/broadinstitute/sting/utils/MathUtils.java | 4 ++-- .../sting/utils/MathUtilsUnitTest.java | 12 ++++++++++++ 5 files changed, 23 insertions(+), 10 deletions(-) diff --git a/build.xml b/build.xml index ce07138b6..9715071cd 100644 --- a/build.xml +++ b/build.xml @@ -955,8 +955,8 @@ - - + + diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java index 0eb35d299..820d58837 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java @@ -221,7 +221,7 @@ public class UnifiedGenotyper extends LocusWalker headerInfo = getHeaderInfo(); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java index 1382306c6..6af3bbd1e 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java @@ -115,11 +115,11 @@ public class UnifiedGenotyperEngine { // --------------------------------------------------------------------------------------------------------- @Requires({"toolkit != null", "UAC != null"}) public UnifiedGenotyperEngine(GenomeAnalysisEngine toolkit, UnifiedArgumentCollection UAC) { - this(toolkit, UAC, Logger.getLogger(UnifiedGenotyperEngine.class), null, null, SampleUtils.getSAMFileSamples(toolkit.getSAMFileHeader())); + this(toolkit, UAC, Logger.getLogger(UnifiedGenotyperEngine.class), null, null, SampleUtils.getSAMFileSamples(toolkit.getSAMFileHeader()), 2*(SampleUtils.getSAMFileSamples(toolkit.getSAMFileHeader()).size())); } - @Requires({"toolkit != null", "UAC != null", "logger != null", "samples != null && samples.size() > 0"}) - public UnifiedGenotyperEngine(GenomeAnalysisEngine toolkit, UnifiedArgumentCollection UAC, Logger logger, PrintStream verboseWriter, VariantAnnotatorEngine engine, Set samples) { + @Requires({"toolkit != null", "UAC != null", "logger != null", "samples != null && samples.size() > 0","N>0"}) + public UnifiedGenotyperEngine(GenomeAnalysisEngine toolkit, UnifiedArgumentCollection UAC, Logger logger, PrintStream verboseWriter, VariantAnnotatorEngine engine, Set samples, int N) { this.BAQEnabledOnCMDLine = toolkit.getArguments().BAQMode != BAQ.CalculationMode.OFF; genomeLocParser = toolkit.getGenomeLocParser(); this.samples = new TreeSet(samples); @@ -130,7 +130,8 @@ public class UnifiedGenotyperEngine { this.verboseWriter = verboseWriter; this.annotationEngine = engine; - N = 2 * this.samples.size(); + this.N = N; + log10AlleleFrequencyPriorsSNPs = new double[UAC.MAX_ALTERNATE_ALLELES][N+1]; log10AlleleFrequencyPriorsIndels = new double[UAC.MAX_ALTERNATE_ALLELES][N+1]; computeAlleleFrequencyPriors(N, log10AlleleFrequencyPriorsSNPs, UAC.heterozygosity); @@ -703,7 +704,7 @@ public class UnifiedGenotyperEngine { for (int i = 0; i < glmClasses.size(); i++) { Class glmClass = glmClasses.get(i); String key = glmClass.getSimpleName().replaceAll("GenotypeLikelihoodsCalculationModel","").toUpperCase(); - System.out.println("KEY:"+key+"\t" + glmClass.getSimpleName()); + //System.out.println("KEY:"+key+"\t" + glmClass.getSimpleName()); try { Object args[] = new Object[]{UAC,logger}; Constructor c = glmClass.getDeclaredConstructor(UnifiedArgumentCollection.class, Logger.class); diff --git a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java index ad4264d4a..34394ff39 100644 --- a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java @@ -1050,7 +1050,7 @@ public class MathUtils { /** * Given two log-probability vectors, compute log of vector product of them: - * in Matlab notation, return log(10.*x'*10.^y) + * in Matlab notation, return log10(10.*x'*10.^y) * @param x vector 1 * @param y vector 2 * @return a double representing log (dotProd(10.^x,10.^y) @@ -1065,7 +1065,7 @@ public class MathUtils { tmpVec[k] = x[k]+y[k]; } - return sumLog10(tmpVec); + return log10sumLog10(tmpVec); diff --git a/public/java/test/org/broadinstitute/sting/utils/MathUtilsUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/MathUtilsUnitTest.java index 482f4da80..9e01eb5ae 100755 --- a/public/java/test/org/broadinstitute/sting/utils/MathUtilsUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/MathUtilsUnitTest.java @@ -284,6 +284,18 @@ public class MathUtilsUnitTest extends BaseTest { Assert.assertTrue(compareDoubleArrays(MathUtils.normalizeFromLog10(new double[] {-1.0, -3.0, -1.0, -2.0}), new double[] {0.1 * 1.0 / 0.211, 0.001 * 1.0 / 0.211, 0.1 * 1.0 / 0.211, 0.01 * 1.0 / 0.211})); } + @Test + public void testDotProduct() { + Assert.assertEquals(MathUtils.dotProduct(new Double[]{-5.0,-3.0,2.0}, new Double[]{6.0,7.0,8.0}),-35.0); + Assert.assertEquals(MathUtils.dotProduct(new Double[]{-5.0}, new Double[]{6.0}),-30.0); + } + + @Test + public void testLogDotProduct() { + Assert.assertEquals(MathUtils.logDotProduct(new double[]{-5.0,-3.0,2.0}, new double[]{6.0,7.0,8.0}),10.0); + Assert.assertEquals(MathUtils.logDotProduct(new double[]{-5.0}, new double[]{6.0}),1.0); + } + /** * Private function used by testNormalizeFromLog10() */