Add group labels to SV plots

Snakefile

@@ -278,15 +278,16 @@ ruleorder: plot_SV_calls_simulated > plot_SV_calls
 rule plot_SV_calls:
     input:
         counts = "counts/{sample}/{windows}.txt.gz",
-        calls  = "sv_calls/{sample}/{windows}.{bpdens}/{method}.txt",
-        complex = "sv_calls/{sample}/{windows}.{bpdens}/{method}.complex.tsv",
+        calls  = "sv_calls/{sample}/{windows}.{bpdens}/{method}_filter{filter}.txt",
+        complex = "sv_calls/{sample}/{windows}.{bpdens}/{method}_filter{filter}.complex.tsv",
         strand = "strand_states/{sample}/{windows}.{bpdens}/final.txt",
         segments = "segmentation2/{sample}/{windows}.{bpdens}.txt",
         scsegments = "segmentation-singlecell/{sample}/{windows}.{bpdens}.txt",
+        grouptrack = "postprocessing/group-table/{sample}/{windows}.{bpdens}/{method}.tsv",
     output:
-        "sv_calls/{sample}/{windows}.{bpdens,selected_j[0-9\\.]+_s[0-9\\.]+}/plots/sv_calls/{method}.{chrom}.pdf"
+        "sv_calls/{sample}/{windows}.{bpdens,selected_j[0-9\\.]+_s[0-9\\.]+}/plots/sv_calls/{method}_filter{filter,(TRUE|FALSE)}.{chrom}.pdf"
     log:
-        "log/plot_SV_calls/{sample}/{windows}.{bpdens}.{method}.{chrom}.log"
+        "log/plot_SV_calls/{sample}/{windows}.{bpdens}.{method}_filter{filter}.{chrom}.log"
     shell:
         """
         Rscript utils/plot-sv-calls.R \
@@ -294,6 +295,7 @@ rule plot_SV_calls:
             singlecellsegments={input.scsegments} \
             strand={input.strand} \
             complex={input.complex} \
+            groups={input.grouptrack} \
             calls={input.calls} \
             {input.counts} \
             {wildcards.chrom} \
@@ -669,13 +671,22 @@ rule postprocessing_filter:
 
 rule postprocessing_merge:
     input: 
-        calls = "postprocessing/filter/{sample}/{window}_fixed_norm.{bpdens,selected_j[0-9\\.]+_s[0-9\\.]+}/simpleCalls_llr{llr}_poppriors{pop_priors,(TRUE|FALSE)}_haplotags{use_haplotags,(TRUE|FALSE)}_gtcutoff{gtcutoff,[0-9\\.]+}_regfactor{regfactor,[0-9]+}.txt"
+        calls = "postprocessing/filter/{sample}/{window}_fixed_norm.{bpdens}/simpleCalls_llr{llr}_poppriors{pop_priors}_haplotags{use_haplotags}_gtcutoff{gtcutoff}_regfactor{regfactor}.txt"
     output: 
         calls = "postprocessing/merge/{sample}/{window}_fixed_norm.{bpdens,selected_j[0-9\\.]+_s[0-9\\.]+}/simpleCalls_llr{llr}_poppriors{pop_priors,(TRUE|FALSE)}_haplotags{use_haplotags,(TRUE|FALSE)}_gtcutoff{gtcutoff,[0-9\\.]+}_regfactor{regfactor,[0-9]+}.txt"
     shell:
         'utils/group_nearby_calls_of_same_AF_and_generate_output_table.pl {input.calls}  > {output.calls}'
 
 
+rule postprocessing_sv_group_table:
+    input: 
+        calls = "postprocessing/merge/{sample}/{window}_fixed_norm.{bpdens}/simpleCalls_llr{llr}_poppriors{pop_priors}_haplotags{use_haplotags}_gtcutoff{gtcutoff}_regfactor{regfactor}.txt"
+    output: 
+        grouptrack = "postprocessing/group-table/{sample}/{window}_fixed_norm.{bpdens,selected_j[0-9\\.]+_s[0-9\\.]+}/simpleCalls_llr{llr}_poppriors{pop_priors,(TRUE|FALSE)}_haplotags{use_haplotags,(TRUE|FALSE)}_gtcutoff{gtcutoff,[0-9\\.]+}_regfactor{regfactor,[0-9]+}.tsv"
+    shell:
+        'utils/create-sv-group-track.py {input.calls}  > {output.grouptrack}'
+
+
 ################################################################################
 # Strand states & phasing                                                      #
 ################################################################################

utils/plot-sv-calls.R

@@ -9,7 +9,7 @@ suppressMessages(library(ggplot2))
 library(scales) %>% invisible
 library(assertthat) %>% invisible
 library(stringr) %>% invisible
-
+library(RColorBrewer) %>% invisible
 
 ################################################################################
 # Settings                                                                     #
@@ -76,6 +76,7 @@ print_usage_and_stop = function(msg = NULL) {
   message("    truth=<file>              Mark the `true`` SVs provided in a table          ")
   message("    strand=<file>             Mark the strand states which calls are based on   ")
   message("    complex=<file>            Mark complex regions given in file                ")
+  message("    groups=<file>             Table with SV call grouping                       ")
   message("    no-none                   Do not hightlight black-listed (i.e. None) bins   ")
   message("                                                                                ")
   message("Generates one plot per chromosome listing all cells below another, separated    ")
@@ -112,7 +113,7 @@ cells_per_page = 8
 show_none  = T
 
 if (length(args)>3) {
-  if (!all(grepl("^(strand|calls|segments|per-page|truth|no-none|complex|singlecellsegments)=?", args[1:(length(args)-3)]))) {
+  if (!all(grepl("^(strand|calls|segments|per-page|truth|no-none|complex|singlecellsegments|groups)=?", args[1:(length(args)-3)]))) {
     print_usage_and_stop("[Error]: Options must be one of `calls`, `segments`, `per-page`, or `truth`") }
   for (op in args[1:(length(args)-3)]) {
     if (grepl("^segments=", op)) f_segments = str_sub(op, 10)
@@ -124,6 +125,7 @@ if (length(args)>3) {
     }
     if (grepl("^strand=", op))   f_strand = str_sub(op, 8)
     if (grepl("^complex=", op))   f_complex = str_sub(op, 9)
+    if (grepl("^groups=", op))   f_groups = str_sub(op, 8)
     if (grepl("^singlecellsegments=", op))   f_scsegments = str_sub(op, 20)
     if (grepl("^no-none$", op)) show_none = F
   }
@@ -249,6 +251,19 @@ if (!is.null(f_complex)) {
 
 }
 
+### Check SV groups file
+if (!is.null(f_groups)) {
+  message(" * Reading SV group file from ", f_groups, "...")
+  groups = fread(f_groups)
+  assert_that("chrom"   %in% colnames(groups),
+              "start"   %in% colnames(groups),
+              "end"     %in% colnames(groups),
+              "group_id"     %in% colnames(groups)) %>% invisible
+  groups[, group_id := paste("SV group", group_id)]
+  groups = groups[chrom == CHROM]
+  message("   --> Found ", nrow(groups), " SV groups in chromosome ", CHROM)
+}
+
 ### Check single cell segmentation file
 if (!is.null(f_scsegments)) {
   message(" * Reading per-cell segmentation regions state file from ", f_scsegments, "...")
@@ -345,6 +360,18 @@ while (i <= n_cells) {
       }
     }
 
+    # Add bars for SV group, if available
+    if (!is.null(f_groups)) {
+      message("   * Adding SV groups")
+      if (nrow(groups) > 0) {
+        plt <- plt +
+          geom_rect(data = groups,
+                    aes(xmin = start, xmax = end, ymin = .85*y_lim, ymax = Inf, fill = group_id))
+          # Add colors for SV classes
+          manual_colors = c(manual_colors, setNames(colorRampPalette(brewer.pal(12,"Set2"))(nrow(groups)),groups$group_id))
+      }
+    }
+
     # Add bars for complex states, if available
     if (!is.null(f_complex)) {
       message("   * Adding complex intervals")