R/genomewide.lsn.plot.R

In GRIN2: Genomic Random Interval (GRIN)

#' Genome-wide Lesion Plot
#'
#' @description
#' Generates a genome-wide lesion plot displaying genomic lesions across
#' chromosomes and patients, together with lesion recurrence and statistical
#' significance from the GRIN analysis.
#'
#' @usage
#' genomewide.lsn.plot(grin.res,
#'                     ordered = FALSE,
#'                     pt.order = NULL,
#'                     lsn.colors = NULL,
#'                     max.log10q = NULL)
#'
#' @param grin.res GRIN results (output from the `grin.stats` function).
#' @param ordered Logical; if `TRUE`, patients are displayed according to the
#' order specified in `pt.order`. If `FALSE` (default), patient IDs are ordered
#' alphabetically.
#' @param pt.order A data frame with two columns: `"ID"` containing patient
#' identifiers matching those in the lesion data, and `"pts.order"` containing
#' the numeric plotting order from 1 to the total number of patients. Required
#' only when `ordered = TRUE`.
#' @param lsn.colors A named vector of colors assigned to lesion types. Names
#' must correspond to lesion types in the GRIN results. If `NULL`, colors are
#' automatically assigned using `default.grin.colors`.
#' @param max.log10q Numeric; optional maximum value of -log10(q-value) displayed
#' in the left panel. 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
#' The function first uses \code{compute.gw.coordinates} to convert
#' chromosome-specific genomic coordinates into continuous genome-wide
#' coordinates used for plotting lesions across the genome.
#'
#' The resulting plot contains three aligned panels. The center panel displays
#' genomic lesions across patients and chromosomes. Each patient occupies one
#' column, and lesions are colored according to lesion type. The left panel
#' displays the -log10(q-value) for each gene-lesion type combination from the
#' GRIN analysis, whereas the right panel displays the corresponding number of
#' affected patients.
#'
#' By default, patients are ordered alphabetically by patient ID. A user-defined
#' patient order can instead be supplied through `pt.order` by setting
#' `ordered = TRUE`.
#'
#' @return
#' Generates a genome-wide lesion plot on the active graphics device and
#' invisibly returns `NULL`. The plot consists of three aligned panels:
#' \itemize{
#'   \item Left panel: Displays the statistical significance of lesion
#'   recurrence for each gene and lesion type as -log10(q-value). Horizontal
#'   lines are colored by lesion type, with longer lines representing greater
#'   statistical significance.
#'   \item Middle panel: Displays the genome-wide distribution of lesions across
#'   chromosomes and patients. Each column represents a patient, chromosomes
#'   are arranged consecutively along the vertical axis, and lesions are colored
#'   according to lesion type.
#'   \item Right panel: Displays the number of patients affected by each lesion
#'   type at each gene. Horizontal lines are colored by lesion type, with longer
#'   lines representing a greater number of affected patients.
#' }
#'
#' @export
#'
#' @importFrom graphics plot rect legend text segments
#'
#' @references
#' Pounds, S., et al. (2013). A genomic random interval model for statistical
#' analysis of genomic lesion data.
#'
#' 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{compute.gw.coordinates}},
#' \code{\link{default.grin.colors}}
#'
#' @examples
#' data(lesion_data)
#' data(hg38_gene_annotation)
#' data(hg38_chrom_size)
#'
#' # Run GRIN analysis
#' grin.results <- grin.stats(lesion_data,
#'                            hg38_gene_annotation,
#'                            hg38_chrom_size)
#'
#' # Generate genome-wide lesion plot with alphabetical patient ordering
#' genomewide.lsn.plot(grin.results,
#'                     max.log10q = 50)
#'
#'
#' # Generate plot using a user-defined patient order, such as grouping patients
#  # by molecular subtype, clinical subgroup, or another user-defined criterion
#' patient.ids <- sort(unique(lesion_data$ID))
#' pt.order <- data.frame(ID = patient.ids,
#'                        pts.order = seq_along(patient.ids))
#'
#' genomewide.lsn.plot(grin.results,
#'                     ordered = TRUE,
#'                     pt.order = pt.order,
#'                     max.log10q = 50)
#'
genomewide.lsn.plot=function(grin.res,        # GRIN results (output of the grin.stats function)
                             ordered=FALSE,   # user can specify a certain patient order in the genome-wide lesion plot by specifying ordered = TRUE and passing pt.order data.frame with the new patient order, otherwise patients will be ordered alphabetically
                             pt.order=NULL,   # data.frame with two columns "ID" that has the patient ID and "pts.order" that has the patient order in numbers
                             lsn.colors=NULL, # Lesion colors (If not provided by the user, colors will be automatically assigned using default.grin.colors function).
                             max.log10q=NULL) # Maximum -log10(q) value for genes in the GRIN results table to be added to the plot. Values above max.log10q will be capped at this value
{
  if (!is.list(grin.res))
    stop("'grin.res' must be GRIN results returned by grin.stats().")

  required.objects=c("lsn.data","lsn.index","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 (!is.logical(ordered) || length(ordered)!=1 || is.na(ordered))
    stop("'ordered' must be either TRUE or FALSE.")

  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$lsn.data)) &&
    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)

  if (ordered == FALSE) {
    grin.res$lsn.data$pts.order=as.numeric(as.factor(grin.res$lsn.data$ID))
    n=max(grin.res$lsn.data$pts.order)
    n.chr=nrow(grin.res$chr.size)
  } else {

    if (is.null(pt.order))
      stop("'pt.order' must be provided when 'ordered = TRUE'.")

    if (!is.data.frame(pt.order))
      stop("'pt.order' must be a data frame.")

    if (!all(c("ID","pts.order") %in% colnames(pt.order)))
      stop("'pt.order' must contain columns named 'ID' and 'pts.order'.")

    if (anyNA(pt.order$ID) ||
        any(pt.order$ID=="") ||
        anyDuplicated(pt.order$ID))
      stop("'pt.order$ID' must contain unique, non-missing patient identifiers.")

    if (!is.numeric(pt.order$pts.order) ||
        anyNA(pt.order$pts.order) ||
        anyDuplicated(pt.order$pts.order))
      stop("'pt.order$pts.order' must contain unique numeric patient-order values.")

    if (!setequal(pt.order$pts.order,seq_len(nrow(pt.order))))
      stop("'pt.order$pts.order' must contain consecutive values from 1 to the number of patients.")

    missing.pts=setdiff(unique(grin.res$lsn.data$ID),pt.order$ID)

    if (length(missing.pts)>0)
      stop(
        "'pt.order' does not contain an ordering value for all patients in the lesion data."
      )

    ordered.pts=pt.order
    n=max(ordered.pts$pts.order)
    lesions=grin.res$lsn.data
    merged.data=merge(ordered.pts,lesions,by="ID",all.y=TRUE)
    grin.res$lsn.data=merged.data

    n.chr=nrow(grin.res$chr.size)
  }

  # set up plotting region

  graphics::plot(c(-0.2,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=F,xlab="",ylab="")

  # background colors for chromosomes

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

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

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

  lsn.types=sort(unique(grin.res$lsn.index$lsn.type))

  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]
  }

  grin.res$lsn.data$lsn.colors=lsn.colors[grin.res$lsn.data$lsn.type]

  grin.res$lsn.data$lsn.size=
    grin.res$lsn.data$x.end-grin.res$lsn.data$x.start

  ord=order(grin.res$lsn.data$lsn.size,decreasing=TRUE)
  grin.res$lsn.data=grin.res$lsn.data[ord,]

  graphics::rect(grin.res$lsn.data$pts.order-1,
                 -grin.res$lsn.data$x.start,
                 grin.res$lsn.data$pts.order,
                 -grin.res$lsn.data$x.end,
                 col=grin.res$lsn.data$lsn.colors,
                 border=grin.res$lsn.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.5)

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

  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 gene-lesion events affecting at least one patient were found to plot.")

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

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

  graphics::segments(1.075*n,
                     -(nsubj.data$x.start+nsubj.data$x.end)/2,
                     1.075*n+0.125*nsubj.data$nsubj/max(nsubj.data$nsubj)*n,
                     col=nsubj.data$lsn.colors)

  # prepare -log10(q-value) data for the left panel

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

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

  if (is.null(max.log10q))
  {
    finite.log10q=q.data$log10q[is.finite(q.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)
  }

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

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

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

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

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

  graphics::text(-(0.075+0.20)*n/2,0,
                 "-log10(q)",cex=0.75,
                 pos=3)

  graphics::text((1.075+1.2)*n/2,0,
                 "Subjects",cex=0.75,
                 pos=3)

  graphics::text(c(-0.075,1.075)*n,
                 -grin.res$chr.size$x.end[n.chr],
                 0,cex=0.75,pos=1)

  graphics::text(c(-0.2,1.2)*n,
                 -grin.res$chr.size$x.end[n.chr],
                 c(max(q.data$log10q),
                   max(nsubj.data$nsubj)),
                 cex=0.75,pos=1)

  invisible(NULL)
}