Jasmine benchmarks seems worse. Back to SURVIVOR for now.

Author: ricardovialle

example/aj_trio/study_case_1/config_nopanel.yaml

@@ -0,0 +1,14 @@
+OUT_FOLDER: "results_study_case_1_nopanel"
+SAMPLE_KEY: "example/aj_trio/study_case_1/sampleKey.txt"
+
+TOOLS:
+  - "manta"
+  - "smoove"
+  - "delly"
+
+# Reference genome files
+REFERENCE_FASTA: "data/ref/human_g1k_v37.fasta"
+REF_BUILD: "37"
+
+# Gencode GTF for SV annotation
+GENCODE_GTF: "data/annotation/gtf/gencode.v38lift37.annotation.nochr.gtf.gz"

workflow/Snakefile

@@ -3,7 +3,7 @@ import glob
 import yaml
 from snakemake.io import *
 
-version="0.3"
+version="0.4"
 
 ###################################################################################################
 ## Check inputs

workflow/envs/jasmine.yaml

@@ -1,6 +1,6 @@
 channels:
   - bioconda
 dependencies:
-  - jasminesv=1.1.0
+  - jasminesv=1.1.2
   - bcftools=1.9
   - samtools=1.9

workflow/rules/benchmark.smk

@@ -1,7 +1,12 @@
 ###################################################################################################
 ## Benchmarking against GiaB
-###################################################################################################
 # Deactivated in the pipeline.
+###################################################################################################
+# Default
+GIAB_VCF = "data/bench/HG002_SVs_Tier1_v0.6.vcf.gz"
+GIAB_BED = "data/bench/HG002_SVs_Tier1_v0.6.22XY.bed"
+
+# Benchmark raw SV callers
 rule benchmark_raw:
 	input:
 		vcf = OUT_FOLDER + "/sv_discovery/{tool}/{sample}/{sample}.{tool}.vcf.gz"
@@ -10,8 +15,8 @@ rule benchmark_raw:
 	conda:
 		SNAKEDIR + "envs/truvari.yaml"		
 	params:
-		base = "data/bench/HG002_SVs_Tier1_v0.6.vcf.gz",
-		bed = "data/bench/HG002_SVs_Tier1_v0.6.22XY.bed",
+		base = GIAB_VCF,
+		bed = GIAB_BED,
 		ref = REFERENCE_FASTA,
 		outdir = OUT_FOLDER + "/bench/{sample}/sv_discovery/{tool}/"
 	shell:
@@ -27,14 +32,15 @@ rule benchmark_raw:
  			-r 2000 \
  			--giabreport; "
 
+# Benchmark each SV caller after genotyping (e.g. Manta + GraphTyper, or Delly + GraphTyper)
 rule benchmark_raw_gt:
 	input:
 		vcf = OUT_FOLDER + "/sv_discovery/{tool}/{sample}/{sample}.{tool}.gt.vcf.gz"
 	output:
 		OUT_FOLDER + "/bench/{sample}/sv_discovery/{tool}.gt/giab_report.txt"
 	params:
-		base = "data/bench/HG002_SVs_Tier1_v0.6.vcf.gz",
-		bed = "data/bench/HG002_SVs_Tier1_v0.6.22XY.bed",
+		base = GIAB_VCF,
+		bed = GIAB_BED,
 		ref = REFERENCE_FASTA,
 		outdir = OUT_FOLDER + "/bench/{sample}/sv_discovery/{tool}.gt/"
 	conda:
@@ -54,14 +60,15 @@ rule benchmark_raw_gt:
  			-r 2000 \
  			--giabreport; "
 
+# Benchmark each SV caller after joint-genotying. This may also include SV reference panels, if specified.
 rule benchmark_gt:
 	input:
 		vcf = OUT_FOLDER + "/sv_genotyping/{sample}/{sample}.vcf.gz"
 	output:
 		OUT_FOLDER + "/bench/{sample}/sv_genotyping/{tool}/giab_report.txt"
 	params:
-		base = "data/bench/HG002_SVs_Tier1_v0.6.vcf.gz",
-		bed = "data/bench/HG002_SVs_Tier1_v0.6.22XY.bed",
+		base = GIAB_VCF,
+		bed = GIAB_BED,
 		ref = REFERENCE_FASTA,
 		outdir = OUT_FOLDER + "/bench/{sample}/sv_genotyping/{tool}/"
 	conda:
@@ -81,19 +88,23 @@ rule benchmark_gt:
  			-r 2000 \
  			--giabreport; "
 
+# Benchmark SVs from the final genotyped and merged call set.
 rule benchmark_merged_gt:
 	input:
 		vcf = OUT_FOLDER + "/merged_cohort/gt_merged.vcf.gz"
 	output:
-		OUT_FOLDER + "/bench/{sample}/merged/giab_report.txt"
+		vcf = OUT_FOLDER + "/bench/{sample}/merged/vcf.gz",
+		report = OUT_FOLDER + "/bench/{sample}/merged/giab_report.txt"
 	params:
-		base = "data/bench/HG002_SVs_Tier1_v0.6.vcf.gz",
-		bed = "data/bench/HG002_SVs_Tier1_v0.6.22XY.bed",
+		base = GIAB_VCF,
+		bed = GIAB_BED,
 		ref = REFERENCE_FASTA,
 		outdir = OUT_FOLDER + "/bench/{sample}/merged/"
 	conda:
 		SNAKEDIR + "envs/truvari.yaml"
 	shell:
+		"bcftools view -c1 -s {wildcards.sample} -Oz -o {output.vcf} {input.vcf}; "
+		"tabix -p vcf {output.vcf}; "
 		"rm -rf {params.outdir}; "
 		"truvari bench \
  			-b {params.base} \
@@ -104,3 +115,5 @@ rule benchmark_merged_gt:
  			-p 0 \
  			-r 2000 \
  			--giabreport; "
+
+

workflow/rules/genotype.smk

@@ -54,7 +54,7 @@ rule merge_samples_01:
 					shell("ls {sv_ref} >> {output.vcf_list}; ")
 
 # Options of using SURVIVOR or Jasmine. Jasmine by default.
-useJAMINE=True
+useJAMINE=False
 if (useJAMINE):
 	checkpoint merge_samples_02:
 		input:
@@ -69,6 +69,7 @@ if (useJAMINE):
 			vcf_list = OUT_FOLDER + "/merged_cohort/vcf.list",
 		shell:
 			"jasmine file_list={params.vcf_list} out_file={output.vcf}; " 
+			"jasmine --dup_to_ins --postprocess_only out_file={output.vcf}; "
 			"bcftools sort -Oz -o {output.vcfgz} {output.vcf}; "
 			"tabix -p vcf {output.vcfgz}; "
 else:
@@ -83,7 +84,7 @@ else:
 			SNAKEDIR + "envs/survivor.yaml"
 		params:
 			vcf_list = OUT_FOLDER + "/merged_cohort/vcf.list",
-			breakpoint_dist = "500",
+			breakpoint_dist = "100",
 			min_num_calls = "1",
 			use_type = "1",
 			use_strand = "1",
@@ -146,14 +147,18 @@ rule genotype_02:
 	output:
 		vcf = OUT_FOLDER + "/sv_genotyping/{sample}/{sample}.vcf",
 		vcfgz = OUT_FOLDER + "/sv_genotyping/{sample}/{sample}.vcf.gz",
-		tbi = OUT_FOLDER + "/sv_genotyping/{sample}/{sample}.vcf.gz.tbi"
+		tbi = OUT_FOLDER + "/sv_genotyping/{sample}/{sample}.vcf.gz.tbi",
+		filt_vcfgz = OUT_FOLDER + "/sv_genotyping/{sample}/{sample}.filt.vcf.gz",
+		filt_tbi = OUT_FOLDER + "/sv_genotyping/{sample}/{sample}.filt.vcf.gz.tbi"
 	conda:
 		SNAKEDIR + "envs/graphtyper.yaml"
 	priority: 0
 	shell:
-		"graphtyper vcf_concatenate {input.vcf} | bcftools sort > {output.vcf}; "
+		"graphtyper vcf_concatenate {input.vcf} | bcftools view --include \"SVMODEL='AGGREGATED'\" | bcftools sort > {output.vcf}; "
 		"bgzip -c {output.vcf} > {output.vcfgz}; "
 		"tabix -p vcf {output.vcfgz}; "
+		"bcftools view -e QUAL==0 -Oz -o {output.filt_vcfgz} {output.vcfgz}; "
+		"tabix -p vcf {output.filt_vcfgz}; "
 
 rule merge_genotyped_samples:
 	input:
@@ -168,11 +173,10 @@ rule merge_genotyped_samples:
 	conda:
 		SNAKEDIR + "envs/graphtyper.yaml"		
 	shell:
-		"graphtyper vcf_merge {input.vcf} --sv | \
-		grep -E -v \"BREAKPOINT|COVERAGE\" | grep -E -v \"VarType=XG|VarType=IG\" > {output.vcf_full}; "
+		"graphtyper vcf_merge {input.vcf} --sv | grep -E -v \"VarType=XG|VarType=IG\" > {output.vcf_full}; "
 		"bcftools sort -Oz -o {output.vcfgz_full} {output.vcf_full}; "
 		"tabix -p vcf {output.vcfgz_full}; "
-		"bcftools view -f PASS {output.vcfgz_full} > {output.vcf}; "
+		"bcftools view -e QUAL==0 {output.vcfgz_full} > {output.vcf}; "
 		"bgzip -c {output.vcf} > {output.vcfgz}; "
 		"tabix -p vcf {output.vcfgz}; "
 
@@ -207,6 +211,6 @@ rule genotype_discovery:
 			--output={params.outdir}; \
 		done; "
 		"graphtyper vcf_concatenate --no_sort {params.outdir}/*/*.vcf.gz | \
-			grep -E -v \"BREAKPOINT|COVERAGE\" | bcftools sort > {output.vcf};" 
+			bcftools view --include \"SVMODEL='AGGREGATED'\" | bcftools view -e QUAL==0 | bcftools sort > {output.vcf};" 
 		"bgzip -c {output.vcf} > {output.vcfgz}; "
 		"tabix -p vcf {output.vcfgz}; "