Skip to content

Commit

Permalink
Merge pull request #1954 from broadinstitute/ms_ExcessHetMaxVal
Browse files Browse the repository at this point in the history
Changes max value output for ExcessHet
  • Loading branch information
akiezun authored Jul 1, 2016
2 parents 3f08383 + 0fe79e0 commit 0cebc03
Show file tree
Hide file tree
Showing 2 changed files with 7 additions and 2 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
*/
public final class ExcessHet extends InfoFieldAnnotation implements StandardAnnotation {
private static final double MIN_NEEDED_VALUE = 1.0E-16;
public static final double PHRED_SCALED_MIN_P_VALUE = -10.0 * Math.log10(MIN_NEEDED_VALUE);

@Override
public Map<String, Object> annotate(final ReferenceContext ref,
Expand Down Expand Up @@ -58,8 +59,9 @@ Pair<Integer, Double> calculateEH(final VariantContext vc, final GenotypesContex
final double pval = exactTest(hetCount, refCount, homCount);

//If the actual phredPval would be infinity we will probably still filter out just a very large number
if (Math.abs(pval - 0.0) < 10e-60) {
return Pair.of(Integer.MAX_VALUE, Double.POSITIVE_INFINITY);
//Since the method does not guarantee precision for any p-value smaller than 1e-16, we can return the phred scaled version
if (pval < 10e-60) {
return Pair.of(sampleCount, PHRED_SCALED_MIN_P_VALUE);
}
final double phredPval = -10.0 * Math.log10(pval);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -182,6 +182,9 @@ public void testAllHetsForLargeCohorts() {
final double hetsValue = new ExcessHet().calculateEH(allHet, allHet.getGenotypes()).getValue();

Assert.assertTrue(Math.abs(singletonValue) < Math.abs(hetsValue), String.format("singleton=%f allHets=%f", singletonValue, hetsValue));

//Since all hets is such an extreme case and the sample size is large here, we know that the p-value should be 0
Assert.assertEquals(hetsValue, ExcessHet.PHRED_SCALED_MIN_P_VALUE, DELTA_PRECISION, String.format("P-value of 0 should be phred scaled to " + ExcessHet.PHRED_SCALED_MIN_P_VALUE));
}

@DataProvider(name = "smallSets")
Expand Down

0 comments on commit 0cebc03

Please sign in to comment.