Permalink
Browse files

Critical bugfix to AFCalcResult affecting UG/HC quality score emissio…

…n thresholds

As reported by Menachem Fromer: a critical bug in AFCalcResult:

Specifically, the implementation:
    public boolean isPolymorphic(final Allele allele, final double log10minPNonRef) {
        return getLog10PosteriorOfAFGt0ForAllele(allele) >= log10minPNonRef;
    }

seems incorrect and should probably be:

getLog10PosteriorOfAFEq0ForAllele(allele) <= log10minPNonRef

The issue here is that the 30 represents a Phred-scaled probability of *error* and it's currently being compared to a log probability of *non-error*.

Instead, we need to require that our probability of error be less than the error threshold.
This bug has only a minor impact on the calls -- hardly any sites change -- which is good.  But the inverted logic effects multi-allelic sites significantly.  Basically you only hit this logic with multiple alleles, and in that case it'\s including extra alt alleles incorrectly, and throwing out good ones.

Change was to create a new function that properly handles thresholds that are PhredScaled quality scores:

    /**
     * Same as #isPolymorphic but takes a phred-scaled quality score as input
     */
    public boolean isPolymorphicPhredScaledQual(final Allele allele, final double minPNonRefPhredScaledQual) {
        if ( minPNonRefPhredScaledQual < 0 ) throw new IllegalArgumentException("phredScaledQual " + minPNonRefPhredScaledQual + " < 0 ");
        final double log10Threshold = Math.log10(QualityUtils.qualToProb(minPNonRefPhredScaledQual));
        return isPolymorphic(allele, log10Threshold);
    }
  • Loading branch information...
Mark DePristo
Mark DePristo committed Nov 28, 2012
1 parent f5b7f54 commit 3cb83c20f229024d4ca65a65746490ae633dacb7
@@ -378,11 +378,9 @@ public VariantCallContext calculateGenotypes(final RefMetaDataTracker tracker, f
if ( alternateAllele.isReference() )
continue;
- // we are non-ref if the probability of being non-ref > the emit confidence.
- // the emit confidence is phred-scaled, say 30 => 10^-3.
- // the posterior AF > 0 is log10: -5 => 10^-5
- // we are non-ref if 10^-5 < 10^-3 => -5 < -3
- final boolean isNonRef = AFresult.isPolymorphic(alternateAllele, UAC.STANDARD_CONFIDENCE_FOR_EMITTING / -10.0);
+ // Compute if the site is considered polymorphic with sufficient confidence relative to our
+ // phred-scaled emission QUAL
+ final boolean isNonRef = AFresult.isPolymorphicPhredScaledQual(alternateAllele, UAC.STANDARD_CONFIDENCE_FOR_EMITTING);
// if the most likely AC is not 0, then this is a good alternate allele to use
if ( isNonRef ) {
@@ -28,6 +28,7 @@
import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.MathUtils;
+import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.variantcontext.Allele;
@@ -234,10 +235,20 @@ public String toString() {
*
* @return true if there's enough confidence (relative to log10minPNonRef) to reject AF == 0
*/
+ @Requires("MathUtils.goodLog10Probability(log10minPNonRef)")
public boolean isPolymorphic(final Allele allele, final double log10minPNonRef) {
return getLog10PosteriorOfAFGt0ForAllele(allele) >= log10minPNonRef;
}
+ /**
+ * Same as #isPolymorphic but takes a phred-scaled quality score as input
+ */
+ public boolean isPolymorphicPhredScaledQual(final Allele allele, final double minPNonRefPhredScaledQual) {
+ if ( minPNonRefPhredScaledQual < 0 ) throw new IllegalArgumentException("phredScaledQual " + minPNonRefPhredScaledQual + " < 0 ");
+ final double log10Threshold = Math.log10(QualityUtils.qualToProb(minPNonRefPhredScaledQual));
+ return isPolymorphic(allele, log10Threshold);
+ }
+
/**
* Are any of the alleles polymorphic w.r.t. #isPolymorphic?
*

0 comments on commit 3cb83c2

Please sign in to comment.