diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java index e16419580..4e63cd29b 100644 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/MultiSampleCaller.java @@ -8,7 +8,6 @@ import org.broadinstitute.sting.gatk.GenomeAnalysisEngine; import org.broadinstitute.sting.gatk.LocusContext; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.walkers.LocusWalker; -import org.broadinstitute.sting.playground.utils.AlleleFrequencyEstimate; import org.broadinstitute.sting.playground.utils.*; import org.broadinstitute.sting.utils.*; import org.broadinstitute.sting.utils.ReadBackedPileup; @@ -167,7 +166,7 @@ public class MultiSampleCaller extends LocusWalker { +public class SingleSampleGenotyper extends LocusWalker { // Control output settings @Argument(fullName = "variants_out", shortName = "varout", doc = "File to which variants should be written", required = true) public File VARIANTS_FILE; @Argument(fullName = "metrics_out", shortName = "metout", doc = "File to which metrics should be written", required = false) public File METRICS_FILE = new File("/dev/stderr"); @@ -131,14 +133,14 @@ public class SingleSampleGenotyper extends LocusWalker onNextBestBasePriors = new HashMap(); @@ -287,7 +295,7 @@ public class GenotypeLikelihoods implements GenotypeGenerator { return s; } - public void ApplyPrior(char ref, double[] allele_likelihoods) { + public void applyPrior(char ref, double[] allele_likelihoods) { int k = 0; for (int i = 0; i < 4; i++) { for (int j = i; j < 4; j++) { @@ -302,7 +310,7 @@ public class GenotypeLikelihoods implements GenotypeGenerator { this.sort(); } - public void ApplyPrior(char ref) { + public void applyPrior(char ref) { for (int i = 0; i < genotypes.length; i++) { if ((genotypes[i].charAt(0) == ref) && (genotypes[i].charAt(1) == ref)) { // hom-ref @@ -420,8 +428,7 @@ public class GenotypeLikelihoods implements GenotypeGenerator { */ @Override public GenotypeCall callGenotypes(RefMetaDataTracker tracker, char ref, ReadBackedPileup pileup) { - //filterQ0Bases(!keepQ0Bases); // Set the filtering / keeping flag - + filterQ0Bases(!keepQ0Bases); // Set the filtering / keeping flag // for calculating the rms of the mapping qualities double squared = 0.0; @@ -436,13 +443,8 @@ public class GenotypeLikelihoods implements GenotypeGenerator { // save off the likelihoods if (likelihoods == null || likelihoods.length == 0) return null; - double lklihoods[] = new double[likelihoods.length]; - - System.arraycopy(likelihoods, 0, lklihoods, 0, likelihoods.length); - - - ApplyPrior(ref); - + // Apply the two calculations + applyPrior(ref); applySecondBaseDistributionPrior(pileup.getBases(), pileup.getSecondaryBasePileup()); // lets setup the locus @@ -450,11 +452,7 @@ public class GenotypeLikelihoods implements GenotypeGenerator { for (int x = 0; x < this.likelihoods.length; x++) { lst.add(new BasicGenotype(pileup.getLocation(),this.genotypes[x],new BayesianConfidenceScore(this.likelihoods[x]))); } - return new SSGGenotypeCall(ref,2,pileup.getLocation(),lst,likelihoods,pileup); - } - //TODO: add this to the above code - boolean discovery = false; - public void setDiscovery(boolean isInDiscoveryMode) { - discovery = isInDiscoveryMode; + return new SSGGenotypeCall(discoveryMode,ref,2,lst,likelihoods,pileup); } + } diff --git a/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java b/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java index d3519c573..05b0cf994 100644 --- a/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java +++ b/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java @@ -22,32 +22,55 @@ import java.util.TreeMap; *

* Class GenotypeCallImpl *

- * The single sample genotypers implementation of the genotype call, plus - * some extra booty for the genotype writers of this world + * The single sample genotypers implementation of the genotype call, which contains + * extra information for the various genotype outputs */ public class SSGGenotypeCall implements GenotypeCall { - // TODO: make SSG into a more robust Genotype call interface // our stored genotype locus private final String mRefBase; private final int mPloidy; - private final GenomeLoc mLoc; private TreeMap mGenotypes = new TreeMap(); private double mLikelihoods[]; private double bestNext = 0; private double bestRef = 0; - private int readDepth; - private double rmsMapping; + private int readDepth = 0; + private double rmsMapping = 0; - public SSGGenotypeCall(char mRefBase, int mPloidy, GenomeLoc mLoc, List genotypes, double likelihoods[], ReadBackedPileup pileup) { + /** + * Generate a single sample genotype object, containing everything we need to represent calls out of a genotyper object + * + * @param discoveryMode are we in discovery mode or genotyping mode + * @param mRefBase the ref base + * @param mPloidy the ploidy + * @param genotypes the genotype + * @param likelihoods the likelihood + * @param pileup the pile-up of reads at the specified locus + */ + public SSGGenotypeCall(boolean discoveryMode, char mRefBase, int mPloidy, List genotypes, double likelihoods[], ReadBackedPileup pileup) { this.mRefBase = String.valueOf(mRefBase).toUpperCase(); this.mPloidy = mPloidy; - this.mLoc = mLoc; if (genotypes.size() < 1) throw new IllegalArgumentException("Genotypes list size must be greater than 0"); - int index = 0; String refStr = Utils.dupString(mRefBase, mPloidy).toUpperCase(); + calculateBestNext(discoveryMode, genotypes, likelihoods, refStr); + + this.readDepth = pileup.getReads().size(); + rmsMapping = Math.sqrt(calculateRMS(pileup) / readDepth); + } + + + /** + * calculate the best and next best, and update the stored confidence scores based on these values + * + * @param discoveryMode are we in discovery mode? + * @param genotypes the genotypes + * @param likelihoods the likelihoods corresponding to the genotypes + * @param refStr the reference string, i.e. reb base[ploidy] + */ + private void calculateBestNext(boolean discoveryMode, List genotypes, double[] likelihoods, String refStr) { + int index = 0; double ref = 0.0; - double best = Double.NEGATIVE_INFINITY; // plus one + double best = Double.NEGATIVE_INFINITY; double next = Double.NEGATIVE_INFINITY; for (Genotype g : genotypes) { if (g.getBases().toUpperCase().equals(refStr)) ref = likelihoods[index]; @@ -61,16 +84,23 @@ public class SSGGenotypeCall implements GenotypeCall { bestRef = Math.abs(best - ref); mLikelihoods = likelihoods; index = 0; - for (Genotype g : genotypes) { - ((BasicGenotype)g).setConfidenceScore( new BayesianConfidenceScore(Math.abs(likelihoods[index] - ref))); - mGenotypes.put(likelihoods[index],g); - index++; - } - this.readDepth = pileup.getReads().size(); - rmsMapping = Math.sqrt(calculateRMS(pileup) / readDepth); + // reset the confidence based on either the discovery mode or the genotype mode + for (Genotype g : genotypes) { + double val = (discoveryMode) ? Math.abs(likelihoods[index] - best) : Math.abs(likelihoods[index] - ref); + ((BasicGenotype) g).setConfidenceScore(new BayesianConfidenceScore(val)); + mGenotypes.put(likelihoods[index], g); + index++; + } } + /** + * calculate the rms , given the read pileup + * + * @param pileup + * + * @return + */ private double calculateRMS(ReadBackedPileup pileup) { double rms = 0.0; for (SAMRecord r : pileup.getReads()) { @@ -174,6 +204,7 @@ public class SSGGenotypeCall implements GenotypeCall { * given the reference, are we a variant? (non-ref) * * @param ref the reference base or bases + * * @return true if we're a variant */ @Override diff --git a/java/src/org/broadinstitute/sting/utils/genotype/glf/GLFWriter.java b/java/src/org/broadinstitute/sting/utils/genotype/glf/GLFWriter.java index 7af77998a..e38e39432 100755 --- a/java/src/org/broadinstitute/sting/utils/genotype/glf/GLFWriter.java +++ b/java/src/org/broadinstitute/sting/utils/genotype/glf/GLFWriter.java @@ -111,7 +111,7 @@ public class GLFWriter implements GenotypeWriter { public void addGenotypeCall(GenotypeCall locus) { SSGGenotypeCall call = (SSGGenotypeCall)locus; LikelihoodObject obj = new LikelihoodObject(call.getLikelihoods(), LikelihoodObject.LIKELIHOOD_TYPE.LOG); - // TODO: fix me aaron + obj.setLikelihoodType(LikelihoodObject.LIKELIHOOD_TYPE.NEGITIVE_LOG); // transform! ... to negitive log likelihoods this.addGenotypeCall(GenomeLocParser.getContigInfo(locus.getLocation().getContig()),(int)locus.getLocation().getStart(),(float)0.0,locus.getReferencebase(),0,obj); }