Detect intergenic variants using SnpEff results

Author: ahmig

workflow/rules/vaf.smk

@@ -198,6 +198,7 @@ rule merge_annotation:
     conda: "../envs/renv.yaml"
     params:
         ref_name = config["ALIGNMENT_REFERENCE"],
+        snpeff_columns = config["ANNOTATION"]["SNPEFF_COLS"],
     input:
         tsv = OUTDIR/"vaf"/"{sample}.masked.prefiltered.tsv",
         annot = OUTDIR/"vaf"/"{sample}.vcf_fields.longer.tsv",

workflow/scripts/merge_annotation.R

@@ -38,7 +38,10 @@ variants <- read_tsv(
 log_info("Reading annotation table")
 annotation <- read_tsv(
   snakemake@input$annot,
-  col_select = c("CHROM", "POS", "REF", "ALT", "VARIANT_NAME"),
+  col_select = c(
+    snakemake@params$snpeff_columns,
+    "VARIANT_NAME"
+  ),
   col_types = list(
     CHROM = col_character(),
     POS = col_integer(),

workflow/scripts/report/NV_description.R

@@ -42,7 +42,7 @@ empty_vcf <- tibble(
   REGION = as.character(NA),
   VARIANT_NAME = as.character(NA),
   ALT_FREQ = as.numeric(NA),
-  GB_FEATURE = as.character(NA),
+  EFFECT = as.character(NA),
   SYNONYMOUS = as.character(NA),
   POS = as.numeric(NA),
   ALT = as.character(NA),
@@ -57,7 +57,7 @@ vcf <- vcf %>%
     REGION,
     VARIANT_NAME,
     ALT_FREQ,
-    GB_FEATURE,
+    EFFECT,
     SYNONYMOUS,
     POS,
     ALT
@@ -93,7 +93,7 @@ vcf <- vcf %>%
       TRUE ~ "SNP"
     ),
     Class = case_when(
-      GB_FEATURE == "Intergenic" ~ "Intergenic",
+      EFFECT == "intergenic_region" ~ "Intergenic",
       TRUE ~ SYNONYMOUS
     ),
     POS = as.numeric(POS)
@@ -104,7 +104,7 @@ vcf <- vcf %>%
       NV_class == "INDEL" ~ str_length(ALT) - 1
     ),
     indel_class = case_when(
-      GB_FEATURE == "Intergenic" ~ "Intergenic",
+      EFFECT == "intergenic_region" ~ "Intergenic",
       NV_class == "INDEL" &
         indel_len %% 3 == 0 ~
         "In frame",
@@ -116,7 +116,7 @@ vcf <- vcf %>%
   ungroup() %>%
   mutate(
     group = case_when(
-      GB_FEATURE == "Intergenic" ~ "Intergenic",
+      EFFECT == "intergenic_region" ~ "Intergenic",
       NV_class == "SNP" ~ Class,
       NV_class == "INDEL" ~ indel_class
     )
@@ -446,7 +446,6 @@ window %>%
 # Heterozygous sites per sample table
 vcf_snp %>%
   filter(ALT_FREQ <= snakemake@params$max_alt_freq) %>%
-  select(!GB_FEATURE) %>%
   left_join(
     metadata,
     by = c("SAMPLE" = "ID")

workflow/scripts/report/evo_plots.R

@@ -25,7 +25,6 @@ vcf <- vcf %>%
     VARIANT_NAME,
     REGION,
     ALT_FREQ,
-    GB_FEATURE,
     SYNONYMOUS,
     POS
   )