style: plot broad & narrow peaks separately

assets/multiqc_config.yaml

@@ -138,10 +138,14 @@ custom_data:
     parent_id: peak_qc
     id: frip_nbases
     section_name: "FRiP vs. number of bases"
-  peak_widths:
+  peak_widths_broad:
     parent_id: peak_qc
-    id: peak_widths
-    section_name: "Distribution of peak widths"
+    id: peak_widths_broad
+    section_name: "Broad peak widths"
+  peak_widths_narrow:
+    parent_id: peak_qc
+    id: peak_widths_narrow
+    section_name: "Narrow peak widths"
   jaccard_heatmap:
     parent_id: peak_qc
     id: jaccard_heatmap
@@ -166,8 +170,10 @@ sp:
     fn: "FRiP_samples*.png"
   frip_nbases:
     fn: "FRiP_nbases*.png"
-  peak_widths:
-    fn: "peak_widths*.png"
+  peak_widths_broad:
+    fn: "peak_widths_broad*.png"
+  peak_widths_narrow:
+    fn: "peak_widths_narrow*.png"
   jaccard_heatmap:
     fn: "jaccard_heatmap_all*.png"
   jaccard_pca_all:

bin/plot_frip.R

@@ -11,7 +11,12 @@ if (length(args) < 1) {
 }
 
 frip_filename <- args[1]
-frip_dat <- read_tsv(frip_filename)
+frip_dat <- read_tsv(frip_filename) %>%
+  mutate(peak_type = case_when(
+    bedtool == "sicer" | bedtool == "macs_broad" ~ "broad",
+    bedtool == "gem" | bedtool == "macs_narrow" ~ "narrow",
+    TRUE ~ NA_character_
+  ))
 
 sample_frip_plot <- frip_dat %>%
   mutate(n_overlap_reads_mil = n_overlap_reads / 10^6) %>%
@@ -23,7 +28,7 @@ sample_frip_plot <- frip_dat %>%
     names_to = "metric"
   ) %>%
   ggplot(aes(value, bedsample, color = bedtool)) +
-  facet_wrap("metric", nrow = 1, scales = "free_x", strip.position = "bottom") +
+  facet_wrap(peak_type ~ metric, scales = "free_x", strip.position = "bottom") +
   geom_jitter(alpha = 0.8, height = 0.1) +
   geom_hline(
     yintercept = seq(1.5, length(unique(frip_dat %>% pull(bedsample))) - 0.5, 1),
@@ -47,6 +52,7 @@ ggsave(filename = "FRiP_samples.png", plot = sample_frip_plot, device = "png", d
 nbasesM_frip_plot <- frip_dat %>%
   ggplot(aes(n_basesM, FRiP, color = bedsample, shape = bedtool)) +
   geom_point(alpha = 0.8) +
+  facet_wrap(~peak_type, scales = "free") +
   guides(
     color = guide_legend(title = "Sample"),
     shape = guide_legend(title = "Peak caller")

bin/plot_peak_widths.R

@@ -2,6 +2,7 @@
 library(dplyr)
 library(ggplot2)
 library(readr)
+library(scales)
 
 set.seed(20230926)
 args <- commandArgs(trailingOnly = TRUE)
@@ -10,30 +11,63 @@ if (length(args) < 1) {
 }
 tsv_filename <- args[1]
 peak_dat <- read_tsv(tsv_filename) %>%
-  mutate(peak_width = chromEnd - chromStart)
-xmin <- peak_dat %>%
-  pull(peak_width) %>%
-  min()
-xmax <- peak_dat %>%
-  pull(peak_width) %>%
-  max()
+  mutate(
+    peak_width = chromEnd - chromStart,
+    peak_type = case_when(
+      tool == "sicer" | tool == "macs_broad" ~ "broad",
+      tool == "gem" | tool == "macs_narrow" ~ "narrow",
+      TRUE ~ NA_character_
+    )
+  )
 
-peak_widths_hist <- peak_dat %>%
-  ggplot(aes(peak_width, fill = tool)) +
-  geom_histogram(alpha = 0.7, position = "identity") +
-  scale_x_log10(
-    limits = c(xmin - 10^2, xmax + 10^2),
-    labels = scales::label_log(digits = 2)
-  ) +
-  scale_fill_viridis_d() +
-  guides(fill = guide_legend(
-    label.position = "bottom",
-    title = "Peak caller",
-    title.position = "top"
-  )) +
-  facet_wrap("sample_id") +
-  labs(x = expression("" * log[10] * " Peak Widths")) +
-  theme_bw() +
-  theme(legend.position = "top")
+tool_colors <- viridisLite::viridis(n = 4, alpha = 0.7, begin = 0, end = 1, direction = 1, option = "D")
+names(tool_colors) <- peak_dat %>%
+  pull(tool) %>%
+  unique()
 
-ggsave(filename = "peak_widths_histogram.png", plot = peak_widths_hist, device = "png", dpi = 300, height = 4, width = 5)
+plot_hist <- function(peak_dat) {
+  xmin <- peak_dat %>%
+    pull(peak_width) %>%
+    min()
+  xmax <- peak_dat %>%
+    pull(peak_width) %>%
+    max()
+  peak_dat %>%
+    ggplot(aes(peak_width, fill = tool)) +
+    geom_histogram(alpha = 0.7, position = "identity") +
+    scale_x_log10(
+      limits = c(xmin - 10^2, xmax + 10^2),
+      labels = label_log(digits = 2)
+    ) +
+    scale_fill_manual(
+      values = tool_colors,
+      breaks = names(tool_colors)
+    ) +
+    guides(fill = guide_legend(
+      label.position = "bottom",
+      title = "Peak caller",
+      title.position = "top"
+    )) +
+    facet_wrap(~sample_id) +
+    labs(x = expression("" * log[10] * " Peak Widths")) +
+    theme_bw() +
+    theme(legend.position = "top")
+}
+
+hist_broad <- peak_dat %>%
+  filter(peak_type == "broad") %>%
+  plot_hist() +
+  scale_y_log10(labels = label_log(digits = 2))
+hist_narrow <- peak_dat %>%
+  filter(peak_type == "narrow") %>%
+  plot_hist() +
+  scale_y_log10(labels = label_log(digits = 2))
+
+ggsave(
+  filename = "peak_widths_broad_histogram.png", plot = hist_broad,
+  device = "png", dpi = 300, height = 4, width = 5
+)
+ggsave(
+  filename = "peak_widths_narrow_histogram.png", plot = hist_narrow,
+  device = "png", dpi = 300, height = 4, width = 5
+)