Recalibration stuff...

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1810 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
chartl 2009-10-12 20:51:39 +00:00
parent caf689821f
commit ee0afba0af
5 changed files with 24 additions and 51 deletions

View File

@ -6,9 +6,7 @@ import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.QualityUtils;
import java.util.List;
import java.util.LinkedList;
import org.broadinstitute.sting.playground.gatk.walkers.Recalibration.LocalMapType;
import net.sf.samtools.SAMRecord;
@ -21,21 +19,24 @@ import net.sf.samtools.SAMRecord;
*/
public class NQSClusteredZScoreWalker extends LocusWalker<LocalMapType, int[][][]> {
static final int WIN_SIDE_SIZE = 5;
static final int Z_SCORE_MAX = 8;
static final int Z_SCORE_MULTIPLIER = 50; // bins are Z_SCORE * (this) rounded to the nearst int
static final int Z_SCORE_MIN = -6;
static final int Z_SCORE_MAX = 6;
static final int Z_SCORE_MULTIPLIER = 20; // bins are Z_SCORE * (this) rounded to the nearst int
static final int MM_OFFSET = 1;
static final int MATCH_OFFSET = 0;
static final int MAX_Q_SCORE = 2 + QualityUtils.MAX_REASONABLE_Q_SCORE;
protected int WINDOW;
protected int Z_SCORE_RANGE;
public void initialize() {
WINDOW = 2*WIN_SIDE_SIZE+1;
Z_SCORE_RANGE = Z_SCORE_MAX - Z_SCORE_MIN;
}
public int[][][] reduceInit() {
int[][][] q = new int[Z_SCORE_MAX*Z_SCORE_MULTIPLIER+1][MAX_Q_SCORE][2];
for ( int i = 0; i < Z_SCORE_MAX*Z_SCORE_MULTIPLIER+1; i ++ ) {
int[][][] q = new int[Z_SCORE_RANGE*Z_SCORE_MULTIPLIER+1][MAX_Q_SCORE][2];
for ( int i = 0; i < Z_SCORE_RANGE*Z_SCORE_MULTIPLIER+1; i ++ ) {
for ( int j = 0; j < MAX_Q_SCORE; j ++ ) {
for ( int k = 0; k < 2; k ++ ) {
q[i][j][k] = 0;
@ -67,7 +68,7 @@ public class NQSClusteredZScoreWalker extends LocusWalker<LocalMapType, int[][][
public void onTraversalDone( int[][][] zScoreBins ) {
out.print( header() );
for ( int i = 0; i < Z_SCORE_MAX*Z_SCORE_MULTIPLIER; i ++ ) {
for ( int i = 0; i < Z_SCORE_RANGE*Z_SCORE_MULTIPLIER; i ++ ) {
for ( int j = 0; j < MAX_Q_SCORE; j ++ ) {
out.print( formatData(zScoreBins[i][j], i, j) );
}
@ -84,12 +85,16 @@ public class NQSClusteredZScoreWalker extends LocusWalker<LocalMapType, int[][][
public int calcZScoreBin( SAMRecord read, int offset ) {
Pair<Double,Double> meanVar = calcWindowMeanVariance(read, offset);
double rawZ = Math.abs((getQScoreAsInt(read, offset) - meanVar.first)/Math.sqrt(meanVar.second));
double rawZ = (getQScoreAsInt(read, offset) - meanVar.first)/Math.sqrt(meanVar.second);
int zBin;
if ( rawZ < Z_SCORE_MAX ) {
zBin = (int) Math.floor(rawZ*Z_SCORE_MULTIPLIER);
if ( rawZ >= Z_SCORE_MAX ) {
zBin = (int) Math.floor(Z_SCORE_RANGE*Z_SCORE_MULTIPLIER);
} else if ( rawZ <= Z_SCORE_MIN ) {
zBin = 0;
} else if ( rawZ > 0 ) {
zBin = (int) Math.floor(rawZ*Z_SCORE_MULTIPLIER) - (int) Math.floor(Z_SCORE_MIN*Z_SCORE_MULTIPLIER);
} else {
zBin = Z_SCORE_MULTIPLIER*Z_SCORE_MAX;
zBin = (int) Math.floor(-Z_SCORE_MIN*Z_SCORE_MULTIPLIER) + (int) Math.floor(rawZ*Z_SCORE_MULTIPLIER);
}
return zBin;
@ -138,10 +143,6 @@ public class NQSClusteredZScoreWalker extends LocusWalker<LocalMapType, int[][][
public String formatData ( int[] matchMismatch, int zScoreBin, int q ) {
String format = "%f\t%d\t%f\t%f\t%d\t%d%n";
for ( int i = 0; i < MAX_Q_SCORE; i ++ ) {
}
double zScore = ( (double) zScoreBin )/Z_SCORE_MULTIPLIER;
int counts = matchMismatch[MATCH_OFFSET] + matchMismatch[MM_OFFSET];
double empMMR = (((double)matchMismatch[MM_OFFSET])/counts);

View File

@ -6,9 +6,7 @@ import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.QualityUtils;
import java.util.List;
import java.util.LinkedList;
import org.broadinstitute.sting.playground.gatk.walkers.Recalibration.LocalMapType;
import net.sf.samtools.SAMRecord;
@ -19,11 +17,11 @@ import net.sf.samtools.SAMRecord;
* Time: 10:54:18 AM
* To change this template use File | Settings | File Templates.
*/
public class NQSCovariantByCountsWalker extends LocusWalker< LocalMapType, int[][][][] > {
public class NQSCovariantByCountsWalker extends LocusWalker<LocalMapType, int[][][][] > {
final static int SMALL_DEVIATION = 5;
final static int LARGE_DEVIATION = 8;
final static int WIN_SIZE_SIDE = 4;
final static int WIN_SIZE_SIDE = 6;
final static int MATCH_OFFSET = 1;
final static int MM_OFFSET = 0;
final static int MAX_Q_SCORE = 2 + QualityUtils.MAX_REASONABLE_Q_SCORE;

View File

@ -5,10 +5,7 @@ import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.utils.QualityUtils;
import java.util.List;
import net.sf.samtools.SAMRecord;
import org.broadinstitute.sting.playground.gatk.walkers.Recalibration.LocalMapType;
/**

View File

@ -1,5 +1,6 @@
package org.broadinstitute.sting.playground.gatk.walkers.Recalibration;
import org.broadinstitute.sting.playground.gatk.walkers.Recalibration.LocalMapType;
import org.broadinstitute.sting.gatk.walkers.LocusWalker;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
@ -7,8 +8,6 @@ import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.utils.QualityUtils;
import net.sf.samtools.SAMRecord;
import java.util.List;
/**
* Created by IntelliJ IDEA.
* User: chartl
@ -29,7 +28,7 @@ public class NQSMismatchCovariantWalker extends LocusWalker<LocalMapType, long[]
public static final String HEADER_FORMAT = "%s %s %s %s %s %s%n";
public long[][][] reduceInit() {
long[][][] mismatchesByNQS = new long[(int)QualityUtils.MAX_QUAL_SCORE+1][this.numNQSGroups()+1][2];
long[][][] mismatchesByNQS = new long[(int) QualityUtils.MAX_QUAL_SCORE+1][this.numNQSGroups()+1][2];
for ( int qualityScore = 0; qualityScore <= (int) QualityUtils.MAX_QUAL_SCORE; qualityScore ++) {
for ( int nqsGroup = 0; nqsGroup <= numNQSGroups(); nqsGroup++ ) {
mismatchesByNQS[qualityScore][nqsGroup][MM_OFFSET] = 0;
@ -124,7 +123,7 @@ public class NQSMismatchCovariantWalker extends LocusWalker<LocalMapType, long[]
protected boolean isDBSNP(RefMetaDataTracker tracker) {
return false;
return false;
// return ( tracker.lookup("dbSNP",null) != null );
}
@ -187,24 +186,3 @@ public class NQSMismatchCovariantWalker extends LocusWalker<LocalMapType, long[]
}
class LocalMapType {
public AlignmentContext context;
public ReferenceContext ref;
public RefMetaDataTracker tracker;
public LocalMapType(AlignmentContext context, ReferenceContext ref, RefMetaDataTracker tracker) {
this.context = context;
this.ref = ref;
this.tracker = tracker;
}
public int numReads() {
return context.numReads();
}
public int qScore(int read) {
return (int) context.getReads().get(read).getBaseQualities()[context.getOffsets().get(read)];
}
}

View File

@ -9,7 +9,6 @@ import org.broadinstitute.sting.utils.StingException;
import org.broadinstitute.sting.utils.Pair;
import java.util.LinkedList;
import java.util.ArrayList;
import java.util.List;
import java.util.ListIterator;
import org.apache.log4j.Logger;