From 0f1c1fa26fede61a098255124ab6148553402900 Mon Sep 17 00:00:00 2001 From: chartl Date: Thu, 24 Feb 2011 19:55:33 +0000 Subject: [PATCH] First general association module. Let the bug fixing begin! git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5307 348d0f76-0448-11de-a6fe-93d51630548a --- .../association/AssociationContext.java | 4 +- .../association/AssociationTestRunner.java | 11 ++- .../association/interfaces/TStatistic.java | 1 + .../modules/LowMappingQuality.java | 44 +++++++++ .../modules/LowMappingQualityAtom.java | 91 +++++++++++++++++++ 5 files changed, 146 insertions(+), 5 deletions(-) create mode 100755 java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQuality.java create mode 100755 java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQualityAtom.java diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationContext.java b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationContext.java index b8b92a9df..95f177fef 100755 --- a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationContext.java +++ b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationContext.java @@ -9,8 +9,8 @@ import java.util.*; * A general context for windowed association */ public abstract class AssociationContext { - private Class clazz; - private List window; + protected Class clazz; + protected List window; public AssociationContext( Class zclaz ) { window = new ArrayList(getWindowSize()); diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationTestRunner.java b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationTestRunner.java index c0c469c25..351ce3cc3 100755 --- a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationTestRunner.java +++ b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationTestRunner.java @@ -2,6 +2,8 @@ package org.broadinstitute.sting.oneoffprojects.walkers.association; import java.util.ArrayList; import java.util.List; + +import cern.jet.random.Normal; import org.broadinstitute.sting.oneoffprojects.walkers.association.interfaces.*; /** @@ -12,6 +14,7 @@ import org.broadinstitute.sting.oneoffprojects.walkers.association.interfaces.*; * To change this template use File | Settings | File Templates. */ public class AssociationTestRunner { + static Normal standardNormal = new Normal(0.0,1.0,null); public static List runTests(AssociationContext context) { List results = new ArrayList(); @@ -20,7 +23,7 @@ public class AssociationTestRunner { } if ( context.getClass().isInstance(ZStatistic.class)) { - results.add(runZ(context)); + results.add(runZ((ZStatistic) context)); } if ( context.getClass().isInstance(FisherExact.class) ) { @@ -34,8 +37,10 @@ public class AssociationTestRunner { return "Test not yet implemented"; } - public static String runZ(AssociationContext context) { - return "Test not yet implemented"; + public static String runZ(ZStatistic context) { + double z = context.getZStatistic(); + double p = z < 0 ? 2*standardNormal.cdf(z) : 2*(1-standardNormal.cdf(z)); + return String.format("Z: %.2f\tP: %.2f",z,p); } public static String runFisherExact(AssociationContext context) { diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/interfaces/TStatistic.java b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/interfaces/TStatistic.java index 613db45f0..a2125004b 100755 --- a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/interfaces/TStatistic.java +++ b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/interfaces/TStatistic.java @@ -9,4 +9,5 @@ package org.broadinstitute.sting.oneoffprojects.walkers.association.interfaces; */ public interface TStatistic { public abstract double getTStatistic(); + public abstract int getDegreesOfFreedom(); } \ No newline at end of file diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQuality.java b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQuality.java new file mode 100755 index 000000000..195475a73 --- /dev/null +++ b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQuality.java @@ -0,0 +1,44 @@ +package org.broadinstitute.sting.oneoffprojects.walkers.association.modules; + +import org.broadinstitute.sting.oneoffprojects.walkers.association.AssociationContext; +import org.broadinstitute.sting.oneoffprojects.walkers.association.AssociationContextAtom; +import org.broadinstitute.sting.oneoffprojects.walkers.association.interfaces.ZStatistic; +import org.broadinstitute.sting.utils.collections.Pair; + +import java.util.Map; + +/** + * Created by IntelliJ IDEA. + * User: chartl + * Date: Feb 24, 2011 + * Time: 1:49:37 PM + * To change this template use File | Settings | File Templates. + */ +public class LowMappingQuality extends AssociationContext implements ZStatistic { + + public LowMappingQuality(Class atomClass ) { + super(atomClass); + } + + public int getWindowSize() { return 20; } + public int slideByValue() { return 5; } + + public double getZStatistic() { + int case_low= 0; + int case_tot = 0; + int ctrl_low = 0; + int ctrl_tot = 0; + for ( LowMappingQualityAtom atom : window ) { + Map> ccc = atom.getCaseControlCounts(); + case_low += ccc.get("case").first; + case_tot += ccc.get("case").second; + ctrl_low += ccc.get("control").first; + ctrl_tot += ccc.get("control").second; + } + double p_case = ( (double) case_low)/case_tot; + double p_ctrl = ((double) ctrl_low)/ctrl_tot; + double p_pool = ( (double) case_low + (double) ctrl_low )/(case_tot+ctrl_tot); + double SE = Math.sqrt( p_pool*(1-p_pool)*( 1/((double)ctrl_tot) + 1/((double)case_tot) )); + return (p_case - p_ctrl)/SE; + } +} diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQualityAtom.java b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQualityAtom.java new file mode 100755 index 000000000..4b1c6db43 --- /dev/null +++ b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQualityAtom.java @@ -0,0 +1,91 @@ +package org.broadinstitute.sting.oneoffprojects.walkers.association.modules; + +import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext; +import org.broadinstitute.sting.gatk.datasources.sample.Sample; +import org.broadinstitute.sting.oneoffprojects.walkers.association.AssociationContextAtom; +import org.broadinstitute.sting.oneoffprojects.walkers.association.MapExtender; +import org.broadinstitute.sting.oneoffprojects.walkers.association.MapHolder; +import org.broadinstitute.sting.utils.collections.Pair; +import org.broadinstitute.sting.utils.pileup.PileupElement; + +import java.util.HashMap; +import java.util.HashSet; +import java.util.Map; + +/** + * Created by IntelliJ IDEA. + * User: chartl + * Date: Feb 24, 2011 + * Time: 1:49:55 PM + * To change this template use File | Settings | File Templates. + */ +public class LowMappingQualityAtom extends AssociationContextAtom { + private static int LOW_THRESH = 5; + private static StratifiedAlignmentContext.StratifiedContextType TYPE = StratifiedAlignmentContext.StratifiedContextType.COMPLETE; + + protected Map> stratifiedCounts; + + public LowMappingQualityAtom(MapExtender e) { + super(e); + stratifiedCounts = new HashMap>(); + for ( Sample s : e.getContext().keySet() ) { + makeCounts(s,e.getPreviousContext(),e.getContext()); + } + } + + public void makeCounts( Sample sample, Map prev, Map cur ) { + int low = 0; + int tot = 0; + HashSet rnames = new HashSet(prev.get(sample).getContext(TYPE).size()); + for ( PileupElement e : prev.get(sample).getContext(TYPE).getBasePileup() ) { + rnames.add(e.getRead().getReadName()); + } + for ( PileupElement e : prev.get(sample).getContext(TYPE).getExtendedEventPileup() ) { + rnames.add(e.getRead().getReadName()); + } + for ( PileupElement e : cur.get(sample).getContext(TYPE).getBasePileup() ) { + if ( ! rnames.contains(e.getRead().getReadName()) ) { + ++tot; + if ( e.getMappingQual() < LOW_THRESH ) { + ++low; + } + } + } + for ( PileupElement e : cur.get(sample).getContext(TYPE).getExtendedEventPileup() ) { + if ( ! rnames.contains(e.getRead().getReadName()) ) { + ++tot; + if ( e.getMappingQual() < LOW_THRESH ) { + ++low; + } + } + } + + stratifiedCounts.put(sample,new Pair(low,tot)); + + } + + public Map> getCounts() { + return stratifiedCounts; + } + + public Map> getCaseControlCounts() { + Pair cases = new Pair(0,0); + Pair controls = new Pair(0,0); + for ( Map.Entry> e : stratifiedCounts.entrySet() ) { + if ( e.getKey().getProperty("cohort").equals("case") ) { + cases.first += e.getValue().first; + cases.second += e.getValue().second; + } + if ( e.getKey().getProperty("cohort").equals("control") ) { + controls.first += e.getValue().first; + controls.second += e.getValue().second; + } + } + + Map> caseControlCounts= new HashMap>(2); + caseControlCounts.put("case",cases); + caseControlCounts.put("control",controls); + return caseControlCounts; + } + +}