zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Told people this worked...forgot to commit! 

-c



git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4306 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
chartl 2010-09-18 03:46:00 +00:00
parent a3f31e5df0
commit 6f6d2eb31f
1 changed files with 13 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
scala/qscript/fullCallingPipeline.q
View File

@ -114,12 +114,14 @@ class fullCallingPipeline extends QScript {
cleanBamFiles :+= realigner.out cleanBamFiles :+= realigner.out
add(targetCreator,realigner,samtoolsindex) // COMMENT THIS NEXT LINE TO SKIP CLEANING
//add(targetCreator,realigner,samtoolsindex)
} }
// actually make calls // actually make calls
endToEnd(uncleanedBase,qscript.bamFiles) endToEnd(uncleanedBase,qscript.bamFiles)
endToEnd(cleanedBase,cleanBamFiles) // COMMENT THIS NEXT LINE TO AVOID CALLING ON CLEANED FILES
//endToEnd(cleanedBase,cleanBamFiles)
} }
def endToEnd(base: String, bamFiles: List[File]) = { def endToEnd(base: String, bamFiles: List[File]) = {
@ -134,7 +136,7 @@ class fullCallingPipeline extends QScript {
snps.min_mapping_quality_score = Some(20) snps.min_mapping_quality_score = Some(20)
snps.min_base_quality_score = Some(20) snps.min_base_quality_score = Some(20)
snps.downsample_to_coverage = Some(200) snps.downsample_to_coverage = Some(200)
snps.annotation :+= "QualByDepthV2" //snps.annotation :+= "QualByDepthV2"
if (qscript.trigger != null) { if (qscript.trigger != null) {
snps.trigger_min_confidence_threshold_for_calling = Some(30) snps.trigger_min_confidence_threshold_for_calling = Some(30)
@ -174,8 +176,8 @@ class fullCallingPipeline extends QScript {
loopNo += 1 loopNo += 1
} }
val mergeIndels = new CombineVariants with CommandLineGATKArgs val mergeIndels = new CombineVariants with CommandLineGATKArgs
mergeIndels.out = new File(qscript.project+".indels.vcf") mergeIndels.out = new TaggedFile(qscript.project+".indels.vcf","vcf")
mergeIndels.genotypemergeoption = Some(org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils.GenotypeMergeType.REQUIRE_UNIQUE) mergeIndels.genotypemergeoption = Some(org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils.GenotypeMergeType.UNIQUIFY)
mergeIndels.priority = priority mergeIndels.priority = priority
mergeIndels.variantmergeoption = Some(org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils.VariantMergeType.UNION) mergeIndels.variantmergeoption = Some(org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils.VariantMergeType.UNION)
mergeIndels.rodBind = indelCallFiles mergeIndels.rodBind = indelCallFiles
@ -193,7 +195,7 @@ class fullCallingPipeline extends QScript {
// 2.a filter on cluster and near indels // 2.a filter on cluster and near indels
val masker = new VariantFiltration with CommandLineGATKArgs val masker = new VariantFiltration with CommandLineGATKArgs
masker.rodBind :+= RodBind("variant", "VCF", annotated.out) masker.variantVCF = annotated.out
masker.rodBind :+= RodBind("mask", "VCF", mergeIndels.out) masker.rodBind :+= RodBind("mask", "VCF", mergeIndels.out)
masker.maskName = "NearIndel" masker.maskName = "NearIndel"
masker.clusterWindowSize = Some(qscript.snpClusterWindow) masker.clusterWindowSize = Some(qscript.snpClusterWindow)
@ -203,10 +205,10 @@ class fullCallingPipeline extends QScript {
// 2.b hand filter with standard filter // 2.b hand filter with standard filter
val handFilter = new VariantFiltration with CommandLineGATKArgs val handFilter = new VariantFiltration with CommandLineGATKArgs
handFilter.rodBind :+= RodBind("variant", "VCF", annotated.out) handFilter.variantVCF = masker.out
handFilter.rodBind :+= RodBind("mask", "VCF", mergeIndels.out) handFilter.rodBind :+= RodBind("mask", "VCF", mergeIndels.out)
handFilter.filterName ++= List("StrandBias","AlleleBalance","QualByDepth","HomopolymerRun") handFilter.filterName ++= List("StrandBias","AlleleBalance","QualByDepth","HomopolymerRun")
handFilter.filterExpression ++= List("\"SB>=0.10\"","\"AB>=0.75\"","QD<5","\"HRun>=4\"") handFilter.filterExpression ++= List("\"SB>=0.10\"","\"AB>=0.75\"","\"QD<5.0\"","\"HRun>=4\"")
handFilter.out = swapExt(annotated.out,".vcf",".handfiltered.vcf") handFilter.out = swapExt(annotated.out,".vcf",".handfiltered.vcf")
@ -215,10 +217,11 @@ class fullCallingPipeline extends QScript {
// todo -- args for resources (properties file) // todo -- args for resources (properties file)
val clusters = new GenerateVariantClusters with CommandLineGATKArgs val clusters = new GenerateVariantClusters with CommandLineGATKArgs
clusters.rodBind :+= RodBind("input", "VCF", masker.out) clusters.rodBind :+= RodBind("input", "VCF", masker.out)
clusters.DBSNP = qscript.dbSNP
val clusters_clusterFile = swapExt(new File(snps.out.getAbsolutePath),".vcf",".cluster") val clusters_clusterFile = swapExt(new File(snps.out.getAbsolutePath),".vcf",".cluster")
clusters.clusterFile = clusters_clusterFile clusters.clusterFile = clusters_clusterFile
clusters.memoryLimit = Some(8) clusters.memoryLimit = Some(8)
clusters.jobQueue = "hugemem" clusters.jobQueue = "gsa"
clusters.use_annotation ++= List("QD", "SB", "HaplotypeScore", "HRun") clusters.use_annotation ++= List("QD", "SB", "HaplotypeScore", "HRun")
@ -226,6 +229,7 @@ class fullCallingPipeline extends QScript {
// 3.ii apply gaussian clusters to the masked vcf // 3.ii apply gaussian clusters to the masked vcf
val recalibrate = new VariantRecalibrator with CommandLineGATKArgs val recalibrate = new VariantRecalibrator with CommandLineGATKArgs
recalibrate.clusterFile = clusters.clusterFile recalibrate.clusterFile = clusters.clusterFile
recalibrate.DBSNP = qscript.dbSNP
recalibrate.rodBind :+= RodBind("input", "VCF", masker.out) recalibrate.rodBind :+= RodBind("input", "VCF", masker.out)
recalibrate.out = swapExt(masker.out,".vcf",".recalibrated.vcf") recalibrate.out = swapExt(masker.out,".vcf",".recalibrated.vcf")
recalibrate.target_titv = qscript.target_titv recalibrate.target_titv = qscript.target_titv