From 0cc634ed5d8d3ec2d06c590ace36f7e3b24400e3 Mon Sep 17 00:00:00 2001 From: aaron Date: Fri, 4 Sep 2009 18:40:43 +0000 Subject: [PATCH] -Renamed rodVariants to RodGeliText -Remove KGenomesSNPROD -Remove rodFLT -Renamed rodGFF to RodGenotypeChipAsGFF -Fixed a problem in SSGenotypeCall -Added basic SSGenotype Test class -Make VCFHeader constructors public git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1536 348d0f76-0448-11de-a6fe-93d51630548a --- .../sting/gatk/contexts/VariantContext.java | 8 +- .../sting/gatk/refdata/KGenomesSNPROD.java | 61 ---------- .../gatk/refdata/ReferenceOrderedData.java | 8 +- .../{rodVariants.java => RodGeliText.java} | 4 +- ...{rodGFF.java => RodGenotypeChipAsGFF.java} | 4 +- .../sting/gatk/refdata/rodFLT.java | 66 ---------- .../gatk/walkers/PileupWithContextWalker.java | 7 +- .../gatk/walkers/filters/RatioFilter.java | 7 +- .../walkers/filters/VECAlleleBalance.java | 4 +- .../gatk/walkers/filters/VECFisherStrand.java | 4 +- .../filters/VECOnOffGenotypeRatio.java | 4 +- .../filters/VariantFiltrationWalker.java | 6 +- .../walkers/genotyper/SSGenotypeCall.java | 3 +- .../gatk/walkers/CoverageEvalWalker.java | 4 +- .../sting/utils/genotype/vcf/VCFHeader.java | 4 +- .../walkers/genotyper/SSGenotypeCallTest.java | 115 ++++++++++++++++++ 16 files changed, 145 insertions(+), 164 deletions(-) delete mode 100755 java/src/org/broadinstitute/sting/gatk/refdata/KGenomesSNPROD.java rename java/src/org/broadinstitute/sting/gatk/refdata/{rodVariants.java => RodGeliText.java} (98%) rename java/src/org/broadinstitute/sting/gatk/refdata/{rodGFF.java => RodGenotypeChipAsGFF.java} (97%) delete mode 100755 java/src/org/broadinstitute/sting/gatk/refdata/rodFLT.java create mode 100644 java/test/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCallTest.java diff --git a/java/src/org/broadinstitute/sting/gatk/contexts/VariantContext.java b/java/src/org/broadinstitute/sting/gatk/contexts/VariantContext.java index 63b87a0ee..438f844b7 100755 --- a/java/src/org/broadinstitute/sting/gatk/contexts/VariantContext.java +++ b/java/src/org/broadinstitute/sting/gatk/contexts/VariantContext.java @@ -25,7 +25,7 @@ package org.broadinstitute.sting.gatk.contexts; -import org.broadinstitute.sting.gatk.refdata.rodVariants; +import org.broadinstitute.sting.gatk.refdata.RodGeliText; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import net.sf.samtools.SAMRecord; @@ -36,11 +36,11 @@ public class VariantContext { private RefMetaDataTracker tracker; private ReferenceContext ref; private AlignmentContext context; - private rodVariants variant; + private RodGeliText variant; private AlignmentContext Q0freeContext = null; public VariantContext(RefMetaDataTracker tracker, ReferenceContext ref, - AlignmentContext context, rodVariants variant) { + AlignmentContext context, RodGeliText variant) { this.tracker = tracker; this.ref = ref; this.context = context; @@ -49,7 +49,7 @@ public class VariantContext { public RefMetaDataTracker getTracker() { return tracker; } public ReferenceContext getReferenceContext() { return ref; } - public rodVariants getVariant() { return variant; } + public RodGeliText getVariant() { return variant; } public AlignmentContext getAlignmentContext(boolean useMQ0Reads) { return (useMQ0Reads ? context : getQ0freeContext()); } diff --git a/java/src/org/broadinstitute/sting/gatk/refdata/KGenomesSNPROD.java b/java/src/org/broadinstitute/sting/gatk/refdata/KGenomesSNPROD.java deleted file mode 100755 index 7fe6a39e5..000000000 --- a/java/src/org/broadinstitute/sting/gatk/refdata/KGenomesSNPROD.java +++ /dev/null @@ -1,61 +0,0 @@ -package org.broadinstitute.sting.gatk.refdata; - -import java.util.*; - -import org.broadinstitute.sting.utils.GenomeLoc; -import org.broadinstitute.sting.utils.GenomeLocParser; - -/** - * loc ref alt EM_alt_freq discovery_likelihood discovery_null discovery_prior discovery_lod EM_N n_ref n_het n_hom - * chr1:1104840 A N 0.000000 -85.341265 -85.341265 0.000000 0.000000 324.000000 162 0 0 - * chr1:1104841 C N 0.000000 -69.937928 -69.937928 0.000000 0.000000 324.000000 162 0 0 - * chr1:1104842 A N 0.000000 -84.816002 -84.816002 0.000000 0.000000 324.000000 162 0 0 - * - */ -public class KGenomesSNPROD extends TabularROD implements SNPCallFromGenotypes { - public KGenomesSNPROD(final String name) { - super(name); - } - - public GenomeLoc getLocation() { - loc = GenomeLocParser.createGenomeLoc(this.get("0"), Long.parseLong(this.get("1"))); - return loc; - } - public String getRefBasesFWD() { return this.get("2"); } - public char getRefSnpFWD() throws IllegalStateException { return getRefBasesFWD().charAt(0); } - public String getAltBasesFWD() { return this.get("3"); } - public char getAltSnpFWD() throws IllegalStateException { - if ( getAltBasesFWD().charAt(0) != getRefSnpFWD() ) - return getAltBasesFWD().charAt(0); - return getAltBasesFWD().charAt(1); - } - public boolean isReference() { return getVariationConfidence() < 0.01; } - public boolean isSNP() { return ! isReference(); } - public boolean isInsertion() { return false; } - public boolean isDeletion() { return false; } - public boolean isIndel() { return false; } - public double getMAF() { return Double.parseDouble(this.get("4")); } - public double getHeterozygosity() { return 2 * getMAF() * (1 - getMAF()); } - public boolean isGenotype() { return false; } - public double getVariationConfidence() { return Double.parseDouble(this.get("8")); } - public double getConsensusConfidence() { return -1; } - public List getGenotype() throws IllegalStateException { throw new IllegalStateException(); } - public int getPloidy() throws IllegalStateException { return 2; } - public boolean isBiallelic() { return true; } - - // SNPCallFromGenotypes interface - public int nIndividuals() { return -1; } - public int nHomRefGenotypes() { return -1; } - public int nHetGenotypes() { return -1; } - public int nHomVarGenotypes() { return -1; } - public List getGenotypes() { return null; } - public int length() { return 1; } - - public String toString() { - StringBuffer sb = new StringBuffer(); - sb.append(getLocation().getContig() + "\t" + getLocation().getStart() + "\t"); - sb.append(getRefSnpFWD() + "\t" + getGenotype().get(0) + "\t"); - sb.append(getMAF() + "\t0\t0\t0\t" + getVariationConfidence()); - return sb.toString(); - } -} \ No newline at end of file diff --git a/java/src/org/broadinstitute/sting/gatk/refdata/ReferenceOrderedData.java b/java/src/org/broadinstitute/sting/gatk/refdata/ReferenceOrderedData.java index 496d5a0db..44650604a 100644 --- a/java/src/org/broadinstitute/sting/gatk/refdata/ReferenceOrderedData.java +++ b/java/src/org/broadinstitute/sting/gatk/refdata/ReferenceOrderedData.java @@ -62,23 +62,21 @@ public class ReferenceOrderedData implements static { // All known ROD types - addModule("GFF", rodGFF.class); + addModule("GFF", RodGenotypeChipAsGFF.class); addModule("dbSNP", rodDbSNP.class); addModule("HapMapAlleleFrequencies", HapMapAlleleFrequenciesROD.class); addModule("SAMPileup", rodSAMPileup.class); addModule("GELI", rodGELI.class); - addModule("FLT", rodFLT.class); addModule("RefSeq", rodRefSeq.class); addModule("Table", TabularROD.class); addModule("PooledEM", PooledEMSNPROD.class); - addModule("1KGSNPs", KGenomesSNPROD.class); addModule("CleanedOutSNP", CleanedOutSNPROD.class); addModule("SangerSNP", SangerSNPROD.class); addModule("SimpleIndel", SimpleIndelROD.class); addModule("PointIndel", PointIndelROD.class); addModule("HapMapGenotype", HapMapGenotypeROD.class); addModule("Intervals", IntervalRod.class); - addModule("Variants", rodVariants.class); + addModule("Variants", RodGeliText.class); addModule("GLF", RodGLF.class); } @@ -232,7 +230,7 @@ public class ReferenceOrderedData implements for (ReferenceOrderedDatum rec : this) { System.out.println(rec.toString()); - rodGFF gff = (rodGFF) rec; + RodGenotypeChipAsGFF gff = (RodGenotypeChipAsGFF) rec; String[] keys = {"LENGTH", "ALT", "FOBARBAR"}; for (String key : keys) { System.out.printf(" -> %s is (%s)%n", key, gff.containsAttribute(key) ? gff.getAttribute(key) : "none"); diff --git a/java/src/org/broadinstitute/sting/gatk/refdata/rodVariants.java b/java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java similarity index 98% rename from java/src/org/broadinstitute/sting/gatk/refdata/rodVariants.java rename to java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java index 0710046d1..b83d510a3 100755 --- a/java/src/org/broadinstitute/sting/gatk/refdata/rodVariants.java +++ b/java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java @@ -32,7 +32,7 @@ import java.io.IOException; import java.util.List; import java.util.Arrays; -public class rodVariants extends BasicReferenceOrderedDatum implements AllelicVariant { +public class RodGeliText extends BasicReferenceOrderedDatum implements AllelicVariant { public enum Genotype { AA, AC, AG, AT, CC, CG, CT, GG, GT, TT } public GenomeLoc loc; public char refBase = 'N'; @@ -43,7 +43,7 @@ public class rodVariants extends BasicReferenceOrderedDatum implements AllelicVa public double lodBtnb; public double[] genotypeLikelihoods = new double[10]; - public rodVariants(final String name) { + public RodGeliText(final String name) { super(name); } diff --git a/java/src/org/broadinstitute/sting/gatk/refdata/rodGFF.java b/java/src/org/broadinstitute/sting/gatk/refdata/RodGenotypeChipAsGFF.java similarity index 97% rename from java/src/org/broadinstitute/sting/gatk/refdata/rodGFF.java rename to java/src/org/broadinstitute/sting/gatk/refdata/RodGenotypeChipAsGFF.java index 5e8ecfacd..3f8107700 100644 --- a/java/src/org/broadinstitute/sting/gatk/refdata/rodGFF.java +++ b/java/src/org/broadinstitute/sting/gatk/refdata/RodGenotypeChipAsGFF.java @@ -16,7 +16,7 @@ import java.util.regex.Pattern; * Time: 10:47:14 AM * To change this template use File | Settings | File Templates. */ -public class rodGFF extends BasicReferenceOrderedDatum implements AllelicVariant { +public class RodGenotypeChipAsGFF extends BasicReferenceOrderedDatum implements AllelicVariant { private String contig, source, feature, strand, frame; private long start, stop; private double score; @@ -27,7 +27,7 @@ public class rodGFF extends BasicReferenceOrderedDatum implements AllelicVariant // Constructors // // ---------------------------------------------------------------------- - public rodGFF(final String name) { + public RodGenotypeChipAsGFF(final String name) { super(name); } diff --git a/java/src/org/broadinstitute/sting/gatk/refdata/rodFLT.java b/java/src/org/broadinstitute/sting/gatk/refdata/rodFLT.java deleted file mode 100755 index 687fd04b8..000000000 --- a/java/src/org/broadinstitute/sting/gatk/refdata/rodFLT.java +++ /dev/null @@ -1,66 +0,0 @@ -package org.broadinstitute.sting.gatk.refdata; - -import java.util.*; - -import org.broadinstitute.sting.utils.*; - -/** - * 1 469 C S 51 52 0.00 44 1 C 0 M - * 1 492 C Y 130 51 0.00 46 1 C 0 T - * 1 20786 G K 5 6 0.00 46 1 T 38 G - * - */ -public class rodFLT extends TabularROD implements SNPCallFromGenotypes { - public rodFLT(final String name) { - super(name); - } - - public GenomeLoc getLocation() { - loc = GenomeLocParser.createGenomeLoc(this.get("0"), Long.parseLong(this.get("1"))); - return loc; - } - public String getRefBasesFWD() { return this.get("2"); } - public char getRefSnpFWD() throws IllegalStateException { return getRefBasesFWD().charAt(0); } - public String getAltBasesFWD() { return new String(BaseUtils.iupacToBases(this.get("3").charAt(0))); } - public char getAltSnpFWD() throws IllegalStateException { - char[] bases = BaseUtils.iupacToBases(this.get("3").charAt(0)); - if ( bases[0] != getRefSnpFWD() ) - return bases[0]; - else - return bases[1]; - } - - public String toString() { - StringBuffer sb = new StringBuffer(); - sb.append(loc.getContig() + "\t" + loc.getStart() + "\t"); - sb.append(getRefSnpFWD() + "\t-1\t-1\t" + getAltBasesFWD()); - for (int i=0; i < 12; i++) - sb.append("\t0"); - return sb.toString(); - } - - public boolean isReference() { return false; } - public boolean isSNP() { return true; } - public boolean isInsertion() { return false; } - public boolean isDeletion() { return false; } - public boolean isIndel() { return false; } - public double getMAF() { return 0.0; } - public double getHeterozygosity() { return 0.0; } - public boolean isGenotype() { return false; } - public double getVariationConfidence() { - double value = Double.parseDouble(this.get("6")); - return (value > 0.0 ? value : 1.0); - } - public double getConsensusConfidence() { return -1; } - public List getGenotype() throws IllegalStateException { return Arrays.asList(getAltBasesFWD()); } - public int getPloidy() throws IllegalStateException { return 2; } - public boolean isBiallelic() { return true; } - public int length() { return 1; } - - // SNPCallFromGenotypes interface - public int nIndividuals() { return -1; } - public int nHomRefGenotypes() { return -1; } - public int nHetGenotypes() { return -1; } - public int nHomVarGenotypes() { return -1; } - public List getGenotypes() { return null; } -} \ No newline at end of file diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/PileupWithContextWalker.java b/java/src/org/broadinstitute/sting/gatk/walkers/PileupWithContextWalker.java index 9088b58da..33182c032 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/PileupWithContextWalker.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/PileupWithContextWalker.java @@ -2,10 +2,7 @@ package org.broadinstitute.sting.gatk.walkers; import org.broadinstitute.sting.gatk.contexts.AlignmentContext; import org.broadinstitute.sting.gatk.contexts.ReferenceContext; -import org.broadinstitute.sting.gatk.refdata.ReferenceOrderedDatum; -import org.broadinstitute.sting.gatk.refdata.rodDbSNP; -import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; -import org.broadinstitute.sting.gatk.refdata.rodGFF; +import org.broadinstitute.sting.gatk.refdata.*; import org.broadinstitute.sting.utils.cmdLine.Argument; import org.broadinstitute.sting.utils.ReadBackedPileup; import org.broadinstitute.sting.utils.Utils; @@ -85,7 +82,7 @@ public class PileupWithContextWalker extends LocusWalker imple } */ - if ( datum != null && (datum instanceof rodGFF)) { + if ( datum != null && (datum instanceof RodGenotypeChipAsGFF)) { rodString += "Hapmap: " + datum.toSimpleString() + " "; } } diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java index 4e54caa0a..470b069ed 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/RatioFilter.java @@ -1,10 +1,9 @@ package org.broadinstitute.sting.gatk.walkers.filters; import org.broadinstitute.sting.gatk.contexts.VariantContext; -import org.broadinstitute.sting.gatk.refdata.rodVariants; +import org.broadinstitute.sting.gatk.refdata.RodGeliText; import org.broadinstitute.sting.utils.*; import org.apache.log4j.Logger; -import cern.jet.math.Arithmetic; public abstract class RatioFilter implements VariantExclusionCriterion { @@ -37,14 +36,14 @@ public abstract class RatioFilter implements VariantExclusionCriterion { highThreshold = threshold; } - protected abstract Pair scoreVariant(char ref, ReadBackedPileup pileup, rodVariants variant); + protected abstract Pair scoreVariant(char ref, ReadBackedPileup pileup, RodGeliText variant); protected abstract boolean excludeHetsOnly(); public boolean useZeroQualityReads() { return false; } public void compute(VariantContextWindow contextWindow) { VariantContext context = contextWindow.getContext(); - rodVariants variant = context.getVariant(); + RodGeliText variant = context.getVariant(); char ref = context.getReferenceContext().getBase(); ReadBackedPileup pileup = new ReadBackedPileup(ref, context.getAlignmentContext(useZeroQualityReads())); diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java index 8966abba7..b008b3715 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECAlleleBalance.java @@ -1,6 +1,6 @@ package org.broadinstitute.sting.gatk.walkers.filters; -import org.broadinstitute.sting.gatk.refdata.rodVariants; +import org.broadinstitute.sting.gatk.refdata.RodGeliText; import org.broadinstitute.sting.utils.ReadBackedPileup; import org.broadinstitute.sting.utils.Utils; import org.broadinstitute.sting.utils.Pair; @@ -30,7 +30,7 @@ public class VECAlleleBalance extends RatioFilter { * Return the count of bases matching the major (first) and minor (second) alleles as a pair. * */ - protected Pair scoreVariant(char ref, ReadBackedPileup pileup, rodVariants variant) { + protected Pair scoreVariant(char ref, ReadBackedPileup pileup, RodGeliText variant) { final String genotype = variant.getBestGenotype(); final String bases = pileup.getBases(); diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECFisherStrand.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECFisherStrand.java index 7d7fbfc99..7bdcd64d9 100644 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECFisherStrand.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECFisherStrand.java @@ -2,7 +2,7 @@ package org.broadinstitute.sting.gatk.walkers.filters; import org.broadinstitute.sting.gatk.contexts.AlignmentContext; import org.broadinstitute.sting.gatk.contexts.VariantContext; -import org.broadinstitute.sting.gatk.refdata.rodVariants; +import org.broadinstitute.sting.gatk.refdata.RodGeliText; import org.broadinstitute.sting.utils.BaseUtils; import net.sf.samtools.SAMRecord; import cern.jet.math.Arithmetic; @@ -23,7 +23,7 @@ public class VECFisherStrand implements VariantExclusionCriterion { public void compute(VariantContextWindow contextWindow) { VariantContext context = contextWindow.getContext(); - rodVariants variant = context.getVariant(); + RodGeliText variant = context.getVariant(); int allele1 = BaseUtils.simpleBaseToBaseIndex(variant.getBestGenotype().charAt(0)); int allele2 = BaseUtils.simpleBaseToBaseIndex(variant.getBestGenotype().charAt(1)); diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java index 594f4d3a0..e25285aa2 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VECOnOffGenotypeRatio.java @@ -1,6 +1,6 @@ package org.broadinstitute.sting.gatk.walkers.filters; -import org.broadinstitute.sting.gatk.refdata.rodVariants; +import org.broadinstitute.sting.gatk.refdata.RodGeliText; import org.broadinstitute.sting.utils.*; public class VECOnOffGenotypeRatio extends RatioFilter { @@ -26,7 +26,7 @@ public class VECOnOffGenotypeRatio extends RatioFilter { * best genotype). On are in the first field, off in the second. * */ - protected Pair scoreVariant(char ref, ReadBackedPileup pileup, rodVariants variant) { + protected Pair scoreVariant(char ref, ReadBackedPileup pileup, RodGeliText variant) { final String genotype = variant.getBestGenotype().toUpperCase(); final String bases = pileup.getBases(); diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VariantFiltrationWalker.java b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VariantFiltrationWalker.java index 274ee48a9..6e579ac92 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/filters/VariantFiltrationWalker.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/filters/VariantFiltrationWalker.java @@ -17,7 +17,7 @@ import java.util.*; * The former modifiesthe likelihoods of each genotype, while the latter makes a decision to include or exclude a * variant outright. At the moment, the variants are expected to be in gelitext format. */ -@Requires(value={DataSource.READS, DataSource.REFERENCE},referenceMetaData=@RMD(name="variant",type=rodVariants.class)) +@Requires(value={DataSource.READS, DataSource.REFERENCE},referenceMetaData=@RMD(name="variant",type= RodGeliText.class)) public class VariantFiltrationWalker extends LocusWalker { @Argument(fullName="variants_out_head", shortName="VOH", doc="File to which modified variants should be written") public String VARIANTS_OUT_HEAD; @Argument(fullName="features", shortName="F", doc="Feature test (optionally with arguments) to apply to genotype posteriors. Syntax: 'testname[:arguments]'", required=false) public String[] FEATURES; @@ -176,7 +176,7 @@ public class VariantFiltrationWalker extends LocusWalker { * @return 1 if the locus was successfully processed, 0 if otherwise */ public Integer map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { - rodVariants variant = (rodVariants) tracker.lookup("variant", null); + RodGeliText variant = (RodGeliText) tracker.lookup("variant", null); // Ignore places where we don't have a variant or where the reference base is ambiguous. if ( variant == null || BaseUtils.simpleBaseToBaseIndex(ref.getBase()) == -1 ) @@ -204,7 +204,7 @@ public class VariantFiltrationWalker extends LocusWalker { VariantContext context = variantContextWindow.getContext(); if ( context == null ) return; - rodVariants variant = context.getVariant(); + RodGeliText variant = context.getVariant(); HashMap exclusionResults = new HashMap(); diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCall.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCall.java index 802fd8895..4ba7c8d9e 100644 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCall.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCall.java @@ -2,7 +2,6 @@ package org.broadinstitute.sting.gatk.walkers.genotyper; import net.sf.samtools.SAMRecord; import org.broadinstitute.sting.utils.GenomeLoc; -import org.broadinstitute.sting.utils.MathUtils; import org.broadinstitute.sting.utils.ReadBackedPileup; import org.broadinstitute.sting.utils.Utils; import org.broadinstitute.sting.utils.genotype.*; @@ -129,7 +128,7 @@ public class SSGenotypeCall implements Genotype, ReadBacked, GenotypesBacked, Li mCompareTo = DiploidGenotype.valueOf(Utils.dupString(this.getReference(),2)); } else { Integer sorted[] = Utils.SortPermutation(mGenotypeLikelihoods.getPosteriors()); - mCompareTo = DiploidGenotype.values()[sorted[DiploidGenotype.values().length - 1]]; + mCompareTo = DiploidGenotype.values()[sorted[DiploidGenotype.values().length - 2]]; } } return mCompareTo; diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/CoverageEvalWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/CoverageEvalWalker.java index 046057a84..3357715f3 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/CoverageEvalWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/CoverageEvalWalker.java @@ -4,7 +4,7 @@ import net.sf.samtools.SAMRecord; import org.broadinstitute.sting.gatk.contexts.AlignmentContext; import org.broadinstitute.sting.gatk.contexts.ReferenceContext; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; -import org.broadinstitute.sting.gatk.refdata.rodGFF; +import org.broadinstitute.sting.gatk.refdata.RodGenotypeChipAsGFF; import org.broadinstitute.sting.gatk.walkers.LocusWalker; import org.broadinstitute.sting.gatk.walkers.genotyper.SSGenotypeCall; import org.broadinstitute.sting.gatk.walkers.genotyper.SingleSampleGenotyper; @@ -52,7 +52,7 @@ public class CoverageEvalWalker extends LocusWalker, String> { public List map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { - rodGFF hapmap_chip = (rodGFF)tracker.lookup("hapmap-chip", null); + RodGenotypeChipAsGFF hapmap_chip = (RodGenotypeChipAsGFF)tracker.lookup("hapmap-chip", null); String hc_genotype; if (hapmap_chip != null) { diff --git a/java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFHeader.java b/java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFHeader.java index 7cc3392d2..34af9b685 100644 --- a/java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFHeader.java +++ b/java/src/org/broadinstitute/sting/utils/genotype/vcf/VCFHeader.java @@ -45,7 +45,7 @@ public class VCFHeader { * * @param metaData the meta data associated with this header */ - protected VCFHeader(Map metaData) { + public VCFHeader(Map metaData) { for (String key : metaData.keySet()) mMetaData.put(key, metaData.get(key)); checkVCFVersion(); } @@ -56,7 +56,7 @@ public class VCFHeader { * @param metaData the meta data associated with this header * @param genotypeSampleNames the genotype format field, and the sample names */ - protected VCFHeader(Map metaData, List genotypeSampleNames) { + public VCFHeader(Map metaData, List genotypeSampleNames) { for (String key : metaData.keySet()) mMetaData.put(key, metaData.get(key)); for (String col : genotypeSampleNames) { if (!col.equals("FORMAT")) diff --git a/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCallTest.java b/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCallTest.java new file mode 100644 index 000000000..0b651306a --- /dev/null +++ b/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/SSGenotypeCallTest.java @@ -0,0 +1,115 @@ +package org.broadinstitute.sting.gatk.walkers.genotyper; + +import net.sf.samtools.SAMFileHeader; +import net.sf.samtools.SAMRecord; +import org.broadinstitute.sting.BaseTest; +import org.broadinstitute.sting.utils.GenomeLoc; +import org.broadinstitute.sting.utils.GenomeLocParser; +import org.broadinstitute.sting.utils.Pair; +import org.broadinstitute.sting.utils.ReadBackedPileup; +import org.broadinstitute.sting.utils.sam.ArtificialSAMUtils; +import org.junit.Assert; +import org.junit.Test; + +import java.util.ArrayList; +import java.util.List; + + +/** + * @author aaron + *

+ * Class SSGenotypeCallTest + *

+ * test the SS Genotype call class + */ +public class SSGenotypeCallTest extends BaseTest { + + // we need a fake GenotypeLikelihoods class + public class GenotypeLikelihoodsImpl extends GenotypeLikelihoods { + + /** + * Must be overridden by concrete subclasses + * + * @param observedBase + * @param chromBase + * @param read + * @param offset + * + * @return + */ + @Override + protected double log10PofTrueBaseGivenMiscall(char observedBase, char chromBase, SAMRecord read, int offset) { + return 0; //To change body of implemented methods use File | Settings | File Templates. + } + + public GenotypeLikelihoodsImpl() { + this.posteriors = new double[10]; + for (int x = 0; x < 10; x++) { + posteriors[x] = -(x); + } + this.likelihoods = new double[10]; + for (int x = 0; x < 10; x++) { + likelihoods[x] = -(x); + } + } + } + + + // make a fake read pile-up + public Pair makePileup() { + List reads = new ArrayList(); + List offset = new ArrayList(); + SAMFileHeader header = ArtificialSAMUtils.createArtificialSamHeader(1, 1, 10000); + for (int x = 0; x < 10; x++) { + reads.add(ArtificialSAMUtils.createArtificialRead(header, "test", 0, 1, 100)); + offset.add(10 - x); + } + GenomeLocParser.setupRefContigOrdering(header.getSequenceDictionary()); + GenomeLoc loc = GenomeLocParser.createGenomeLoc("chr1", 1, 10); + return new Pair(new ReadBackedPileup(loc, 'A', reads, offset), loc); + } + + + @Test + public void testBestVrsRefSame() { + Pair myPair = makePileup(); + SSGenotypeCall call = new SSGenotypeCall(true, myPair.second, 'A', new GenotypeLikelihoodsImpl(), myPair.first); + Assert.assertEquals(0, call.getNegLog10PError(), 0.0000001); + } + + @Test + public void testBestVrsRef2() { + Pair myPair = makePileup(); + SSGenotypeCall call2 = new SSGenotypeCall(true, myPair.second, 'T', new GenotypeLikelihoodsImpl(), myPair.first); + Assert.assertEquals(9, call2.getNegLog10PError(), 0.0000001); + } + + @Test + public void testBestVrsRef3() { + Pair myPair = makePileup(); + SSGenotypeCall call3 = new SSGenotypeCall(true, myPair.second, 'C', new GenotypeLikelihoodsImpl(), myPair.first); + Assert.assertEquals(4, call3.getNegLog10PError(), 0.0000001); + } + + + @Test + public void testBestVrsNextSame() { + Pair myPair = makePileup(); + SSGenotypeCall call = new SSGenotypeCall(false, myPair.second, 'A', new GenotypeLikelihoodsImpl(), myPair.first); + Assert.assertEquals(1, call.getNegLog10PError(), 0.0000001); + } + + @Test + public void testBestVrsNext2() { + Pair myPair = makePileup(); + SSGenotypeCall call2 = new SSGenotypeCall(false, myPair.second, 'A', new GenotypeLikelihoodsImpl(), myPair.first); + Assert.assertEquals(1, call2.getNegLog10PError(), 0.0000001); + } + + @Test + public void testBestVrsNext3() { + Pair myPair = makePileup(); + SSGenotypeCall call3 = new SSGenotypeCall(false, myPair.second, 'C', new GenotypeLikelihoodsImpl(), myPair.first); + Assert.assertEquals(1, call3.getNegLog10PError(), 0.0000001); + } +}