Merge pull request #228 from mskcc/feature/somatic_mutation_pon

Implementation of panel of normal (PoN) annotation of somatic SNVs/indels
mskcc · May 9, 2019 · fa1f570 · fa1f570
2 parents 173c9a9 + 1332ece
commit fa1f570
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 8 deletions.
diff --git a/conf/references.config b/conf/references.config
@@ -57,6 +57,10 @@ params {
       mapabilityBlacklist = "${params.reference_base}/mskcc-igenomes/grch37/annotation/wgEncodeDacMapabilityConsensusExcludable.bed.gz"
       mapabilityBlacklistIndex = "${mapabilityBlacklist}.tbi"
       isoforms = "${params.reference_base}/mskcc-igenomes/grch37/annotation/isoforms"
+      exomePoN = "${params.reference_base}/mskcc-igenomes/grch37/annotation/wes.pon.vcf.gz"
+      exomePoNIndex = "${exomePoN}.tbi"
+      wgsPoN = "${params.reference_base}/mskcc-igenomes/grch37/annotation/wgs.pon.vcf.gz"
+      wgsPoNIndex = "${wgsPoN}.tbi"
     } 
     'GRCh38' {
       acLoci           = "${params.genome_base}/Annotation/ASCAT/1000G_phase3_GRCh38_maf0.3.loci"

diff --git a/somatic.nf b/somatic.nf
@@ -367,14 +367,14 @@ process RunStrelka2 {
   when: 'manta' in tools && 'strelka2' in tools
 
   script:
-  options = "" 
-  if(params.assayType == "exome") options = "--exome"
-
+  options = ""
   intervals = wgsIntervals
   if(params.assayType == "exome") {
+    options = "--exome"
     if(target == 'agilent') intervals = agilentTargets
     if(target == 'idt') intervals = idtTargets
-  }
+   }
+
   """
   configureStrelkaSomaticWorkflow.py \
     ${options} \
@@ -500,6 +500,12 @@ process CombineChannel {
       referenceMap.mapabilityBlacklist,
       referenceMap.mapabilityBlacklistIndex
     ])
+    set file(exomePoN), file(wgsPoN), file(exomePoNIndex), file(wgsPoNIndex) from Channel.value([
+      referenceMap.exomePoN,
+      referenceMap.wgsPoN,
+      referenceMap.exomePoNIndex,
+      referenceMap.wgsPoNIndex
+    ])
 
   output:
     file("${idTumor}.union.pass.vcf") into vcfMergedOutput
@@ -508,10 +514,15 @@ process CombineChannel {
 
   script:
   isec_dir = "${idTumor}.isec"
+  pon = wgsPoN
+  if(params.exome) {
+    pon = wesPon
+  }
   """
   echo -e "##INFO=<ID=MuTect2,Number=0,Type=Flag,Description=\"Variant was called by MuTect2\">\n##INFO=<ID=Strelka2,Number=0,Type=Flag,Description=\"Variant was called by Strelka2\">" > vcf.header
   echo -e '##INFO=<ID=RepeatMasker,Number=1,Type=String,Description="RepeatMasker">' > vcf.rm.header
   echo -e '##INFO=<ID=EncodeDacMapability,Number=1,Type=String,Description="EncodeDacMapability">' > vcf.map.header
+  echo -e '##INFO=<ID=PoN,Number=1,Type=Integer,Description="Count in panel of normals">' > vcf.pon.header
 
   bcftools isec \
     --output-type z \
@@ -562,14 +573,22 @@ process CombineChannel {
     --header-lines vcf.map.header \
     --annotations ${mapabilityBlacklist} \
     --columns CHROM,FROM,TO,EncodeDacMapability \
-    --output-type v \
-    --output ${idTumor}.union.vcf
+    --output-type z \
+    --output ${idTumor}.union.vcf.gz
+
+  bcftools annotate \
+    --header-lines vcf.pon.header \
+    --annotations ${pon} \
+    --columns PoN:=AC \
+    --output-type z \
+    --output ${idTumor}.union.pon.vcf.gz \
+    ${idTumor}.union.vcf.gz
 
   bcftools filter \
     --include 'FILTER=\"PASS\"' \
     --output-type v \
     --output ${idTumor}.union.pass.vcf \
-    ${idTumor}.union.vcf
+    ${idTumor}.union.pon.vcf.gz
   """
 }
 
@@ -611,7 +630,7 @@ process RunVcf2Maf {
     --vcf-normal-id ${idNormal} \
     --input-vcf ${vcfMerged} \
     --ref-fasta ${genomeFile} \
-    --retain-info MuTect2,Strelka2,RepeatMasker,EncodeDacMapability \
+    --retain-info MuTect2,Strelka2,RepeatMasker,EncodeDacMapability,PoN \
     --custom-enst ${isoforms} \
     --output-maf ${outfile} \
     --filter-vcf 0