diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/BeagleOutputToVCFWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/BeagleOutputToVCFWalker.java index 680f0b4e4..d54efee0f 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/BeagleOutputToVCFWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/BeagleOutputToVCFWalker.java @@ -192,6 +192,7 @@ public class BeagleOutputToVCFWalker extends RodWalker { String markerKey = phasedLine[1]; HashMap haplotypePairs = new HashMap(); + System.out.println(markerKey); j = 2; for (String sample : samples) { @@ -235,30 +236,8 @@ public class BeagleOutputToVCFWalker extends RodWalker { return allThere; } - /* - private class MyFileReader { - private FileReader reader; - private String fileName; - public MyFileReader(String fileName) { - try{ - reader = new FileReader(fileName); - } catch ( FileNotFoundException e) { - throw new StingException("Could not find required input file: " + fileName); - } - - this.fileName = fileName; - - } - public String GetNextLine() { - String line; - int x = reader.read(); - line. - //reader.read() - } - } */ - private class BeagleFileReader { private String headerString; private BufferedReader reader; diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/ProduceBeagleInputWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/ProduceBeagleInputWalker.java index 4f13fa6c5..d83181963 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/ProduceBeagleInputWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/ProduceBeagleInputWalker.java @@ -90,9 +90,6 @@ public class ProduceBeagleInputWalker extends RodWalker { return 0; } - // Get Reference base for this site: will be output to screen, not directly used by Beagle but rather by output analysis tools - // char re = (char)ref.getBase(); - // output marker ID to Beagle input file beagleWriter.print(String.format("%s ",vc_eval.getLocation().toString())); @@ -117,10 +114,21 @@ public class ProduceBeagleInputWalker extends RodWalker { if (genotype.isCalled() && genotype.hasAttribute(VCFGenotypeRecord.GENOTYPE_LIKELIHOODS_KEY)) { String[] glArray = genotype.getAttributeAsString(VCFGenotypeRecord.GENOTYPE_LIKELIHOODS_KEY).split(","); + Double maxLikelihood = -100.0; + ArrayList likeArray = new ArrayList(); + for (String gl : glArray) { - Double d_gl = 100*Math.pow(10, Double.valueOf(gl)); - beagleWriter.print(String.format("%5.2f ",d_gl)); - } + // need to normalize likelihoods to avoid precision loss. In worst case, if all 3 log-likelihoods are too + // small, we could end up with linear likelihoods of form 0.00 0.00 0.00 which will mess up imputation. + Double dg = Double.valueOf(gl); + if (dg> maxLikelihood) + maxLikelihood = dg; + + likeArray.add(dg); + } + + for (Double likeVal: likeArray) + beagleWriter.print(String.format("%5.4f ",Math.pow(10, likeVal-maxLikelihood))); } else beagleWriter.print("0.33 0.33 0.33 "); // write 1/3 likelihoods for uncalled genotypes.