diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyMetricsWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyMetricsWalker.java index 5fe59a4fd..654592a5f 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyMetricsWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyMetricsWalker.java @@ -61,9 +61,9 @@ public class AlleleFrequencyMetricsWalker extends BasicLociWalker= LOD_cutoff) { num_loci_confident += 1; } + if (Math.abs(alleleFreq.lodVsRef) >= LOD_cutoff) { num_loci_confident += 1; } - if (alleleFreq.qstar > 0.0 && alleleFreq.LOD >= LOD_cutoff) + if (alleleFreq.qstar > 0.0 && alleleFreq.lodVsRef >= LOD_cutoff) { // Confident variant. @@ -75,20 +75,7 @@ public class AlleleFrequencyMetricsWalker extends BasicLociWalker= LOD_cutoff && has_hapmap_chip_genotype) { - //confident_calls_at_hapmap_chip_sites++; + if (has_hapmap_chip_genotype) { // convert hapmap call to mixture of ref/nonref String hapmap_genotype = hapmap_chip_genotype.getFeature(); @@ -107,43 +94,39 @@ public class AlleleFrequencyMetricsWalker extends BasicLociWalker= LOD_cutoff) { + + // Calculate genotyping performance - did we get the correct genotype of the N+1 choices? + if (hapmap_q != -1 && hapmap_q == alleleFreq.qstar) { + hapmap_genotype_correct++; + }else{ + hapmap_genotype_incorrect++; + System.out.printf(" INCORRECT GENOTYPE Bases: %s", AlleleFrequencyWalker.getBases(context)); + //AlleleFrequencyWalker.print_base_qual_matrix(AlleleFrequencyWalker.getOneBaseQuals(context)); + } } - // Now calculate ref / var performance - did we correctly classify the site as - // reference or variant without regard to genotype; i.e. het/hom "miscalls" don't matter here - boolean hapmap_var = hapmap_q != 0.0; - boolean called_var = qstar != 0.0; - if (hapmap_var != called_var) { - hapmap_refvar_incorrect++; - }else{ - hapmap_refvar_correct++; + if (alleleFreq.lodVsRef >= LOD_cutoff || -1 * alleleFreq.lodVsRef >= LOD_cutoff) { + + // Now calculate ref / var performance - did we correctly classify the site as + // reference or variant without regard to genotype; i.e. het/hom "miscalls" don't matter here + boolean hapmap_var = hapmap_q != 0.0; + boolean called_var = alleleFreq.qstar != 0.0; + if (hapmap_q != -1 && hapmap_var != called_var) { + hapmap_refvar_incorrect++; + }else{ + hapmap_refvar_correct++; + System.out.printf(" INCORRECT REFVAR CALL Bases: %s\n", AlleleFrequencyWalker.getBases(context)); + } } + System.out.print("\n"); } - /*String hapmap_genotype = hapmap_chip_genotype.getFeature(); - String called_genotype = alleleFreq.asString(); - System.out.format("HAPMAP %s %s %.2f\n", hapmap_genotype, called_genotype, alleleFreq.LOD); - // There is a hapmap site here, is it correct? - if (hapmap_genotype.equals(called_genotype)) - { - hapmap_tp++; - } else { - hapmap_fp++; - } */ - return alleleFreq; } @@ -190,7 +173,7 @@ public class AlleleFrequencyMetricsWalker extends BasicLociWalker= 5) || (alleleFreq.LOD <= -5)) { System.out.print(alleleFreq.asTabularString()); } + //if ((alleleFreq.lodVsRef >= 5) || (alleleFreq.lodVsRef <= -5)) { System.out.print(alleleFreq.asTabularString()); } System.out.print(alleleFreq.asTabularString()); if (this.num_loci_total % 1000 == 0) { printMetrics(); } diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java index 6e78f1d60..425628da8 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java @@ -8,6 +8,7 @@ import org.broadinstitute.sting.playground.utils.AlleleFrequencyEstimate; import net.sf.samtools.SAMRecord; import java.util.List; +import java.util.Arrays; public class AlleleFrequencyWalker extends BasicLociWalker { @@ -159,13 +160,14 @@ public class AlleleFrequencyWalker extends BasicLociWalker best_posterior) - { - best_qstar = qstar; - best_qhat = q; - best_posterior = posterior; + if (posterior > bestMixtures[qstar_N].posterior) + bestMixtures[qstar_N] = new MixtureLikelihood(posterior, qstar, q); - best_pDq = pDq; - best_pqG = pqG; - best_pG = pG; - } //System.out.format("%.2f %.2f %5.2f %5.2f %5.2f %5.2f\n", q, qstar, pDq, pqG, pG, posterior); } } - double posterior_null_hyp = P_D_q(bases, quals, 0.0, refnum, altnum) + P_q_G(bases, N, 0.0, epsilon, 0) + P_G(N, 0); - double LOD = best_posterior - posterior_null_hyp; + // First reverset sort NONREF mixtures according to highest posterior probabililty + Arrays.sort(bestMixtures, 1, N+1); + + // Calculate Lod of any variant call versus the reference call + // Answers how confident are we in the best variant (nonref) mixture versus the null hypothesis + // reference mixture - qhat = qstar = 0.0 + double lodVarVsRef = bestMixtures[1].posterior - posterior_null_hyp; + + // Now reverse sort ALL mixtures according to highest posterior probability + Arrays.sort(bestMixtures); + + // Calculate Lod of the mixture versus other possible + // Answers how confident are we in the best mixture versus the next best mixture + double lodBestVsNextBest = bestMixtures[0].posterior - bestMixtures[1].posterior; AlleleFrequencyEstimate alleleFreq = new AlleleFrequencyEstimate(location, num2nuc[refnum], num2nuc[altnum], N, - best_qhat, - best_qstar, - LOD, + bestMixtures[0].qhat, + bestMixtures[0].qstar, + lodVarVsRef, + lodBestVsNextBest, depth); return alleleFreq; } @@ -213,6 +221,12 @@ public class AlleleFrequencyWalker extends BasicLociWalker= 5) || (alleleFreq.LOD <= -5)) { System.out.print(alleleFreq.asTabularString()); } + if ((alleleFreq.lodVsRef >= 5) || (alleleFreq.lodVsRef <= -5)) { System.out.print(alleleFreq.asTabularString()); } return 0; } @@ -393,7 +407,7 @@ public class AlleleFrequencyWalker extends BasicLociWalker