zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix a minor rounding bug and putz around with fractional counts in the pooled caller. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@472 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
jmaguire 2009-04-19 17:52:24 +00:00
parent af6788fa3d
commit bcba1ff424
3 changed files with 31 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java
View File

@ -16,6 +16,7 @@ import java.util.List;
import java.util.Arrays; import java.util.Arrays;
import java.util.Random; import java.util.Random;
import java.io.PrintStream; import java.io.PrintStream;
import java.io.DataInputStream;
public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate, String>// implements AllelicVariant public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate, String>// implements AllelicVariant
{ {
@ -144,25 +145,6 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
return alleleFreq; return alleleFreq;
} }
/*
static public String getBases (LocusContext context)
{
// Convert bases to CharArray
int numReads = context.getReads().size(); //numReads();
char[] bases = new char[numReads];
//int refnum = nuc2num[ref];
List<SAMRecord> reads = context.getReads();
List<Integer> offsets = context.getOffsets();
for (int i =0; i < numReads; i++ ) {
SAMRecord read = reads.get(i);
int offset = offsets.get(i);
bases[i] = read.getReadString().charAt(offset);
}
return new String(bases);
}
*/
/* /*
static public String[] getIndels(LocusContext context) static public String[] getIndels(LocusContext context)
{ {
@ -325,6 +307,8 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
bestMixtures[0] = new MixtureLikelihood(posterior_null_hyp, 0.0, 0.0); bestMixtures[0] = new MixtureLikelihood(posterior_null_hyp, 0.0, 0.0);
posteriors[0] = posterior_null_hyp; posteriors[0] = posterior_null_hyp;
//System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n", location, 0.0, 0.0, P_D_q(bases, quals, 0.0, refnum, altnum), P_q_G(bases, N, 0.0, 0, 0), P_G(N, 0), prior_alt_frequency, posterior_null_hyp, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases));
assert(! Double.isNaN(P_D_q(bases, quals, 0.0, refnum, altnum))); assert(! Double.isNaN(P_D_q(bases, quals, 0.0, refnum, altnum)));
assert(! Double.isNaN(P_q_G(bases, N, 0.0, 0, 0))); assert(! Double.isNaN(P_q_G(bases, N, 0.0, 0, 0)));
assert(! Double.isNaN(P_G(N, 0))); assert(! Double.isNaN(P_G(N, 0)));
@ -332,8 +316,6 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
assert(! Double.isNaN(posteriors[0])); assert(! Double.isNaN(posteriors[0]));
assert(! Double.isInfinite(posteriors[0])); assert(! Double.isInfinite(posteriors[0]));
//System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n", location, 0.0, 0.0, P_D_q(bases, quals, 0.0, refnum, altnum), P_q_G(bases, N, 0.0, 0, 0), P_G(N, 0), posterior_null_hyp, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases));
// qstar - true allele balances // qstar - true allele balances
//for (qstar = epsilon + ((1.0 - 2*epsilon)/N), qstar_N = 1; qstar <= 1.0; qstar += (1.0 - 2*epsilon)/N, qstar_N++) //for (qstar = epsilon + ((1.0 - 2*epsilon)/N), qstar_N = 1; qstar <= 1.0; qstar += (1.0 - 2*epsilon)/N, qstar_N++)
for (qstar_N = 1; qstar_N <= N; qstar_N += 1) for (qstar_N = 1; qstar_N <= N; qstar_N += 1)
@ -347,8 +329,8 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
bestMixtures[qstar_N] = new MixtureLikelihood(posterior, qstar, q); bestMixtures[qstar_N] = new MixtureLikelihood(posterior, qstar, q);
posteriors[qstar_N] = posterior; posteriors[qstar_N] = posterior;
//System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n", //System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n",
// location, q, qstar, pDq, pqG, pG, posterior, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases)); location, q, qstar, pDq, pqG, pG, prior_alt_frequency, posterior, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases));
} }
} }
@ -512,7 +494,10 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
{ {
// For small n and the edges, compute it directly. // For small n and the edges, compute it directly.
ans = Math.log10((double)nchoosek(n, k)) + Math.log10(Math.pow(p, k)) + Math.log10(Math.pow(1-p, n-k)); ans = Math.log10((double)nchoosek(n, k)) + Math.log10(Math.pow(p, k)) + Math.log10(Math.pow(1-p, n-k));
//System.out.printf("DBG1: %d %d %f %f\n", k, n, p, ans); //System.out.printf("DBG1: %d %d %f %f %f %f %f\n",
// k, n, p,
// nchoosek(n,k), Math.pow(p,k), Math.pow(1-p, n-k),
// ans);
} }
else else
{ {
@ -537,7 +522,7 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
ans = ans_1; ans = ans_1;
} }
//System.out.printf("DBG3: %f\n", ans); //System.out.printf("DBG3: %d %d %f %f\n", n, k, p, ans);
return ans; return ans;
} }
@ -555,7 +540,8 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
for (long i = 1; i <= k; i++) for (long i = 1; i <= k; i++)
accum = accum * (n-k+i) / i; accum = accum * (n-k+i) / i;
return accum + 0.5; // avoid rounding error //return accum + 0.5; // avoid rounding error
return accum; // avoid rounding error
/* /*
long m = n - k; long m = n - k;
@ -664,6 +650,9 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
private double prior_alt_frequency = -1.0; private double prior_alt_frequency = -1.0;
public void setAlleleFrequencyPrior(double frequency) public void setAlleleFrequencyPrior(double frequency)
{ {
assert(! Double.isNaN(frequency) );
assert(! Double.isInfinite(frequency));
if (frequency == 0) { frequency = 1e-5; } // how many epsilons do we need!? This is worrisome.
this.prior_alt_frequency = frequency; this.prior_alt_frequency = frequency;
} }
@ -755,6 +744,7 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
AlleleFrequencyWalker.binomialProb(1500,2965,0.508065); AlleleFrequencyWalker.binomialProb(1500,2965,0.508065);
AlleleFrequencyWalker.binomialProb(0,197,0.5); AlleleFrequencyWalker.binomialProb(0,197,0.5);
AlleleFrequencyWalker.binomialProb(13,197,0.5); AlleleFrequencyWalker.binomialProb(13,197,0.5);
AlleleFrequencyWalker.binomialProb(0, 7, 1e-3);
System.exit(0); System.exit(0);
{ {

9
java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCallingExperiment.java
View File

@ -183,9 +183,10 @@ public class PoolCallingExperiment extends LocusWalker<AlleleFrequencyEstimate,
{ {
for (int j = 0; j <= shallow_calls[i].N; j++) for (int j = 0; j <= shallow_calls[i].N; j++)
{ {
//System.out.printf("DBG3: %d %f %d\n", j, shallow_calls[i].posteriors[j], shallow_calls[i].N); if (Double.isInfinite(shallow_calls[i].posteriors[j])) { shallow_calls[i].posteriors[j] = -10000; }
EM_sum += shallow_calls[i].posteriors[j] * (double)j; System.out.printf("DBG3: %d %f %d\n", j, shallow_calls[i].posteriors[j], shallow_calls[i].N);
EM_sum += shallow_calls[i].N; EM_sum += Math.pow(10,shallow_calls[i].posteriors[j]) * (double)j;
EM_N += shallow_calls[i].N;
} }
} }
@ -195,6 +196,8 @@ public class PoolCallingExperiment extends LocusWalker<AlleleFrequencyEstimate,
shallow_calls_fraction_correct = correct_shallow_calls / total_shallow_calls; shallow_calls_fraction_correct = correct_shallow_calls / total_shallow_calls;
trajectory[iterations+1] = EM_alt_freq; trajectory[iterations+1] = EM_alt_freq;
likelihood_trajectory[iterations+1] = likelihood/(double)total_shallow_calls; likelihood_trajectory[iterations+1] = likelihood/(double)total_shallow_calls;
System.out.printf("DBGTRAJ %f %f %f %f\n", EM_sum, EM_N, trajectory[iterations], trajectory[iterations+1]);
} }
// 7. Compare to estimation from the pool. // 7. Compare to estimation from the pool.

9
java/src/org/broadinstitute/sting/playground/utils/AlleleFrequencyEstimate.java
View File

@ -46,7 +46,14 @@ public class AlleleFrequencyEstimate {
if( Double.isNaN(pBest)) { System.out.printf("pBest is NaN\n"); } if( Double.isNaN(pBest)) { System.out.printf("pBest is NaN\n"); }
if( Double.isNaN(pRef)) { System.out.printf("pRef is NaN\n"); } if( Double.isNaN(pRef)) { System.out.printf("pRef is NaN\n"); }
if( Double.isInfinite(lodVsRef)) { System.out.printf("lodVsRef is Infinite\n"); } if( Double.isInfinite(lodVsRef))
{
System.out.printf("lodVsRef is Infinite: %c %s\n", ref, bases);
for (int i = 0; i < posteriors.length; i++)
{
System.out.printf("POSTERIOR %d %f\n", i, posteriors[i]);
}
}
if( Double.isInfinite(lodVsNextBest)) { System.out.printf("lodVsNextBest is Infinite\n"); } if( Double.isInfinite(lodVsNextBest)) { System.out.printf("lodVsNextBest is Infinite\n"); }
if( Double.isInfinite(qhat)) { System.out.printf("qhat is Infinite\n"); } if( Double.isInfinite(qhat)) { System.out.printf("qhat is Infinite\n"); }
if( Double.isInfinite(qstar)) { System.out.printf("qstar is Infinite\n"); } if( Double.isInfinite(qstar)) { System.out.printf("qstar is Infinite\n"); }