R/plot_region.R

In Shians/NanoMethViz: Visualise methylation data from Oxford Nanopore sequencing

plot_region_impl <- function(
        x,
        chr,
        start,
        end,
        anno_regions = NULL,
        binary_threshold = NULL,
        avg_method = c("mean", "median"),
        spaghetti = FALSE,
        heatmap = TRUE,
        heatmap_subsample = 50,
        smoothing_window = 2000,
        gene_anno = TRUE,
        window_prop = 0,
        palette = ggplot2::scale_colour_brewer(palette = "Set1"),
        line_size = 1,
        mod_scale = c(0, 1),
        span = NULL
) {
    sample_anno <- samples(x)
    chr <- as.character(chr)
    start <- as.numeric(start)
    end <- as.numeric(end)

    if (length(window_prop) == 1) {
        window_prop <- c(window_prop, window_prop)
    }

    feature_width <- end - start
    window_left <- feature_width * window_prop[1]
    window_right <- feature_width * window_prop[2]

    # query data
    methy_data <- query_methy(
        x,
        chr,
        floor(start - window_left * 1.1),
        ceiling(end + window_right * 1.1),
        simplify = TRUE
    )

    if (nrow(methy_data) == 0) {
        warning("no methylation data in region, returning empty plot")
        return(ggplot() + theme_void())
    }

    methy_data <- methy_data %>%
        dplyr::select(-"strand") %>%
        tibble::as_tibble()


    # setup base plot
    title <- glue::glue("{chr}:{start}-{end}")
    xlim <- round(c(start - window_left, end + window_right))
    p1 <- plot_methylation_data(
            methy_data = methy_data,
            start = start,
            end = end,
            chr = chr,
            title = title,
            anno_regions = anno_regions,
            binary_threshold = binary_threshold,
            avg_method = avg_method,
            spaghetti = spaghetti,
            sample_anno = sample_anno,
            smoothing_window = smoothing_window,
            palette_col = palette,
            line_size = line_size,
            mod_scale = mod_scale
        ) +
        ggplot2::coord_cartesian(xlim = xlim, expand = FALSE) +
        ggplot2::labs(x = "Position", y = "Mean Modification Probability")

    p_out <- p1

    # if exon anno exists, append it to plot
    if (gene_anno && nrow(exons(x)) != 0) {
        exons_anno <- query_exons_region(x, chr = chr, start = start, end = end)

        p2 <- plot_gene_annotation(exons_anno, xlim[1], xlim[2]) +
            ggplot2::coord_cartesian(xlim = xlim, expand = FALSE) +
            ggplot2::scale_x_continuous(labels = scales::label_number(scale_cut = scales::cut_si("b")))

        anno_height <- attr(p2, "plot_height")
        heights <- c(1, 0.075 * anno_height)

        p_out <- p1 / p2 + patchwork::plot_layout(heights = heights)
    }

    # if heatmap requested, append it to plot
    if (heatmap) {
        p_heatmap <- plot_region_heatmap(x, chr, start, end, window_prop = window_prop, subsample = heatmap_subsample) +
            ggplot2::coord_cartesian(
                xlim = xlim
            )

        p_out <- stack_plots(p_out, ggrastr::rasterise(p_heatmap, dpi = 300))
    }

    p_out
}

#' @rdname plot_region
#'
#' @param anno_regions the data.frame of regions to be annotated.
#' @param binary_threshold the modification probability such that calls with
#'   modification probability above the threshold are set to 1 and probabilities
#'   equal to or below the threshold are set to 0.
#' @param avg_method the average method for pre-smoothing at each genomic position.
#'   Data is pre-smoothed at each genomic position before the smoothed aggregate line
#'   is generated for performance reasons. The default is "mean" which corresponds
#'   to the average methylation fraction. The alternative "median" option is
#'   closer to an average within the more common methylation state.
#' @param spaghetti whether or not individual reads should be shown.
#' @param heatmap whether or not read-methylation heatmap should be shown.
#' @param heatmap_subsample how many packed rows of reads to subsample to.
#' @param smoothing_window the window size for smoothing the trend line.
#' @param gene_anno whether to show gene annotation.
#' @param window_prop the size of flanking region to plot. Can be a vector of two
#'   values for left and right window size. Values indicate proportion of gene
#'   length.
#' @param palette the ggplot colour palette used for groups.
#' @param line_size the size of the lines.
#' @param mod_scale the scale range for modification probabilities. Default c(0, 1), set to "auto" for automatic
#'   limits.
#' @param span DEPRECATED, use smoothing_window instead. Will be removed in next version.
#'
#' @details
#' This function plots the methylation data for a given region. The main trendline plot shows the average methylation
#' probability across the region. The heatmap plot shows the methylation probability for each read across the region.
#' The gene annotation plot shows the exons of the region. In the heatmap, each row represents one or more
#' non-overlapping reads where the coloured segments represent the methylation probability at each position. Data along
#' a read is connected by a grey line. The gene annotation plot shows the isoforms and exons of genes within the region,
#' with arrows indicating the direction of transcription.
#'
#' Since V3.0.0 NanoMethViz has changed the smoothing strategy from a loess smoothing to a weighted moving average. This
#' is because the loess smoothing was too computationally expensive for large datasets and had a span parameter that was
#' difficult to tune. The new smoothing strategy is controlled by the smoothing_window argument.
#'
#' @examples
#' nmr <- load_example_nanomethresult()
#' plot_region(nmr, "chr7", 6703892, 6730431)
#'
#' @export
setMethod("plot_region",
    signature(x = "NanoMethResult", chr = "character", start = "numeric", end = "numeric"),
    plot_region_impl
)

#' @rdname plot_region
#' @export
setMethod("plot_region",
    signature(x = "ModBamResult", chr = "character", start = "numeric", end = "numeric"),
    plot_region_impl
)

#' @rdname plot_region
#' @export
setMethod("plot_region",
    signature(x = "NanoMethResult", chr = "factor", start = "numeric", end = "numeric"),
    plot_region_impl
)

#' @rdname plot_region
#' @export
setMethod("plot_region",
    signature(x = "ModBamResult", chr = "factor", start = "numeric", end = "numeric"),
    plot_region_impl
)