From 69de3e51bf60f509d20471098619430c7129db06 Mon Sep 17 00:00:00 2001 From: ebanks Date: Fri, 12 Nov 2010 08:31:40 +0000 Subject: [PATCH] Better precision for the calculated AF value. Now looks at the total number of samples to determine how much precision is necessary. Also, changing default min BQ used for calling in UGv2 to Q17. git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4655 348d0f76-0448-11de-a6fe-93d51630548a --- .../variantcontext/VariantContextUtils.java | 29 +++++++++++-------- .../genotyper/UnifiedArgumentCollection.java | 2 +- .../UnifiedGenotyperIntegrationTest.java | 2 +- 3 files changed, 19 insertions(+), 14 deletions(-) diff --git a/java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/VariantContextUtils.java b/java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/VariantContextUtils.java index c54f10aaa..230d1b892 100755 --- a/java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/VariantContextUtils.java +++ b/java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/VariantContextUtils.java @@ -112,31 +112,36 @@ public class VariantContextUtils { // if there are alternate alleles, record the relevant tags if ( vc.getAlternateAlleles().size() > 0 ) { - ArrayList alleleFreqs = new ArrayList(); + ArrayList alleleFreqs = new ArrayList(); ArrayList alleleCounts = new ArrayList(); + double totalChromosomes = (double)vc.getChromosomeCount(); for ( Allele allele : vc.getAlternateAlleles() ) { - alleleCounts.add(vc.getChromosomeCount(allele)); - alleleFreqs.add((double)vc.getChromosomeCount(allele) / (double)vc.getChromosomeCount()); + int altChromosomes = vc.getChromosomeCount(allele); + alleleCounts.add(altChromosomes); + String freq = String.format(makePrecisionFormatStringFromDenominatorValue(totalChromosomes), ((double)altChromosomes / totalChromosomes)); + alleleFreqs.add(freq); } attributes.put(VCFConstants.ALLELE_COUNT_KEY, alleleCounts); attributes.put(VCFConstants.ALLELE_FREQUENCY_KEY, alleleFreqs); } - -// // otherwise, remove them if present and requested -// else if ( removeStaleValues ) { -// if ( attributes.containsKey(VCFConstants.ALLELE_COUNT_KEY) ) -// attributes.remove(VCFConstants.ALLELE_COUNT_KEY); -// if ( attributes.containsKey(VCFConstants.ALLELE_FREQUENCY_KEY) ) -// attributes.remove(VCFConstants.ALLELE_FREQUENCY_KEY); -// } - // otherwise, set them to 0 else { attributes.put(VCFConstants.ALLELE_COUNT_KEY, 0); attributes.put(VCFConstants.ALLELE_FREQUENCY_KEY, 0.0); } } + private static String makePrecisionFormatStringFromDenominatorValue(double maxValue) { + int precision = 1; + + while ( maxValue > 1 ) { + precision++; + maxValue /= 10.0; + } + + return "%." + precision + "f"; + } + /** * A simple but common wrapper for matching VariantContext objects using JEXL expressions */ diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedArgumentCollection.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedArgumentCollection.java index 1f66afe83..8478b2d65 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedArgumentCollection.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedArgumentCollection.java @@ -79,7 +79,7 @@ public class UnifiedArgumentCollection { // control the various parameters to be used @Argument(fullName = "min_base_quality_score", shortName = "mbq", doc = "Minimum base quality required to consider a base for calling", required = false) - public int MIN_BASE_QUALTY_SCORE = 20; + public int MIN_BASE_QUALTY_SCORE = 17; @Argument(fullName = "min_mapping_quality_score", shortName = "mmq", doc = "Minimum read mapping quality required to consider a read for calling", required = false) public int MIN_MAPPING_QUALTY_SCORE = 20; 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 3d737500f..776f61362 100755 --- a/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java +++ b/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java @@ -25,7 +25,7 @@ public class public void testMultiSamplePilot1Joint() { WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec( baseCommand + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10,022,000-10,025,000", 1, - Arrays.asList("9870dd089d0f76c8dd3e74bc9f2ba0f0")); + Arrays.asList("1dde074afe990fd38941de18bd35e188")); executeTest("testMultiSamplePilot1 - Joint Estimate", spec); }