make heatmap versions scaled based on both counts and relative abundance. Also make the cutoff zone a bit more opaque and scale it to the total length of the primers

Author: cjfields

bin/MergeCheck_Plot.R

@@ -2,16 +2,22 @@
 suppressPackageStartupMessages(library(tidyverse))
 suppressPackageStartupMessages(library(optparse))
 
-# A very simple plot script for generating a 'heatmap' checking overlaps. 
+# A very simple plot script for generating a 'heatmap' checking overlaps.
 option_list = list(
-    # make_option(c("--forward_clip"), type="character", default=NULL, help="Forward (5') trim"),
-    # make_option(c("--reverse_clip"), type="character", default=NULL, help="Reverse (5') trim"),
-    # make_option(c("--minMergedLen"), type="numeric", default=0, help="cpus"),
+  make_option("--forward",
+    type    = "character",
+    default = "",
+    help    = "Forward (5') primer"),
+  make_option("--reverse", 
+    type    = "character",
+    default = "",
+    help    = "Reverse (5') primer")
 )
 
-opt <- parse_args(OptionParser(option_list=option_list))
-len_files <- list.files(".", 
-                        pattern = "*.lengthstats.txt", 
+opt <- parse_args(OptionParser(option_list = option_list))
+
+len_files <- list.files(".",
+                        pattern    = "*.lengthstats.txt",
                         full.names = TRUE)
 
 lens_tmp <- lapply(len_files,
@@ -20,50 +26,69 @@ lens_tmp <- lapply(len_files,
 
 names(lens_tmp) <- gsub("\\S+/(\\S+).lengthstats.txt", "\\1", len_files)
 
-# bind all the data, then group by Sample, add in relative abundance per Sample and binning info, then group by Sample + Bin and summarize counts and RelAb per bin
-# It's a bit of a hack but it generally works; however it's not perfect
-lens_all <- bind_rows(lens_tmp, .id = "Sample") %>% 
-  group_by(Sample) %>% 
-  mutate(RelAb=Count/sum(Count),
-         Bin = cut(Length, 
-                   seq(min(Length), 
-                       max(Length),5), 
-                   include.lowest = TRUE)) %>% 
-  group_by(Sample, Bin) %>% 
-  mutate(ReadCountPerBin=sum(Count),
-         RelAbPerBin=sum(RelAb)) 
+# bind all the data,
+# group by Sample,
+# add in relative abundance per Sample and binning info,
+# group by Sample + Bin,
+# summarize counts and RelAb per bin.
+lens_all <- bind_rows(lens_tmp, .id = "Sample") |>
+  group_by(Sample) |>
+  mutate(RelAb = Count / sum(Count),
+         Bin   = cut(Length,
+                     seq(min(Length),
+                         max(Length), 5),
+                     include.lowest = TRUE)) |>
+  group_by(Sample, Bin) |>
+  mutate(ReadCountPerBin = sum(Count),
+         RelAbPerBin     = sum(RelAb))
 
 lens_all$Sample <- factor(lens_all$Sample)
 
 # This will become settable, but essentially anything 50nt or less is not kept
-cutoff <- 50
+cutoff <- nchar(opt$forward) + nchar(opt$reverse)
+
+if (cutoff == 0) {
+  cutoff <- 50
+}
 
-gg <- lens_all |> ggplot(aes(x=Length, y=Sample, fill=ReadCountPerBin)) + 
+cat("Cutoff is ", cutoff)
+
+gg <- lens_all |>
+  ggplot(aes(x = Length, y = Sample, fill = ReadCountPerBin)) +
   geom_tile(stat = "identity") +
-  scale_fill_viridis_c(option="plasma", direction = -1) +
+  scale_fill_viridis_c(option = "plasma", direction = -1) +
   annotate("rect",
-           xmin = 0,
-           xmax = cutoff,
-           ymin = 0.5,
-           ymax = Inf,
-           alpha=0.3, fill="blue") +
+           xmin  = 0,
+           xmax  = cutoff,
+           ymin  = 0.5,
+           ymax  = Inf,
+           alpha = 0.1,
+           fill  = "blue") +
   theme_minimal()
 
-if (nlevels(lens_all$Sample) > 40) {
-    gg <- gg + theme(axis.text.y=element_blank())
+if (nlevels(lens_all$Sample) > 50) {
+  gg <- gg + theme(axis.text.y = element_blank())
 }
 
-# if(fprimer > 0) {
-#   gg <- gg+ geom_vline(xintercept=fprimer, color = "red", alpha = 0.5)
-# }
+ggsave("MergedCheck_heatmap.counts.pdf", gg)
 
-# if(rprimer > 0) {
-#   gg <- gg+ geom_vline(xintercept=rprimer, color = "black", alpha = 0.5)
-# }
+gg2 <- lens_all |>
+  ggplot(aes(x = Length, y = Sample, fill = RelAbPerBin)) +
+  geom_tile(stat = "identity") +
+  scale_fill_viridis_c(option = "plasma", direction = -1) +
+  annotate("rect",
+           xmin  = 0,
+           xmax  = cutoff,
+           ymin  = 0.5,
+           ymax  = Inf,
+           alpha = 0.1,
+           fill  = "blue") +
+  theme_minimal()
+
+if (nlevels(lens_all$Sample) > 50) {
+  gg2 <- gg2 + theme(axis.text.y = element_blank())
+}
 
-# if(maxsizeprimers > 0) {
-#   gg <- gg+ geom_vline(xintercept=maxsizeprimers, color = "green", alpha = 0.5)
-# }
+ggsave("MergedCheck_heatmap.RelAb.pdf", gg2)
 
-ggsave("MergedCheck_heatmap.pdf")
-saveRDS(gg, "MergedCheck_heatmap.RDS")
+saveRDS(lens_all, "stats.RDS")

modules/local/overlapheatmap.nf

@@ -8,8 +8,8 @@ process OVERLAP_HEATMAP {
     path(merged_tables)
 
     output:
-    path("MergedCheck_heatmap.pdf"), emit: overlap_check_pdf
-    path("MergedCheck_heatmap.RDS"), emit: overlap_check_rds
+    path("MergedCheck_heatmap*.pdf"),     emit: overlap_check_pdf
+    path("stats.RDS"), emit: overlap_check_rds
     // path "versions.yml"           , emit: versions
 
     when:
@@ -18,12 +18,13 @@ process OVERLAP_HEATMAP {
     script:
     def args = task.ext.args ?: ''
     """
-    MergeCheck_Plot.R
+    MergeCheck_Plot.R --forward ${params.for_primer} --reverse ${params.rev_primer}
     """
 
     stub:
     def args = task.ext.args ?: ''
 
     """
+    touch 
     """
 }