Added a module to test for reference mismatch associations, and a self-normalized/self-normalizing version.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5408 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
chartl 2011-03-08 20:01:28 +00:00
parent 63e1625cc5
commit da88c29b6e
2 changed files with 139 additions and 0 deletions

View File

@ -0,0 +1,91 @@
package org.broadinstitute.sting.oneoffprojects.walkers.association.modules;
import org.broadinstitute.sting.gatk.datasources.sample.Sample;
import org.broadinstitute.sting.oneoffprojects.walkers.association.MapExtender;
import org.broadinstitute.sting.oneoffprojects.walkers.association.RegionalAssociationWalker;
import org.broadinstitute.sting.oneoffprojects.walkers.association.statistics.casecontrol.UStatistic;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import java.util.*;
/**
* Created by IntelliJ IDEA.
* User: chartl
* Date: 3/8/11
* Time: 2:32 PM
* To change this template use File | Settings | File Templates.
*/
public class MismatchRate extends UStatistic {
private Map<Sample,Object> sampleStats = new HashMap<Sample,Object>();
private int currentRefBase = 0;
public void init(RegionalAssociationWalker walker) {
Set<Sample> samples = walker.getSamples();
for ( Sample s : samples ) {
if ( s.hasProperty("mismatch_rate.mean") && s.hasProperty("mismatch_rate.std") ) {
double mn = (Double) s.getProperty("mismatch_rate.mean");
double std = (Double) s.getProperty("mismatch_rate.std");
sampleStats.put(s,new Pair<Double,Double>(mn,std));
} else {
sampleStats.put(s,new MathUtils.RunningAverage());
}
}
}
@Override
public Map<Sample,Object> mapLocus(MapExtender extender) {
currentRefBase = BaseUtils.simpleBaseToBaseIndex(extender.getReferenceContext().getBase());
Map<Sample,ReadBackedPileup> pileups = extender.getReadFilteredPileup();
Map<Sample,Object> maps = new HashMap<Sample,Object>(pileups.size());
for ( Map.Entry<Sample,ReadBackedPileup> samPileup : pileups.entrySet() ) {
maps.put(samPileup.getKey(),map(samPileup.getKey(),samPileup.getValue()));
}
return maps;
}
public Collection<Number> map(Sample sample, ReadBackedPileup pileup) {
Object stats = sampleStats.get(sample);
double mn;
double std;
if ( stats instanceof Pair ) {
mn = ((Pair<Double,Double>)stats).first;
std = ((Pair<Double,Double>)stats).second;
} else {
MathUtils.RunningAverage ra = (MathUtils.RunningAverage) stats;
mn = ra.mean();
std = ra.stddev();
if ( std <= 0.0 ) {
std = 1.0;
}
ra.add(pileup.size());
}
return Arrays.asList((Number) ((calcMMR(pileup) - mn) / std));
}
public double calcMMR(ReadBackedPileup rbp) {
int[] counts = rbp.getBaseCounts();
int total = 0;
int nonref = 0;
for ( int base : new int[]{0,1,2,3} ) {
total += counts[base];
if ( base != currentRefBase ) {
nonref += counts[base];
}
}
return ((double)nonref)/total;
}
// note: this is to satisfy the interface, and is never called due to override
public Collection<Number> map(ReadBackedPileup pileup) { return null; }
public int getWindowSize() { return 100; }
public int slideByValue() { return 10; }
public boolean usePreviouslySeenReads() { return true; }
}

View File

@ -0,0 +1,48 @@
package org.broadinstitute.sting.oneoffprojects.walkers.association.modules;
import org.broadinstitute.sting.gatk.datasources.sample.Sample;
import org.broadinstitute.sting.oneoffprojects.walkers.association.MapExtender;
import org.broadinstitute.sting.oneoffprojects.walkers.association.statistics.casecontrol.ZStatistic;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import java.util.Map;
/**
* Created by IntelliJ IDEA.
* User: chartl
* Date: 3/8/11
* Time: 1:17 PM
* To change this template use File | Settings | File Templates.
*/
public class ReferenceMismatches extends ZStatistic {
final static int[] BASE_INDEX = {0,1,2,3};
int currentRefBase = 0;
@Override
public Map<Sample,Object> mapLocus(MapExtender extender) {
currentRefBase = BaseUtils.simpleBaseToBaseIndex(extender.getReferenceContext().getBase());
return super.mapLocus(extender);
}
public Pair<Number,Number> map(ReadBackedPileup rbp) {
int[] counts = rbp.getBaseCounts();
int total = 0;
int nonref = 0;
for ( int base : BASE_INDEX ) {
total += counts[base];
if ( base != currentRefBase ) {
nonref += counts[base];
}
}
return new Pair<Number,Number>(nonref,total);
}
public int getWindowSize() { return 100; }
public int slideByValue() { return 25; }
public boolean usePreviouslySeenReads() { return true; }
}