From 50c870cfce7d0117f89abbece8bdfcdc42ed0d18 Mon Sep 17 00:00:00 2001 From: carneiro Date: Wed, 16 Feb 2011 17:22:30 +0000 Subject: [PATCH] Data Processing Pipeline: local indel realignment, mark duplicates and BQSR. Done. Pacbio pipeline: now all pacbio bams have baq annotated in so running UG is uber fast. Methods pipeline: minor cosmetic changes. git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5253 348d0f76-0448-11de-a6fe-93d51630548a --- .../MethodsDevelopmentCallingPipeline.scala | 2 +- .../oneoffs/carneiro/dataProcessing.scala | 39 +++++++++++++------ .../qscript/oneoffs/carneiro/pbCalling.scala | 6 +-- 3 files changed, 29 insertions(+), 18 deletions(-) diff --git a/scala/qscript/core/MethodsDevelopmentCallingPipeline.scala b/scala/qscript/core/MethodsDevelopmentCallingPipeline.scala index 928c6666d..f8c60b55c 100755 --- a/scala/qscript/core/MethodsDevelopmentCallingPipeline.scala +++ b/scala/qscript/core/MethodsDevelopmentCallingPipeline.scala @@ -1,4 +1,4 @@ -import org.broadinstitute.sting.queue.extensions.gatk.CommandLineGATK +import org.broadinstitute.sting.gatk.CommandLineGATK import org.broadinstitute.sting.queue.extensions.gatk._ import org.broadinstitute.sting.queue.QScript import org.broadinstitute.sting.gatk.phonehome.GATKRunReport diff --git a/scala/qscript/oneoffs/carneiro/dataProcessing.scala b/scala/qscript/oneoffs/carneiro/dataProcessing.scala index 05173738e..69ae3887c 100755 --- a/scala/qscript/oneoffs/carneiro/dataProcessing.scala +++ b/scala/qscript/oneoffs/carneiro/dataProcessing.scala @@ -1,7 +1,7 @@ - import org.broadinstitute.sting.queue.extensions.gatk._ import org.broadinstitute.sting.queue.extensions.picard.PicardBamJarFunction import org.broadinstitute.sting.queue.QScript +import org.broadinstitute.sting.queue.extensions.samtools.SamtoolsIndexFunction class dataProcessing extends QScript { @@ -31,12 +31,16 @@ class dataProcessing extends QScript { def script = { + // Files generated by the pipeline val baseName: String = swapExt(qscript.inputBam, ".bam", "").toString() - def cleanedBam: String = baseName + ".cleaned.bam" - def dedupedBam: String = baseName + ".cleaned.dedup.bam" - def metricsFile: String = swapExt(qscript.inputBam, "bam", "metrics").toString() - def targetIntervals: String = baseName + ".indel.intervals" + val cleanedBam: String = baseName + ".cleaned.bam" + val dedupedBam: String = baseName + ".cleaned.dedup.bam" + val metricsFile: String = swapExt(qscript.inputBam, "bam", "metrics").toString() + val targetIntervals: String = baseName + ".indel.intervals" + val recalFile: String = baseName + ".recal.csv" + val recalBam: String = baseName + ".cleaned.dedup.recal.bam" + // Reference sequence, dbsnps and RODs used by the pipeline val reference: File = new File("/humgen/1kg/reference/human_g1k_v37.fasta") val dbSNP: File = new File("/humgen/gsa-hpprojects/GATK/data/dbsnp_132_b37.leftAligned.vcf") val dindelPilotCalls: String = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Unvalidated/1kg.pilot_release.merged.indels.sites.hg19.vcf" @@ -44,16 +48,16 @@ class dataProcessing extends QScript { val dindelASNCalls: String = "/humgen/1kg/DCC/ftp/technical/working/20110126_dindel_august/ASN.dindel_august_release_merged_pilot1.20110126.sites.vcf.gz" val dindelEURCalls: String = "/humgen/1kg/DCC/ftp/technical/working/20110126_dindel_august/EUR.dindel_august_release_merged_pilot1.20110126.sites.vcf.gz" + // General arguments to all programs trait CommandLineGATKArgs extends CommandLineGATK { this.jarFile = qscript.gatkJar this.reference_sequence = reference - this.memoryLimit = Some(2) + this.memoryLimit = Some(4) this.jobTempDir = qscript.outputTmpDir } val target = new RealignerTargetCreator with CommandLineGATKArgs - target.memoryLimit = Some(4) target.input_file :+= qscript.inputBam target.out = new File(targetIntervals) target.mismatchFraction = Some(0.0) @@ -69,7 +73,6 @@ class dataProcessing extends QScript { // 2.) Clean without SW val clean = new IndelRealigner with CommandLineGATKArgs - clean.memoryLimit = Some(4) clean.input_file :+= qscript.inputBam clean.targetIntervals = new File(targetIntervals) clean.out = new File(cleanedBam) @@ -88,9 +91,9 @@ class dataProcessing extends QScript { // 3.) Mark Duplicates val dedup = new PicardBamJarFunction{ @Input(doc="cleaned bam") var clean: File = new File(cleanedBam) - @Output(doc="deduped bam") var dedup: File = new File(dedupedBam) + @Output(doc="deduped bam") var deduped: File = new File(dedupedBam) override def inputBams = List(clean) - override def outputBam = dedup + override def outputBam = deduped override def commandLine = super.commandLine + " M=" + metricsFile sortOrder = null } @@ -98,9 +101,21 @@ class dataProcessing extends QScript { dedup.jarFile = qscript.dedupJar dedup.jobName = baseName + ".dedup" - // val cov = new CountCovariates with CommandLineGATKArgs + val index = new SamtoolsIndexFunction + index.bamFile = new File(dedupedBam) + index.analysisName = baseName + ".index" + val cov = new CountCovariates with CommandLineGATKArgs + cov.rodBind :+= RodBind("dbsnp", "VCF", dbSNP) + cov.covariate ++= List("ReadGroupCovariate", "QualityScoreCovariate", "CycleCovariate", "DinucCovariate") + cov.input_file :+= new File(dedupedBam) + cov.recal_file = new File(recalFile) - add(target, clean, dedup) + val recal = new TableRecalibration with CommandLineGATKArgs + recal.input_file :+= new File (dedupedBam) + recal.recal_file = new File(recalFile) + recal.out = new File(recalBam) + + add(target, clean, dedup, index, cov, recal) } } \ No newline at end of file diff --git a/scala/qscript/oneoffs/carneiro/pbCalling.scala b/scala/qscript/oneoffs/carneiro/pbCalling.scala index 1d31216db..31e525b55 100755 --- a/scala/qscript/oneoffs/carneiro/pbCalling.scala +++ b/scala/qscript/oneoffs/carneiro/pbCalling.scala @@ -15,10 +15,6 @@ class pbCalling extends QScript { @Argument(shortName="dataset", doc="selects the datasets to run. If not provided, all datasets will be used", required=false) var datasets: List[String] = Nil - @Argument(shortName="noBAQ", doc="turns off BAQ calculation", required=false) - var noBAQ: Boolean = false - - class Target( val baseName: String, @@ -169,7 +165,7 @@ class pbCalling extends QScript { this.stand_emit_conf = Some( if ( t.isLowpass ) { 4.0 } else { 30.0 } ) this.input_file :+= t.bamList this.out = t.rawVCF - this.baq = Some( if (noBAQ) {org.broadinstitute.sting.utils.baq.BAQ.CalculationMode.OFF} else {org.broadinstitute.sting.utils.baq.BAQ.CalculationMode.RECALCULATE}) + this.baq = Some(org.broadinstitute.sting.utils.baq.BAQ.CalculationMode.CALCULATE_AS_NECESSARY) this.analysisName = t.name + "_UG" if (t.dbsnpFile.endsWith(".rod")) this.DBSNP = new File(t.dbsnpFile)