R/drugInteractions.R

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

#' Drug-Gene Interactions
#' @description Checks for drug-gene interactions and druggable categories
#' @references Griffith, M., Griffith, O. L., Coffman, A. C., Weible, J. V., McMichael, J. F., Spies, N. C., et. al,. 2013. DGIdb - Mining the druggable genome. Nature Methods.
#' @details This function takes a list of genes and checks for known/reported drug-gene interactions or Druggable categories.
#' All gene-drug interactions and drug claims are compiled from Drug Gene Interaction Databse. See reference for details and cite it if you use this function.
#' @param maf an \code{\link{MAF}} object generated by \code{\link{read.maf}}
#' @param top Top number genes to check for. Default 20
#' @param genes Manually specify gene list
#' @param plotType Can be bar, pie Default bar plot.
#' @param drugs Check for known/reported drugs. Default FALSE
#' @param fontSize Default 0.8
#' @export
#' @examples
#' laml.maf <- system.file("extdata", "tcga_laml.maf.gz", package = "maftools")
#' laml <- read.maf(maf = laml.maf)
#' drugInteractions(maf = laml)

drugInteractions = function(maf, top = 20, genes = NULL, plotType = "bar", drugs = FALSE, fontSize = 0.8){

  if(drugs){
    db = system.file('extdata', 'drugs.tsv.gz', package = 'maftools')
  }else{
    db = system.file('extdata', 'categories.tsv.gz', package = 'maftools')
  }

  if(Sys.info()[['sysname']] == 'Windows'){
    db.gz = gzfile(description = db, open = 'r')
    suppressWarnings(db <- data.table::as.data.table(read.csv(file = db.gz, header = TRUE, sep = '\t', stringsAsFactors = FALSE)))
    close(db.gz)
  } else{
    db = suppressWarnings(data.table::fread(cmd = paste('zcat <', db), sep = '\t', stringsAsFactors = FALSE, verbose = FALSE, data.table = TRUE, showProgress = TRUE, header = TRUE, fill = TRUE))
  }

  if(is.null(genes)){
    genes = unique(getGeneSummary(x = maf)[1:top,Hugo_Symbol])
  }

  if(drugs){
    db = db[Gene %in% genes]
    if(nrow(db) == 0){
      stop("No claimed drugs found for given genes.")
    }

    ndb = db[order(Gene)][,.N,Gene][order(N, decreasing = TRUE)]

    cat("Number of claimed drugs for given genes:\n")
    print(ndb)
    return(db)
  }else{
    db = db[Gene %in% genes]
    if(nrow(db) == 0){
      stop("No druggable categories found for given genes.")
    }

    db = db[Gene %in% genes]
    ndb = db[,.N,category][order(N, decreasing = TRUE)]
    dbg = db[,paste(Gene, collapse = ", "), category][,.(category, V1)]

    dbg$V1 = unlist(lapply(strsplit(dbg$V1, split = ", "), function(x){
      x = unlist(x)
      x = x[!is.na(x)]
      x = names(sort(table(x), decreasing = TRUE))
      if(length(x) >=5){
        return(paste(x[1:5], collapse = ","))
      }else{
        return(paste(x, collapse = ","))
      }
    }))

    ndb = merge(ndb, dbg)[order(N, decreasing = TRUE)]
    ndb[,label := paste0(category, " [", V1, "]")]

    if(plotType == "pie"){
      par(mar = c(2, 5, 2, 5))
      pie(x = ndb$N, labels = ndb$category, clockwise = TRUE,
          border = NA, lty = 1, col = heat.colors(n = nrow(ndb), alpha = 0.8), cex = fontSize)
    }else if(plotType == "bar"){
      par(mar = c(2, 2, 3, 4))
      b = barplot(height = ndb$N, names.arg = NA, las = 2,
                  col = heat.colors(n = nrow(ndb), alpha = 0.8),
                  border = NA, axes = FALSE, horiz = TRUE)
      #mtext(text = ndb$category, side = 2, at = b, las = 2, adj = 0, cex = 0.5)
      text(labels = ndb$label, las = 2, adj = 0, cex = fontSize, x = 0.1, y = b, font = 2)
      axis(side = 1, at = c(0, max(ndb$N)))
      mtext(text = "# Genes", side = 1, cex = 1, line = 0.5)
      title(main = "Druggable categories", adj = 0)
    }
    return(db)
  }
}