zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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
This commit is contained in:
chartl 2011-02-24 19:55:33 +00:00
parent 292b421113
commit 0f1c1fa26f
5 changed files with 146 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationContext.java
View File

@ -9,8 +9,8 @@ import java.util.*;
* A general context for windowed association
*/
public abstract class AssociationContext<X extends AssociationContextAtom> {
private Class<? extends AssociationContextAtom> clazz;
private List<X> window;
protected Class<? extends AssociationContextAtom> clazz;
protected List<X> window;
public AssociationContext( Class<X> zclaz ) {
window = new ArrayList<X>(getWindowSize());

11
java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/AssociationTestRunner.java
View File

@ -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<String> runTests(AssociationContext context) {
List<String> results = new ArrayList<String>();
@ -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) {

1
java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/interfaces/TStatistic.java
View File

@ -9,4 +9,5 @@ package org.broadinstitute.sting.oneoffprojects.walkers.association.interfaces;
*/
public interface TStatistic {
public abstract double getTStatistic();
public abstract int getDegreesOfFreedom();
}

44
java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQuality.java 100755
View File

@ -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<LowMappingQualityAtom> implements ZStatistic {
public LowMappingQuality(Class<LowMappingQualityAtom> 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<String, Pair<Integer,Integer>> 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;
}
}

91
java/src/org/broadinstitute/sting/oneoffprojects/walkers/association/modules/LowMappingQualityAtom.java 100755
View File

@ -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<Sample, Pair<Integer,Integer>> stratifiedCounts;
public LowMappingQualityAtom(MapExtender e) {
super(e);
stratifiedCounts = new HashMap<Sample,Pair<Integer,Integer>>();
for ( Sample s : e.getContext().keySet() ) {
makeCounts(s,e.getPreviousContext(),e.getContext());
}
}
public void makeCounts( Sample sample, Map<Sample, StratifiedAlignmentContext> prev, Map<Sample,StratifiedAlignmentContext> cur ) {
int low = 0;
int tot = 0;
HashSet<String> rnames = new HashSet<String>(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<Integer,Integer>(low,tot));
}
public Map<Sample,Pair<Integer,Integer>> getCounts() {
return stratifiedCounts;
}
public Map<String,Pair<Integer,Integer>> getCaseControlCounts() {
Pair<Integer,Integer> cases = new Pair<Integer,Integer>(0,0);
Pair<Integer,Integer> controls = new Pair<Integer,Integer>(0,0);
for ( Map.Entry<Sample,Pair<Integer,Integer>> 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<String,Pair<Integer,Integer>> caseControlCounts= new HashMap<String,Pair<Integer,Integer>>(2);
caseControlCounts.put("case",cases);
caseControlCounts.put("control",controls);
return caseControlCounts;
}
}