Expand plotValsBySeqContext

DESCRIPTION

@@ -55,6 +55,7 @@ Suggests:
     ComplexUpset,
     e1071,
     edgeR,
+    GGally,
     knitr,
     limma,
     MatrixGenerics,

NAMESPACE

@@ -150,6 +150,7 @@ importFrom(cli,cli_warn)
 importFrom(data.table,fread)
 importFrom(dplyr,across)
 importFrom(dplyr,all_of)
+importFrom(dplyr,any_of)
 importFrom(dplyr,arrange)
 importFrom(dplyr,between)
 importFrom(dplyr,bind_cols)
@@ -197,6 +198,7 @@ importFrom(ggplot2,geom_vline)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,labs)
 importFrom(ggplot2,margin)
+importFrom(ggplot2,position_dodge)
 importFrom(ggplot2,scale_color_gradient)
 importFrom(ggplot2,scale_color_manual)
 importFrom(ggplot2,scale_fill_gradientn)

R/plotValsBySeqContext.R

@@ -19,8 +19,10 @@
 #'     \code{facetBy} must be \code{"sample"}), \code{"sample_union"} (
 #'     in which case the selection of sequence contexts will be done for each
 #'     sample, and the union of all selected contexts will be shown in the
-#'     plot) or \code{"overall"} (in which case the selection of sequence
-#'     contexts will be made without considering the sample information). Note
+#'     plot), \code{"overall"} (in which case the selection of sequence
+#'     contexts will be made without considering the sample information) or
+#'     \code{"sample_var"} (in which case the variance across samples will be
+#'     used to select/order contexts). Note
 #'     that for \code{"sample_union"}, the plots may contain more than
 #'     \code{topN + bottomN} contexts. For \code{"sample_union"} and
 #'     \code{"overall"}, the order of the contexts in the plots will be
@@ -29,10 +31,19 @@
 #'     Must be either \code{NULL} (in which case no facetting is done, and
 #'     values are aggregated across all samples in \code{se}) or \code{"sample"}
 #'     (in which case values are aggregated within each sample, and plots are
-#'     facetted accordingly).
+#'     facetted accordingly). Ignored if \code{plotType = "pairs"}.
+#' @param fillBy Character scalar indicating how to fill the bars or violins.
+#'     Must be either \code{NULL} (in which case a single value should be
+#'     specified to \code{"fillColors"} and used for all bars/violins) or
+#'     \code{"sample"} (in which case bars or violins will be split and filled
+#'     by sample). Ignored if \code{plotType = "pairs"}.
+#' @param fillColors Either a (preferably named) character vector defining the
+#'     color to use for each sample if \code{fillBy = "sample"}, or a single
+#'     color to use for all violins/bars.
 #' @param plotType Character scalar indicating what type of plot to create.
 #'     Should be one of \code{"violin"} (note that the violins will be plotted
-#'     with \code{scale="width"}), \code{"bar"} and \code{"errorbar"}.
+#'     with \code{scale="width"}), \code{"bar"}, \code{"errorbar"} or
+#'     \code{"pairs"}.
 #' @param topN,bottomN Numeric scalars determining the number of sequence
 #'     contexts to include in the plot. The \code{topN} contexts with the
 #'     highest average values and the \code{bottomN} contexts with the lowest
@@ -67,21 +78,24 @@
 #'                      topN = 3, bottomN = 3, facetBy = "sample")
 #'
 #' @importFrom dplyr group_by mutate ungroup arrange desc select distinct
-#'     bind_rows slice_max slice_min filter summarize inner_join
+#'     bind_rows slice_max slice_min filter summarize inner_join any_of across
 #' @importFrom SummarizedExperiment rowData assay colnames assayNames
 #' @importFrom ggplot2 ggplot aes geom_violin theme_bw geom_col geom_linerange
-#'     geom_errorbar coord_flip
-#' @importFrom SparseArray rowSums is_nonna
+#'     geom_errorbar coord_flip position_dodge scale_fill_manual
+#' @importFrom SparseArray rowSums is_nonna nnawhich nnavals
 #' @importFrom tidytext reorder_within scale_x_reordered
 #' @importFrom stats sd
 #' @importFrom rlang .data
+#' @importFrom cli cli_warn
 #'
 plotValsBySeqContext <- function(se,
                                  seqContextColumn = "sequenceContext",
                                  assayName = "mod_prob",
                                  aggregation = "mean",
                                  selectContextsBy = "sample_union",
                                  facetBy = NULL,
+                                 fillBy = NULL,
+                                 fillColors = "grey",
                                  plotType = "violin",
                                  topN = 10,
                                  bottomN = 10,
@@ -95,17 +109,30 @@ plotValsBySeqContext <- function(se,
                   validValues = c("none", "mean"))
     .assertScalar(x = facetBy, type = "character", validValues = "sample",
                   allowNULL = TRUE)
-    if (is.null(facetBy)) {
+    .assertScalar(x = plotType, type = "character",
+                  validValues = c("violin", "bar", "errorbar", "pairs"))
+    if (plotType == "pairs") {
+        .assertPackagesAvailable("GGally")
+    }
+    if (is.null(facetBy) || plotType == "pairs") {
         .assertScalar(x = selectContextsBy, type = "character",
-                      validValues = c("sample_union", "overall"))
+                      validValues = c("sample_union", "overall", "sample_var"))
     } else {
         .assertScalar(x = selectContextsBy, type = "character",
-                      validValues = c("sample", "sample_union", "overall"))
+                      validValues = c("sample", "sample_union", "overall",
+                                      "sample_var"))
+    }
+    if (plotType == "pairs" && ncol(se) < 2) {
+        cli_abort("{.arg se} must have at least two samples for a pairs plot")
     }
+    .assertScalar(x = fillBy, type = "character", validValues = "sample",
+                  allowNULL = TRUE)
+    .assertVector(x = fillColors, type = "character", rngLen = c(1, Inf))
     .assertScalar(x = assayName, type = "character",
                   validValues = .getReadLevelAssayNames(se))
-    .assertScalar(x = plotType, type = "character",
-                  validValues = c("violin", "bar", "errorbar"))
+    if (selectContextsBy %in% c("sample_var") && ncol(se) < 2) {
+        cli_abort("{.arg se} must have at least two samples if {.arg selectContextsBy} is {.val {selectContextsBy}}")
+    }
     .assertScalar(x = topN, type = "numeric", rngIncl = c(0, Inf))
     .assertScalar(x = bottomN, type = "numeric", rngIncl = c(0, Inf))
     .assertScalar(x = flipCoord, type = "logical")
@@ -133,9 +160,10 @@ plotValsBySeqContext <- function(se,
     }
 
     # ... sample-wise ones (required if either facetBy = "sample" or
-    #     selectContextsBy = "sample" or "sample_union")
+    #     selectContextsBy = "sample", "sample_union" or "sample_var")
     if ((!is.null(facetBy) && facetBy == "sample") ||
-        selectContextsBy %in% c("sample", "sample_union")) {
+        (!is.null(fillBy) && fillBy == "sample") ||
+        selectContextsBy %in% c("sample", "sample_union", "sample_var")) {
         if (aggregation == "mean") {
             dfSample <- do.call(bind_rows, lapply(colnames(se), function(cn) {
                 a <- assay(se, assayName)[[cn]]
@@ -157,30 +185,39 @@ plotValsBySeqContext <- function(se,
     }
 
     # choose data frames to use for selection/plotting
-    if (selectContextsBy %in% c("sample", "sample_union")) {
+    if (selectContextsBy %in% c("sample", "sample_union", "sample_var")) {
         dfSel <- dfSample
     } else {
         dfSel <- dfGlobal
     }
 
-    if (is.null(facetBy)) {
+    if (is.null(facetBy) && is.null(fillBy)) {
         dfPlot <- dfGlobal
     } else {
         dfPlot <- dfSample
     }
 
     dfPlot$flipCoord <- flipCoord
+    dfPlot$selectContextsBy <- selectContextsBy
 
     # calculate mean/sd for each context and select top/bottom ones to include
     dfSelSum <- dfSel |>
         group_by(.data$seqContext, .data$sample) |>
         summarize(valsMean = mean(.data$vals),
                   .groups = "drop")
+    if (selectContextsBy == "sample_var") {
+        dfSelSum <- dfSelSum |>
+            group_by(.data$seqContext) |>
+            mutate(valsMeanVar = var(.data$valsMean))
+        selCol <- "valsMeanVar"
+    } else {
+        selCol <- "valsMean"
+    }
     dfSelSum <- bind_rows(
         dfSelSum |> group_by(sample) |>
-            slice_max(.data$valsMean, n = topN, with_ties = FALSE),
+            slice_max(.data[[selCol]], n = topN, with_ties = FALSE),
         dfSelSum |> group_by(sample) |>
-            slice_min(.data$valsMean, n = bottomN, with_ties = FALSE)
+            slice_min(.data[[selCol]], n = bottomN, with_ties = FALSE)
     ) |>
         ungroup() |>
         distinct()
@@ -197,62 +234,137 @@ plotValsBySeqContext <- function(se,
     } else {
         contextsToKeep <- unique(dfSelSum$seqContext)
         dfPlot <- dfPlot |>
-            dplyr::filter(seqContext %in% contextsToKeep) |>
+            dplyr::filter(.data$seqContext %in% contextsToKeep) |>
             mutate(orderCol = "overall")
+        if (selectContextsBy == "sample_var") {
+            dfPlot <- dfPlot |>
+                left_join(dfSelSum |>
+                              dplyr::select(c("seqContext", "valsMeanVar")) |>
+                              distinct(),
+                          by = "seqContext")
+        }
     }
 
     # plot
-    if (plotType == "violin") {
-        dfPlot <- dfPlot |>
-            mutate(seqContext = reorder_within(.data$seqContext, by = ifelse(
-                .data$flipCoord, .data$vals, -.data$vals),
-                within = orderCol, fun = mean))
-        gg <- ggplot(dfPlot, aes(x = .data$seqContext, y = .data$vals))  +
-            geom_violin(scale = "width")
-    } else if (plotType %in% c("bar", "errorbar")) {
+    if (plotType == "pairs") {
         dfPlot <- dfPlot |>
             group_by(.data$seqContext, .data$sample, .data$orderCol,
                      .data$flipCoord) |>
             summarize(valsMean = mean(.data$vals),
                       valsSd = sd(.data$vals),
                       .groups = "drop") |>
-            mutate(seqContext = reorder_within(.data$seqContext, by = ifelse(
-                .data$flipCoord, .data$valsMean, -.data$valsMean),
-                within = orderCol, fun = mean))
-        gg <- ggplot(dfPlot,
-                     aes(x = .data$seqContext, y = .data$valsMean)) +
-            geom_col()
-        if (plotType == "errorbar") {
+            dplyr::select(c("seqContext", "sample", "valsMean")) |>
+            tidyr::pivot_wider(names_from = "sample", values_from = "valsMean")
+        GGally::ggpairs(dfPlot, columns = setdiff(colnames(dfPlot), "seqContext")) +
+            theme_bw()
+    } else {
+        if (plotType == "violin") {
+            dfPlot <- dfPlot |>
+                mutate(seqContext = reorder_within(
+                    .data$seqContext,
+                    by = ifelse(flipCoord, 1, -1) *
+                        ifelse(.data$selectContextsBy == "sample_var",
+                               .data$valsMeanVar, .data$vals),
+                    within = .data$orderCol,
+                    fun = mean
+                ))
+            gg <- ggplot(dfPlot, aes(x = .data$seqContext, y = .data$vals))
+            if (is.null(fillBy)) {
+                gg <- gg +
+                    geom_violin(scale = "width", fill = fillColors[1],
+                                quantiles = 0.5, quantile.linetype = "solid")
+            } else {
+                gg <- gg +
+                    geom_violin(scale = "width", aes(fill = .data[[fillBy]]),
+                                quantiles = 0.5, quantile.linetype = "solid")
+                if (length(fillColors) >= length(unique(dfPlot[[fillBy]]))) {
+                    gg <- gg +
+                        scale_fill_manual(values = fillColors)
+                } else {
+                    cli_warn("Not enough colors - using defaults")
+                }
+            }
+        } else if (plotType %in% c("bar", "errorbar")) {
+            dfPlot <- dfPlot |>
+                group_by(across(any_of(c("seqContext", "sample", "orderCol",
+                                         "flipCoord", "valsMeanVar",
+                                         "selectContextsBy")))) |>
+                summarize(valsMean = mean(.data$vals),
+                          valsSd = sd(.data$vals),
+                          .groups = "drop") |>
+                mutate(seqContext = reorder_within(
+                    .data$seqContext,
+                    by = ifelse(flipCoord, 1, -1) *
+                        ifelse(.data$selectContextsBy == "sample_var",
+                               .data$valsMeanVar, .data$valsMean),
+                    within = .data$orderCol,
+                    fun = mean
+                ))
+            gg <- ggplot(dfPlot,
+                         aes(x = .data$seqContext, y = .data$valsMean))
+            if (is.null(fillBy)) {
+                gg <- gg +
+                    geom_col(fill = fillColors[1])
+            } else {
+                gg <- gg +
+                    geom_col(aes(fill = .data[[fillBy]]),
+                             position = position_dodge())
+                if (length(fillColors) >= length(unique(dfPlot[[fillBy]]))) {
+                    gg <- gg +
+                        scale_fill_manual(values = fillColors)
+                } else {
+                    cli_warn("Not enough colors - using defaults")
+                }
+            }
+            if (plotType == "errorbar") {
+                if (is.null(fillBy)) {
+                    gg <- gg +
+                        geom_linerange(
+                            aes(ymin = .data$valsMean,
+                                ymax = .data$valsMean + .data$valsSd)
+                        ) +
+                        geom_errorbar(
+                            aes(ymin = .data$valsMean + .data$valsSd,
+                                ymax = .data$valsMean + .data$valsSd),
+                            width = 0.2
+                        )
+                } else {
+                    gg <- gg +
+                        geom_linerange(
+                            aes(group = .data[[fillBy]],
+                                ymin = .data$valsMean,
+                                ymax = .data$valsMean + .data$valsSd),
+                            position = position_dodge(width = 0.9)
+                        ) +
+                        geom_errorbar(
+                            aes(group = .data[[fillBy]],
+                                ymin = .data$valsMean + .data$valsSd,
+                                ymax = .data$valsMean + .data$valsSd),
+                            width = 0.2,
+                            position = position_dodge(width = 0.9)
+                        )
+                }
+            }
+        }
+        gg <- gg +
+            theme_bw() +
+            labs(x = "Sequence context",
+                 y = paste0(yAxisLabel, ifelse(plotType == "violin", "",
+                                               ifelse(plotType == "bar", " (mean)",
+                                                      " (mean + sd)"))))
+        if (!is.null(facetBy)) {
+            gg <- gg + facet_wrap(~ sample,
+                                  scales = ifelse(flipCoord, "free_y", "free_x"))
+        }
+
+        if (flipCoord) {
+            gg <- gg +
+                coord_flip()
+        } else {
             gg <- gg +
-                geom_linerange(
-                    aes(ymin = .data$valsMean,
-                        ymax = .data$valsMean + .data$valsSd)
-                ) +
-                geom_errorbar(
-                    aes(ymin = .data$valsMean + .data$valsSd,
-                        ymax = .data$valsMean + .data$valsSd),
-                    width = 0.2
-                )
+                theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
         }
-    }
-    gg <- gg +
-        theme_bw() +
-        labs(x = "Sequence context",
-             y = paste0(yAxisLabel, ifelse(plotType == "violin", "",
-                                           ifelse(plotType == "bar", " (mean)",
-                                                  " (mean + sd)"))))
-    if (!is.null(facetBy)) {
-        gg <- gg + facet_wrap(~ sample,
-                              scales = ifelse(flipCoord, "free_y", "free_x"))
-    }
 
-    if (flipCoord) {
-        gg <- gg +
-            coord_flip()
-    } else {
-        gg <- gg +
-            theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
+        gg + scale_x_reordered()
     }
-
-    gg + scale_x_reordered()
 }