From e7c823b4bf591d8eaa3b3213ba8e2db8b56a210f Mon Sep 17 00:00:00 2001 From: David Benjamin Date: Fri, 29 Jun 2018 12:06:53 -0400 Subject: [PATCH 1/3] put the epsilon is isPolymorphic for GVCF mode --- .../walkers/genotyper/afcalc/AFCalculationResult.java | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AFCalculationResult.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AFCalculationResult.java index 3810079af17..a23cd91a256 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AFCalculationResult.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AFCalculationResult.java @@ -21,6 +21,12 @@ public final class AFCalculationResult { private static final int AF1p = 1; private static final int LOG_10_ARRAY_SIZES = 2; + // In GVCF mode the STANDARD_CONFIDENCE_FOR_CALLING is 0 by default, and it's nice having this easily-interpretable + // threshold that says "call anything with any evidence at all." The problem is that *everything* has at least some evidence, + // so this would end up putting every site, or at least too many sites, in the gvcf. Thus this parameter is in place to say + // that "0" really means "epsilon." + private static final double EPSILON = 1.0e-10; + private final double[] log10LikelihoodsOfAC; private final double[] log10PriorsOfAC; private final double[] log10PosteriorsOfAC; @@ -206,7 +212,7 @@ public String toString() { */ public boolean isPolymorphic(final Allele allele, final double log10minPNonRef) { Utils.nonNull(allele); - return getLog10PosteriorOfAFEq0ForAllele(allele) < log10minPNonRef; + return getLog10PosteriorOfAFEq0ForAllele(allele) + EPSILON < log10minPNonRef; } /** From e7f1e4186e13fdca9c6b9b422b493118556373ab Mon Sep 17 00:00:00 2001 From: David Benjamin Date: Fri, 29 Jun 2018 12:23:01 -0400 Subject: [PATCH 2/3] updated tests after checking that all changes were at QUAL=0 sites as expected --- .../expected.testGVCFMode.gatk4.alleleSpecific.g.vcf | 8 ++++---- .../haplotypecaller/expected.testGVCFMode.gatk4.g.vcf | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.alleleSpecific.g.vcf b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.alleleSpecific.g.vcf index a7ac4db918f..f3f2e6551be 100644 --- a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.alleleSpecific.g.vcf +++ b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.alleleSpecific.g.vcf @@ -227,7 +227,7 @@ 20 10004219 . C . . END=10004221 GT:DP:GQ:MIN_DP:PL 0/0:52:72:51:0,72,1080 20 10004222 . C CA, 12.96 . AS_RAW_BaseQRankSum=|0.3,1|NaN;AS_RAW_MQ=121951.00|23791.00|0.00;AS_RAW_MQRankSum=|-2.1,1|NaN;AS_RAW_ReadPosRankSum=|-0.1,1|NaN;AS_SB_TABLE=18,13|4,4|0,0;BaseQRankSum=0.335;ClippingRankSum=0.000;DP=55;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=-2.042;RAW_MQ=198329.00;ReadPosRankSum=-0.052 GT:AD:DP:GQ:PL:SB 0/1:31,8,0:39:50:50,0,923,153,952,1105:18,13,4,4 20 10004223 . A AG, 47.73 . AS_RAW_BaseQRankSum=|-0.5,1|NaN;AS_RAW_MQ=108277.00|39441.00|0.00;AS_RAW_MQRankSum=|1.8,1|NaN;AS_RAW_ReadPosRankSum=|2.2,1|NaN;AS_SB_TABLE=16,20|5,6|0,0;BaseQRankSum=-0.480;ClippingRankSum=0.000;DP=63;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=1.886;RAW_MQ=214369.00;ReadPosRankSum=2.234 GT:AD:DP:GQ:PL:SB 0/1:36,11,0:47:85:85,0,982,193,1017,1210:16,20,5,6 -20 10004224 . A G, 0 . AS_RAW_BaseQRankSum=|-0.4,1|NaN;AS_RAW_MQ=152618.00|10800.00|0.00;AS_RAW_MQRankSum=|0.6,1|NaN;AS_RAW_ReadPosRankSum=|1.7,1|NaN;AS_SB_TABLE=20,28|2,1|0,0;BaseQRankSum=-0.343;ClippingRankSum=0.000;DP=67;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=0.609;RAW_MQ=227310.00;ReadPosRankSum=1.728 GT:AD:DP:GQ:PL:SB 0/0:48,3,0:51:95:0,95,1492,143,1500,1549:20,28,2,1 +20 10004224 . A . . END=10004224 GT:DP:GQ:MIN_DP:PL 0/0:70:0:70:0,0,450 20 10004225 . A . . END=10004240 GT:DP:GQ:MIN_DP:PL 0/0:70:99:65:0,114,1710 20 10004241 . G . . END=10004241 GT:DP:GQ:MIN_DP:PL 0/0:66:0:66:0,0,1965 20 10004242 . A . . END=10004350 GT:DP:GQ:MIN_DP:PL 0/0:73:99:58:0,120,1800 @@ -1210,7 +1210,7 @@ 20 10074807 . C . . END=10075042 GT:DP:GQ:MIN_DP:PL 0/0:68:99:56:0,120,1800 20 10075043 . T C, 2405.77 . AS_RAW_BaseQRankSum=||;AS_RAW_MQ=0.00|216841.00|0.00;AS_RAW_MQRankSum=||;AS_RAW_ReadPosRankSum=||;AS_SB_TABLE=0,0|36,25|0,0;DP=62;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=220441.00 GT:AD:DP:GQ:PL:SB 1/1:0,61,0:61:99:2434,184,0,2434,184,2434:0,0,36,25 20 10075044 . C . . END=10075141 GT:DP:GQ:MIN_DP:PL 0/0:55:99:51:0,120,1800 -20 10075142 . G A, 0 . AS_RAW_BaseQRankSum=|-2.5,1|NaN;AS_RAW_MQ=211318.00|11641.00|0.00;AS_RAW_MQRankSum=|-1.4,1|NaN;AS_RAW_ReadPosRankSum=|-2.4,1|NaN;AS_SB_TABLE=36,25|0,4|0,0;BaseQRankSum=-2.480;ClippingRankSum=0.000;DP=65;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=-1.359;RAW_MQ=222959.00;ReadPosRankSum=-2.348 GT:AD:DP:GQ:PL:SB 0/0:61,4,0:65:92:0,92,2165,183,2177,2268:36,25,0,4 +20 10075142 . G . . END=10075142 GT:DP:GQ:MIN_DP:PL 0/0:65:91:65:0,91,2199 20 10075143 . A . . END=10075167 GT:DP:GQ:MIN_DP:PL 0/0:81:99:68:0,120,1800 20 10075168 . C T, 3612.77 . AS_RAW_BaseQRankSum=||;AS_RAW_MQ=0.00|309981.00|0.00;AS_RAW_MQRankSum=||;AS_RAW_ReadPosRankSum=||;AS_SB_TABLE=0,0|38,50|0,0;DP=88;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=309981.00 GT:AD:DP:GQ:PL:SB 1/1:0,88,0:88:99:3641,264,0,3641,264,3641:0,0,38,50 20 10075169 . G . . END=10075488 GT:DP:GQ:MIN_DP:PL 0/0:67:99:49:0,99,1485 @@ -1249,9 +1249,9 @@ 20 10081751 . T . . END=10081799 GT:DP:GQ:MIN_DP:PL 0/0:77:99:68:0,120,1800 20 10081800 . C T, 1109.77 . AS_RAW_BaseQRankSum=|0.2,1|NaN;AS_RAW_MQ=126000.00|118800.00|0.00;AS_RAW_MQRankSum=|0.0,1|NaN;AS_RAW_ReadPosRankSum=|-1.1,1|NaN;AS_SB_TABLE=26,9|17,16|0,0;BaseQRankSum=0.241;ClippingRankSum=0.000;DP=68;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.000;RAW_MQ=244800.00;ReadPosRankSum=-1.082 GT:AD:DP:GQ:PL:SB 0/1:35,33,0:68:99:1138,0,1222,1244,1322,2565:26,9,17,16 20 10081801 . A . . END=10082884 GT:DP:GQ:MIN_DP:PL 0/0:76:99:43:0,111,1665 -20 10082885 . A C, 0 . AS_RAW_BaseQRankSum=|-2.3,1|NaN;AS_RAW_MQ=147600.00|7200.00|0.00;AS_RAW_MQRankSum=|0.0,1|NaN;AS_RAW_ReadPosRankSum=|-2.4,1|NaN;AS_SB_TABLE=18,23|0,2|0,0;BaseQRankSum=-2.266;ClippingRankSum=0.000;DP=43;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=0.000;RAW_MQ=154800.00;ReadPosRankSum=-2.350 GT:AD:DP:GQ:PGT:PID:PL:SB 0/0:41,2,0:43:93:0|1:10082885_A_C:0,93,1525,123,1531,1561:18,23,0,2 +20 10082885 . A . . END=10082885 GT:DP:GQ:MIN_DP:PL 0/0:43:83:43:0,83,1521 20 10082886 . C . . END=10082891 GT:DP:GQ:MIN_DP:PL 0/0:43:99:42:0,99,1485 -20 10082892 . C T, 1712.77 . AS_RAW_BaseQRankSum=||;AS_RAW_MQ=0.00|155641.00|0.00;AS_RAW_MQRankSum=||;AS_RAW_ReadPosRankSum=||;AS_SB_TABLE=0,0|20,24|0,0;DP=44;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=155641.00 GT:AD:DP:GQ:PGT:PID:PL:SB 1/1:0,44,0:44:99:0|1:10082885_A_C:1741,132,0,1741,132,1741:0,0,20,24 +20 10082892 . C T, 1712.77 . AS_RAW_BaseQRankSum=||;AS_RAW_MQ=0.00|155641.00|0.00;AS_RAW_MQRankSum=||;AS_RAW_ReadPosRankSum=||;AS_SB_TABLE=0,0|20,24|0,0;DP=44;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=155641.00 GT:AD:DP:GQ:PL:SB 1/1:0,44,0:44:99:1741,132,0,1741,132,1741:0,0,20,24 20 10082893 . C . . END=10082901 GT:DP:GQ:MIN_DP:PL 0/0:41:99:38:0,99,1485 20 10082902 . C . . END=10082903 GT:DP:GQ:MIN_DP:PL 0/0:38:96:38:0,96,1440 20 10082904 . C . . END=10082905 GT:DP:GQ:MIN_DP:PL 0/0:40:99:39:0,102,1530 diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.g.vcf b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.g.vcf index 18ca377ecdb..38214c214fb 100644 --- a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.g.vcf +++ b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.g.vcf @@ -220,7 +220,7 @@ 20 10004219 . C . . END=10004221 GT:DP:GQ:MIN_DP:PL 0/0:52:72:51:0,72,1080 20 10004222 . C CA, 12.96 . BaseQRankSum=0.335;ClippingRankSum=0.000;DP=55;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=-2.042;RAW_MQ=198329.00;ReadPosRankSum=-0.052 GT:AD:DP:GQ:PL:SB 0/1:31,8,0:39:50:50,0,923,153,952,1105:18,13,4,4 20 10004223 . A AG, 47.73 . BaseQRankSum=-0.480;ClippingRankSum=0.000;DP=63;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=1.886;RAW_MQ=214369.00;ReadPosRankSum=2.234 GT:AD:DP:GQ:PL:SB 0/1:36,11,0:47:85:85,0,982,193,1017,1210:16,20,5,6 -20 10004224 . A G, 0 . BaseQRankSum=-0.343;ClippingRankSum=0.000;DP=67;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=0.609;RAW_MQ=227310.00;ReadPosRankSum=1.728 GT:AD:DP:GQ:PL:SB 0/0:48,3,0:51:95:0,95,1492,143,1500,1549:20,28,2,1 +20 10004224 . A . . END=10004224 GT:DP:GQ:MIN_DP:PL 0/0:70:0:70:0,0,450 20 10004225 . A . . END=10004240 GT:DP:GQ:MIN_DP:PL 0/0:70:99:65:0,114,1710 20 10004241 . G . . END=10004241 GT:DP:GQ:MIN_DP:PL 0/0:66:0:66:0,0,1965 20 10004242 . A . . END=10004350 GT:DP:GQ:MIN_DP:PL 0/0:73:99:58:0,120,1800 @@ -1203,7 +1203,7 @@ 20 10074807 . C . . END=10075042 GT:DP:GQ:MIN_DP:PL 0/0:68:99:56:0,120,1800 20 10075043 . T C, 2405.77 . DP=62;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=220441.00 GT:AD:DP:GQ:PL:SB 1/1:0,61,0:61:99:2434,184,0,2434,184,2434:0,0,36,25 20 10075044 . C . . END=10075141 GT:DP:GQ:MIN_DP:PL 0/0:55:99:51:0,120,1800 -20 10075142 . G A, 0 . BaseQRankSum=-2.480;ClippingRankSum=0.000;DP=65;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=-1.359;RAW_MQ=222959.00;ReadPosRankSum=-2.348 GT:AD:DP:GQ:PL:SB 0/0:61,4,0:65:92:0,92,2165,183,2177,2268:36,25,0,4 +20 10075142 . G . . END=10075142 GT:DP:GQ:MIN_DP:PL 0/0:65:91:65:0,91,2199 20 10075143 . A . . END=10075167 GT:DP:GQ:MIN_DP:PL 0/0:81:99:68:0,120,1800 20 10075168 . C T, 3612.77 . DP=88;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=309981.00 GT:AD:DP:GQ:PL:SB 1/1:0,88,0:88:99:3641,264,0,3641,264,3641:0,0,38,50 20 10075169 . G . . END=10075488 GT:DP:GQ:MIN_DP:PL 0/0:67:99:49:0,99,1485 @@ -1242,9 +1242,9 @@ 20 10081751 . T . . END=10081799 GT:DP:GQ:MIN_DP:PL 0/0:77:99:68:0,120,1800 20 10081800 . C T, 1109.77 . BaseQRankSum=0.241;ClippingRankSum=0.000;DP=68;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.000;RAW_MQ=244800.00;ReadPosRankSum=-1.082 GT:AD:DP:GQ:PL:SB 0/1:35,33,0:68:99:1138,0,1222,1244,1322,2565:26,9,17,16 20 10081801 . A . . END=10082884 GT:DP:GQ:MIN_DP:PL 0/0:76:99:43:0,111,1665 -20 10082885 . A C, 0 . BaseQRankSum=-2.266;ClippingRankSum=0.000;DP=43;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=0.000;RAW_MQ=154800.00;ReadPosRankSum=-2.350 GT:AD:DP:GQ:PGT:PID:PL:SB 0/0:41,2,0:43:93:0|1:10082885_A_C:0,93,1525,123,1531,1561:18,23,0,2 +20 10082885 . A . . END=10082885 GT:DP:GQ:MIN_DP:PL 0/0:43:83:43:0,83,1521 20 10082886 . C . . END=10082891 GT:DP:GQ:MIN_DP:PL 0/0:43:99:42:0,99,1485 -20 10082892 . C T, 1712.77 . DP=44;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=155641.00 GT:AD:DP:GQ:PGT:PID:PL:SB 1/1:0,44,0:44:99:0|1:10082885_A_C:1741,132,0,1741,132,1741:0,0,20,24 +20 10082892 . C T, 1712.77 . DP=44;ExcessHet=3.0103;MLEAC=2,0;MLEAF=1.00,0.00;RAW_MQ=155641.00 GT:AD:DP:GQ:PL:SB 1/1:0,44,0:44:99:1741,132,0,1741,132,1741:0,0,20,24 20 10082893 . C . . END=10082901 GT:DP:GQ:MIN_DP:PL 0/0:41:99:38:0,99,1485 20 10082902 . C . . END=10082903 GT:DP:GQ:MIN_DP:PL 0/0:38:96:38:0,96,1440 20 10082904 . C . . END=10082905 GT:DP:GQ:MIN_DP:PL 0/0:40:99:39:0,102,1530 From dd34dff2b6c91f6abbd8e1d2f6a5e0b6281a3bb0 Mon Sep 17 00:00:00 2001 From: David Benjamin Date: Fri, 29 Jun 2018 12:41:59 -0400 Subject: [PATCH 3/3] edited another test vcf --- ...S.15PCT.20.10100000-10150000.gatk3.8-1-1-gdde23f56a6.g.vcf | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.CEUTrio.HiSeq.WGS.b37.NA12878.CONTAMINATED.WITH.HCC1143.NORMALS.15PCT.20.10100000-10150000.gatk3.8-1-1-gdde23f56a6.g.vcf b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.CEUTrio.HiSeq.WGS.b37.NA12878.CONTAMINATED.WITH.HCC1143.NORMALS.15PCT.20.10100000-10150000.gatk3.8-1-1-gdde23f56a6.g.vcf index 1dd89aae24b..8fffa2774bf 100644 --- a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.CEUTrio.HiSeq.WGS.b37.NA12878.CONTAMINATED.WITH.HCC1143.NORMALS.15PCT.20.10100000-10150000.gatk3.8-1-1-gdde23f56a6.g.vcf +++ b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.CEUTrio.HiSeq.WGS.b37.NA12878.CONTAMINATED.WITH.HCC1143.NORMALS.15PCT.20.10100000-10150000.gatk3.8-1-1-gdde23f56a6.g.vcf @@ -494,7 +494,7 @@ 20 10132679 . T . . END=10132679 GT:DP:GQ:MIN_DP:PL 0/0:63:87:63:0,87,1305 20 10132680 . G . . END=10132680 GT:DP:GQ:MIN_DP:PL 0/0:64:90:64:0,90,1350 20 10132681 . T . . END=10132681 GT:DP:GQ:MIN_DP:PL 0/0:63:87:63:0,87,1305 -20 10132682 . ATGTGTGTGTGTGTG A,ATGTGTGTGTGTG, 0 . BaseQRankSum=-1.153;ClippingRankSum=0.000;DP=72;ExcessHet=3.0103;MLEAC=0,0,0;MLEAF=0.00,0.00,0.00;MQRankSum=2.033;RAW_MQ=187276.00;ReadPosRankSum=-1.246 GT:AD:DP:GQ:PL:SB 0/0:20,0,3,0:23:81:0,81,1219,85,714,1074,113,1009,938,1152:10,10,0,3 +20 10132682 . A . . END=10132696 GT:DP:GQ:MIN_DP:PL 0/0:63:0:52:0,0,1240 20 10132697 . T . . END=10132697 GT:DP:GQ:MIN_DP:PL 0/0:52:0:52:0,0,1312 20 10132698 . G . . END=10132704 GT:DP:GQ:MIN_DP:PL 0/0:51:60:48:0,60,900 20 10132705 . T . . END=10132705 GT:DP:GQ:MIN_DP:PL 0/0:50:32:50:0,32,1691 @@ -591,7 +591,7 @@ 20 10140159 . A . . END=10141258 GT:DP:GQ:MIN_DP:PL 0/0:81:99:41:0,114,1540 20 10141259 . G . . END=10141259 GT:DP:GQ:MIN_DP:PL 0/0:98:0:98:0,0,2798 20 10141260 . C . . END=10141439 GT:DP:GQ:MIN_DP:PL 0/0:76:99:67:0,120,1800 -20 10141440 . T C, 0 . BaseQRankSum=-3.862;ClippingRankSum=0.000;DP=75;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=-4.968;RAW_MQ=253436.00;ReadPosRankSum=-0.008 GT:AD:DP:GQ:PL:SB 0/0:60,13,0:73:91:0,91,2184,180,2193,2281:30,30,8,5 +20 10141440 . T . . END=10141440 GT:DP:GQ:MIN_DP:PL 0/0:75:0:75:0,0,1964 20 10141441 . G . . END=10141691 GT:DP:GQ:MIN_DP:PL 0/0:74:99:59:0,120,1800 20 10141692 . A G, 0 . BaseQRankSum=-4.974;ClippingRankSum=0.000;DP=72;ExcessHet=3.0103;MLEAC=0,0;MLEAF=0.00,0.00;MQRankSum=1.855;RAW_MQ=229044.00;ReadPosRankSum=-0.506 GT:AD:DP:GQ:PGT:PID:PL:SB 0/0:58,14,0:72:55:0|1:10141692_A_G:0,55,3517,181,3526,3652:35,23,8,6 20 10141693 . A . . END=10141699 GT:DP:GQ:MIN_DP:PL 0/0:76:99:74:0,120,1800