From ddb5045e1406bd7bd4d0da97b7c0bf5525ca4997 Mon Sep 17 00:00:00 2001 From: Ryan Poplin Date: Fri, 19 Aug 2011 19:29:51 -0400 Subject: [PATCH] Updating the methods development calling pipeline for the new rod binding syntax and the new best practices. --- .../VariantDataManager.java | 11 +- .../VariantRecalibrator.java | 9 +- .../MethodsDevelopmentCallingPipeline.scala | 109 +++++++++--------- 3 files changed, 74 insertions(+), 55 deletions(-) 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 cb4d94332..429becfc7 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 @@ -233,13 +233,15 @@ public class VariantDataManager { } 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) { + final List> training, final List> truth, final List> known, final List> badSites, final List> resource) { datum.isKnown = false; datum.atTruthSite = false; datum.atTrainingSite = false; datum.atAntiTrainingSite = false; datum.prior = 2.0; + //BUGBUG: need to clean this up + for( final RodBinding rod : training ) { for( final VariantContext trainVC : tracker.getValues(rod, genomeLoc) ) { if( isValidVariant( evalVC, trainVC, TRUST_ALL_POLYMORPHIC ) ) { @@ -264,6 +266,13 @@ public class VariantDataManager { } } } + for( final RodBinding rod : resource ) { + for( final VariantContext trainVC : tracker.getValues(rod, genomeLoc) ) { + if( isValidVariant( evalVC, trainVC, TRUST_ALL_POLYMORPHIC ) ) { + 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 ) { 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 517c2362a..d34a13427 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 @@ -138,6 +138,12 @@ public class VariantRecalibrator extends RodWalker> badSites = Collections.emptyList(); + /** + * Any set of sites for which you would like to apply a prior probability but for which you don't want to use as training, truth, or known sites. + */ + @Input(fullName="resource", shortName = "resource", doc="A list of sites for which to apply a prior probability of being correct but which aren't used by the algorithm", required=false) + public List> resource = Collections.emptyList(); + ///////////////////////////// // Outputs ///////////////////////////// @@ -226,6 +232,7 @@ public class VariantRecalibrator extends RodWalker rod : allInputBindings ) { try { rodToPriorMap.put(rod.getName(), (rod.getTags().containsKey("prior") ? Double.parseDouble(rod.getTags().getValue("prior")) : 0.0) ); @@ -263,7 +270,7 @@ public class VariantRecalibrator extends RodWalker new Target("NA12878.HiSeq", hg18, dbSNP_hg18_129, hapmap_hg18, "/humgen/gsa-hpprojects/dev/depristo/oneOffProjects/1000GenomesProcessingPaper/wgs.v13/HiSeq.WGS.cleaned.indels.10.mask", new File("/humgen/gsa-hpprojects/NA12878Collection/bams/NA12878.HiSeq.WGS.bwa.cleaned.recal.bam"), new File("/home/radon01/depristo/work/oneOffProjects/1000GenomesProcessingPaper/wgs.v13/HiSeq.WGS.cleaned.ug.snpfiltered.indelfiltered.vcf"), - "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg18.intervals", 2.07, 99.0, !lowPass), - "HiSeq19" -> new Target("NA12878.HiSeq19", hg19, dbSNP_b37_129, hapmap_b37, indelMask_b37, + "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg18.intervals", 2.17, 99.0, !lowPass, !exome, 1), + "HiSeq19" -> new Target("NA12878.HiSeq19", hg19, dbSNP_b37, hapmap_b37, indelMask_b37, new File("/humgen/gsa-hpprojects/NA12878Collection/bams/NA12878.HiSeq.WGS.bwa.cleaned.recal.hg19.bam"), new File("/humgen/gsa-hpprojects/dev/carneiro/hiseq19/analysis/snps/NA12878.HiSeq19.filtered.vcf"), - "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg19.intervals", 2.3, 99.0, !lowPass), - "GA2hg19" -> new Target("NA12878.GA2.hg19", hg19, dbSNP_b37_129, hapmap_b37, indelMask_b37, + "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg19.intervals", 2.17, 99.0, !lowPass, !exome, 1), + "GA2hg19" -> new Target("NA12878.GA2.hg19", hg19, dbSNP_b37, hapmap_b37, indelMask_b37, new File("/humgen/gsa-hpprojects/NA12878Collection/bams/NA12878.GA2.WGS.bwa.cleaned.hg19.bam"), new File("/humgen/gsa-hpprojects/dev/carneiro/hiseq19/analysis/snps/NA12878.GA2.hg19.filtered.vcf"), - "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg19.intervals", 2.3, 99.0, !lowPass), + "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg19.intervals", 2.17, 99.0, !lowPass, !exome, 1), "WEx" -> new Target("NA12878.WEx", hg18, dbSNP_hg18_129, hapmap_hg18, "/humgen/gsa-hpprojects/dev/depristo/oneOffProjects/1000GenomesProcessingPaper/wgs.v13/GA2.WEx.cleaned.indels.10.mask", new File("/humgen/gsa-hpprojects/NA12878Collection/bams/NA12878.WEx.cleaned.recal.bam"), new File("/home/radon01/depristo/work/oneOffProjects/1000GenomesProcessingPaper/wgs.v13/GA2.WEx.cleaned.ug.snpfiltered.indelfiltered.vcf"), - "/seq/references/HybSelOligos/whole_exome_agilent_1.1_refseq_plus_3_boosters/whole_exome_agilent_1.1_refseq_plus_3_boosters.targets.interval_list", 2.6, 97.0, !lowPass), - "WExTrio" -> new Target("CEUTrio.WEx", hg19, dbSNP_b37_129, hapmap_b37, indelMask_b37, + "/seq/references/HybSelOligos/whole_exome_agilent_1.1_refseq_plus_3_boosters/whole_exome_agilent_1.1_refseq_plus_3_boosters.targets.interval_list", 2.6, 97.0, !lowPass, exome, 1), + "WExTrio" -> new Target("CEUTrio.WEx", hg19, dbSNP_b37, hapmap_b37, indelMask_b37, new File("/humgen/gsa-hpprojects/NA12878Collection/bams/CEUTrio.HiSeq.WEx.bwa.cleaned.recal.bam"), new File("/humgen/gsa-hpprojects/dev/carneiro/trio/analysis/snps/CEUTrio.WEx.filtered.vcf"), - "/seq/references/HybSelOligos/whole_exome_agilent_1.1_refseq_plus_3_boosters/whole_exome_agilent_1.1_refseq_plus_3_boosters.Homo_sapiens_assembly19.targets.interval_list", 2.6, 97.0, !lowPass), - "WGSTrio" -> new Target("CEUTrio.WGS", hg19, dbSNP_b37_129, hapmap_b37, indelMask_b37, + "/seq/references/HybSelOligos/whole_exome_agilent_1.1_refseq_plus_3_boosters/whole_exome_agilent_1.1_refseq_plus_3_boosters.Homo_sapiens_assembly19.targets.interval_list", 2.6, 97.0, !lowPass, exome, 3), + "WGSTrio" -> new Target("CEUTrio.WGS", hg19, dbSNP_b37, hapmap_b37, indelMask_b37, new File("/humgen/gsa-hpprojects/NA12878Collection/bams/CEUTrio.HiSeq.WGS.bwa.cleaned.recal.bam"), new File("/humgen/gsa-hpprojects/dev/carneiro/trio/analysis/snps/CEUTrio.WEx.filtered.vcf"), // ** THIS GOLD STANDARD NEEDS TO BE CORRECTED ** - "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg19.intervals", 2.3, 99.0, !lowPass), + "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.hg19.intervals", 2.3, 99.0, !lowPass, !exome, 3), "FIN" -> new Target("FIN", b37, dbSNP_b37, hapmap_b37, indelMask_b37, new File("/humgen/1kg/processing/pipeline_test_bams/FIN.79sample.Nov2010.chr20.bam"), new File("/humgen/gsa-hpprojects/dev/data/AugChr20Calls_v4_3state/ALL.august.v4.chr20.filtered.vcf"), // ** THIS GOLD STANDARD NEEDS TO BE CORRECTED ** - "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.chr20.hg19.intervals", 2.3, 99.0, lowPass), + "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.chr20.hg19.intervals", 2.3, 99.0, lowPass, !exome, 79), "TGPWExGdA" -> new Target("1000G.WEx.GdA", b37, dbSNP_b37, hapmap_b37, indelMask_b37, new File("/humgen/1kg/processing/pipeline_test_bams/Barcoded_1000G_WEx_Reduced_Plate_1.cleaned.list"), // BUGBUG: reduce from 60 to 20 people new File("/humgen/gsa-scr1/delangel/NewUG/calls/AugustRelease.filtered_Q50_QD5.0_SB0.0.allSamples.SNPs_hg19.WEx_UG_newUG_MQC.vcf"), // ** THIS GOLD STANDARD NEEDS TO BE CORRECTED ** - "/seq/references/HybSelOligos/whole_exome_agilent_1.1_refseq_plus_3_boosters/whole_exome_agilent_1.1_refseq_plus_3_boosters.Homo_sapiens_assembly19.targets.interval_list", 2.6, 99.0, !lowPass), + "/seq/references/HybSelOligos/whole_exome_agilent_1.1_refseq_plus_3_boosters/whole_exome_agilent_1.1_refseq_plus_3_boosters.Homo_sapiens_assembly19.targets.interval_list", 2.6, 99.0, !lowPass, exome, 96), "LowPassN60" -> new Target("lowpass.N60", b36, dbSNP_b36, hapmap_b36, indelMask_b36, new File("/humgen/1kg/analysis/bamsForDataProcessingPapers/lowpass_b36/lowpass.chr20.cleaned.matefixed.bam"), // the bam list to call from new File("/home/radon01/depristo/work/oneOffProjects/VQSRCutByNRS/lowpass.N60.chr20.filtered.vcf"), // the gold standard VCF file to run through the VQSR - "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.chr20.b36.intervals", 2.3, 99.0, lowPass), // chunked interval list to use with Queue's scatter/gather functionality + "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.chr20.b36.intervals", 2.3, 99.0, lowPass, !exome, 60), // chunked interval list to use with Queue's scatter/gather functionality "LowPassEUR363Nov" -> new Target("EUR.nov2010", b37, dbSNP_b37, hapmap_b37, indelMask_b37, new File("/humgen/1kg/processing/pipeline_test_bams/EUR.363sample.Nov2010.chr20.bam"), new File("/humgen/gsa-hpprojects/dev/data/AugChr20Calls_v4_3state/ALL.august.v4.chr20.filtered.vcf"), // ** THIS GOLD STANDARD NEEDS TO BE CORRECTED ** - "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.chr20.hg19.intervals", 2.3, 99.0, lowPass) + "/humgen/1kg/processing/pipeline_test_bams/whole_genome_chunked.chr20.hg19.intervals", 2.3, 99.0, lowPass, !exome, 363) ) @@ -187,22 +194,18 @@ class MethodsDevelopmentCallingPipeline extends QScript { } def bai(bam: File) = new File(bam + ".bai") - val FiltersToIgnore = List("DPFilter", "ABFilter", "ESPStandard", "QualByDepth", "StrandBias", "HomopolymerRun") // 1.) Unified Genotyper Base class GenotyperBase (t: Target) extends UnifiedGenotyper with UNIVERSAL_GATK_ARGS { this.memoryLimit = 3 this.reference_sequence = t.reference this.intervalsString ++= List(t.intervals) - this.scatterCount = 63 // the smallest interval list has 63 intervals, one for each Mb on chr20 + this.scatterCount = 120 this.dcov = if ( t.isLowpass ) { 50 } else { 250 } this.stand_call_conf = if ( t.isLowpass ) { 4.0 } else { 30.0 } this.stand_emit_conf = if ( t.isLowpass ) { 4.0 } else { 30.0 } this.input_file :+= t.bamList - if (t.dbsnpFile.endsWith(".rod")) - this.DBSNP = new File(t.dbsnpFile) - else if (t.dbsnpFile.endsWith(".vcf")) - this.rodBind :+= RodBind("dbsnp", "VCF", t.dbsnpFile) + this.D = new File(t.dbsnpFile) } // 1a.) Call SNPs with UG @@ -216,7 +219,6 @@ class MethodsDevelopmentCallingPipeline extends QScript { this.baq = if (noBAQ) {org.broadinstitute.sting.utils.baq.BAQ.CalculationMode.OFF} else {org.broadinstitute.sting.utils.baq.BAQ.CalculationMode.CALCULATE_AS_NECESSARY} this.analysisName = t.name + "_UGs" this.jobName = queueLogDir + t.name + ".snpcall" - this.A ++= List("FisherStrand") } // 1b.) Call Indels with UG @@ -234,15 +236,14 @@ class MethodsDevelopmentCallingPipeline extends QScript { this.reference_sequence = t.reference this.intervalsString ++= List(t.intervals) this.scatterCount = 10 - this.filterName ++= List("HARD_TO_VALIDATE") - this.filterExpression ++= List("\"MQ0 >= 4 && (MQ0 / (1.0 * DP)) > 0.1\"") - this.variantVCF = t.rawIndelVCF + this.V = t.rawIndelVCF this.out = t.filteredIndelVCF - this.filterName ++= List("LowQual", "StrandBias", "QualByDepth", "HomopolymerRun") - if (t.isLowpass) - this.filterExpression ++= List("\"QUAL<30.0\"", "\"SB>=-1.0\"", "\"QD<1.0\"", "\"HRun>=15\"") - else - this.filterExpression ++= List("\"QUAL<50.0\"", "\"SB>=-1.0\"", "\"QD<5.0\"", "\"HRun>=15\"") + this.filterName ++= List("IndelQD", "IndelReadPosRankSum", "IndelFS") + this.filterExpression ++= List("\"QD < 2.0\"", "\"ReadPosRankSum < -20.0\"", "\"FS > 200.0\"") + if (t.nSamples >= 10) { + this.filterName ++= List("IndelInbreedingCoeff") + this.filterExpression ++= List("\"InbreedingCoeff < -0.8\"") + } this.analysisName = t.name + "_VF" this.jobName = queueLogDir + t.name + ".indelfilter" } @@ -252,17 +253,22 @@ class MethodsDevelopmentCallingPipeline extends QScript { this.memoryLimit = 4 this.reference_sequence = t.reference this.intervalsString ++= List(t.intervals) - this.rodBind :+= RodBind("input", "VCF", if ( goldStandard ) { t.goldStandard_VCF } else { t.rawVCF } ) - this.rodBind :+= RodBind("hapmap", "VCF", t.hapmapFile, "known=false,training=true,truth=true,prior=15.0") - if( t.hapmapFile.contains("b37") ) - this.rodBind :+= RodBind("omni", "VCF", omni_b37, "known=false,training=true,truth=true,prior=12.0") - else if( t.hapmapFile.contains("b36") ) - this.rodBind :+= RodBind("omni", "VCF", omni_b36, "known=false,training=true,truth=true,prior=12.0") - if (t.dbsnpFile.endsWith(".rod")) - this.rodBind :+= RodBind("dbsnp", "DBSNP", t.dbsnpFile, "known=true,training=false,truth=false,prior=10.0") - else if (t.dbsnpFile.endsWith(".vcf")) - this.rodBind :+= RodBind("dbsnp", "VCF", t.dbsnpFile, "known=true,training=false,truth=false,prior=10.0") - this.use_annotation ++= List("QD", "HaplotypeScore", "MQRankSum", "ReadPosRankSum", "HRun", "FS") + this.input :+= ( if ( goldStandard ) { t.goldStandard_VCF } else { t.rawVCF } ) + this.training :+= new TaggedFile( t.hapmapFile, "prior=15.0") + this.truth :+= new TaggedFile( t.hapmapFile, "prior=15.0") + this.training :+= new TaggedFile( omni_b37, "prior=12.0") + this.truth :+= new TaggedFile( omni_b37, "prior=12.0") + this.training :+= new TaggedFile( training_1000G, "prior=10.0" ) + this.badSites :+= new TaggedFile( badSites_1000G, "prior=2.0" ) + this.known :+= new TaggedFile( t.dbsnpFile, "prior=2.0" ) + this.resource :+= new TaggedFile( projectConsensus_1000G, "prior=10.0" ) + this.use_annotation ++= List("QD", "HaplotypeScore", "MQRankSum", "ReadPosRankSum", "MQ", "FS") + if(t.nSamples >= 10) { + this.use_annotation ++= List("InbreedingCoeff") + } + if(!t.isExome) { + this.use_annotation ++= List("DP") + } this.tranches_file = if ( goldStandard ) { t.goldStandardTranchesFile } else { t.tranchesFile } this.recal_file = if ( goldStandard ) { t.goldStandardRecalFile } else { t.recalFile } this.allPoly = true @@ -277,7 +283,7 @@ class MethodsDevelopmentCallingPipeline extends QScript { this.memoryLimit = 4 this.reference_sequence = t.reference this.intervalsString ++= List(t.intervals) - this.rodBind :+= RodBind("input", "VCF", if ( goldStandard ) { t.goldStandard_VCF } else { t.rawVCF } ) + this.input :+= ( if ( goldStandard ) { t.goldStandard_VCF } else { t.rawVCF } ) this.tranches_file = if ( goldStandard ) { t.goldStandardTranchesFile } else { t.tranchesFile} this.recal_file = if ( goldStandard ) { t.goldStandardRecalFile } else { t.recalFile } this.ts_filter_level = t.trancheTarget @@ -290,19 +296,16 @@ class MethodsDevelopmentCallingPipeline extends QScript { class EvalBase(t: Target) extends VariantEval with UNIVERSAL_GATK_ARGS { this.memoryLimit = 3 this.reference_sequence = t.reference - this.rodBind :+= RodBind("comphapmap", "VCF", t.hapmapFile) + this.comp :+= new TaggedFile(t.hapmapFile, "hapmap" ) this.intervalsString ++= List(t.intervals) - if (t.dbsnpFile.endsWith(".rod")) - this.DBSNP = new File(t.dbsnpFile) - else if (t.dbsnpFile.endsWith(".vcf")) - this.rodBind :+= RodBind("dbsnp", "VCF", t.dbsnpFile) + this.D = new File(t.dbsnpFile) this.sample = samples } // 5a.) SNP Evaluation (OPTIONAL) based on the cut vcf class snpEvaluation(t: Target) extends EvalBase(t) { - if (t.reference == b37 || t.reference == hg19) this.rodBind :+= RodBind("compomni", "VCF", omni_b37) - this.rodBind :+= RodBind("eval", "VCF", t.recalibratedVCF ) + if (t.reference == b37 || t.reference == hg19) this.comp :+= new TaggedFile( omni_b37, "omni" ) + this.eval :+= t.recalibratedVCF this.out = t.evalFile this.analysisName = t.name + "_VEs" this.jobName = queueLogDir + t.name + ".snp.eval" @@ -310,7 +313,7 @@ class MethodsDevelopmentCallingPipeline extends QScript { // 5b.) Indel Evaluation (OPTIONAL) class indelEvaluation(t: Target) extends EvalBase(t) { - this.rodBind :+= RodBind("eval", "VCF", t.filteredIndelVCF) + this.eval :+= t.filteredIndelVCF this.evalModule :+= "IndelStatistics" this.out = t.evalIndelFile this.analysisName = t.name + "_VEi"