feat: stratification of VAFs

.test/config/config.yaml

@@ -13,6 +13,9 @@ variant-calls:
     path: "resources/variants/seqc2-somatic/all.truth.format-added.vcf.gz"
     benchmark: "seqc2-somatic-ea"
     tumor_sample_name: "truth"
+    vaf-field: 
+      - INFO
+      - TVAF
 
 custom-benchmarks:
   my-custom-benchmark: # custom benchmark name, choose freely (no whitespace allowed)

config/config.yaml

@@ -17,6 +17,7 @@ variant-calls:
     # Uncomment and point to file containing contig name replacements
     # as needed for 'bcftools annotate --rename-chrs'
     # rename-contigs: path/to/rename-contigs.txt
+    vaf-field: FORMAT/AF # needs to be checked with bcftools view -h
 
 # Optional custom benchmarks
 custom-benchmarks:

workflow/resources/datavzrd/precision-recall-config.yte.yaml

@@ -67,47 +67,108 @@ views:
 
   results-table:
     dataset: results
-    desc: |
-      Precision and recall are calculated by matching variants between each callset 
-      and truth, stratified by coverage categories. The matching ignores genotype
-      differences. Instead, genotype mismatches are displayed in the "genotype mismatch rate" 
-      column.
-    page-size: 12
-    render-table:
-      columns:
-        callset:
-          plot:
-            heatmap:
-              scale: ordinal
-              color-scheme: category20
-        precision:
-          precision: 3
-          plot:
-            ticks:
-              scale: linear
-        recall:
-          precision: 3
-          plot:
-            ticks:
-              scale: linear
-        genotype_mismatch_rate:
-          ?if params.somatic:
-            display-mode: hidden
-          ?else:
-            display-mode: normal
-          precision: 3
-          plot:
-            ticks:
-              scale: linear
-        coverage:
-          plot:
-            heatmap:
-              scale: ordinal
-              domain:
-                - low
-                - medium
-                - high
-              range:
-                - "#c6dbef"
-                - "#9ecae1"
-                - "#6baed6"
+    ?if params.vaf:
+      desc: |
+        Precision and recall are calculated by matching variants between each callset 
+        and truth, stratified by coverage categories. Stratified by VAF.
+      page-size: 12
+      render-table:
+        columns:
+          callset:
+            plot:
+              heatmap:
+                scale: ordinal
+                color-scheme: category20
+          precision:
+            precision: 3
+            plot:
+              ticks:
+                scale: linear
+          recall:
+            precision: 3
+            plot:
+              ticks:
+                scale: linear
+          coverage:
+            plot:
+              heatmap:
+                scale: ordinal
+                domain:
+                  - low
+                  - medium
+                  - high
+                range:
+                  - "#c6dbef"
+                  - "#9ecae1"
+                  - "#6baed6"
+          vaf:
+            plot:
+             heatmap:
+                scale: linear
+                domain:
+                  - 0.1
+                  - 0.2
+                  - 0.3
+                  - 0.4
+                  - 0.5
+                  - 0.6
+                  - 0.7
+                  - 0.8
+                  - 0.9
+                  - 1.0
+                range:
+                  - "#c0e6baff"
+                  - "#abdda5ff"
+                  - "#94d391ff"
+                  - "#7bc77dff"
+                  - "#60ba6cff"
+                  - "#46ab5eff"
+                  - "#329a51ff"
+                  - "#208943ff"
+                  - "#0e7735ff"
+                  - "#1a833fff"
+    ?else:
+      desc: |
+        Precision and recall are calculated by matching variants between each callset 
+        and truth, stratified by coverage categories. The matching ignores genotype
+        differences. Instead, genotype mismatches are displayed in the "genotype mismatch rate" 
+        column.
+      page-size: 12
+      render-table:
+        columns:
+          callset:
+            plot:
+              heatmap:
+                scale: ordinal
+                color-scheme: category20
+          precision:
+            precision: 3
+            plot:
+              ticks:
+                scale: linear
+          recall:
+            precision: 3
+            plot:
+              ticks:
+                scale: linear
+          genotype_mismatch_rate:
+            ?if params.somatic:
+              display-mode: hidden
+            ?else:
+              display-mode: normal
+            precision: 3
+            plot:
+              ticks:
+                scale: linear
+          coverage:
+            plot:
+              heatmap:
+                scale: ordinal
+                domain:
+                  - low
+                  - medium
+                  - high
+                range:
+                  - "#c6dbef"
+                  - "#9ecae1"
+                  - "#6baed6"

workflow/resources/presets.yaml

@@ -15,6 +15,9 @@ benchmarks:
     bam-url: https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/seqc/Somatic_Mutation_WG/data/WES/WES_EA_T_1.bwa.dedup.bam
     target-regions: https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/seqc/Somatic_Mutation_WG/technical/reference_genome/Exome_Target_bed/S07604624_Covered_human_all_v6_plus_UTR.liftover.to.hg38.bed6.gz
     grch37: false
+    vaf-field: 
+      - INFO # either FORMAT or INFO
+      - TVAF # name of tumor variant allele frequency
 
   imgag-somatic-5perc:
     genome: na12878-somatic

workflow/rules/common.smk

@@ -378,6 +378,24 @@ def get_somatic_flag(wildcards):
     return somatic_flag
 
 
+def get_vaf_fields(wildcards):
+    vaf_callset = config["variant-calls"][wildcards.callset].get("vaf-field")
+
+    benchmark = config["variant-calls"][wildcards.callset]["benchmark"]
+    vaf_benchmark = benchmarks[benchmark].get("vaf-field")
+
+    # can return (None, None) if param not set
+    return (vaf_callset, vaf_benchmark)
+
+
+def get_vaf_status(wildcards):
+    vaf_benchmark = benchmarks[wildcards.benchmark].get("vaf-field")
+    if vaf_benchmark is None:
+        return False
+    else:
+        return True
+
+
 def get_collect_stratifications_input(wildcards):
     import json
 

workflow/rules/eval.smk

@@ -169,10 +169,14 @@ rule calc_precision_recall:
         calls="results/vcfeval/{callset}/{cov}/output.vcf.gz",
         idx="results/vcfeval/{callset}/{cov}/output.vcf.gz.tbi",
         common_src=common_src,
+        truth=get_stratified_truth(),
+        query="results/stratified-variants/{callset}/{cov}.vcf.gz",
     output:
         snvs="results/precision-recall/callsets/{callset}/{cov}.{vartype}.tsv",
     log:
         "logs/calc-precision-recall/{callset}/{cov}/{vartype}.log",
+    params:
+        vaf_fields=get_vaf_fields,
     conda:
         "../envs/pysam.yaml"
     script:
@@ -206,6 +210,7 @@ rule collect_precision_recall:
     params:
         callsets=lambda w: get_benchmark_callsets(w.benchmark),
         labels=get_collect_precision_recall_labels,
+        vaf=get_vaf_status,
     log:
         "logs/collect-precision-recall/{benchmark}/{vartype}.log",
     conda:
@@ -222,6 +227,7 @@ rule render_precision_recall_report_config:
         ),
     params:
         somatic=get_somatic_status,
+        vaf=get_vaf_status,
     output:
         "results/datavzrd-config/precision-recall/{benchmark}/{vartype}.config.yaml",
     log: