|
| 1 | +#### libraries #### |
| 2 | +library(data.table) |
| 3 | +library(tidyverse) |
| 4 | +library(patchwork) |
| 5 | +library(ggplot2) |
| 6 | +library(ggnewscale) |
| 7 | +library(stringr) |
| 8 | +# for interactive plotting |
| 9 | +library(plotly) |
| 10 | +library(htmlwidgets) |
| 11 | +# set plot base size |
| 12 | +theme_set(theme_minimal(base_size = 6)) |
| 13 | + |
| 14 | +#### configs #### |
| 15 | +# input |
| 16 | +sample_annotation_path <- snakemake@input[["sample_annotation"]] |
| 17 | +sample_annotation_w_QC_path <- snakemake@input[["sample_annotation_w_QC"]] |
| 18 | + |
| 19 | +# output |
| 20 | +sample_annotation_plot_path <- snakemake@output[["sample_annotation_plot"]] |
| 21 | +sample_annotation_html_path <- snakemake@output[["sample_annotation_html"]] |
| 22 | + |
| 23 | +#### load & prepare data #### |
| 24 | +# load data |
| 25 | +sample_annotation <- data.table::fread(file.path(sample_annotation_path), header = TRUE) |
| 26 | +sample_annotation <- data.frame(sample_annotation[!duplicated(sample_annotation[[1]]), ], row.names = 1, check.names = FALSE) |
| 27 | + |
| 28 | +anno <- data.frame(fread(file.path(sample_annotation_w_QC_path), header=TRUE), row.names=1) |
| 29 | + |
| 30 | +# determine QC (pipeline provided) columns |
| 31 | +names(sample_annotation) <- gsub(" +", "_", names(sample_annotation)) # replace empty space ` ` with underscore `_` |
| 32 | +qc_cols <- setdiff(names(anno), names(sample_annotation)) |
| 33 | + |
| 34 | +# determine metadata (user provided) columns by removing non-numeric columns that are unique for each row (e.g., bam_file) |
| 35 | +sample_annotation <- sample_annotation %>% select(where(\(.x){ n <- dplyr::n_distinct(na.omit(.x)); (is.numeric(.x) || n < nrow(sample_annotation)) })) |
| 36 | +meta_cols <- names(sample_annotation) |
| 37 | + |
| 38 | +# drop columns with no variation (e.g., read_type) |
| 39 | +has_var <- function(x) length(unique(na.omit(x))) > 1 |
| 40 | +qc_cols <- keep(qc_cols, ~ has_var(anno[[.x]])) |
| 41 | +meta_cols <- keep(meta_cols, ~ has_var(anno[[.x]])) |
| 42 | + |
| 43 | +# collapse "duplicates" (ie redundant columns) only among non-numeric |
| 44 | +num_cols <- meta_cols[vapply(anno[meta_cols], is.numeric, logical(1))] |
| 45 | +cat_cols <- setdiff(meta_cols, num_cols) |
| 46 | +sig <- vapply(cat_cols, |
| 47 | + \(col) paste(match(anno[[col]], unique(anno[[col]])), collapse = "|"), |
| 48 | + character(1)) |
| 49 | +meta_cols <- c(cat_cols[!duplicated(sig)], num_cols) |
| 50 | + |
| 51 | +#### Z-score & cluster QC data #### |
| 52 | +qc_mat <- anno |> select(all_of(qc_cols)) |> scale() |> as.matrix() |
| 53 | + |
| 54 | +row_ord <- hclust(dist(qc_mat))$order |
| 55 | +col_ord <- hclust(dist(t(qc_mat)))$order |
| 56 | + |
| 57 | +qc_long <- as_tibble(qc_mat[row_ord, col_ord], rownames = "sample") |> |
| 58 | + pivot_longer(-sample, names_to = "metric", values_to = "z") |
| 59 | + |
| 60 | +#### prepare metadata for plotting #### |
| 61 | +meta_long <- anno[row_ord, meta_cols] %>% |
| 62 | + mutate(sample = rownames(anno)[row_ord]) %>% |
| 63 | + mutate(across(-sample, as.character)) %>% |
| 64 | + pivot_longer(-sample, names_to = "meta", values_to = "value") %>% |
| 65 | + group_by(meta) %>% |
| 66 | + mutate(num_val = suppressWarnings(as.numeric(value)), |
| 67 | + type = if (all(!is.na(num_val))) "numeric" else "factor", |
| 68 | + col = if (type[1] == "numeric") { |
| 69 | + scales::col_numeric("plasma", |
| 70 | + domain = range(num_val, na.rm = TRUE))(num_val) |
| 71 | + } else { |
| 72 | + pal <- scales::hue_pal(l = 65)(n_distinct(value)) |
| 73 | + setNames(pal, sort(unique(value)))[value] |
| 74 | + } |
| 75 | + ) %>% |
| 76 | + ungroup() %>% |
| 77 | + select(-num_val) |
| 78 | + |
| 79 | +#### embed ALL metadata into a tooltip string of interactive plot #### |
| 80 | +meta_txt <- anno[row_ord, meta_cols] %>% mutate(across(everything(), as.character)) |
| 81 | +meta_txt <- pmap_chr(meta_txt, \(...) { |
| 82 | + vals <- c(...) |
| 83 | + paste(paste(names(vals), vals, sep = ": "), collapse = "<br>") |
| 84 | + }) |
| 85 | +names(meta_txt) <- rownames(anno)[row_ord] |
| 86 | + |
| 87 | +#### plot heatmaps #### |
| 88 | + |
| 89 | +#### QC heatmap #### |
| 90 | + |
| 91 | +# add (un-scaled) metric values (raw) for tooltip of interactive plot |
| 92 | +qc_raw_long <- anno[row_ord, qc_cols] %>% |
| 93 | + mutate(sample = rownames(anno)[row_ord]) %>% |
| 94 | + pivot_longer(-sample, names_to = "metric", values_to = "raw") |
| 95 | + |
| 96 | +qc_long <- qc_long %>% left_join(qc_raw_long, by = c("sample", "metric")) |
| 97 | + |
| 98 | +# keep the clustered order in the plot and add hover-tooltip |
| 99 | +qc_long <- qc_long %>% |
| 100 | + mutate(sample = factor(sample, levels = rownames(qc_mat)[row_ord]), # row-order |
| 101 | + metric = factor(metric, levels = colnames(qc_mat)[col_ord]), # col-order |
| 102 | + hover = paste0("Sample: ", sample, |
| 103 | + "<br>Metric: ", metric, |
| 104 | + "<br>Value: ", signif(raw, 4), |
| 105 | + "<br>", meta_txt[as.character(sample)]) |
| 106 | + ) |
| 107 | + |
| 108 | +# plot |
| 109 | +p_qc <- ggplot(qc_long, aes(x = metric, |
| 110 | + y = sample, |
| 111 | + fill = z, |
| 112 | + text = hover)) + |
| 113 | + geom_tile() + |
| 114 | + scale_x_discrete(limits = colnames(qc_mat)[col_ord]) + # enforce col order |
| 115 | + scale_y_discrete(limits = rownames(qc_mat)[row_ord]) + # enforce row order |
| 116 | + scale_fill_gradient2(low = "blue", mid = "white", high = "red", name = "z-score", |
| 117 | + guide = guide_colourbar(barheight = 2, # thinner |
| 118 | + barwidth = 0.15)) + |
| 119 | + labs(x = NULL, y = NULL, title = "QC metrics (scaled)") + |
| 120 | + theme(axis.text.x = element_text(angle = 45, hjust = 1), |
| 121 | + panel.grid = element_blank()) |
| 122 | + |
| 123 | +#### metadata heatmap as "annotation" #### |
| 124 | +p_meta <- NULL |
| 125 | +if(length(meta_cols) > 0){ |
| 126 | + # order columns (x) exactly like the QC heatmap |
| 127 | + meta_long <- meta_long %>% mutate(sample = factor(sample, levels = rownames(qc_mat)[row_ord])) # row-order |
| 128 | + meta_levels <- unique(meta_long$meta) |
| 129 | + meta_long <- meta_long %>% mutate(meta = factor(meta, levels = meta_levels)) |
| 130 | + |
| 131 | + p_meta <- ggplot() + |
| 132 | + scale_y_discrete(limits = levels(qc_long$sample)) + |
| 133 | + scale_x_discrete(limits = meta_levels) + |
| 134 | + labs(x = NULL, y = NULL, title = "Metadata") + |
| 135 | + theme(axis.text.x = element_text(angle = 45, hjust = 1), |
| 136 | + axis.text.y = element_blank(), |
| 137 | + axis.title = element_blank(), |
| 138 | + panel.grid = element_blank()) |
| 139 | + |
| 140 | + for (v in meta_levels) { |
| 141 | + dat <- dplyr::filter(meta_long, meta == v) |
| 142 | + |
| 143 | + p_meta <- p_meta + ggnewscale::new_scale_fill() # reset “fill” for this column |
| 144 | + |
| 145 | + if (dat$type[1] == "numeric") { # continuous legend |
| 146 | + p_meta <- p_meta + |
| 147 | + geom_tile(data = dat, aes(x = meta, y = sample, fill = as.numeric(value)), colour = "grey60", linewidth = 0.1) + |
| 148 | + scale_fill_viridis_c(name = v, option = "plasma", |
| 149 | + guide = guide_colourbar(barheight = 2, # thinner |
| 150 | + barwidth = 0.15)) |
| 151 | + } else { # categorical legend |
| 152 | + pal <- setNames(dat$col, dat$value) |
| 153 | + |
| 154 | + # reduce legend in case of more than 10 levels |
| 155 | + max_items <- min(10, length(unique(dat$value))) |
| 156 | + all_levels <- unique(names(pal)) |
| 157 | + show_levels <- all_levels[1:max_items] |
| 158 | + |
| 159 | + p_meta <- p_meta + |
| 160 | + geom_tile(data = dat, aes(x = meta, y = sample, fill = value), colour = "grey60", linewidth = 0.1) + |
| 161 | + scale_fill_manual(values = pal, |
| 162 | + # name = v, |
| 163 | + breaks = show_levels, |
| 164 | + guide = guide_legend(keywidth = 0.25, |
| 165 | + keyheight = 0.4, |
| 166 | + ncol=1, |
| 167 | + byrow = TRUE, |
| 168 | + title = ifelse( |
| 169 | + length(all_levels) <= max_items, |
| 170 | + v, |
| 171 | + paste0(v, " (showing ", max_items, "/", length(all_levels), ")") |
| 172 | + ) |
| 173 | + ) |
| 174 | + ) |
| 175 | + } |
| 176 | + } |
| 177 | +} |
| 178 | + |
| 179 | +#### combine and save plots #### |
| 180 | +p_combined <- if (is.null(p_meta)) p_qc else (p_qc | p_meta) + plot_layout(widths = c(length(qc_cols), length(meta_cols)), guides = "collect") & theme(legend.position = "right") |
| 181 | + |
| 182 | +# determine sizes |
| 183 | +n_rows <- nrow(qc_mat) |
| 184 | +n_cols <- length(qc_cols) + length(meta_cols) |
| 185 | +max_row_label <- max(nchar(rownames(anno))) |
| 186 | +max_col_label <- max(nchar(c(qc_cols, meta_cols))) |
| 187 | + |
| 188 | +height_in <- n_rows * 0.08 + max_col_label * 0.05 + 1 |
| 189 | +width_in <- n_cols * 0.10 + max_row_label * 0.05 + 2 |
| 190 | + |
| 191 | +# options(repr.plot.width = width_in, repr.plot.height = height_in) |
| 192 | +# p_combined |
| 193 | + |
| 194 | +ggsave(sample_annotation_plot_path, plot = p_combined, width = width_in, height = height_in, units = "in", dpi = 300) |
| 195 | + |
| 196 | +#### interactive plot #### |
| 197 | +# determine sizes in pixels |
| 198 | +width_px <- round((length(qc_cols) * 0.10 + max_row_label * 0.05 + 2) * 96) |
| 199 | +height_px <- round(height_in * 96) |
| 200 | + |
| 201 | +p_qc_interactive <- plotly::ggplotly(p_qc, tooltip = "text", width = width_px, height = height_px) |
| 202 | + |
| 203 | +# p_qc_interactive |
| 204 | + |
| 205 | +htmlwidgets::saveWidget(p_qc_interactive, sample_annotation_html_path, selfcontained = TRUE, title = "Sample annotation") |
0 commit comments