Decouple dN/dS plots from data

Author: ahmig

config/design_plots.R

@@ -93,11 +93,11 @@ DNDS_LABELS <- c(
 )
 
 DNDS_COLORS <- c(
-  dn = "#E53E47",
-  ds = "#2C47F5"
+  dN = "#E53E47",
+  dS = "#2C47F5"
 )
 
 DNDS_SHAPES <- c(
-  dn = 2,
-  ds = 4
+  dN = 2,
+  dS = 4
 )

workflow/rules/evolution.smk

@@ -14,3 +14,17 @@ rule N_S_sites:
         LOGDIR / "N_S_sites" / "log.txt"
     script:
         "../scripts/N_S_sites_from_fasta.py"
+
+
+rule dnds_data:
+    conda: "../envs/renv.yaml"
+    input: 
+        n_s_sites = OUTDIR/f"{OUTPUT_NAME}.ancestor.N_S.sites.csv",
+        variants =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
+        metadata = config["METADATA"]
+    output:
+        table = report(REPORT_DIR_TABLES/"dnds.csv")
+    log:
+        LOGDIR / "dnds_data" / "log.txt"
+    script:
+        "../scripts/dnds_data.R"

workflow/rules/report.smk

@@ -225,22 +225,19 @@ rule time_signal_plot:
         "../scripts/report/time_signal_plot.R"
 
 
-rule evo_plots:
+rule dnds_plots:
     conda: "../envs/renv.yaml"
     params: 
         design = config["PLOTS"]
     input: 
-        n_s_sites = OUTDIR/f"{OUTPUT_NAME}.ancestor.N_S.sites.csv",
-        vcf =  OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
-        metadata = config["METADATA"]
+        table = REPORT_DIR_TABLES/"dnds.csv",
     output:
-        plot = report(REPORT_DIR_PLOTS/"figure_11.png"),
-        plot_omega = report(REPORT_DIR_PLOTS/"figure_12.png"),
-        table = report(REPORT_DIR_TABLES/"figure_11.csv")
+        plot_dn_ds = report(REPORT_DIR_PLOTS/"dn_and_ds.png"),
+        plot_omega = report(REPORT_DIR_PLOTS/"dnds.png"),
     log:
         LOGDIR / "evo_plots" / "log.txt"
     script:
-        "../scripts/report/evo_plots.R"
+        "../scripts/report/dnds_plots.R"
 
 
 rule snp_plots:
@@ -292,8 +289,8 @@ rule report:
         tree       = report(REPORT_DIR_PLOTS/"allele_freq_tree.png"),
         temest     = report(REPORT_DIR_PLOTS/"time_signal.png"),
         heat_table = report(REPORT_DIR_TABLES/"heatmap.csv"),
-        evo        = report(REPORT_DIR_PLOTS/"figure_11.png"),
-        omega_plot = report(REPORT_DIR_PLOTS/"figure_12.png"),
+        evo        = report(REPORT_DIR_PLOTS/"dn_and_ds.png"),
+        omega_plot = report(REPORT_DIR_PLOTS/"dnds.png"),
         freyja_ts  = OUTDIR/"demixing"/"freyja_data"/"last_barcode_update.txt",
         value      = REPORT_DIR_TABLES/"diversity.json",
         stats_lm   = REPORT_DIR_TABLES/"time_signal.json",

workflow/scripts/dnds_data.R

@@ -0,0 +1,68 @@
+#!/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)
+
+# Read inputs
+log_info("Reading variants table")
+variants <- read_delim(
+    snakemake@input[["variants"]],
+    col_select = c(
+        "SAMPLE",
+        "VARIANT_NAME",
+        "REGION",
+        "ALT_FREQ",
+        "SYNONYMOUS",
+        "POS"
+    )
+)
+
+log_info("Reading metadata table")
+metadata <- read_delim(snakemake@input[["metadata"]]) %>%
+    mutate(
+        interval = as.numeric(
+            as.Date(CollectionDate) - min(as.Date(CollectionDate))
+        )
+    ) %>%
+    select(ID, interval) %>%
+    rename(SAMPLE = ID)
+
+log_debug("Adding metadata to variants table")
+variants <- left_join(variants, metadata)
+
+log_info("Reading N/S sites")
+N_S_position <- read_delim(snakemake@input[["n_s_sites"]])
+
+log_info("Computing dN/dS over time (NG86)")
+dn.ds <- variants %>%
+    group_by(SAMPLE, SYNONYMOUS) %>%
+    summarise(
+        Freq = sum(ALT_FREQ, na.rm = TRUE)
+    ) %>%
+    pivot_wider(
+        names_from = SYNONYMOUS,
+        values_from = Freq,
+        values_fill = 0
+    ) %>%
+    transmute(
+        dn = No / sum(N_S_position$N),
+        ds = Yes / sum(N_S_position$S)
+    ) %>%
+    ungroup() %>%
+    left_join(unique(select(variants, SAMPLE, interval))) %>%
+    transmute(
+        sample = SAMPLE,
+        day = interval,
+        dN = dn,
+        dS = ds,
+        w = dn / ds
+    )
+
+log_info("Writing results")
+write_csv(dn.ds, snakemake@output[["table"]])

workflow/scripts/report/dnds_plots.R

@@ -0,0 +1,85 @@
+#!/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)
+
+# Import file with plots style
+source(snakemake@params[["design"]])
+
+# Read inputs
+log_info("Reading plot data")
+plot.df <- read_delim(snakemake@input$table) %>%
+  filter(w != Inf) %>%
+  pivot_longer(
+    c("dN", "dS", "w"),
+    values_to = "value",
+    names_to = "ratio"
+  )
+
+log_info("Plotting dN and dS")
+p.dn.ds <- plot.df %>%
+  filter(ratio != "w") %>%
+  ggplot() +
+  aes(
+    x = day,
+    y = value,
+    color = ratio,
+    shape = ratio
+  ) +
+  geom_point(size = 3) +
+  geom_line() +
+  scale_color_manual(
+    name = "",
+    labels = DNDS_LABELS,
+    values = DNDS_COLORS
+  ) +
+  scale_shape_manual(
+    name = "",
+    values = DNDS_SHAPES,
+    labels = DNDS_LABELS
+  ) +
+  labs(
+    y = "Substitution rate",
+    x = "Days since the initial sampling",
+    color = ""
+  )
+
+ggsave(
+  filename = snakemake@output$plot_dn_ds,
+  plot = p.dn.ds,
+  width = 159.2,
+  height = 119.4,
+  units = "mm",
+  dpi = 250
+)
+
+log_info("Plotting w (dN/dS)")
+p.omega <- plot.df %>%
+  filter(ratio == "w") %>%
+  ggplot() +
+  aes(
+    x = day,
+    y = value,
+  ) +
+  geom_point(color = "black") +
+  geom_line(color = "black") +
+  labs(
+    y = "w (dN/dS)",
+    x = "Days since the initial sampling",
+    color = ""
+  )
+
+ggsave(
+  filename = snakemake@output$plot_omega,
+  plot = p.omega,
+  width = 159.2,
+  height = 119.4,
+  units = "mm",
+  dpi = 250
+)