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

Write N for the alt allele when we're hom-ref. 

Stop EM loop when we've converged (likelihood[t-1] == likelihood[t]).


git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@737 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
jmaguire 2009-05-17 13:58:11 +00:00
parent bd53bc18f9
commit 94e324b844
1 changed files with 17 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java
View File

@ -114,6 +114,7 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
double[] trajectory = new double[MAX_ITERATIONS + 1]; trajectory[0] = EM_alt_freq; double[] trajectory = new double[MAX_ITERATIONS + 1]; trajectory[0] = EM_alt_freq;
double[] likelihood_trajectory = new double[MAX_ITERATIONS + 1]; likelihood_trajectory[0] = 0.0; double[] likelihood_trajectory = new double[MAX_ITERATIONS + 1]; likelihood_trajectory[0] = 0.0;
boolean is_a_snp = false; boolean is_a_snp = false;
char alt = 'N';
for (int iterations = 0; iterations < MAX_ITERATIONS; iterations++) for (int iterations = 0; iterations < MAX_ITERATIONS; iterations++)
{ {
// 6. Re-call from shallow coverage using the estimated frequency as a prior, // 6. Re-call from shallow coverage using the estimated frequency as a prior,
@ -131,7 +132,9 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
calls[i] = callers.get(i).map(tracker, ref, contexts[i]); calls[i] = callers.get(i).map(tracker, ref, contexts[i]);
String genotype = calls[i].genotype(); String genotype = calls[i].genotype();
likelihood += calls[i].posterior(); if (calls[i].alt != 'N') { alt = calls[i].alt; }
likelihood += calls[i].pBest;
if (! FRACTIONAL_COUNTS) if (! FRACTIONAL_COUNTS)
{ {
@ -156,9 +159,18 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
} }
EM_alt_freq = EM_sum / EM_N; EM_alt_freq = EM_sum / EM_N;
trajectory[iterations+1] = EM_alt_freq; trajectory[iterations+1] = EM_alt_freq;
likelihood_trajectory[iterations+1] = likelihood/(double)EM_N; likelihood_trajectory[iterations+1] = likelihood;
//System.out.printf("DBGTRAJ %f %f %f %f\n", EM_sum, EM_N, trajectory[iterations], trajectory[iterations+1]); if (likelihood_trajectory[iterations] == likelihood_trajectory[iterations+1]) { break; }
//System.out.printf("DBGTRAJ %s %f %f %f %f %f %f\n",
// context.getLocation(),
// EM_sum,
// EM_N,
// trajectory[iterations],
// trajectory[iterations+1],
// likelihood_trajectory[iterations],
// likelihood_trajectory[iterations+1]);
} }
// 7. Output some statistics. // 7. Output some statistics.
@ -171,7 +183,8 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
//System.out.printf("DBG %f %f %c %s\n", calls[i].pBest, calls[i].pRef, ref, calls[i].bases); //System.out.printf("DBG %f %f %c %s\n", calls[i].pBest, calls[i].pRef, ref, calls[i].bases);
} }
double discovery_lod = discovery_posterior - discovery_null; double discovery_lod = discovery_posterior - discovery_null;
discovery_output_file.printf("%s %f %f %f %f\n", context.getLocation(), EM_alt_freq, discovery_posterior, discovery_null, discovery_lod); if (discovery_lod == 0) { alt = 'N'; }
discovery_output_file.printf("%s %c %c %f %f %f %f\n", context.getLocation(), ref, alt, EM_alt_freq, discovery_posterior, discovery_null, discovery_lod);
for (int i = 0; i < sample_names.size(); i++) for (int i = 0; i < sample_names.size(); i++)
{ {