Filter out errors in SnpEff

Also replaces empty values ("" in config YAML) with NA in R
and adds code comments

Author: ahmig

config/config.yaml

@@ -35,7 +35,9 @@ ANNOTATION:
     FEATUREID: "ANN[*].FEATUREID"
     HGVS_P: "ANN[*].HGVS_P"
     HGVS_C: "ANN[*].HGVS_C"
+    ERRORS: "ANN[*].ERRORS"
   FILTER_INCLUDE:  # see: https://pcingola.github.io/SnpEff/adds/VCFannotationformat_v1.0.pdf
+    ERRORS: [""]
   #  IMPACT: [HIGH, MODERATE, LOW]
   FILTER_EXCLUDE:
     EFFECT: [upstream_gene_variant, downstream_gene_variant]

workflow/scripts/format_vcf_fields_longer.R

@@ -7,22 +7,45 @@ sink(log, type = "output")
 
 library(tidyverse)
 
+# Replace "" values with NA in R filter list
+# Snakemake passes filters like: list(ERRORS = c(""))
+empty.to.na <- function(x) {
+    x[x == ""] <- NA
+    x
+}
+filter.include <- lapply(snakemake@params$filter_include, empty.to.na)
+filter.exclude <- lapply(snakemake@params$filter_exclude, empty.to.na)
+
+# Process input table
 read_tsv(snakemake@input$tsv) %>%
+    # Separate <sep>-delimited "...[*]..." columns (e.g. ANN[*].EFFECT)
     separate_longer_delim(contains("[*]"), delim = snakemake@params$sep) %>%
+    # Replace empty ("") fields with NA
     mutate(across(contains("[*]"), ~ na_if(., ""))) %>%
+    # Rename "...[*]..." columns using the provided lookup via Snakemake config
     rename(all_of(unlist(snakemake@params$colnames_mapping))) %>%
+    # Apply dynamic filters from the Snakemake config:
+    # map2 pairs column names (.x) with value vectors (.y) and builds boolean expressions.
+    # Inside the expr call, !! injects a single value into each expression.
+    # The resulting list of expressions is spliced with !!! so each becomes its
+    # own condition as if written directly inside the filter call.
     filter(
+        # Keep variants that include the required values in each defined field (e.g. empty ERRORS)
         !!!map2(
-            names(snakemake@params$filter_include),
-            snakemake@params$filter_include,
+            names(filter.include),
+            filter.include,
             ~ expr(.data[[!!.x]] %in% !!.y)
         ),
+        # Keep variants that exclude the required values in each defined field  (e.g. EFFECT != "upstream_gene_variant")
         !!!map2(
-            names(snakemake@params$filter_exclude),
-            snakemake@params$filter_exclude,
+            names(filter.exclude),
+            filter.exclude,
             ~ expr(!(.data[[!!.x]] %in% !!.y))
         )
     ) %>%
+    # Keep unique rows
     distinct() %>%
+    # Assign variant name using the pattern defined via Snakemake config
     mutate(VARIANT_NAME = str_glue(snakemake@params$variant_name_pattern)) %>%
+    # Write output file
     write_tsv(snakemake@output$tsv)