diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/ApplyVariantClustersWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/ApplyVariantClustersWalker.java index 469a84a0e..c561d7325 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/ApplyVariantClustersWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/ApplyVariantClustersWalker.java @@ -180,7 +180,7 @@ public class ApplyVariantClustersWalker extends RodWalker 0.1 ) { // only use the known prior if the value is specified (meaning not equal to zero) @@ -193,7 +193,7 @@ public class ApplyVariantClustersWalker extends RodWalker annotationMap, final double qualityScore ) { + public static double decodeAnnotation( final String annotationKey, final VariantContext vc ) { + double value = 0.0; + if( annotationKey.equals("AB") && !vc.getAttributes().containsKey(annotationKey) ) { + value = (0.5 - 0.005) + (0.01 * rand.nextDouble()); // HomVar calls don't have an allele balance + } else if( annotationKey.equals("QUAL") ) { + value = vc.getPhredScaledQual(); + } else { + try { + value = Double.parseDouble( (String)vc.getAttribute( annotationKey, "0.0" ) ); + if( Double.isInfinite(value) ) { + value = ( value > 0 ? 1.0 : -1.0 ) * INFINITE_ANNOTATION_VALUE; + } + } catch( NumberFormatException e ) { + // do nothing, default value is 0.0 + } + } + return value; + } + + public final double evaluateVariant( final VariantContext vc ) { final double[] pVarInCluster = new double[numGaussians]; final double[] annotations = new double[dataManager.numAnnotations]; for( int jjj = 0; jjj < dataManager.numAnnotations; jjj++ ) { - double value = 0.0; - final String annotationKey = dataManager.annotationKeys.get(jjj); - if( annotationKey.equals("QUAL") ) { - value = qualityScore; - } else if( annotationKey.equals("AB") && !annotationMap.containsKey(annotationKey) ) { - value = (0.5 - 0.005) + (0.01 * Math.random()); // HomVar calls don't have an allele balance - } else { - try { - final Object stringValue = annotationMap.get( annotationKey ); - if( stringValue != null ) { - value = Double.parseDouble( stringValue.toString() ); - if( Double.isInfinite(value) ) { - value = ( value > 0 ? 1.0 : -1.0 ) * INFINITE_ANNOTATION_VALUE; - } - } - } catch( NumberFormatException e ) { - // do nothing, default value is 0.0 - } - } - + final double value = decodeAnnotation( dataManager.annotationKeys.get(jjj), vc ); annotations[jjj] = (value - dataManager.meanVector[jjj]) / dataManager.varianceVector[jjj]; } diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/VariantOptimizer.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/VariantOptimizer.java index 545afcd87..83d079582 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/VariantOptimizer.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variantoptimizer/VariantOptimizer.java @@ -129,7 +129,7 @@ public class VariantOptimizer extends RodWalker return mapList; } - double annotationValues[] = new double[annotationKeys.size()]; + final double annotationValues[] = new double[annotationKeys.size()]; // todo -- do we really need to support multiple tracks -- logic is cleaner without this case -- what's the use case? for( final VariantContext vc : tracker.getAllVariantContexts(ref, null, context.getLocation(), false, false) ) { @@ -137,23 +137,7 @@ public class VariantOptimizer extends RodWalker if( !vc.isFiltered() || IGNORE_ALL_INPUT_FILTERS || (ignoreInputFilterSet != null && ignoreInputFilterSet.containsAll(vc.getFilters())) ) { int iii = 0; for( final String key : annotationKeys ) { - - double value = 0.0; - if( key.equals("AB") && !vc.getAttributes().containsKey(key) ) { - value = (0.5 - 0.005) + (0.01 * Math.random()); // HomVar calls don't have an allele balance - } else if( key.equals("QUAL") ) { - value = vc.getPhredScaledQual(); - } else { - try { - value = Double.parseDouble( (String)vc.getAttribute( key, "0.0" ) ); - if( Double.isInfinite(value) ) { - value = ( value > 0 ? 1.0 : -1.0 ) * INFINITE_ANNOTATION_VALUE; - } - } catch( NumberFormatException e ) { - // do nothing, default value is 0.0 - } - } - annotationValues[iii++] = value; + annotationValues[iii++] = VariantGaussianMixtureModel.decodeAnnotation( key, vc ); } final VariantDatum variantDatum = new VariantDatum(); @@ -234,7 +218,6 @@ public class VariantOptimizer extends RodWalker final Process p = Runtime.getRuntime().exec( rScriptCommandLine ); p.waitFor(); } catch (InterruptedException e) { - e.printStackTrace(); throw new StingException(e.getMessage()); } catch ( IOException e ) { throw new StingException( "Unable to execute RScript command: " + rScriptCommandLine );