@@ -1,13 +1,20 @@

 #' Plot aggregate regions

 #'

-#' @param x the NanoMethResults object.

+#' @param x the NanoMethResult object.

 #' @param regions a table of regions or GRanges, or a list of such objects. The table of regions must contain chr, start and end columns.

 #' @param groups if 'features' is a list, a vector of characters of the same length as the list containing names for each member.

 #' @param flank the number of flanking bases to add to each side of each region.

+#' @param span the span for loess smoothing.

 #'

 #' @return

 #' @export

-plot_aggregate_regions <- function(x, regions, groups = NULL, flank = 2000, stranded = TRUE) {

+plot_agg_regions <- function(

+    x,

+    regions,

+    groups = NULL,

+    flank = 2000,

+    stranded = TRUE

+) {

     is_df_or_granges <- function(x) {

         is.data.frame(x) || is(x, "GRanges")

     }

@@ -63,7 +70,7 @@ plot_aggregate_regions <- function(x, regions, groups = NULL, flank = 2000, stra

         ggplot2::theme_minimal()

     # take binned means

-    grid_size <- 2^10

+    grid_size <- 2^12

     binned_pos <- seq(-1.1/3, 1 + 1.1/3, length.out = grid_size+ 2)[2:(1 + grid_size)]

     methy_data <- methy_data %>%

         mutate(interval = cut(rel_pos, breaks = grid_size)) %>%

@@ -72,11 +79,9 @@ plot_aggregate_regions <- function(x, regions, groups = NULL, flank = 2000, stra

         ungroup() %>%

         mutate(rel_pos = binned_pos)

-    span <- 0.10

-

     if (!is.null(groups)) {

         p <- p +

-            .geom_smooth(methy_data, sspan = span, group = TRUE) +

+            .agg_geom_smooth(methy_data, span = span, group = TRUE) +

             ggplot2::scale_colour_brewer(palette = "Set1") +

             ggplot2::coord_cartesian(clip = "off")

     } else {

@@ -100,12 +105,114 @@ plot_aggregate_regions <- function(x, regions, groups = NULL, flank = 2000, stra

         ggplot2::ylab("Average Methylation Probability")

 }

+#' Plot aggregate regions with grouped samples

+#'

+#' @param x the NanoMethResult object.

+#' @param regions a table of regions or GRanges.

+#' @param flank the number of flanking bases to add to each side of each region.

+#' @param span the span for loess smoothing.

+#'

+#' @return

+#' @export

+plot_agg_regions_sample_grouped <- function(

+    x,

+    regions,

+    flank = 2000,

+    stranded = TRUE,

+    span = 0.05

+) {

+    is_df_or_granges <- function(x) {

+        is.data.frame(x) || is(x, "GRanges")

+    }

+

+    assert_that(is_df_or_granges(regions) || is.list(regions))

+

+    # convert GRanges to tibble

+    if (is(regions, "GRanges")) {

+        regions <- tibble::as_tibble(regions) %>%

+            dplyr::rename(chr = "seqnames")

+    }

+

+    # convert data.frame regions to single element list

+    if (!is.null(dim(regions))) {

+        regions <- list(regions)

+    }

+

+    # query methylation data

+    methy_data <- purrr::map(

+        regions,

+        function(features) {

+            .get_features_with_rel_pos(

+                    features,

+                    methy = methy(x),

+                    flank = flank,

+                    stranded = stranded,

+                    by_sample = TRUE) %>%

+                dplyr::bind_rows()

+        }

+    )

+

+    methy_data <- methy_data %>%

+        dplyr::bind_rows()

+

+    methy_data <- methy_data %>%

+        dplyr::group_by(sample, .data$rel_pos) %>%

+        dplyr::summarise(methy_prob = mean(.data$methy_prob)) %>%

+        dplyr::ungroup()

+

+    # set up plot

+    p <- ggplot2::ggplot() +

+        ggplot2::ylim(c(0, 1)) +

+        ggplot2::theme_minimal()

+

+    # take binned means

+    grid_size <- 2^12

+    binned_pos_df <- tibble::tibble(

+        interval = levels(cut(methy_data$rel_pos, breaks = grid_size)),

+        binned_pos = seq(-1.1/3, 1 + 1.1/3, length.out = grid_size+ 2)[2:(1 + grid_size)]

+    )

+

+    methy_data <- methy_data %>%

+        mutate(interval = cut(rel_pos, breaks = grid_size)) %>%

+        group_by(sample, interval) %>%

+        summarise(methy_prob = mean(methy_prob)) %>%

+        ungroup() %>%

+        left_join(binned_pos_df, by = "interval") %>%

+        mutate(rel_pos = binned_pos) %>%

+        left_join(samples(x), by = "sample")

+

+    p <- p +

+        ggplot2::stat_smooth(

+            aes(x = .data$rel_pos, y = .data$methy_prob, group = haplotype, col = haplotype),

+            method = "loess",

+            span = span,

+            na.rm = TRUE,

+            se = FALSE,

+            data = methy_data) +

+        ggplot2::scale_colour_brewer(palette = "Set1") +

+        ggplot2::coord_cartesian(clip = "off")

+

+    kb_marker <- round(flank / 1000, 1)

+

+    labels <- c(glue::glue("-{kb_marker}kb"), "start", "end", glue::glue("+{kb_marker}kb"))

+

+    p +

+        geom_vline(xintercept = 0, linetype = "dashed", color = "grey80") +

+        geom_vline(xintercept = 1, linetype = "dashed", color = "grey80") +

+        ggplot2::scale_x_continuous(

+            name = "Relative Position",

+            breaks = c(-.33, 0, 1, 1.33),

+            limits = c(-0.33, 1.33),

+            labels = labels) +

+        ggplot2::ylab("Average Methylation Probability")

+}

+

 .split_rows <- function(x) {

     split(x, 1:nrow(x))

 }

 #' @importFrom purrr map2

-.get_features_with_rel_pos <- function(features, methy, flank, stranded, feature_ids = NULL) {

+.get_features_with_rel_pos <- function(features, methy, flank, stranded, by_sample = FALSE, feature_ids = NULL) {

     .scale_to_feature <- function(x, feature, stranded = stranded) {

         within <- x >= feature$start & x <= feature$end

         upstream <- x < feature$start

@@ -128,6 +235,10 @@ plot_aggregate_regions <- function(x, regions, groups = NULL, flank = 2000, stra

         features$end + flank * 1.10,

         simplify = FALSE)

+    empty <- purrr::map_lgl(out, function(x) nrow(x) == 0)

+    out <- out[!empty]

+    features <- features[!empty, ]

+

     out <- purrr::map2(out, .split_rows(features), function(x, y) {

         if (nrow(x) == 0) {

             return(tibble::add_column(x, rel_pos = numeric()))

@@ -135,10 +246,17 @@ plot_aggregate_regions <- function(x, regions, groups = NULL, flank = 2000, stra

         x$rel_pos <- .scale_to_feature(x$pos, y, stranded = stranded)

-        x %>%

-            group_by(rel_pos) %>%

-            summarise(methy_prob = mean(e1071::sigmoid(statistic))) %>%

-            ungroup()

+        if (by_sample) {

+            x %>%

+                group_by(sample, rel_pos) %>%

+                summarise(methy_prob = mean(e1071::sigmoid(statistic))) %>%

+                ungroup()

+        } else {

+            x %>%

+                group_by(rel_pos) %>%

+                summarise(methy_prob = mean(e1071::sigmoid(statistic))) %>%

+                ungroup()

+        }

     })

     if (!is.null(feature_ids)) {

@@ -149,7 +267,7 @@ plot_aggregate_regions <- function(x, regions, groups = NULL, flank = 2000, stra

     out

 }

-.geom_smooth <- function(data, span, group) {

+.agg_geom_smooth <- function(data, span, group) {

     if (group) {

         ggplot2::stat_smooth(

             aes(x = .data$rel_pos,