R/genomewide.log10q.plot.R

In GRIN2: Genomic Random Interval (GRIN)

#' Genome-wide -log10(q-value) Plot
#'
#' @description
#' Generates a genome-wide plot of -log10(q-values) for each annotated gene or
#' lesion boundary evaluated by GRIN. Statistical significance can be displayed
#' for one or more selected lesion types.
#'
#' @usage
#' genomewide.log10q.plot(grin.res,
#'                        lsn.grps,
#'                        lsn.colors = NULL,
#'                        max.log10q = NULL)
#'
#' @param grin.res GRIN results object (output from `grin.stats`) generated
#' using either gene annotation or lesion boundaries as marker input.
#' @param lsn.grps A character vector specifying the lesion type(s) to include
#' in the plot.
#' @param lsn.colors A named vector of colors corresponding to the selected
#' lesion types. If `NULL`, colors are automatically assigned using
#' `default.grin.colors`.
#' @param max.log10q Numeric; optional maximum value of -log10(q-value)
#' displayed on the plot. Values greater than `max.log10q` are capped at this
#' value. If `NULL`, the plotting limit is determined automatically from the
#' observed finite -log10(q-values).
#'
#' @details
#' This function displays the genome-wide statistical significance of lesions
#' affecting annotated genomic markers. Depending on the marker data supplied
#' to `grin.stats`, these markers may represent genes or lesion boundaries.
#'
#' The function first adds continuous genome-wide plotting coordinates using
#' `compute.gw.coordinates` when these coordinates are not already present in
#' `grin.res`.
#'
#' Chromosomes are arranged consecutively along the vertical axis. For each
#' selected lesion type, a horizontal line is drawn at the genomic position of
#' each affected marker. Line length represents the corresponding
#' -log10(q-value), with longer lines indicating greater statistical
#' significance, and line color identifies the lesion type.
#'
#' @return
#' Generates a genome-wide -log10(q-value) plot on the active graphics device
#' and invisibly returns `NULL`. Chromosomes are displayed along the vertical
#' axis, while the horizontal axis represents -log10(q-value). Each horizontal
#' line corresponds to an affected gene or lesion boundary, with line length
#' representing statistical significance and line color indicating lesion type.
#'
#' @export
#'
#' @importFrom graphics plot rect legend text segments
#'
#' @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{grin.stats}},
#' \code{\link{grin.lsn.boundaries}},
#' \code{\link{genomewide.lsn.plot}},
#' \code{\link{compute.gw.coordinates}},
#' \code{\link{default.grin.colors}}
#'
#' @examples
#' data(lesion_data)
#' data(hg38_gene_annotation)
#' data(hg38_chrom_size)
#'
#' # Use lesion boundaries as genomic markers for gain lesions
#' gain <- lesion_data[lesion_data$lsn.type == "gain", ]
#' lsn.bound.gain <- grin.lsn.boundaries(gain,
#'                                       hg38_chrom_size)
#'
#' GRIN.results.gain.bound <- grin.stats(gain,
#'                                       lsn.bound.gain,
#'                                       hg38_chrom_size)
#'
#' # Plot genome-wide significance of gain lesion boundaries
#' genomewide.log10q.plot(GRIN.results.gain.bound,
#'                        lsn.grps = "gain",
#'                        lsn.colors = c("gain" = "red"),
#'                        max.log10q = 10)
#'
#' # Gene annotation can also be used as the marker input to grin.stats instead
#' # of lesion boundaries.
#' # Multiple lesion types can be displayed together by including their names
#' # in lsn.grps.
#'
genomewide.log10q.plot=function(grin.res,        # GRIN results (output of the grin.stats function)
                                lsn.grps,        # selected lesion groups to be added to the plot
                                lsn.colors=NULL, # Lesion colors
                                max.log10q=NULL) # Maximum -log10(q) value to be added to the plot
{

  if (!is.list(grin.res))
    stop("'grin.res' must be GRIN results returned by grin.stats().")

  required.objects=c("lsn.data","gene.hits","chr.size")

  if (!all(required.objects %in% names(grin.res)))
    stop(
      "'grin.res' must contain the following components: ",
      paste(required.objects,collapse=", "),
      "."
    )

  if (missing(lsn.grps) ||
      is.null(lsn.grps) ||
      length(lsn.grps)==0 ||
      anyNA(lsn.grps) ||
      any(lsn.grps==""))
    stop("'lsn.grps' must contain at least one valid lesion type.")

  lsn.grps=unique(as.character(lsn.grps))

  required.nsubj=paste0("nsubj.",lsn.grps)
  required.qval=paste0("q.nsubj.",lsn.grps)

  if (!all(required.nsubj %in% colnames(grin.res$gene.hits)) ||
      !all(required.qval %in% colnames(grin.res$gene.hits)))
    stop(
      "Selected lesion types in 'lsn.grps' were not found in the GRIN results."
    )

  if (!is.null(max.log10q))
  {
    if (!is.numeric(max.log10q) ||
        length(max.log10q)!=1 ||
        is.na(max.log10q) ||
        !is.finite(max.log10q) ||
        max.log10q<=0)
      stop("'max.log10q' must be a single positive numeric value or NULL.")
  }

  gw.coords.available=
    all(c("x.start","x.end") %in% colnames(grin.res$gene.hits)) &&
    all(c("x.start","x.end") %in% colnames(grin.res$chr.size))

  if (!gw.coords.available)
    grin.res=compute.gw.coordinates(grin.res)

  grin.res$lsn.data$x.ID=as.numeric(as.factor(grin.res$lsn.data$ID))
  n=max(grin.res$lsn.data$x.ID)
  n.chr=nrow(grin.res$chr.size)

  # set up plotting region

  graphics::plot(c(-0.05,1.2)*n,
                 c(0.03*grin.res$chr.size$x.end[n.chr],
                   -1.1*grin.res$chr.size$x.end[n.chr]),
                 type="n",axes=FALSE,xlab="",ylab="")

  # background colors for chromosomes

  graphics::rect(0*n,-grin.res$chr.size$x.start,
                 1*n,-grin.res$chr.size$x.end,
                 col=c("lightgray","gray"),
                 border=NA)

  lsn.types=lsn.grps

  if (is.null(lsn.colors))
  {
    lsn.colors=default.grin.colors(lsn.types)
  } else {

    if (is.null(names(lsn.colors)))
      stop("'lsn.colors' must be a named vector with names corresponding to lesion types.")

    missing.colors=setdiff(lsn.types,names(lsn.colors))

    if (length(missing.colors)>0)
      stop(
        "'lsn.colors' does not contain colors for the following lesion types: ",
        paste(missing.colors,collapse=", "),
        "."
      )

    lsn.colors=lsn.colors[lsn.types]
  }

  nsubj.mtx=unlist(grin.res$gene.hits[,paste0("nsubj.",lsn.types)])
  qval.mtx=unlist(grin.res$gene.hits[,paste0("q.nsubj.",lsn.types)])

  nsubj.data=cbind.data.frame(
    gene=grin.res$gene.hits$gene,
    x.start=grin.res$gene.hits$x.start,
    x.end=grin.res$gene.hits$x.end,
    nsubj=nsubj.mtx,
    log10q=-log10(qval.mtx),
    lsn.type=rep(lsn.types,each=nrow(grin.res$gene.hits))
  )

  nsubj.data=nsubj.data[nsubj.data$nsubj>0,]

  if (nrow(nsubj.data)==0)
    stop("No markers affected by the selected lesion types were found to plot.")

  nsubj.data$lsn.colors=lsn.colors[nsubj.data$lsn.type]

  nsubj.data=nsubj.data[!is.na(nsubj.data$log10q),]

  if (nrow(nsubj.data)==0)
    stop("No non-missing q-values were found to plot.")

  if (is.null(max.log10q))
  {
    finite.log10q=nsubj.data$log10q[is.finite(nsubj.data$log10q)]

    if (length(finite.log10q)==0)
      stop(
        "All plotted q-values are zero. Please provide a finite value for 'max.log10q'."
      )

    max.log10q=max(finite.log10q)
  }

  nsubj.data$log10q[nsubj.data$log10q>max.log10q]=max.log10q

  ord=order(nsubj.data$log10q,decreasing=TRUE)
  nsubj.data=nsubj.data[ord,]

  log10q.scale=max(nsubj.data$log10q)

  if (log10q.scale==0)
    log10q.scale=1

  graphics::segments(0*n,
                     -(nsubj.data$x.start+nsubj.data$x.end)/2,
                     0*n+1*nsubj.data$log10q/log10q.scale*n,
                     col=nsubj.data$lsn.colors)

  graphics::text(c(n,0)[(1:n.chr)%%2+1],
                 pos=c(4,2)[(1:n.chr)%%2+1],
                 -(grin.res$chr.size$x.start+grin.res$chr.size$x.end)/2,
                 grin.res$chr.size$chrom,
                 cex=0.75)

  graphics::legend(n/2,-1.07*grin.res$chr.size$x.end[n.chr],
                   fill=lsn.colors,cex=0.9,
                   legend=names(lsn.colors),
                   xjust=0.5,
                   ncol=length(lsn.colors),
                   border=NA,bty="n")

  graphics::text(1*n/2,0,
                 "-log10(q)",cex=0.95,
                 pos=3)

  graphics::text(c(0,0.25,0.5,0.75,1)*n,
                 -grin.res$chr.size$x.end[n.chr],
                 c(0,
                   max.log10q/4,
                   max.log10q/2,
                   round(max.log10q*0.75,1),
                   max.log10q),
                 cex=0.75,pos=1)

  invisible(NULL)
}