R/grin.stats.lsn.plot.R

In GRIN2: Genomic Random Interval (GRIN)

#' GRIN Lesion and Statistics Plot
#'
#' @description
#' Generates a regional plot showing genomic lesions overlapping a selected
#' feature together with its GRIN statistics. This function is particularly
#' designed for regulatory and other genomic features that do not have
#' transcript or exon structure to display.
#'
#' @usage
#' grin.stats.lsn.plot(grin.res,
#'                     feature = NULL,
#'                     lsn.clrs = NULL,
#'                     expand = 0.0005)
#'
#' @param grin.res GRIN results for genes, regulatory elements, or other
#' genomic features, as returned by the `grin.stats` function.
#' @param feature Character string specifying the feature of interest. This is
#' typically an Ensembl regulatory feature ID or another regulatory-region
#' identifier used as marker input to `grin.stats`. An Ensembl gene ID can also
#' be provided.
#' @param lsn.clrs A named vector of colors assigned to lesion types. If
#' `NULL`, colors are automatically assigned using `default.grin.colors`.
#' @param expand Numeric; proportion of the feature length added upstream and
#' downstream to define the plotting region. Default is `0.0005`. Setting
#' `expand = 0` restricts the plot to the exact feature boundaries.
#'
#' @details
#' This function provides a regional visualization of lesions and GRIN
#' statistics without displaying transcript or exon structure. It is therefore
#' particularly useful for regulatory features, such as promoters, enhancers,
#' CTCF binding sites, open chromatin regions, or other genomic elements for
#' which a transcript panel is not applicable.
#'
#' Genes can also be plotted with this function when only the lesion
#' distribution and GRIN statistics are needed. However, for gene-centered
#' visualization that includes transcript and exon structure,
#' `lsn.transcripts.plot` is generally more informative.
#'
#' The plot contains a regional lesion panel followed by statistical summaries:
#' \itemize{
#'   \item The upper panel displays individual genomic lesions overlapping the
#'   selected feature. Each row represents one overlapping lesion event, and
#'   lesions are colored according to lesion type. Dashed vertical lines mark
#'   the start and end coordinates of the selected feature.
#'   \item The first statistical summary reports, for each lesion type, the
#'   number of affected subjects together with the corresponding
#'   \eqn{-\log_{10}(p)} and \eqn{-\log_{10}(q)} values from the GRIN
#'   lesion-recurrence analysis.
#'   \item The second statistical summary reports the constellation
#'   \eqn{-\log_{10}(p)} and \eqn{-\log_{10}(q)} values across lesion orders
#'   evaluated by GRIN.
#' }
#'
#' @return
#' Generates a regional lesion and GRIN statistics plot on the active graphics
#' device and invisibly returns `NULL`. No transcript or exon annotation panel
#' is included.
#'
#' @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{lsn.transcripts.plot}},
#' \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)
#'
#' # Plot lesions and GRIN statistics for a gene without displaying
#' # transcript or exon structure
#' grin.stats.lsn.plot(grin.results,
#'                     feature = "ENSG00000184937")
#'
#' # This function is particularly useful when regulatory elements or other
#' # non-transcript genomic features are used as marker input to grin.stats.
#' # For gene-centered plots that include transcript and exon annotation,
#' # use lsn.transcripts.plot().
#'
grin.stats.lsn.plot=function(grin.res,          # GRIN results for regulatory features (output of the grin.stats function)
                             feature=NULL,      # Feature Ensembl ID from Ensembl regulatory build or other regulatory annotation. An Ensembl ID of a gene can be provided as well.
                             lsn.clrs=NULL,     # Specified colors per lesion types. If not specified, colors will be automatically assigned using default.grin.colors function
                             expand=0.0005)     # Controls the proportion of the feature length added upstream and downstream to the plotting region
{

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

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

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

  if (is.null(feature) ||
      !is.character(feature) ||
      length(feature)!=1 ||
      is.na(feature) ||
      feature=="")
    stop("Exactly one valid feature must be specified.")

  if (!is.numeric(expand) ||
      length(expand)!=1 ||
      is.na(expand) ||
      !is.finite(expand) ||
      expand<0)
    stop("'expand' must be a single non-negative numeric value.")

  # Find the requested regulatory feature

  gene.data=grin.res[["gene.data"]]
  gene.mtch=which(gene.data[,"gene"]==feature)

  if (length(gene.mtch)==0)
    stop(paste0(feature," not found in gene.data."))

  if (length(gene.mtch)>1)
    stop(paste0("Multiple matches of ",feature," found in gene.data."))

  gene.chr=gene.data[gene.mtch,"chrom"]
  gene.start=gene.data[gene.mtch,"loc.start"]
  gene.end=gene.data[gene.mtch,"loc.end"]
  gene.size=(gene.end-gene.start)+1

  # Find lesions on the same chromosome as the feature

  lsn.dset=order.index.lsn.data(grin.res[["lsn.data"]])
  lsn.data=lsn.dset$lsn.data

  lsn.types=unique(as.character(lsn.data$lsn.type))

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

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

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

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

    lsn.clrs=lsn.clrs[lsn.types]
  }

  lsn.index=lsn.dset$lsn.index
  lsn.ind.mtch=which(lsn.index$chrom==gene.chr)

  if (length(lsn.ind.mtch)==0)
    stop(paste0("No lesions were found on chromosome ",gene.chr,"."))

  lsn.chr.rows=NULL

  for (i in lsn.ind.mtch)
  {
    blk.rows=(lsn.index$row.start[i]:lsn.index$row.end[i])
    lsn.chr.rows=c(lsn.chr.rows,blk.rows)
  }

  lsn.chr.rows=unlist(lsn.chr.rows)

  lsn.chr.data=lsn.data[lsn.chr.rows,]

  if (any(lsn.chr.data$chrom!=gene.chr))
    stop(paste0("Error in finding lesions on the same chromosome as feature ",feature,"."))

  # Find lesions that overlap the regulatory feature

  ov.rows=which(
    (lsn.chr.data$loc.start<=gene.end) &
      (lsn.chr.data$loc.end>=gene.start)
  )

  if (length(ov.rows)==0)
    stop(paste0("No lesions overlap feature ",feature,"."))

  lsn.gene=lsn.chr.data[ov.rows,]

  lsn.gene$size=lsn.gene$loc.end-
    lsn.gene$loc.start+1

  lsn.gene=lsn.gene[
    order(lsn.gene$lsn.type,lsn.gene$size),
  ]

  # define plotting data

  x.start=gene.start-expand*gene.size
  x.end=gene.end+expand*gene.size

  lsn.gene$index=1:nrow(lsn.gene)
  lsn.gene$subj.num=as.numeric(as.factor(lsn.gene$index))
  lsn.gene$lsn.clr=lsn.clrs[lsn.gene$lsn.type]
  lsn.gene$type.num=as.numeric(as.factor(lsn.gene$lsn.type))

  n.type=max(lsn.gene$type.num)
  n.subj=max(lsn.gene$subj.num)

  lsn.gene$y0=-lsn.gene$subj.num+(lsn.gene$type.num-1)/n.type
  lsn.gene$y1=-lsn.gene$subj.num+lsn.gene$type.num/n.type

  lsn.gene$x0=pmax(lsn.gene$loc.start,x.start)
  lsn.gene$x1=pmin(lsn.gene$loc.end,x.end)

  graphics::plot(
    c(x.start-0.055*(x.end-x.start),
      x.end+0.13*(x.end-x.start)),
    c(+0.05,-1.75)*n.subj,
    type="n",
    main="",
    xlab="",
    ylab="",
    axes=FALSE
  )

  graphics::rect(
    x.start,
    -(1:n.subj),
    x.end,
    -(1:n.subj)+1,
    col=c("snow","gainsboro")[1+(1:n.subj)%%2],
    border=NA
  )

  graphics::segments(
    c(gene.start,gene.end),
    rep(-n.subj,2),
    c(gene.start,gene.end),
    rep(0,2),
    col="darkgray"
  )

  graphics::rect(
    lsn.gene$x0,
    lsn.gene$y0,
    lsn.gene$x1,
    lsn.gene$y1,
    col=lsn.gene$lsn.clr,
    border=lsn.gene$lsn.clr
  )

  graphics::segments(
    c(gene.start,gene.end),
    rep(-n.subj,2),
    c(gene.start,gene.end),
    rep(0,2),
    col="darkgray",
    lty=2
  )

  graphics::text(
    c(gene.start,gene.end),
    -n.subj,
    c(gene.start,gene.end),
    pos=1,
    cex=0.75
  )

  graphics::text(
    (gene.start+gene.end)/2,
    0,
    feature,
    pos=3,
    cex=1
  )

  lgd=graphics::legend(
    (x.start+x.end)/2,
    -1.10*n.subj,
    fill=lsn.clrs,
    legend=names(lsn.clrs),
    ncol=length(lsn.clrs),
    xjust=0.4,
    border=NA,
    cex=0.7,
    bty="n"
  )

  graphics::text(
    lgd$text$x[1]-0.085*diff(range(lgd$text$x)),
    -c(1.25,1.35,1.45)*n.subj,
    c("n","-log10p","-log10q"),
    pos=2,
    cex=0.7,
    font=2
  )

  gene.stats=grin.res[["gene.hits"]]
  stat.mtch=which(gene.stats$gene==feature)

  if (length(stat.mtch)==0)
    stop(paste0(feature," not found in gene.hits."))

  if (length(stat.mtch)>1)
    stop(paste0("Multiple matches of ",feature," found in gene.hits."))

  gene.stats=gene.stats[stat.mtch,]

  graphics::text(
    lgd$text$x,
    -1.25*n.subj,
    gene.stats[,paste0("nsubj.",names(lsn.clrs))],
    cex=0.7
  )

  graphics::text(
    lgd$text$x,
    -1.35*n.subj,
    round(-log10(
      gene.stats[,paste0("p.nsubj.",names(lsn.clrs))]
    ),2),
    cex=0.7
  )

  graphics::text(
    lgd$text$x,
    -1.45*n.subj,
    round(-log10(
      gene.stats[,paste0("q.nsubj.",names(lsn.clrs))]
    ),2),
    cex=0.7
  )

  lgd2=graphics::legend(
    (x.start+x.end)/2,
    -1.5*n.subj,
    legend=paste0("const",1:length(lsn.types),".typ"),
    ncol=length(lsn.types),
    xjust=0.42,
    border=NA,
    cex=0.7,
    bty="n"
  )

  graphics::text(
    lgd2$text$x[1]-0.0005*diff(range(lgd2$text$x)),
    -c(1.65,1.75)*n.subj,
    c("const.p","const.q"),
    pos=2,
    cex=0.7,
    font=2
  )

  graphics::text(
    lgd2$text$x,
    -1.65*n.subj,
    round(-log10(
      gene.stats[,paste0("p",1:length(lsn.types),".nsubj")]
    ),2),
    cex=0.7,
    pos=4
  )

  graphics::text(
    lgd2$text$x,
    -1.75*n.subj,
    round(-log10(
      gene.stats[,paste0("q",1:length(lsn.types),".nsubj")]
    ),2),
    cex=0.7,
    pos=4
  )

  invisible(NULL)
}