Remove greps in clean_vcf_part1 which could match sample IDs to variant IDs (#157)

cwhelan · web-flow · commit ed9ba7a69cde · 2021-05-06T20:25:41.000-04:00
* make script use bcftools instead of greps to calculate sex-specific median values
* udpate sv-pipeline docker
diff --git a/input_values/dockers.json b/input_values/dockers.json
@@ -13,7 +13,7 @@
   "sv_base_docker" : "us.gcr.io/broad-dsde-methods/gatk-sv/sv-base:mw-gnomad-02-6a66c96",
   "sv_base_mini_docker" : "us.gcr.io/broad-dsde-methods/gatk-sv/sv-base-mini:rlc_posthoc_filtering_cnv_mcnv_compatability_9a8561",
   "sv_pipeline_base_docker" : "us.gcr.io/broad-dsde-methods/gatk-sv/sv-pipeline-base:rlc_posthoc_filtering_cnv_mcnv_compatability_9a8561",
-  "sv_pipeline_docker" : "us.gcr.io/broad-dsde-methods/gatk-sv/sv-pipeline:cw_cleanvcf1_ev_script_7cf6798",
+  "sv_pipeline_docker" : "us.gcr.io/broad-dsde-methods/cwhelan/sv-pipeline:cw_clean_vcf_part1_rm_allosome_greps_8348649",
   "sv_pipeline_qc_docker" : "us.gcr.io/broad-dsde-methods/markw/sv-pipeline-qc:mw-xz-fixes-7cbffee",
   "sv_pipeline_rdtest_docker" : "us.gcr.io/broad-dsde-methods/gatk-sv/sv-pipeline-rdtest:mw-gnomad-02-6a66c96",
   "wham_docker" : "us.gcr.io/broad-dsde-methods/wham:8645aa",
diff --git a/src/sv-pipeline/04_variant_resolution/scripts/clean_vcf_part1.sh b/src/sv-pipeline/04_variant_resolution/scripts/clean_vcf_part1.sh
@@ -19,6 +19,7 @@ BCFTOOLS=/usr/local/bin/bcftools
 vcf=$1
 backgroundlist=$2
 famfile=$3
+allosome_fai=$4
 
 ##get sampleids from VCF##
 zcat $vcf \
@@ -58,6 +59,7 @@ zcat convertsvtype.vcf.gz \
   |awk 'NR==FNR{inFileA[$1]=$2; next} {if ($3 in inFileA && $1!~"#") $6=inFileA[$3]; print }' OFS='\t' vargq.persample - \
   |bgzip \
    >cleaninfo.vcf.gz
+   tabix -p vcf cleaninfo.vcf.gz
 
    
 ##fix sex chr if necessary##
@@ -83,21 +85,20 @@ awk '{if ($5==2) print $2}' $famfile \
    if [ $(cat clean.bed.ids.txt|wc -l) -gt 0 ]
    then
 
-    zcat cleaninfo.vcf.gz \
-     |awk '{if ($1!~"#" && ($1~"X" || $1~"Y")) $1=$3;print}' OFS="\t"\
-     |vcftools --vcf - --stdout --extract-FORMAT-info RD_CN \
-     |awk -F"\t" 'NR==1{for (i=3;i<=NF;i++) header[i]=$i} NR>1{for(j=3;j<=NF;j++) print $1"\t"header[j] "\t" $j }' \
-     |fgrep -wf clean.bed.ids.txt \
-     |awk '{if ($3!=".") print}' \
-     |gzip \
-      >RD_CN.sexcheck.FORMAT.gz
-     
-    zcat RD_CN.sexcheck.FORMAT.gz|fgrep -wf male.txt| Rscript -e 'd<-read.table("stdin")' \
+  awk '{print $1"\t0\t"$2}' < ${allosome_fai} > allosomes.list
+   ${BCFTOOLS} query -R allosomes.list -S male.txt -i 'ID=@clean.bed.ids.txt' -f '[%ID\t%SAMPLE\t%RD_CN\n]' cleaninfo.vcf.gz \
+      | awk '{if ($3!=".") print}' \
+      | gzip > RD_CN.sexcheck.FORMAT.male.gz
+
+   ${BCFTOOLS} query -R allosomes.list -S female.txt -i 'ID=@clean.bed.ids.txt' -f '[%ID\t%SAMPLE\t%RD_CN\n]' cleaninfo.vcf.gz \
+      | awk '{if ($3!=".") print}' \
+      | gzip > RD_CN.sexcheck.FORMAT.female.gz
+
+    zcat RD_CN.sexcheck.FORMAT.male.gz| Rscript -e 'd<-read.table("stdin")' \
     -e 'x<-tapply(d[,3],d[,1],median)' \
     -e 'write.table(x,"male.median.value.pervar.txt",col.names=FALSE,quote=FALSE,sep = "\t")'
 
-
-     zcat RD_CN.sexcheck.FORMAT.gz|fgrep -wf female.txt| Rscript -e 'd<-read.table("stdin")' \
+     zcat RD_CN.sexcheck.FORMAT.female.gz| Rscript -e 'd<-read.table("stdin")' \
      -e 'x<-tapply(d[,3],d[,1],median)' \
      -e 'write.table(x,"female.median.value.pervar.txt",col.names=FALSE,quote=FALSE,sep = "\t")'
     fi
@@ -242,7 +243,6 @@ fi
 
 
   tabix -p vcf combinedsex.vcf.gz
-  tabix -p vcf cleaninfo.vcf.gz
 
 zcat combinedsex.vcf.gz|awk '{if ($1!~"#") print $3}'>modified.ids.txt
 
diff --git a/test_input_templates/module0506/Module0506.json.tmpl b/test_input_templates/module0506/Module0506.json.tmpl
@@ -1,6 +1,7 @@
 {
   "Module0506.bin_exclude": {{ reference_resources.bin_exclude | tojson }},
   "Module0506.contig_list": {{ reference_resources.primary_contigs_fai | tojson }},
+  "Module0506.allosome_fai": {{ reference_resources.allosome_file | tojson }},
   "Module0506.cytobands": {{ reference_resources.cytobands | tojson }},
   "Module0506.mei_bed": {{ reference_resources.mei_bed | tojson }},
   "Module0506.pe_exclude_list": {{ reference_resources.pesr_exclude_list | tojson }},
diff --git a/test_input_templates/module0506/Module0506Test.json.tmpl b/test_input_templates/module0506/Module0506Test.json.tmpl
@@ -120,6 +120,7 @@
 
   "Module0506Test.Module0506.bin_exclude": {{ reference_resources.bin_exclude | tojson }},
   "Module0506Test.Module0506.contig_list": {{ reference_resources.primary_contigs_fai | tojson }},
+  "Module0506Test.Module0506.allosome_fai": {{ reference_resources.allosome_file | tojson }},
   "Module0506Test.Module0506.cytobands": {{ reference_resources.cytobands | tojson }},
   "Module0506Test.Module0506.mei_bed": {{ reference_resources.mei_bed | tojson }},
   "Module0506Test.Module0506.pe_exclude_list": {{ reference_resources.pesr_exclude_list | tojson }},
diff --git a/wdl/CleanVcf.wdl b/wdl/CleanVcf.wdl
@@ -9,6 +9,7 @@ workflow CleanVcf {
     String contig
     File background_list
     File ped_file
+    File allosome_fai
     String prefix
     Int max_shards_per_chrom_step1
     File bothsides_pass_list
@@ -60,6 +61,7 @@ workflow CleanVcf {
         ped_file=ped_file,
         sv_pipeline_docker=sv_pipeline_docker,
         bothsides_pass_list=bothsides_pass_list,
+        allosome_fai=allosome_fai,
         runtime_attr_override=runtime_override_clean_vcf_1a
     }
   }
@@ -203,6 +205,7 @@ task CleanVcf1a {
     File ped_file
     String sv_pipeline_docker
     File bothsides_pass_list
+    File allosome_fai
     RuntimeAttr? runtime_attr_override
   }
 
@@ -235,7 +238,7 @@ task CleanVcf1a {
   command <<<
     set -eu -o pipefail
     
-    /opt/sv-pipeline/04_variant_resolution/scripts/clean_vcf_part1.sh ~{vcf} ~{background_list} ~{ped_file}
+    /opt/sv-pipeline/04_variant_resolution/scripts/clean_vcf_part1.sh ~{vcf} ~{background_list} ~{ped_file} ~{allosome_fai}
     /opt/sv-pipeline/04_variant_resolution/scripts/add_bothsides_support_filter.py \
       --bgzip \
       --outfile int.w_bothsides.vcf.gz \
diff --git a/wdl/GATKSVPipelineBatch.wdl b/wdl/GATKSVPipelineBatch.wdl
@@ -251,6 +251,7 @@ workflow GATKSVPipelineBatch {
       pesr_vcfs=[Module04.genotyped_pesr_vcf],
       depth_vcfs=Module04b.regenotyped_depth_vcfs,
       contig_list=primary_contigs_fai,
+      allosome_fai=allosome_file,
       ref_dict=reference_dict,
       disc_files=[GATKSVPipelinePhase1.merged_PE],
       disc_files_index=[GATKSVPipelinePhase1.merged_PE_index],
diff --git a/wdl/GATKSVPipelineSingleSample.wdl b/wdl/GATKSVPipelineSingleSample.wdl
@@ -871,6 +871,7 @@ workflow GATKSVPipelineSingleSample {
       pesr_vcfs=[ConvertCNVsWithoutDepthSupportToBNDs.out_vcf],
       depth_vcfs=[Module04.genotyped_depth_vcf],
       contig_list=primary_contigs_fai,
+      allosome_fai=allosome_file,
       ref_dict=reference_dict,
 
       merge_complex_genotype_vcfs = true,
diff --git a/wdl/Module0506.wdl b/wdl/Module0506.wdl
@@ -31,6 +31,7 @@ workflow Module0506 {
 
     File bin_exclude
     File contig_list
+    File allosome_fai
     Int max_shards_per_chrom
     Int min_variants_per_shard
     File cytobands
@@ -270,6 +271,7 @@ workflow Module0506 {
       complex_resolve_background_fail_lists=Module0506ComplexResolve.complex_resolve_background_fail_lists,
       merged_ped_file=ped_file,
       contig_list=contig_list,
+      allosome_fai=allosome_fai,
       max_shards_per_chrom=max_shards_per_chrom,
       max_shards_per_chrom_clean_vcf_step1=max_shards_per_chrom_clean_vcf_step1,
       min_records_per_shard_clean_vcf_step1=min_records_per_shard_clean_vcf_step1,
diff --git a/wdl/Module0506Clean.wdl b/wdl/Module0506Clean.wdl
@@ -13,6 +13,7 @@ workflow Module0506Clean {
     File merged_ped_file
 
     File contig_list
+    File allosome_fai
     Int max_shards_per_chrom
     Int max_shards_per_chrom_clean_vcf_step1
     Int min_records_per_shard_clean_vcf_step1
@@ -57,6 +58,7 @@ workflow Module0506Clean {
         background_list=complex_resolve_background_fail_lists[i],
         ped_file=merged_ped_file,
         bothsides_pass_list=complex_resolve_bothside_pass_lists[i],
+        allosome_fai=allosome_fai,
         prefix=cohort_name,
         max_shards_per_chrom_step1=max_shards_per_chrom_clean_vcf_step1,
         min_records_per_shard_step1=min_records_per_shard_clean_vcf_step1,

Original file line number	Diff line number	Diff line change
`@@ -1,6 +1,7 @@`
`1`	`1`	`{`
`2`	`2`	`"Module0506.bin_exclude": {{ reference_resources.bin_exclude \| tojson }},`
`3`	`3`	`"Module0506.contig_list": {{ reference_resources.primary_contigs_fai \| tojson }},`
	`4`	`+ "Module0506.allosome_fai": {{ reference_resources.allosome_file \| tojson }},`
`4`	`5`	`"Module0506.cytobands": {{ reference_resources.cytobands \| tojson }},`
`5`	`6`	`"Module0506.mei_bed": {{ reference_resources.mei_bed \| tojson }},`
`6`	`7`	`"Module0506.pe_exclude_list": {{ reference_resources.pesr_exclude_list \| tojson }},`