Format extracted VCF fields from the SnpEff annotation

Author: ahmig

config/config.yaml

@@ -20,6 +20,26 @@ VC:
   IVAR_QUALITY: 20
   IVAR_FREQ: 0.05
   IVAR_DEPTH: 30
+ANNOTATION:
+  SNPEFF_COLS:
+    CHROM: CHROM
+    POS: POS
+    REF: REF
+    ALT: ALT
+    EFFECT: "ANN[*].EFFECT"
+    IMPACT: "ANN[*].IMPACT"
+    BIOTYPE: "ANN[*].BIOTYPE"
+    GENE: "ANN[*].GENE"
+    GENEID: "ANN[*].GENEID"
+    FEATURE: "ANN[*].FEATURE"
+    FEATUREID: "ANN[*].FEATUREID"
+    HGVS_P: "ANN[*].HGVS_P"
+    HGVS_C: "ANN[*].HGVS_C"
+  FILTER_INCLUDE:  # see: https://pcingola.github.io/SnpEff/adds/VCFannotationformat_v1.0.pdf
+  #  IMPACT: [HIGH, MODERATE, LOW]
+  FILTER_EXCLUDE:
+  #  EFFECT: [upstream_gene_variant, downstream_gene_variant]
+  VARIANT_NAME_PATTERN: "{GENE}:{coalesce(HGVS_P, HGVS_C)}"  # coalesce finds the first non-missing element
 DEMIX:
   PATHOGEN: "SARS-CoV-2"
   MIN_QUALITY: 20

workflow/rules/vaf.smk

@@ -151,12 +151,8 @@ rule extract_vcf_fields:
     threads: 1
     conda: "../envs/snpeff.yaml"
     params:
-        extract_columns = [
-            "CHROM", "POS", "REF", "ALT",
-            '"ANN[*].EFFECT"', '"ANN[*].IMPACT"', '"ANN[*].BIOTYPE"',
-            '"ANN[*].GENE"', '"ANN[*].GENEID"', '"ANN[*].FEATURE"', '"ANN[*].FEATUREID"', '"ANN[*].HGVS_P"', '"ANN[*].HGVS_C"'
-        ],
-        sep = ","
+        extract_columns = [f"'{col}'" for col in config["ANNOTATION"]["SNPEFF_COLS"].values()],
+        sep = ",",
     input:
         vcf = OUTDIR/"vaf"/"{sample}.annotated.vcf"
     output:
@@ -170,7 +166,10 @@ rule extract_vcf_fields:
 rule format_vcf_fields_longer:
     conda: "../envs/renv.yaml"
     params:
-        sep = ","
+        colnames_mapping = config["ANNOTATION"]["SNPEFF_COLS"],
+        filter_include = config["ANNOTATION"]["FILTER_INCLUDE"],
+        variant_name_pattern = config["ANNOTATION"]["VARIANT_NAME_PATTERN"],
+        sep = ",",
     input:
         tsv = OUTDIR/"vaf"/"{sample}.vcf_fields.tsv"
     output:

workflow/scripts/format_vcf_fields_longer.R

@@ -9,4 +9,19 @@ library(tidyverse)
 
 read_tsv(snakemake@input$tsv) %>%
     separate_rows(contains("[*]"), sep = snakemake@params$sep) %>%
+    rename(all_of(unlist(snakemake@params$colnames_mapping))) %>%
+    filter(
+        !!!map2(
+            names(snakemake@params$filter_include),
+            snakemake@params$filter_include,
+            ~ expr(.data[[!!.x]] %in% !!.y)
+        ),
+        !!!map2(
+            names(snakemake@params$filter_exclude),
+            snakemake@params$filter_exclude,
+            ~ expr(!(.data[[!!.x]] %in% !!.y))
+        )
+    ) %>%
+    distinct() %>%
+    mutate(VARIANT_NAME = str_glue(snakemake@params$variant_name_pattern)) %>%
     write_tsv(snakemake@output$tsv)