refactor: calculate AF correlation matrix before rendering report

Author: ahmig

template.qmd

@@ -65,6 +65,7 @@ library(tidyr)
 library(dplyr)
 library(heatmaply)
 library(readr)
+library(tibble)
 
 # Diversity
 div_values <- fromJSON(params$div_value)
@@ -89,10 +90,9 @@ table <- read.csv(params$table)
 n.samples <- table %>% pull(Sample) %>% unique() %>% length()
 
 # Heatmap
-heatmap.df <- read.csv(params$heat_tab)
-row.names(heatmap.df) <- heatmap.df$`...1`
-heatmap.df[1] <- NULL
-cor.mat <- cor(heatmap.df, method = params$cor_method, use = "pairwise.complete.obs")
+cor.mat <- read_csv(params$heat_tab) %>%
+  column_to_rownames("...1") %>%
+  as.matrix()
 if (params$cor_method == "pearson") {
   cor.method.name <- "Pearson's"
 } else if (params$cor_method == "spearman") {

workflow/rules/report.smk

@@ -394,15 +394,19 @@ rule af_trajectory_panel_plot:
 
 rule pairwise_trajectory_correlation_data:
     conda: "../envs/renv.yaml"
+    params:
+        cor_method = config["COR"]["METHOD"],
+        cor_use = "pairwise.complete.obs",
     input:
         variants = OUTDIR/f"{OUTPUT_NAME}.variants.all_sites.tsv",
         metadata = config["METADATA"],
     output:
-        table = report(REPORT_DIR_TABLES/"pairwise_trajectory_correlation.csv"),
+        table = REPORT_DIR_TABLES/"pairwise_trajectory_frequency_data.csv",
+        matrix = report(REPORT_DIR_TABLES/"pairwise_trajectory_correlation_matrix.csv"),
     log:
         LOGDIR / "pairwise_trajectory_correlation_data" / "log.txt"
     script:
-        "../scripts/pairwise_trajectory_correlation_data.R"
+        "../scripts/report/pairwise_trajectory_correlation_data.R"
 
 
 rule summary_table:
@@ -436,7 +440,7 @@ rule report:
         temest     = report(REPORT_DIR_PLOTS/"time_signal.png"),
         evo        = report(REPORT_DIR_PLOTS/"dn_and_ds.png"),
         omega_plot = report(REPORT_DIR_PLOTS/"dnds.png"),
-        heat_table = report(REPORT_DIR_TABLES/"pairwise_trajectory_correlation.csv"),
+        heat_table = report(REPORT_DIR_TABLES/"pairwise_trajectory_correlation_matrix.csv"),
         freyja_ts  = OUTDIR/"demixing"/"freyja_data"/"last_barcode_update.txt",
         value      = REPORT_DIR_TABLES/"diversity.json",
         stats_lm   = REPORT_DIR_TABLES/"time_signal.json",

…s/pairwise_trajectory_correlation_data.R …t/pairwise_trajectory_correlation_data.Rworkflow/scripts/pairwise_trajectory_correlation_data.R renamed to workflow/scripts/report/pairwise_trajectory_correlation_data.R workflow/scripts/pairwise_trajectory_correlation_data.R renamed to workflow/scripts/report/pairwise_trajectory_correlation_data.R

@@ -23,6 +23,7 @@ date_order <- read_csv(snakemake@input[["metadata"]]) %>%
   pull(ID) %>%
   unique()
 
+log_info("Formatting variants")
 all_variants_wider <- variants %>%
   select(SAMPLE, VARIANT_NAME, ALT_FREQ) %>%
   pivot_wider(
@@ -35,5 +36,15 @@ all_variants_wider <- variants %>%
   rename_with(~ str_replace(., "^([^|]+)\\|.*$", "\\1(...)"), -SAMPLE) %>%
   column_to_rownames(var = "SAMPLE")
 
-log_info("Saving table")
+log_info("Saving table of frequencies")
 write.csv(all_variants_wider, snakemake@output[["table"]])
+
+log_info("Calculating correlations")
+cor.mat <- cor(
+  all_variants_wider,
+  method = snakemake@params$cor_method,
+  use = snakemake@params$cor_use
+)
+
+log_info("Saving correlation matrix")
+write.csv(cor.mat, snakemake@output[["matrix"]])