From 32389dc0a93db49be858ea1ee1b2e34b284a7fda Mon Sep 17 00:00:00 2001 From: ebanks Date: Fri, 14 May 2010 19:17:46 +0000 Subject: [PATCH] Fixed GQ estimate when chosen genotype isn't the most likely according to the GLs. git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@3362 348d0f76-0448-11de-a6fe-93d51630548a --- .../DiploidGenotypeCalculationModel.java | 23 ++++++++++++++----- .../UnifiedGenotyperIntegrationTest.java | 4 ++-- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java index 37fecd640..db90022bc 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java @@ -289,12 +289,23 @@ public class DiploidGenotypeCalculationModel extends JointEstimateGenotypeCalcul int genotype = indexes[index]; double score; - if ( genotype == GenotypeType.REF.ordinal() ) - score = matrix[index][GenotypeType.REF.ordinal()] - Math.max(matrix[index][GenotypeType.HET.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]); - else if ( genotype == GenotypeType.HET.ordinal() ) - score = matrix[index][GenotypeType.HET.ordinal()] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]); - else // ( genotype == GenotypeType.HOM.ordinal() ) - score = matrix[index][GenotypeType.HOM.ordinal()] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HET.ordinal()]); + + int maxEntry = MathUtils.maxElementIndex(matrix[index]); + // if the max value is for the most likely genotype, we can compute next vs. next best + if ( genotype == maxEntry ) { + if ( genotype == GenotypeType.REF.ordinal() ) + score = matrix[index][genotype] - Math.max(matrix[index][GenotypeType.HET.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]); + else if ( genotype == GenotypeType.HET.ordinal() ) + score = matrix[index][genotype] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HOM.ordinal()]); + else // ( genotype == GenotypeType.HOM.ordinal() ) + score = matrix[index][genotype] - Math.max(matrix[index][GenotypeType.REF.ordinal()], matrix[index][GenotypeType.HET.ordinal()]); + } + // otherwise, we need to calculate the probability of the genotype + else { + double[] normalized = MathUtils.normalizeFromLog10(matrix[index]); + double chosenGenotype = normalized[genotype]; + score = -1.0 * Math.log10(1.0 - chosenGenotype); + } samplesToGenotypes.put(sample, new Pair(genotype, Math.abs(score))); index++; diff --git a/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java b/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java index 52735164a..d812a60af 100755 --- a/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java +++ b/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java @@ -77,8 +77,8 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest { @Test public void testParameter() { HashMap e = new HashMap(); - e.put( "-genotype", "310d3f8ac2a1100d82fe6241e201e1d8" ); - e.put( "-all_bases", "470196b3f9a2ed622383e25552972a3e" ); + e.put( "-genotype", "03e228232eeb0a32ea3791daec29248c" ); + e.put( "-all_bases", "fe894f335e048d89fccd7a7b39a105a1" ); e.put( "--min_base_quality_score 26", "a1142c69f1951a05d5a33215c69133af" ); e.put( "--min_mapping_quality_score 26", "796d1ea790481a1a257609279f8076bb" ); e.put( "--max_mismatches_in_40bp_window 5", "379bd56d4402e974127e0677ff51d040" );