From 20bffe0430a6fe79b8e9cc558bd3a55e25c91f37 Mon Sep 17 00:00:00 2001 From: Laurent Francioli Date: Wed, 30 Nov 2011 14:46:38 +0100 Subject: [PATCH] Adapted for the new version of MendelianViolation --- .../TransmissionDisequilibriumTest.java | 57 ++++++++----------- 1 file changed, 23 insertions(+), 34 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/TransmissionDisequilibriumTest.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/TransmissionDisequilibriumTest.java index 3de179365..6cc8923e8 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/TransmissionDisequilibriumTest.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/TransmissionDisequilibriumTest.java @@ -12,10 +12,8 @@ import org.broadinstitute.sting.utils.MendelianViolation; import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineType; import org.broadinstitute.sting.utils.codecs.vcf.VCFInfoHeaderLine; import org.broadinstitute.sting.utils.exceptions.UserException; -import org.broadinstitute.sting.utils.text.XReadLines; import org.broadinstitute.sting.utils.variantcontext.VariantContext; -import java.io.FileNotFoundException; import java.util.*; /** @@ -26,42 +24,33 @@ import java.util.*; public class TransmissionDisequilibriumTest extends InfoFieldAnnotation implements ExperimentalAnnotation { - private Set fullMVSet = null; + private Set trios = null; private final static int REF = 0; private final static int HET = 1; private final static int HOM = 2; public Map annotate(RefMetaDataTracker tracker, AnnotatorCompatibleWalker walker, ReferenceContext ref, Map stratifiedContexts, VariantContext vc) { - if ( fullMVSet == null ) { - fullMVSet = new HashSet(); - + if ( trios == null ) { if ( walker instanceof VariantAnnotator ) { - final Map> families = ((VariantAnnotator) walker).getSampleDB().getFamilies(); - for( final Set family : families.values() ) { - for( final Sample sample : family ) { - if( sample.getParents().size() == 2 && family.containsAll(sample.getParents()) ) { // only works with trios for now - fullMVSet.add( new MendelianViolation(sample, 0.0) ); - } - } - } + trios = ((VariantAnnotator) walker).getSampleDB().getChildrenWithParents(); } else { throw new UserException("Transmission disequilibrium test annotation can only be used from the Variant Annotator and requires a valid ped file be passed in."); } } final Map toRet = new HashMap(1); - final HashSet mvsToTest = new HashSet(); + final HashSet triosToTest = new HashSet(); - for( final MendelianViolation mv : fullMVSet ) { - final boolean hasAppropriateGenotypes = vc.hasGenotype(mv.getSampleChild()) && vc.getGenotype(mv.getSampleChild()).hasLikelihoods() && - vc.hasGenotype(mv.getSampleDad()) && vc.getGenotype(mv.getSampleDad()).hasLikelihoods() && - vc.hasGenotype(mv.getSampleMom()) && vc.getGenotype(mv.getSampleMom()).hasLikelihoods(); + for( final Sample child : trios) { + final boolean hasAppropriateGenotypes = vc.hasGenotype(child.getID()) && vc.getGenotype(child.getID()).hasLikelihoods() && + vc.hasGenotype(child.getPaternalID()) && vc.getGenotype(child.getPaternalID()).hasLikelihoods() && + vc.hasGenotype(child.getMaternalID()) && vc.getGenotype(child.getMaternalID()).hasLikelihoods(); if ( hasAppropriateGenotypes ) { - mvsToTest.add(mv); + triosToTest.add(child); } } - toRet.put("TDT", calculateTDT( vc, mvsToTest )); + toRet.put("TDT", calculateTDT( vc, triosToTest )); return toRet; } @@ -72,27 +61,27 @@ public class TransmissionDisequilibriumTest extends InfoFieldAnnotation implemen public List getDescriptions() { return Arrays.asList(new VCFInfoHeaderLine("TDT", 1, VCFHeaderLineType.Float, "Test statistic from Wittkowski transmission disequilibrium test.")); } // Following derivation in http://en.wikipedia.org/wiki/Transmission_disequilibrium_test#A_modified_version_of_the_TDT - private double calculateTDT( final VariantContext vc, final Set mvsToTest ) { + private double calculateTDT( final VariantContext vc, final Set triosToTest ) { - final double nABGivenABandBB = calculateNChildren(vc, mvsToTest, HET, HET, HOM); - final double nBBGivenABandBB = calculateNChildren(vc, mvsToTest, HOM, HET, HOM); - final double nAAGivenABandAB = calculateNChildren(vc, mvsToTest, REF, HET, HET); - final double nBBGivenABandAB = calculateNChildren(vc, mvsToTest, HOM, HET, HET); - final double nAAGivenAAandAB = calculateNChildren(vc, mvsToTest, REF, REF, HET); - final double nABGivenAAandAB = calculateNChildren(vc, mvsToTest, HET, REF, HET); + final double nABGivenABandBB = calculateNChildren(vc, triosToTest, HET, HET, HOM); + final double nBBGivenABandBB = calculateNChildren(vc, triosToTest, HOM, HET, HOM); + final double nAAGivenABandAB = calculateNChildren(vc, triosToTest, REF, HET, HET); + final double nBBGivenABandAB = calculateNChildren(vc, triosToTest, HOM, HET, HET); + final double nAAGivenAAandAB = calculateNChildren(vc, triosToTest, REF, REF, HET); + final double nABGivenAAandAB = calculateNChildren(vc, triosToTest, HET, REF, HET); final double numer = (nABGivenABandBB - nBBGivenABandBB) + 2.0 * (nAAGivenABandAB - nBBGivenABandAB) + (nAAGivenAAandAB - nABGivenAAandAB); final double denom = (nABGivenABandBB + nBBGivenABandBB) + 4.0 * (nAAGivenABandAB + nBBGivenABandAB) + (nAAGivenAAandAB + nABGivenAAandAB); return (numer * numer) / denom; } - private double calculateNChildren( final VariantContext vc, final Set mvsToTest, final int childIdx, final int momIdx, final int dadIdx ) { - final double likelihoodVector[] = new double[mvsToTest.size() * 2]; + private double calculateNChildren( final VariantContext vc, final Set triosToTest, final int childIdx, final int momIdx, final int dadIdx ) { + final double likelihoodVector[] = new double[triosToTest.size() * 2]; int iii = 0; - for( final MendelianViolation mv : mvsToTest ) { - final double[] momGL = vc.getGenotype(mv.getSampleMom()).getLikelihoods().getAsVector(); - final double[] dadGL = vc.getGenotype(mv.getSampleDad()).getLikelihoods().getAsVector(); - final double[] childGL = vc.getGenotype(mv.getSampleChild()).getLikelihoods().getAsVector(); + for( final Sample child : triosToTest ) { + final double[] momGL = vc.getGenotype(child.getMaternalID()).getLikelihoods().getAsVector(); + final double[] dadGL = vc.getGenotype(child.getPaternalID()).getLikelihoods().getAsVector(); + final double[] childGL = vc.getGenotype(child.getID()).getLikelihoods().getAsVector(); likelihoodVector[iii++] = momGL[momIdx] + dadGL[dadIdx] + childGL[childIdx]; likelihoodVector[iii++] = momGL[dadIdx] + dadGL[momIdx] + childGL[childIdx]; }