R/top.alex.waterfall.plots.R

In GRIN2: Genomic Random Interval (GRIN)

#' Generate Waterfall Plots for Significant Genes
#'
#' @description
#' Generates waterfall plots for genes showing significant associations between
#' lesion status and gene expression based on Kruskal-Wallis (KW) test results.
#' Genes with KW q-values below a user-specified threshold are selected for
#' plotting.
#'
#' @usage
#' top.alex.waterfall.plots(out.dir,
#'                          alex.data,
#'                          alex.kw.results,
#'                          q,
#'                          lsn.data)
#'
#' @param out.dir Character string specifying the directory where the waterfall
#' plot PDF files will be saved. The directory must already exist.
#' @param alex.data A list returned by `alex.prep.lsn.expr` containing the
#' matched expression and lesion data used for lesion-expression association
#' analyses.
#' @param alex.kw.results Kruskal-Wallis lesion-expression association results
#' returned by `KW.hit.express`.
#' @param q Numeric; KW q-value threshold used to select genes for plotting.
#' Genes with `q.KW < q` are included.
#' @param lsn.data Lesion data in GRIN-compatible format, as used in
#' `alex.prep.lsn.expr`.
#'
#' @details
#' The function selects genes from `alex.kw.results` with non-missing
#' `q.KW` values below the specified `q` threshold and with available gene
#' symbols. Duplicate gene symbols are represented once.
#'
#' For each selected gene, `alex.waterfall.prep` prepares the corresponding
#' lesion and expression data and `alex.waterfall.plot` generates a waterfall
#' plot showing gene expression across patients according to lesion status.
#'
#' Each gene is saved as a separate PDF file named
#' `<gene>_waterfall_plot.pdf` in `out.dir`.
#'
#' @return
#' Invisibly returns `NULL`. The function creates one PDF waterfall plot for
#' each selected gene in the specified output directory. If no genes meet the
#' specified q-value threshold, no PDF files are created.
#'
#' @export
#'
#' @importFrom grDevices pdf dev.off
#'
#' @references
#' Cao, X., Elsayed, A. H., & Pounds, S. B. (2023). Statistical Methods
#' Inspired by Challenges in Pediatric Cancer Multi-omics.
#'
#' @author
#' Abdelrahman Elsayed \email{abdelrahman.elsayed@stjude.org} and
#' Stanley Pounds \email{stanley.pounds@stjude.org}
#'
#' @seealso \code{\link{alex.prep.lsn.expr}},
#' \code{\link{KW.hit.express}},
#' \code{\link{alex.waterfall.prep}},
#' \code{\link{alex.waterfall.plot}}
#'
#' @examples
#' data(expr_data)
#' data(lesion_data)
#' data(hg38_gene_annotation)
#'
#' # 1) Prepare matched expression and lesion matrices
#' alex.data <- alex.prep.lsn.expr(expr_data,
#'                                 lesion_data,
#'                                 hg38_gene_annotation,
#'                                 min.expr = 5,
#'                                 min.pts.lsn = 5)
#'
#' # 2) Run Kruskal-Wallis lesion-expression association analysis
#' alex.kw.results <- KW.hit.express(alex.data,
#'                                   hg38_gene_annotation,
#'                                   min.grp.size = 5)
#'
#' # 3) Create a temporary output directory and generate waterfall plots
#' resultsFolder <- file.path(tempdir(), "temp.out")
#' dir.create(resultsFolder, showWarnings = FALSE, recursive = TRUE)
#'
#' top.alex.waterfall.plots(out.dir = resultsFolder,
#'                          alex.data = alex.data,
#'                          alex.kw.results = alex.kw.results,
#'                          q = 1e-15,
#'                          lsn.data = lesion_data)
#'
#' # Clean up
#' unlink(resultsFolder, recursive = TRUE)
#'
top.alex.waterfall.plots=function(out.dir,            # path to the folder to add waterfall plots
                                  alex.data,          # output of the alex.prep.lsn.expr function
                                  alex.kw.results,    # Kruskal-Wallis test results (output of the KW.hit.express function)
                                  q,                  # q-value threshold for genes to be plotted
                                  lsn.data)           # Lesion data in a GRIN compatible format
{

  if (!dir.exists(out.dir))
    stop("'out.dir' does not exist. Please provide an existing output directory.")

  if (!is.numeric(q) || length(q)!=1 || is.na(q) || q<0 || q>1)
    stop("'q' must be a single numeric value between 0 and 1.")

  # extract genes below specified q-value for waterfall plots

  KW.results=alex.kw.results[!is.na(alex.kw.results$q.KW),]
  selected.genes=KW.results[KW.results$q.KW<q,]
  selected.genes=selected.genes[!is.na(selected.genes$gene.name),]
  selected.genes=selected.genes[selected.genes$gene.name!="",]
  selected.genes=selected.genes[!duplicated(selected.genes[,"gene.name"]),]
  top.genes=selected.genes$gene.name

  for (i in seq_along(top.genes))
  {
    grDevices::pdf(file.path(out.dir,
                             paste0(top.genes[i],"_waterfall_plot.pdf")),
                   height=8,width=10)

    tryCatch({
      temp=alex.waterfall.prep(alex.data,alex.kw.results,top.genes[i],lsn.data)
      plot.gene=alex.waterfall.plot(temp,lsn.data)
    },
    finally=grDevices::dev.off())
  }

  invisible(NULL)
}