Get variant name for plots from SnpEff annotation

Author: ahmig

workflow/rules/distances.smk

@@ -5,7 +5,7 @@ rule weighted_distances:
         samples = expand("{sample}", sample = iter_samples()),
         mask_class = ["mask"]
     input:
-        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv",
+        tsv = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         vcf = lambda wildcards: select_problematic_vcf(),
         ancestor = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta",
         reference = OUTDIR/"reference.fasta"

workflow/rules/report.smk

@@ -30,7 +30,7 @@ rule demix_plot:
 rule heatmap_plot_data:
     conda: "../envs/renv.yaml"
     input:
-        vcf =  OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv",
+        vcf =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         metadata = config["METADATA"]
     output:
         table = report(REPORT_DIR_TABLES/"heatmap.csv")
@@ -46,7 +46,7 @@ rule window:
         window = config["WINDOW"]["WIDTH"],
         step = config["WINDOW"]["STEP"]
     input:
-        vcf = OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv",
+        vcf = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         gb = OUTDIR/"reference.gb",
         features = config["FEATURES_JSON"]
     output:
@@ -105,7 +105,7 @@ rule general_NV_description:
         coordinates = config["COORDINATES_JSON"],
         regions = config["PLOT_GENOME_REGIONS"],
         window = OUTDIR/f"{OUTPUT_NAME}.window.csv",
-        vcf =  OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv",
+        vcf =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         metadata = config["METADATA"]
     output:
         fig = report(REPORT_DIR_PLOTS/"figure_5a.png"),
@@ -230,7 +230,7 @@ rule evo_plots:
         design = config["PLOTS"]
     input: 
         N_S = OUTDIR/f"{OUTPUT_NAME}.ancestor.N_S.sites.csv",
-        vcf =  OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv",
+        vcf =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         metadata = config["METADATA"]
     output:
         plot = report(REPORT_DIR_PLOTS/"figure_11.png"),
@@ -249,7 +249,7 @@ rule snp_plots:
         cor_method = config["COR"]["METHOD"],
         cor_exact = config["COR"]["EXACT"]
     input:
-        vcf =  OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv",
+        vcf =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         metadata = config["METADATA"]
     output:
         pseudovolcano = report(REPORT_DIR_PLOTS/"figure_6.png"),

workflow/rules/vaf.smk

@@ -193,26 +193,26 @@ rule compile_vcf_fields_longer:
         "../scripts/compile_vcf_fields_longer.R"
 
 
-rule vcf_to_tsv:
+rule merge_annotation:
     threads: 1
     conda: "../envs/renv.yaml"
     input:
-        ann_vcf = OUTDIR/"vaf"/"{sample}.annotated.vcf",
-        pre_tsv = OUTDIR/"vaf"/"{sample}.masked.prefiltered.tsv"
+        tsv = OUTDIR/"vaf"/"{sample}.masked.tsv",
+        annot = OUTDIR/"vaf"/"{sample}.vcf_fields.longer.tsv",
     output:
-        tsv = OUTDIR/"vaf"/"{sample}.masked.filtered.tsv"
+        tsv = OUTDIR/"vaf"/"{sample}.variants.tsv"
     log:
-        LOGDIR / "vcf_to_tsv" / "{sample}.log.txt"
+        LOGDIR / "merge_annotation" / "{sample}.log.txt"
     script:
-        "../scripts/vcf_to_tsv.R"
+        "../scripts/merge_annotation.R"
 
 
 rule compile_variants:
     threads: 1
     conda: "../envs/renv.yaml"
-    input: expand(OUTDIR/"vaf"/"{sample}.masked.filtered.tsv", sample=iter_samples())
+    input: expand(OUTDIR/"vaf"/"{sample}.variants.tsv", sample=iter_samples())
     output:
-        tsv = OUTDIR/f"{OUTPUT_NAME}.masked.filtered.tsv"
+        tsv = OUTDIR/f"{OUTPUT_NAME}.variants.tsv"
     log:
         LOGDIR / "compile_variants" / "log.txt"
     script:

workflow/scripts/compile_variants.R

@@ -14,7 +14,7 @@ tables <- lapply(
     snakemake@input,
     function(path) {
         read_tsv(path) %>%
-            mutate(SAMPLE = sub('\\.masked.filtered.tsv$', '', basename(path)))
+            mutate(SAMPLE = sub('\\.variants.tsv$', '', basename(path)))
     }
 )
 

workflow/scripts/merge_annotation.R

@@ -0,0 +1,30 @@
+#!/usr/bin/env Rscript
+
+# Write stdout and stderr to log file
+log <- file(snakemake@log[[1]], open = "wt")
+sink(log, type = "message")
+sink(log, type = "output")
+
+library(tidyverse)
+library(logger)
+log_threshold(INFO)
+
+log_info("Reading tables")
+variants <- read_tsv(snakemake@input$tsv)
+annotation <- read_tsv(
+  snakemake@input$annot,
+  col_select = c("CHROM", "POS", "REF", "ALT", "VARIANT_NAME")
+)
+
+log_info("Merging tables")
+merged <- left_join(
+  variants,
+  annotation,
+  by = c("CHROM", "POS", "REF", "ALT")
+)
+
+log_info("Saving results")
+write_tsv(
+  tsv,
+  snakemake@output$tsv
+)

workflow/scripts/report/NV_description.R

@@ -40,7 +40,7 @@ date_order <- metadata %>%
 empty_vcf <- tibble(
   SAMPLE = date_order,
   REGION = as.character(NA),
-  variant = as.character(NA),
+  VARIANT_NAME = as.character(NA),
   ALT_FREQ = as.numeric(NA),
   GB_FEATURE = as.character(NA),
   synonimous = as.character(NA),
@@ -55,7 +55,7 @@ vcf <- vcf %>%
   dplyr::select(
     SAMPLE,
     REGION,
-    variant,
+    VARIANT_NAME,
     ALT_FREQ,
     GB_FEATURE,
     synonimous,
@@ -414,14 +414,14 @@ vcf %>%
   select(
     SAMPLE,
     POS,
-    variant,
+    VARIANT_NAME,
     ALT_FREQ,
     NV_class,
     group
   ) %>%
   rename(
     sample = SAMPLE,
-    Variant = variant,
+    Variant = VARIANT_NAME,
     Class = group
   ) %>%
   filter(ALT_FREQ > 0) %>%
@@ -464,11 +464,11 @@ vcf_snp %>%
 # Stats for reporting
 n_indels <- vcf %>%
   filter(NV_class == "INDEL") %>%
-  pull(variant) %>%
+  pull(VARIANT_NAME) %>%
   unique() %>%
   length()
 
-n_snv <- length(unique(vcf$variant)) - n_indels
+n_snv <- length(unique(vcf$VARIANT_NAME)) - n_indels
 model <- lm(n ~ CollectionDate, data = figur_SNP_table)
 
 # Calculate correlation, if possible

workflow/scripts/report/evo_plots.R

@@ -22,7 +22,7 @@ N_S_position <- read_delim(snakemake@input[["N_S"]])
 vcf <- vcf %>%
   dplyr::select(
     SAMPLE,
-    variant,
+    VARIANT_NAME,
     REGION,
     ALT_FREQ,
     GB_FEATURE,

workflow/scripts/report/heatmap.R

@@ -19,11 +19,11 @@ date_order <- read_csv(snakemake@input[["metadata"]]) %>%
 
 # Create SNP variable and select useful variables from vcf
 vcf <- vcf %>%
-  dplyr::select(variant, SAMPLE, ALT_FREQ)
+  dplyr::select(VARIANT_NAME, SAMPLE, ALT_FREQ)
 
 vcf <- vcf %>%
   pivot_wider(
-    names_from = variant,
+    names_from = VARIANT_NAME,
     values_from = ALT_FREQ,
     values_fill = 0,
     values_fn = sum

workflow/scripts/report/snp_plots.R

@@ -43,7 +43,7 @@ date_order <- metadata %>%
 # Simplify features names and create SNP variable
 vcf <- vcf %>%
   dplyr::select(
-    variant,
+    VARIANT_NAME,
     SAMPLE,
     REGION,
     ALT_FREQ,
@@ -52,7 +52,7 @@ vcf <- vcf %>%
 
 # Fill positions without alt frequency with 0
 vcf <- vcf %>%
-  complete(nesting(variant, POS), SAMPLE, fill = list(ALT_FREQ = 0))
+  complete(nesting(VARIANT_NAME, POS), SAMPLE, fill = list(ALT_FREQ = 0))
 
 # Join variants file and metadata file
 vcf <- left_join(vcf, data, by = c("SAMPLE" = "ID"))
@@ -72,7 +72,7 @@ vcf <- arrange(
 # Get list with all different polymorphisms
 SNPs <- pull(
   vcf,
-  variant
+  VARIANT_NAME
 ) %>%
   unique()
 
@@ -89,7 +89,7 @@ cor.df.fill <- lapply(
   function(snp) {
     df <- filter(
       vcf,
-      variant == snp
+      VARIANT_NAME == snp
     )
 
     test <- cor.test(
@@ -161,14 +161,14 @@ sign <- filter(
 # SNPs which are in positions with more than one alternative allele
 dup <- vcf %>%
   select(
-    variant,
+    VARIANT_NAME,
     POS
   ) %>%
   unique() %>%
   group_by(POS) %>%
   filter(n() > 1) %>%
   ungroup() %>%
-  pull(variant) %>%
+  pull(VARIANT_NAME) %>%
   unique()
 
 subset <- c(sign, dup) %>%
@@ -179,12 +179,12 @@ plot.height <- ceiling(length(subset) / 4) * 42
 
 log_info("Plotting SNPs trends in time")
 panel <- vcf %>%
-  filter(variant %in% subset) %>%
+  filter(VARIANT_NAME %in% subset) %>%
   ggplot() +
   aes(
     x = interval,
     y = ALT_FREQ,
-    color = variant
+    color = VARIANT_NAME
   ) +
   scale_color_manual(values = sample(color, length(subset))) +
   geom_point() +
@@ -234,11 +234,11 @@ cor.df.fill %>%
 
 log_info("Saving SNPs trends table")
 vcf %>%
-  filter(variant %in% subset) %>%
+  filter(VARIANT_NAME %in% subset) %>%
   transmute(
     sample = SAMPLE,
     POS = POS,
-    NV = variant,
+    NV = VARIANT_NAME,
     ALT_FREQ = ALT_FREQ,
     DaysSinceFirst = interval
   ) %>%