R/alex.boxplots.R

In GRIN2: Genomic Random Interval (GRIN)

#' Generate Box Plots of Gene Expression by Lesion Groups
#'
#' @description
#' Generates box plots of gene expression levels stratified by lesion groups
#' for genes selected according to a user-specified q-value threshold from the
#' Kruskal-Wallis lesion-expression association results.
#'
#' @usage
#' alex.boxplots(
#'   out.dir,
#'   alex.data,
#'   alex.kw.results,
#'   q,
#'   gene.annotation
#' )
#'
#' @param out.dir Path to an existing directory where the resulting PDF files
#' containing the box plots will be saved. One PDF file is generated for each
#' selected gene.
#'
#' @param alex.data Output from the \code{\link{alex.prep.lsn.expr}} function.
#' A list containing \code{"alex.expr"} (expression data), \code{"alex.lsn"}
#' (lesion-group data), and \code{"alex.row.mtch"} (matched Ensembl gene IDs).
#' The expression and lesion data contain genes in rows and subjects in columns,
#' with matching row and column order.
#'
#' @param alex.kw.results Kruskal-Wallis lesion-expression association results
#' returned by \code{\link{KW.hit.express}}. The data frame must contain
#' \code{"gene"} and \code{"q.KW"}.
#'
#' @param q Numeric q-value threshold. Genes with \code{q.KW < q} are included
#' in the output box plots.
#'
#' @param gene.annotation A gene annotation data frame containing at least
#' \code{"gene"} (Ensembl gene ID) and \code{"gene.name"} (gene symbol).
#' When a gene symbol is unavailable, the Ensembl gene ID is used to label
#' the corresponding plot.
#'
#' @details
#' For each gene with a Kruskal-Wallis q value below the specified threshold,
#' expression values are grouped according to lesion status. Subjects are
#' displayed as individual points together with a box plot showing the
#' expression distribution within each lesion group.
#'
#' Lesion groups are ordered according to their expression values using
#' \code{\link[forcats]{fct_reorder}}. One PDF file is generated for each
#' selected gene and saved in \code{out.dir}.
#'
#' @return
#' Generates one PDF file for each selected gene in \code{out.dir}. Each file
#' contains a box plot showing gene expression across lesion groups.
#'
#' @export
#'
#' @importFrom ggplot2 aes geom_boxplot geom_jitter position_jitter xlab ylab theme_bw scale_fill_discrete guide_legend theme element_text ggplot
#' @importFrom forcats fct_reorder
#' @importFrom grDevices dev.off pdf
#'
#' @author
#' Abdelrahman Elsayed \email{abdelrahman.elsayed@stjude.org},
#' Stanley Pounds \email{stanley.pounds@stjude.org}
#'
#' @references
#' Cao, X., Elsayed, A. H., & Pounds, S. B. (2023). Statistical Methods Inspired
#' by Challenges in Pediatric Cancer Multi-omics.
#'
#' @seealso
#' \code{\link{alex.prep.lsn.expr}},
#' \code{\link{KW.hit.express}}
#'
#' @examples
#' data(expr_data)
#' data(lesion_data)
#' data(hg38_gene_annotation)
#'
#' # Prepare expression and lesion data
#' alex.data <- alex.prep.lsn.expr(expr_data,
#'                                 lesion_data,
#'                                 hg38_gene_annotation,
#'                                 min.expr = 5,
#'                                 min.pts.lsn = 5)
#'
#' # Run Kruskal-Wallis test
#' alex.kw.results <- KW.hit.express(alex.data,
#'                                   hg38_gene_annotation,
#'                                   min.grp.size = 5)
#'
#' # Generate box plots for significant genes
#' dir.create(resultsFolder <- file.path(tempdir(), "temp.out"),
#'            showWarnings = FALSE)
#'
#' alex.boxplots(out.dir = resultsFolder,
#'               alex.data = alex.data,
#'               alex.kw.results = alex.kw.results,
#'               q = 1e-15,
#'               gene.annotation = hg38_gene_annotation)
#'
#' unlink(resultsFolder, recursive = TRUE)
alex.boxplots=function(out.dir,            # Path to the folder where the boxplots of selected genes based on the specified q value of the KW results table will be added
                       alex.data,          # output of the alex.prep.lsn.expr function (list of three data tables "alex.expr" with expression data, "alex.lsn" with overlapped gene lesion data and row.mtch)
                       alex.kw.results,    # ALEX Kruskal-Wallis test results (output of the KW.hit.express function)
                       q,                  # minimum q value for a gene to be included in output PDF file of box plots
                       gene.annotation)    # gene annotation data

{

  # Validate input data

  if (!dir.exists(out.dir))
    stop("out.dir must be an existing directory.")

  if (!is.list(alex.data))
    stop("alex.data must be the output from alex.prep.lsn.expr().")

  required.alex.objects=c("alex.expr","alex.lsn","alex.row.mtch")
  if (!all(required.alex.objects %in% names(alex.data)))
    stop("alex.data must contain: alex.expr, alex.lsn, and alex.row.mtch.")

  if (!is.data.frame(alex.kw.results))
    stop("alex.kw.results must be a data frame.")

  if (!all(c("gene","q.KW") %in% colnames(alex.kw.results)))
    stop("alex.kw.results must contain the columns: gene and q.KW.")

  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.")

  if (!is.data.frame(gene.annotation))
    stop("gene.annotation must be a data frame.")

  if (!all(c("gene","gene.name") %in% colnames(gene.annotation)))
    stop("gene.annotation must contain the columns: gene and gene.name.")

  # To extract genes below specified q.value for boxplots:

  KW.results=alex.kw.results[!is.na(alex.kw.results$q.KW),]
  selected.genes=KW.results[KW.results$q.KW<q,]

  if (nrow(selected.genes)==0)
    stop("No genes have q.KW below the specified q-value threshold.")

  selected.IDS=selected.genes$gene
  selected.lsns=alex.data$alex.lsn[rownames(alex.data$alex.lsn) %in% selected.IDS,,drop=FALSE]
  selected.lsns=t(selected.lsns)
  selected.expr=alex.data$alex.expr[rownames(alex.data$alex.expr) %in% selected.IDS,,drop=FALSE]
  selected.expr=t(selected.expr)

  # To make sure that both lesion and expression matrices have same patients and genes order

  check.rownames=all(rownames(selected.lsns)==rownames(selected.expr))
  if (check.rownames==FALSE)
    stop("Gene-lesion matrix patient IDs must match expression matrix patient IDs.")

  check.colnames=all(colnames(selected.lsns)==colnames(selected.expr))
  if (check.colnames==FALSE)
    stop("Gene-lesion matrix gene IDs must match expression matrix gene IDs.")

  # add ensembl.ID to gene.name column if empty

  #gene.annotation=gene.annotation %>% mutate_all(na_if,"")
  gene.annotation$gene.name <- ifelse(is.na(gene.annotation$gene.name) | gene.annotation$gene.name=="",gene.annotation$gene, gene.annotation$gene.name)

  n=ncol(selected.lsns)

  for (i in 1:n)
  {
    ensembl.id=colnames(selected.lsns)[i]
    ensembl.df=as.data.frame(ensembl.id)
    colnames(ensembl.df) <- ("gene")
    genes.df.merged=merge(gene.annotation,ensembl.df,by="gene", all.y=TRUE)
    gene.names=as.character(genes.df.merged$gene.name)
    lsn.box=unlist(selected.lsns[,i])
    expr.box=unlist(selected.expr[,i])
    df=as.data.frame(cbind(lsn.box, expr.box))
    df$expr.box=as.numeric(df$expr.box)

    {
      grDevices::pdf(file.path(out.dir,paste0(gene.names,"_boxplot.pdf")),
                     height=8,width=10)
      p=ggplot2::ggplot(df, ggplot2::aes(x = forcats::fct_reorder(lsn.box, expr.box, .desc =TRUE), y = expr.box)) +
        ggplot2::geom_boxplot(ggplot2::aes(fill = forcats::fct_reorder(lsn.box, expr.box,  .desc =TRUE))) +
        ggplot2::geom_jitter(position=ggplot2::position_jitter(0.02)) + # Increasing the number from 0.02 increase the space between the dots but add some extra dots if number of patients is small
        ggplot2::xlab(paste0(gene.names, '_lsn'))+
        ggplot2::ylab(paste0(gene.names, '_Expression'))+
        ggplot2::theme_bw(base_size = 12) +
        ggplot2::scale_fill_discrete(guide = ggplot2::guide_legend(title = "Lesion")) +
        ggplot2::theme(axis.text=ggplot2::element_text(size=10), axis.title=ggplot2::element_text(size=12,face="bold"))+
        ggplot2::theme(legend.text=ggplot2::element_text(size=8))
      print(p)
      grDevices::dev.off()
    }

  }
}