diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/MendelianViolationEvaluator.java b/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/MendelianViolationEvaluator.java index 6d73310c4..6f2f89fcf 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/MendelianViolationEvaluator.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/MendelianViolationEvaluator.java @@ -9,6 +9,7 @@ import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.walkers.varianteval.VariantEvalWalker; import org.broadinstitute.sting.gatk.walkers.varianteval.tags.Analysis; import org.broadinstitute.sting.gatk.walkers.varianteval.tags.DataPoint; +import org.broadinstitute.sting.utils.MendelianViolation; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; import java.util.Arrays; @@ -47,7 +48,6 @@ public class MendelianViolationEvaluator extends VariantEvaluator { @DataPoint(description = "Number of mendelian variants found") long nVariants; - @DataPoint(description = "Number of mendelian violations found") long nViolations; @@ -60,47 +60,17 @@ public class MendelianViolationEvaluator extends VariantEvaluator { @DataPoint(description = "number of child hom variant calls where the parent was hom ref") long KidHomVar_ParentHomRef; - TrioStructure trio; - double mendelianViolationQualThreshold; + MendelianViolation mv; - private static Pattern FAMILY_PATTERN = Pattern.compile("(.*)\\+(.*)=(.*)"); - - public static class TrioStructure { - public String mom, dad, child; - } - - public static TrioStructure parseTrioDescription(String family) { - Matcher m = FAMILY_PATTERN.matcher(family); - if (m.matches()) { - TrioStructure trio = new TrioStructure(); - //System.out.printf("Found a family pattern: %s%n", parent.FAMILY_STRUCTURE); - trio.mom = m.group(1); - trio.dad = m.group(2); - trio.child = m.group(3); - return trio; - } else { - throw new IllegalArgumentException("Malformatted family structure string: " + family + " required format is mom+dad=child"); - } - } - - // todo: fix public void initialize(VariantEvalWalker walker) { - trio = parseTrioDescription(walker.getFamilyStructure()); - mendelianViolationQualThreshold = walker.getMendelianViolationQualThreshold(); + mv = new MendelianViolation(walker.getFamilyStructure(), walker.getMendelianViolationQualThreshold()); } - public boolean enabled() { //return getVEWalker().FAMILY_STRUCTURE != null; return true; } - private double getQThreshold() { - //return getVEWalker().MENDELIAN_VIOLATION_QUAL_THRESHOLD / 10; // we aren't 10x scaled in the GATK a la phred - return mendelianViolationQualThreshold / 10; // we aren't 10x scaled in the GATK a la phred - //return 0.0; - } - public String getName() { return "mendelian_violations"; } @@ -111,19 +81,14 @@ public class MendelianViolationEvaluator extends VariantEvaluator { public String update1(VariantContext vc, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { if (vc.isBiallelic() && vc.hasGenotypes()) { // todo -- currently limited to biallelic loci - Genotype momG = vc.getGenotype(trio.mom); - Genotype dadG = vc.getGenotype(trio.dad); - Genotype childG = vc.getGenotype(trio.child); - - if (includeGenotype(momG) && includeGenotype(dadG) && includeGenotype(childG)) { + if (mv.setAlleles(vc)) { nVariants++; - if (momG == null || dadG == null || childG == null) - throw new IllegalArgumentException(String.format("VariantContext didn't contain genotypes for expected trio members: mom=%s dad=%s child=%s", trio.mom, trio.dad, trio.child)); + Genotype momG = vc.getGenotype(mv.getSampleMom()); + Genotype dadG = vc.getGenotype(mv.getSampleDad()); + Genotype childG = vc.getGenotype(mv.getSampleChild()); - // all genotypes are good, so let's see if child is a violation - - if (isViolation(vc, momG, dadG, childG)) { + if (mv.isViolation()) { nViolations++; String label; @@ -151,6 +116,42 @@ public class MendelianViolationEvaluator extends VariantEvaluator { return null; // we don't capture any intersting sites } + +/* + private double getQThreshold() { + //return getVEWalker().MENDELIAN_VIOLATION_QUAL_THRESHOLD / 10; // we aren't 10x scaled in the GATK a la phred + return mendelianViolationQualThreshold / 10; // we aren't 10x scaled in the GATK a la phred + //return 0.0; + } + + TrioStructure trio; + double mendelianViolationQualThreshold; + + private static Pattern FAMILY_PATTERN = Pattern.compile("(.*)\\+(.*)=(.*)"); + + public static class TrioStructure { + public String mom, dad, child; + } + + public static TrioStructure parseTrioDescription(String family) { + Matcher m = FAMILY_PATTERN.matcher(family); + if (m.matches()) { + TrioStructure trio = new TrioStructure(); + //System.out.printf("Found a family pattern: %s%n", parent.FAMILY_STRUCTURE); + trio.mom = m.group(1); + trio.dad = m.group(2); + trio.child = m.group(3); + return trio; + } else { + throw new IllegalArgumentException("Malformatted family structure string: " + family + " required format is mom+dad=child"); + } + } + + public void initialize(VariantEvalWalker walker) { + trio = parseTrioDescription(walker.getFamilyStructure()); + mendelianViolationQualThreshold = walker.getMendelianViolationQualThreshold(); + } + private boolean includeGenotype(Genotype g) { return g.getNegLog10PError() > getQThreshold() && g.isCalled(); } @@ -181,4 +182,9 @@ public class MendelianViolationEvaluator extends VariantEvaluator { return true; } + + +*/ + + } diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/SelectVariants.java b/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/SelectVariants.java index 4925a1708..3f6cda8f7 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/SelectVariants.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/SelectVariants.java @@ -70,8 +70,7 @@ public class SelectVariants extends RodWalker { @Argument(fullName="discordance", shortName = "disc", doc="Output variants that were not called from a ROD comparison track. Use -disc ROD_NAME", required=false) private String discordanceRodName = ""; - @Deprecated - @Argument(fullName="family", shortName="fam", doc="USE YAML FILE INSTEAD (-SM) !!! string formatted as dad+mom=child where these parameters determine which sample names are examined", required=false) + @Argument(fullName="family_structure", shortName="family", doc="USE YAML FILE INSTEAD (-SM) !!! string formatted as dad+mom=child where these parameters determine which sample names are examined", required=false) private String FAMILY_STRUCTURE = ""; @Argument(fullName="mendelianViolation", shortName="mv", doc="output mendelian violation sites only. Sample metadata information will be taken from YAML file (passed with -SM)", required=false) @@ -112,8 +111,10 @@ public class SelectVariants extends RodWalker { if (MENDELIAN_VIOLATIONS) mv = new MendelianViolation(getToolkit(), MENDELIAN_VIOLATION_QUAL_THRESHOLD); - else if (!FAMILY_STRUCTURE.isEmpty()) + else if (!FAMILY_STRUCTURE.isEmpty()) { mv = new MendelianViolation(FAMILY_STRUCTURE, MENDELIAN_VIOLATION_QUAL_THRESHOLD); + MENDELIAN_VIOLATIONS = true; + } // Initialize VCF header Set headerLines = VCFUtils.smartMergeHeaders(vcfRods.values(), logger); diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/MendelianViolationClassifier.java b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/MendelianViolationClassifier.java index c82be962c..b2ea35249 100644 --- a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/MendelianViolationClassifier.java +++ b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/MendelianViolationClassifier.java @@ -18,7 +18,6 @@ import org.broadinstitute.sting.gatk.walkers.LocusWalker; import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedArgumentCollection; import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine; import org.broadinstitute.sting.gatk.walkers.genotyper.VariantCallContext; -import org.broadinstitute.sting.gatk.walkers.varianteval.evaluators.MendelianViolationEvaluator; import org.broadinstitute.sting.utils.*; import org.broadinstitute.sting.utils.collections.Pair; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; @@ -58,16 +57,18 @@ public class MendelianViolationClassifier extends LocusWalker homozygousRegions; public HashMap homozygousRegionCounts; public HashMap regionKeys; + public org.broadinstitute.sting.utils.MendelianViolation mvObject; public ExtendedTrioStructure(String family) { - MendelianViolationEvaluator.TrioStructure struct = MendelianViolationEvaluator.parseTrioDescription(family); - this.child = struct.child; - this.mom = struct.mom; - this.dad = struct.dad; + mvObject = new org.broadinstitute.sting.utils.MendelianViolation(family, 0); + this.child = mvObject.getSampleChild(); + this.mom = mvObject.getSampleMom(); + this.dad = mvObject.getSampleDad(); homozygousRegions = new HashMap(3); homozygousRegionCounts = new HashMap(3); homozygousRegions.put(child,null); @@ -414,7 +415,7 @@ public class MendelianViolationClassifier extends LocusWalker