R/oncoPathway.R

In PoisonAlien/maftools: Summarize, Analyze and Visualize MAF Files

#' Enrichment of known oncogenic or custom pathways
#' @description Checks for enrichment of known or custom pathways
#' @references
#' Sanchez-Vega F, Mina M, Armenia J, Chatila WK, Luna A, La KC, Dimitriadoy S, Liu DL, Kantheti HS, Saghafinia S et al. 2018. Oncogenic Signaling Pathways in The Cancer Genome Atlas. Cell 173: 321-337 e310
#' Bailey, Matthew H et al. “Comprehensive Characterization of Cancer Driver Genes and Mutations.” Cell vol. 173,2 (2018): 371-385.e18.
#'
#' @details
#' Oncogenic signalling and SMG pathways are derived from TCGA cohorts. See references for details.
#' @param maf an \code{\link{MAF}} object generated by \code{\link{read.maf}}
#' @param pathdb Either `sigpw` or `smgbp`,  `sigpw` uses known oncogenic signalling pathways (Sanchez/Vega et al) whereas `smgbp` uses pan cancer significantly mutated genes classified according to biological process (Bailey et al). Default \code{smgbp}
#' @param pathways Can be a two column data.frame/tsv-file with pathway-name and genes involved in them. Default `NULL`. This argument is mutually exclusive with \code{pathdb}
#' @param plotType Can be 'treemap' or 'bar'. Set NA to suppress plotting. Default NA
#' @param fontSize Default 1
#' @param panelWidths Default c(2, 4, 4)
#' @param col Default #f39c12
#' @return fraction of altered pathway. attr genes contain pathway contents
#' @seealso \code{\link{plotPathways}}
#' @export
#' @examples
#' laml.maf <- system.file("extdata", "tcga_laml.maf.gz", package = "maftools")
#' laml <- read.maf(maf = laml.maf)
#' pathways(maf = laml)

pathways = function(maf, pathdb = "sigpw", pathways = NULL, fontSize = 1, panelWidths = c(2, 4, 4), plotType = NA, col = "#f39c12"){


  if(is.null(pathways)){
    pathways = match.arg(arg = pathdb, choices = c("sigpw", "smgbp"))
  }else{
    if(is.data.frame(pathways)){
      colnames(pathways)[1:2] = c("Pathway", "Gene")
    } else if(file.exists(pathways)){
      pathways = data.table::fread(file = pathways)
      colnames(pathways)[1:2] = c("Pathway", "Gene")
    }
  }

  altered_pws = get_pw_summary(maf = maf, pathways = pathways)
  altered_pws = altered_pws[!n_affected_genes %in% 0][order(Fraction_mutated_samples)]

  if(nrow(altered_pws) == 0){
    stop("None of the provided pathways are altered!")
  }

  #return(altered_pws)
  if(!is.na(plotType)){
    plotType = match.arg(arg = plotType, choices = c("treemap", "bar"))

    if(plotType == "treemap"){
      x <- 0  # X-coordinate of the origin
      y <- 0  # Y-coordinate of the origin
      dx <- 100  # Full width of the treemap
      dy <- 100  # Full height of the treemap

      sizes_norm <- altered_pws$Mutated_samples / sum(altered_pws$Mutated_samples) * dx * dy  # Normalize the sizes

      # Compute the treemap rectangles
      rectangles <- squarify(sizes_norm, x, y, dx, dy)

      # Function to draw a rectangle
      draw_rect <- function(x, y, dx, dy, col) {
        rect(x, y, x + dx, y + dy, col = col, border = "#34495e")
      }

      # Function to draw text inside the rectangle
      draw_text <- function(x, y, dx, dy, label, fs = 1) {
        text(x + dx / 2, y + dy / 2, label, cex = fs, xpd = TRUE)
      }

      draw_text_sub <- function(x, y, dx, dy, label, fs = 1) {
        text(x + dx/2, y + dy / 2, label, cex = fs * 0.8, pos = 1, xpd = TRUE)
      }

      # Plot the treemap using base plotting functions
      #par(mar = c(1, 1, 1, 1))
      plot(NA, xlim = c(0, dx), ylim = c(0, dy), xlab = "", ylab = "", xaxt = "n", yaxt = "n", frame.plot = FALSE)

      sub_text = paste0("N:", altered_pws$Mutated_samples,  " [", round(altered_pws$Fraction_mutated_samples * 100, 2), "%]")
      for (i in seq_along(rectangles)) {
        rect_info <- rectangles[[i]]
        draw_rect(rect_info$x, rect_info$y, rect_info$dx, rect_info$dy, col)
        draw_text(rect_info$x, rect_info$y, rect_info$dx, rect_info$dy, altered_pws$Pathway[i], fs = fontSize) # Use original sizes
        draw_text_sub(rect_info$x, rect_info$y, rect_info$dx, rect_info$dy, sub_text[i], fs = fontSize) # Use original sizes
      }
    }else if(plotType == "bar"){
      nsamps = as.numeric(maf@summary[ID == "Samples", summary])

      lo = layout(mat = matrix(c(1, 2, 3), ncol = 3, nrow = 1), widths = panelWidths)
      par(mar = c(2, 2, 2, 0))
      plot(NA, xlim = c(0, 1), ylim = c(0, nrow(altered_pws)), axes = FALSE)
      text(x = 1, y = seq(0.5, nrow(altered_pws), by = 1), labels = altered_pws$Pathway, adj = 1, xpd = TRUE, cex = fontSize)
      title(main = "Pathway", adj = 0)

      par(mar = c(2, 0, 2, 1), xpd = TRUE)
      plot(NA, xlim = c(0, 1.2), ylim = c(0, nrow(altered_pws)), axes = FALSE)
      rect(xleft = 0, ybottom = seq(0.1, nrow(altered_pws), by = 1), xright = 1, ytop = seq(0.2, nrow(altered_pws), by = 1)+0.7, col = '#ecf0f1', border = "white")
      rect(xleft = 0, ybottom = seq(0.1, nrow(altered_pws), by = 1), xright = altered_pws$fraction_affected, ytop = seq(0.2, nrow(altered_pws), by = 1)+0.7, col = col, border = "white")
      text(x = 1.05, y = seq(0.5, nrow(altered_pws), by = 1), labels = paste0(altered_pws$n_affected_genes, "/", altered_pws$N), adj = 0, cex = fontSize)
      axis(side = 1, at = seq(0, 1, 0.25), line = -1, cex.axis = fontSize)
      title(main = "Fraction of pathway affected", adj = 0)

      par(mar = c(2, 0, 2, 1), xpd = TRUE)
      plot(NA, xlim = c(0, 1.2), ylim = c(0, nrow(altered_pws)), axes = FALSE)
      rect(xleft = 0, ybottom = seq(0.1, nrow(altered_pws), by = 1), xright = 1, ytop = seq(0.2, nrow(altered_pws), by = 1)+0.7, col = '#ecf0f1', border = "white")
      rect(xleft = 0, ybottom = seq(0.1, nrow(altered_pws), by = 1), xright = altered_pws$Fraction_mutated_samples, ytop = seq(0.2, nrow(altered_pws), by = 1)+0.7, col = col, border = "white")
      text(x = 1.05, y = seq(0.5, nrow(altered_pws), by = 1), labels = paste0(altered_pws$Mutated_samples, "/", nsamps), adj = 0, cex = fontSize)
      axis(side = 1, at = seq(0, 1, 0.25), line = -1, cex.axis = fontSize)
      title(main = "Fraction of samples affected")
    }
  }

  altered_pws[order(-Fraction_mutated_samples)]
}