Plot differentially methylated regions

plot_dmr.R

# # MIT License
#
# Copyright (c) 2024 Mickaël Canouil
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# library(ggplot2)
# library(ggbio)
# library(Homo.sapiens)
# library(GenomicRanges)
# library(IRanges)
# library(scales)
# library(purrr)
# library(patchwork)
draw_dmr <- function(.dmr, .dmp, base_colour = viridis_pal(begin = 0.5, end = 0.5)(1)) {
  pval_trans_inf <- function() {
    neglog10_breaks <- function(n = 5) {
      function(x) {
        rng <- -log(range(x, na.rm = TRUE), base = 10)
        min <- 0
        max <- floor(rng[1])
        if (max == min) {
          10^-min
        } else {
          by <- floor((max - min) / n) + 1
          10^-seq(min, max, by = by)
        }
      }
    }
    trans_new(
      name = "pval",
      transform = function(x) {
        x[x == Inf] <- 0
        x[x == -Inf] <- 1
        -log(x, 10)
      },
      inverse = function(x) 10^-x,
      breaks = neglog10_breaks(),
      format = function(x) {
        parse(
          text = gsub("e", " %*% 10^", gsub("1e+00", "1", scientific_format()(x), fixed = TRUE))
        )
      },
      domain = c(0, 1)
    )
  }

  dmr_range <- as.data.frame(.dmr)[, c("start", "end")]
  .dmp <- .dmp[.dmp[["cpg_pos"]] >= dmr_range[["start"]] & .dmp[["cpg_pos"]] <= dmr_range[["end"]], ]

  if (!is.na(.dmr[["overlapping.genes"]])) {
    p_gene <- autoplot(
      object = Homo.sapiens,
      which = GRanges(
        seqnames = .dmr[["seqnames"]], 
        ranges = IRanges(start = .dmr[["start"]], end = .dmr[["end"]], width = .dmr[["width"]]), 
        strand = .dmr[["strand"]]
      ),
      colour = base_colour,
      fill = base_colour
    )@ggplot +
      scale_y_continuous(expand = expansion(0.15), position = "right") +
      facet_grid(rows = vars("Gene"), switch = "y") +
      theme(plot.margin = unit(x = c(5.5, 5.5, 0, 5.5), units = "pt"))
  }

  p_estimate <- ggplot(
    data = .dmp,
    mapping = aes(x = cpg_pos, y = estimate, ymin = estimate - se, ymax = estimate + se)
  ) +
    geom_hline(yintercept = 0, linetype = 2, size = 0.2)  +
    geom_pointrange(size = 0.2, show.legend = FALSE, colour = base_colour) +
    geom_line(colour = base_colour) +
    scale_y_continuous(position = "right") +
    facet_grid(rows = vars("Estimate (M-value)"), switch = "y") +
    if (!is.na(.dmr[["overlapping.genes"]])) {
      theme(plot.margin = unit(x = c(0, 5.5, 0, 5.5), units = "pt"))
    } else {
      theme(plot.margin = unit(x = c(5.5, 5.5, 0, 5.5), units = "pt"))
    }

  p_pval <- ggplot(
    data = distinct(.dmp, CpG, cpg_chr, cpg_pos, pvalue), #cpg_strand
    mapping = aes(x = cpg_pos, y = pvalue)
  ) +
    geom_point(colour = base_colour) +
    annotate(
      geom = "rect",
      xmin = -Inf, xmax = Inf, ymin = 1, ymax = 0.05,
      fill = "firebrick2", alpha = 0.2, colour = NA
    ) +
    geom_hline(yintercept = 0.05, linetype = 2, colour = "firebrick2") +
    scale_y_continuous(
      position = "right",
      trans = pval_trans_inf(),
      expand = expansion(mult = c(0.002, 0.05))
    ) +
    facet_grid(rows = vars("P-value"), switch = "y") +
    theme(plot.margin = unit(x = c(0, 5.5, 0, 5.5), units = "pt"))

  p_density <- ggplot(
    data = distinct(.dmp, CpG, cpg_chr, cpg_pos), #cpg_strand,
    mapping = aes(x = cpg_pos)
  ) +
    geom_density(stat = "density", colour = base_colour, fill = base_colour, alpha = 0.6) +
    scale_y_continuous(
      expand = expansion(mult = c(0, 0.05)),
      position = "right",
      labels = function(x) {
        parse(
          text = gsub("e", " %*% 10^", gsub("0e+00", "0", gsub("1e+00", "1", scientific_format()(x), fixed = TRUE)))
        )
      }
    ) +
    facet_grid(rows = vars("Density of CpG"), switch = "y") +
    theme(plot.margin = unit(x = c(0, 5.5, 5.5, 5.5), units = "pt"))

  if (!is.na(.dmr[["overlapping.genes"]]) & length(p_gene$layers) != 0) {
    lp <- list(p_gene, p_estimate, p_pval, p_density)
  } else {
    lp <- list(p_estimate, p_pval, p_density)
  }

  lp <- map2(
    .x = lp,
    .y = as.list(c(rep(FALSE, length(lp) - 1), TRUE)),
    .range = range(unlist(map(lp, ~ layer_scales(.x)$x$range$range))),
    .f = function(p, x_axis, .range) {
      out <- p +
        scale_x_continuous(
          labels = comma_format(scale = 1e-6, suffix = " Mb", accuracy = 1e-3),
          limits = c(.range)
        ) +
        theme_light() +
        theme(
          axis.title = element_blank(),
          legend.title = element_blank(),
          strip.background = element_rect(fill = NA),
          strip.text = element_text(colour = "grey30")
        )
        if (!x_axis) {
          out <- out +
            theme(
              axis.text.x = element_blank(),
              axis.ticks.x = element_blank(),
              axis.title.x = element_blank()
            )
        }
      out
    }
  )

  lp <- map(lp, `+`,
    geom_rect(
      data = dmr_range,
      mapping = aes(
        xmin = start,
        xmax = end,
        ymin = -Inf,
        ymax = Inf
      ),
      inherit.aes = FALSE,
      fill = "dodgerblue",
      colour = NA,
      alpha = 0.1
    )
  )
  if (!is.na(.dmr[["overlapping.genes"]]) & length(p_gene$layers) != 0) {
    lp[[1]] <- lp[[1]] + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
      )
  }
  
  if (is.na(.dmr[["overlapping.genes"]])) {
    wrap_plots(lp, ncol = 1)
  } else {
    wrap_plots(lp, ncol = 1) + 
      plot_annotation(title = gsub("((?:[^,]*,){3}).*", "\\1 ...", .dmr[["overlapping.genes"]], perl = TRUE))
  }
}