Decouple demix plot from data

Fix plot design source

Update sample sorting

Update demix table column names

Fix column names

Update temp variable name

Update structure

Author: ahmig

workflow/rules/demix.smk

@@ -1,4 +1,4 @@
-rule demix_update:
+rule demix_barcode_update:
     threads: 1
     shadow: "shallow"
     conda:
@@ -14,7 +14,7 @@ rule demix_update:
         pathogens = OUTDIR/"demixing"/"freyja_data"/"pathogen_config.yml",
         usher_barcodes = OUTDIR/"demixing"/"freyja_data"/"usher_barcodes.feather"
     log:
-        LOGDIR / "demix_update" / "log.txt"
+        LOGDIR / "demix_barcode_update" / "log.txt"
     shell:
         "mkdir -p {output.folder:q} && "
         "freyja update --outdir {output.folder:q} --pathogen {params.pathogen:q} >{log} 2>&1"
@@ -63,7 +63,7 @@ rule demix:
         relaxed_mrca_thresh = config["DEMIX"]["RELAXED_MRCA_THRESH"],
         solver = config["DEMIX"]["SOLVER"],
     output:
-        demix_file = OUTDIR/"demixing"/"{sample}/{sample}_demixed.tsv"
+        demix_file = OUTDIR/"demixing"/"samples"/"{sample}/{sample}_demixed.tsv"
     log:
         LOGDIR / "demix" / "{sample}.log.txt"
     shell:
@@ -91,15 +91,15 @@ rule demix:
         ">{log} 2>&1"
 
 
-rule summarise_demixing:
+rule summarise_demix:
     threads: 1
     conda: "../envs/renv.yaml"
     shadow: "shallow"
     input:
-        tables = expand(OUTDIR/"demixing"/"{sample}/{sample}_demixed.tsv", sample=iter_samples())
+        tables = expand(OUTDIR/"demixing"/"samples"/"{sample}/{sample}_demixed.tsv", sample=iter_samples())
     output:
-        summary_df = report(OUTDIR/"summary_freyja_demixing.csv")
+        summary_df = report(OUTDIR/"demixing"/"summary.csv")
     log:
-        LOGDIR / "summarise_demixing" / "log.txt"
+        LOGDIR / "summarise_demix" / "log.txt"
     script: 
-        "../scripts/summary_demixing.R"
+        "../scripts/summarise_demix.R"

workflow/rules/report.smk

@@ -49,16 +49,27 @@ rule diversity:
         "../scripts/report/diversity_plot.R"
 
 
+rule freyja_plot_data:
+    conda: "../envs/renv.yaml"
+    input:
+        summary_demixing = OUTDIR/"demixing"/"summary.csv",
+        metadata = config["METADATA"]
+    output:
+        data = report(REPORT_DIR_TABLES/"freyja.csv")
+    log:
+        LOGDIR / "freyja_plot_data" / "log.txt"
+    script:
+        "../scripts/report/freyja_plot_data.R"
+
+
 rule freyja_plot:
     conda: "../envs/renv.yaml"
     params:
         design = config["PLOTS"]
     input:
-        summary_demixing =  OUTDIR/"summary_freyja_demixing.csv",
-        metadata = config["METADATA"]
+        data = REPORT_DIR_TABLES/"freyja.csv"
     output:
-        fig = report(REPORT_DIR_PLOTS/"figure_1.png"),
-        table = report(REPORT_DIR_TABLES/"figure_1.csv")
+        plot = report(REPORT_DIR_PLOTS/"freyja.png")
     log:
         LOGDIR / "freyja_plot" / "log.txt"
     script:
@@ -175,7 +186,7 @@ rule report:
     shadow: "shallow"
     input:
         qmd        = Path(config["REPORT_QMD"]).resolve(),
-        freyja     = report(REPORT_DIR_PLOTS/"figure_1.png"),
+        freyja     = report(REPORT_DIR_PLOTS/"freyja.png"),
         tree_ml    = report(REPORT_DIR_PLOTS/"figure_2.png"),
         diversity  = report(REPORT_DIR_PLOTS/"figure_3.png"),
         fig_cor    = report(REPORT_DIR_PLOTS/"figure_4.png"),

workflow/scripts/report/freyja_plot.R

@@ -13,43 +13,17 @@ log_threshold(INFO)
 source(snakemake@params[["design"]])
 
 # Read inputs
-demix <- read_csv(snakemake@input[["summary_demixing"]])
-
-# DATA PROCESSING
-log_info("Obtaining main lineages")
-main_lineages <- demix %>%
-  group_by(sample) %>%
-  top_n(1, abundances) %>%
-  ungroup() %>%
-  pull(lineages) %>%
-  unique()
-
-# Obtain sample names ordered by CollectionDate
-metadata <- read_csv(snakemake@input[["metadata"]])
-date_order <- metadata %>%
-  arrange(CollectionDate) %>%
-  filter(ID %in% demix$sample) %>%
-  pull(ID) %>%
-  unique()
-
-# PLOT
-log_info("Plotting summary demixing")
-demix_plot <- demix %>%
-  mutate(
-    lineages = case_when(
-      lineages %in% main_lineages ~ lineages
-    )
-  ) %>%
-  group_by(lineages, sample) %>%
+log_info("Plotting")
+p <- read_csv(snakemake@input$data) %>%
   mutate(
-    abundances = sum(abundances)
+    sample = reorder(sample, CollectionDate),
+    lineage = ifelse(is_main, lineage, NA)
   ) %>%
-  unique() %>%
   ggplot() +
   aes(
-    x = factor(sample, date_order),
-    y = as.numeric(abundances),
-    fill = lineages
+    x = sample,
+    y = abundance,
+    fill = lineage
   ) +
   scale_fill_viridis_d(
     na.value = "gray50",
@@ -68,34 +42,12 @@ demix_plot <- demix %>%
     fill = "Lineage"
   )
 
+log_info("Saving plot")
 ggsave(
-  filename = snakemake@output[["fig"]],
-  plot = demix_plot,
+  filename = snakemake@output[["plot"]],
+  plot = p,
   width = 159.2,
   height = 119.4,
   units = "mm",
   dpi = 250
 )
-
-
-# PLOT TABLES
-log_info("Saving plot table")
-demix %>%
-  mutate(
-    lineages = case_when(
-      lineages %in% main_lineages ~ lineages,
-      TRUE ~ "Other"
-    )
-  ) %>%
-  group_by(sample, lineages) %>%
-  summarise(abundances = sum(abundances)) %>%
-  ungroup() %>%
-  left_join(
-    select(
-      metadata,
-      ID,
-      CollectionDate
-    ),
-    by = c("sample" = "ID")
-  ) %>%
-  write.csv(snakemake@output[["table"]], row.names = FALSE)

workflow/scripts/report/freyja_plot_data.R

@@ -0,0 +1,56 @@
+#!/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
+demix <- read_csv(snakemake@input[["summary_demixing"]])
+
+# Data processing
+log_info("Obtaining main lineages")
+main_lineages <- demix %>%
+  group_by(sample) %>%
+  top_n(1, abundance) %>%
+  ungroup() %>%
+  pull(lineage) %>%
+  unique()
+
+# Obtain sample names ordered by CollectionDate
+log_info("Sorting dates")
+metadata <- read_csv(snakemake@input[["metadata"]])
+date_order <- metadata %>%
+  arrange(CollectionDate) %>%
+  filter(ID %in% demix$sample) %>%
+  pull(ID) %>%
+  unique()
+
+# Build plot data
+log_info("Building plot data")
+plot.data <- demix %>%
+  group_by(lineage, sample) %>%
+  summarize(abundance = sum(abundance)) %>%
+  ungroup() %>%
+  left_join(
+    select(
+      metadata,
+      ID,
+      CollectionDate
+    ),
+    by = c("sample" = "ID")
+  ) %>%
+  mutate(
+    is_main = lineage %in% main_lineages
+  )
+
+# Save plot data
+log_info("Saving plot data")
+write_csv(
+  plot.data,
+  snakemake@output[["data"]]
+)

workflow/scripts/summary_demixing.R workflow/scripts/summarise_demix.Rworkflow/scripts/summary_demixing.R renamed to workflow/scripts/summarise_demix.R

@@ -11,8 +11,8 @@ log_threshold(INFO)
 
 # Empty dataframe to be filled with data
 demix <- data.frame(
-  "lineages" = NA,
-  "abundances" = NA,
+  "lineage" = NA,
+  "abundance" = NA,
   "sample" = NA
 ) %>%
   filter(!is.na(sample))
@@ -23,21 +23,25 @@ lapply(
   function(tsv_file) {
     read_tsv(
       tsv_file,
-      col_names = c("variable", "valor"),
+      col_names = c("variable", "value"),
       show_col_types = FALSE
     ) %>%
       filter(
-        row_number() %in% c(3, 4)
+        variable %in% c("lineages", "abundances")
       ) %>%
       pivot_wider(
         names_from = variable,
-        values_from = valor
+        values_from = value
       ) %>%
       separate_rows(
         lineages,
         abundances,
         sep = " "
       ) %>%
+      rename(
+        lineage = lineages,
+        abundance = abundances
+      ) %>%
       mutate(
         sample = str_extract(
           basename(tsv_file),