diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyRecalibration.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyRecalibration.java index 74b7b8e7d..abe27e483 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyRecalibration.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyRecalibration.java @@ -28,9 +28,9 @@ package org.broadinstitute.sting.gatk.walkers.variantrecalibration; import org.broadinstitute.sting.commandline.Argument; import org.broadinstitute.sting.commandline.Input; import org.broadinstitute.sting.commandline.Output; +import org.broadinstitute.sting.commandline.RodBinding; import org.broadinstitute.sting.gatk.contexts.AlignmentContext; import org.broadinstitute.sting.gatk.contexts.ReferenceContext; -import org.broadinstitute.sting.gatk.datasources.rmd.ReferenceOrderedDataSource; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.walkers.RodWalker; import org.broadinstitute.sting.utils.SampleUtils; @@ -56,6 +56,11 @@ public class ApplyRecalibration extends RodWalker { ///////////////////////////// // Inputs ///////////////////////////// + /** + * The raw input variants to be recalibrated. + */ + @Input(fullName="input", shortName = "input", doc="The raw input variants to be recalibrated", required=true) + public List> input; @Input(fullName="recal_file", shortName="recalFile", doc="The output recal file used by ApplyRecalibration", required=true) private File RECAL_FILE; @Input(fullName="tranches_file", shortName="tranchesFile", doc="The input tranches file describing where to cut the data", required=true) @@ -101,17 +106,8 @@ public class ApplyRecalibration extends RodWalker { } Collections.reverse(tranches); // this algorithm wants the tranches ordered from best (lowest truth sensitivity) to worst (highest truth sensitivity) - for( final ReferenceOrderedDataSource d : this.getToolkit().getRodDataSources() ) { - if( d.getName().startsWith("input") ) { - inputNames.add(d.getName()); - logger.info("Found input variant track with name " + d.getName()); - } else { - logger.info("Not evaluating ROD binding " + d.getName()); - } - } - - if( inputNames.size() == 0 ) { - throw new UserException.BadInput( "No input variant tracks found. Input variant binding names must begin with 'input'." ); + for( final RodBinding rod : input ) { + inputNames.add( rod.getName() ); } if( IGNORE_INPUT_FILTERS != null ) { @@ -168,7 +164,7 @@ public class ApplyRecalibration extends RodWalker { return 1; } - for( VariantContext vc : tracker.getValues(VariantContext.class, inputNames, context.getLocation()) ) { + for( VariantContext vc : tracker.getValues(input, context.getLocation()) ) { if( vc != null ) { if( VariantRecalibrator.checkRecalibrationMode( vc, MODE ) && (vc.isNotFiltered() || ignoreInputFilterSet.containsAll(vc.getFilters())) ) { String filterString = null; diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java index 8687b5796..08026a45e 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java @@ -26,10 +26,10 @@ package org.broadinstitute.sting.gatk.walkers.variantrecalibration; import org.apache.log4j.Logger; +import org.broadinstitute.sting.commandline.RodBinding; import org.broadinstitute.sting.gatk.GenomeAnalysisEngine; -import org.broadinstitute.sting.gatk.contexts.AlignmentContext; -import org.broadinstitute.sting.gatk.contexts.ReferenceContext; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; +import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.MathUtils; import org.broadinstitute.sting.utils.collections.ExpandingArrayList; import org.broadinstitute.sting.utils.exceptions.UserException; @@ -38,6 +38,7 @@ import org.broadinstitute.sting.utils.variantcontext.VariantContext; import java.io.PrintStream; import java.util.ArrayList; import java.util.Collections; +import java.util.HashMap; import java.util.List; /** @@ -105,31 +106,6 @@ public class VariantDataManager { } } - public void addTrainingSet( final TrainingSet trainingSet ) { - trainingSets.add( trainingSet ); - } - - public boolean checkHasTrainingSet() { - for( final TrainingSet trainingSet : trainingSets ) { - if( trainingSet.isTraining ) { return true; } - } - return false; - } - - public boolean checkHasTruthSet() { - for( final TrainingSet trainingSet : trainingSets ) { - if( trainingSet.isTruth ) { return true; } - } - return false; - } - - public boolean checkHasKnownSet() { - for( final TrainingSet trainingSet : trainingSets ) { - if( trainingSet.isKnown ) { return true; } - } - return false; - } - public ExpandingArrayList getTrainingData() { final ExpandingArrayList trainingData = new ExpandingArrayList(); for( final VariantDatum datum : data ) { @@ -240,13 +216,14 @@ public class VariantDataManager { if( jitter && annotationKey.equalsIgnoreCase("HRUN") ) { // Integer valued annotations must be jittered a bit to work in this GMM value += -0.25 + 0.5 * GenomeAnalysisEngine.getRandomGenerator().nextDouble(); } - if (vc.isIndel() && annotationKey.equalsIgnoreCase("QD")) { - // normalize QD by event length for indel case - int eventLength = Math.abs(vc.getAlternateAllele(0).getBaseString().length() - vc.getReference().getBaseString().length()); // ignore multi-allelic complication here for now - if (eventLength > 0) // sanity check - value /= (double)eventLength; - } + if (vc.isIndel() && annotationKey.equalsIgnoreCase("QD")) { + // normalize QD by event length for indel case + int eventLength = Math.abs(vc.getAlternateAllele(0).getBaseString().length() - vc.getReference().getBaseString().length()); // ignore multi-allelic complication here for now + if (eventLength > 0) { // sanity check + value /= (double)eventLength; + } + } if( jitter && annotationKey.equalsIgnoreCase("HaplotypeScore") && MathUtils.compareDoubles(value, 0.0, 0.0001) == 0 ) { value = -0.2 + 0.4*GenomeAnalysisEngine.getRandomGenerator().nextDouble(); } if( jitter && annotationKey.equalsIgnoreCase("FS") && MathUtils.compareDoubles(value, 0.0, 0.001) == 0 ) { value = -0.2 + 0.4*GenomeAnalysisEngine.getRandomGenerator().nextDouble(); } @@ -257,30 +234,44 @@ public class VariantDataManager { return value; } - public void parseTrainingSets( final RefMetaDataTracker tracker, final ReferenceContext ref, final AlignmentContext context, final VariantContext evalVC, final VariantDatum datum, final boolean TRUST_ALL_POLYMORPHIC ) { + public void parseTrainingSets( final RefMetaDataTracker tracker, final GenomeLoc genomeLoc, final VariantContext evalVC, final VariantDatum datum, final boolean TRUST_ALL_POLYMORPHIC, final HashMap rodToPriorMap, + final List> training, final List> truth, final List> known, final List> badSites) { datum.isKnown = false; datum.atTruthSite = false; datum.atTrainingSite = false; datum.atAntiTrainingSite = false; datum.prior = 2.0; - datum.consensusCount = 0; - for( final TrainingSet trainingSet : trainingSets ) { - for( final VariantContext trainVC : tracker.getValues(VariantContext.class, trainingSet.name, ref.getLocus()) ) { - if( trainVC != null && trainVC.isNotFiltered() && trainVC.isVariant() && - ((evalVC.isSNP() && trainVC.isSNP()) || ((evalVC.isIndel()||evalVC.isMixed()) && (trainVC.isIndel()||trainVC.isMixed()))) && - (TRUST_ALL_POLYMORPHIC || !trainVC.hasGenotypes() || trainVC.isPolymorphic()) ) { - - datum.isKnown = datum.isKnown || trainingSet.isKnown; - datum.atTruthSite = datum.atTruthSite || trainingSet.isTruth; - datum.atTrainingSite = datum.atTrainingSite || trainingSet.isTraining; - datum.prior = Math.max( datum.prior, trainingSet.prior ); - datum.consensusCount += ( trainingSet.isConsensus ? 1 : 0 ); + for( final RodBinding rod : training ) { + for( final VariantContext trainVC : tracker.getValues(rod, genomeLoc) ) { + if( isValidVariant( evalVC, trainVC, TRUST_ALL_POLYMORPHIC ) ) { + datum.atTrainingSite = true; + datum.prior = Math.max( datum.prior, (rodToPriorMap.containsKey(rod.getName()) ? rodToPriorMap.get(rod.getName()) : 0.0) ); } + } + } + for( final RodBinding rod : truth ) { + for( final VariantContext trainVC : tracker.getValues(rod, genomeLoc) ) { + if( isValidVariant( evalVC, trainVC, TRUST_ALL_POLYMORPHIC ) ) { + datum.atTruthSite = true; + datum.prior = Math.max( datum.prior, (rodToPriorMap.containsKey(rod.getName()) ? rodToPriorMap.get(rod.getName()) : 0.0) ); + } + } + } + for( final RodBinding rod : known ) { + for( final VariantContext trainVC : tracker.getValues(rod, genomeLoc) ) { + if( isValidVariant( evalVC, trainVC, TRUST_ALL_POLYMORPHIC ) ) { + datum.isKnown = true; + datum.prior = Math.max( datum.prior, (rodToPriorMap.containsKey(rod.getName()) ? rodToPriorMap.get(rod.getName()) : 0.0) ); + } + } + } + for( final RodBinding rod : badSites ) { + for( final VariantContext trainVC : tracker.getValues(rod, genomeLoc) ) { if( trainVC != null ) { - datum.atAntiTrainingSite = datum.atAntiTrainingSite || trainingSet.isAntiTraining; + datum.atAntiTrainingSite = true; + datum.prior = Math.max( datum.prior, (rodToPriorMap.containsKey(rod.getName()) ? rodToPriorMap.get(rod.getName()) : 0.0) ); } - } } } @@ -292,4 +283,10 @@ public class VariantDataManager { (datum.worstAnnotation != -1 ? annotationKeys.get(datum.worstAnnotation) : "NULL"))); } } + + private boolean isValidVariant( final VariantContext evalVC, final VariantContext trainVC, final boolean TRUST_ALL_POLYMORPHIC) { + return trainVC != null && trainVC.isNotFiltered() && trainVC.isVariant() && + ((evalVC.isSNP() && trainVC.isSNP()) || ((evalVC.isIndel()||evalVC.isMixed()) && (trainVC.isIndel()||trainVC.isMixed()))) && + (TRUST_ALL_POLYMORPHIC || !trainVC.hasGenotypes() || trainVC.isPolymorphic()); + } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java index 7bd7ea46d..da9da936b 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java @@ -25,13 +25,9 @@ package org.broadinstitute.sting.gatk.walkers.variantrecalibration; -import org.broadinstitute.sting.commandline.Argument; -import org.broadinstitute.sting.commandline.ArgumentCollection; -import org.broadinstitute.sting.commandline.Hidden; -import org.broadinstitute.sting.commandline.Output; +import org.broadinstitute.sting.commandline.*; import org.broadinstitute.sting.gatk.contexts.AlignmentContext; import org.broadinstitute.sting.gatk.contexts.ReferenceContext; -import org.broadinstitute.sting.gatk.datasources.rmd.ReferenceOrderedDataSource; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.walkers.RodWalker; import org.broadinstitute.sting.gatk.walkers.TreeReducible; @@ -57,11 +53,51 @@ import java.util.*; public class VariantRecalibrator extends RodWalker, ExpandingArrayList> implements TreeReducible> { - public static final String VQS_LOD_KEY = "VQSLOD"; - public static final String CULPRIT_KEY = "culprit"; + public static final String VQS_LOD_KEY = "VQSLOD"; // Log odds ratio of being a true variant versus being false under the trained gaussian mixture model + public static final String CULPRIT_KEY = "culprit"; // The annotation which was the worst performing in the Gaussian mixture model, likely the reason why the variant was filtered out @ArgumentCollection private VariantRecalibratorArgumentCollection VRAC = new VariantRecalibratorArgumentCollection(); + ///////////////////////////// + // Inputs + ///////////////////////////// + /** + * The raw input variants to be recalibrated. + */ + @Input(fullName="input", shortName = "input", doc="The raw input variants to be recalibrated", required=true) + public List> input; + /** + * A list of training variants used to train the Gaussian mixture model. + * + * Input variants which are found to overlap with these training sites are used to build the Gaussian mixture model. + */ + @Input(fullName="training", shortName = "training", doc="A list of training variants used to train the Gaussian mixture model", required=true) + public List> training; + /** + * A list of true variants to be used when deciding the truth sensitivity cut of the final callset. + * + * When deciding where to set the cutoff in VQSLOD sensitivity to these truth sites is used. + * Typically one might want to say I dropped my threshold until I got back 99% of HapMap sites, for example. + */ + @Input(fullName="truth", shortName = "truth", doc="A list of true variants to be used when deciding the truth sensitivity cut of the final callset", required=true) + public List> truth; + /** + * A list of known variants to be used for metric comparison purposes. + * + * The known / novel status of a variant isn't used by the algorithm itself and is only used for reporting / display purposes. + * The output metrics are stratified by known status in order to aid in comparisons with other call sets. + */ + @Input(fullName="known", shortName = "known", doc="A list of known variants to be used for metric comparison purposes", required=false) + public List> known = Collections.emptyList(); + /** + * A list of known bad variants used to supplement training the negative model. + * + * In addition to using the worst 3% of variants as compared to the Gaussian mixture model, we can also supplement the list + * with a database of known bad variants. Maybe these are loci which are frequently filtered out in many projects (centromere, for example). + */ + @Input(fullName="badSites", shortName = "badSites", doc="A list of known bad variants used to supplement training the negative model", required=false) + public List> badSites = Collections.emptyList(); + ///////////////////////////// // Outputs ///////////////////////////// @@ -96,9 +132,9 @@ public class VariantRecalibrator extends RodWalker ignoreInputFilterSet = new TreeSet(); - private final Set inputNames = new HashSet(); private final VariantRecalibratorEngine engine = new VariantRecalibratorEngine( VRAC ); + private final HashMap rodToPriorMap = new HashMap(); //--------------------------------------------------------------------------------------------------------------- // @@ -123,31 +159,24 @@ public class VariantRecalibrator extends RodWalker> allInputBindings = new ArrayList>(); + allInputBindings.addAll(truth); + allInputBindings.addAll(training); + allInputBindings.addAll(known); + allInputBindings.addAll(badSites); + for( final RodBinding rod : allInputBindings ) { + try { + rodToPriorMap.put(rod.getName(), (rod.getTags().containsKey("prior") ? Double.parseDouble(rod.getTags().getValue("prior")) : 0.0) ); + } catch( NumberFormatException e ) { + throw new UserException.BadInput("Bad rod binding syntax. Prior key-value tag detected but isn't parsable. Expecting something like -training:prior=12.0 my.set.vcf"); + } + } } //--------------------------------------------------------------------------------------------------------------- @@ -163,10 +192,12 @@ public class VariantRecalibrator extends RodWalker