R/lsn.transcripts.plot.R

In GRIN2: Genomic Random Interval (GRIN)

#' Plot Genomic Lesions with Transcript Annotations
#'
#' @description
#' Creates an integrated genomic-region plot combining patient-level genomic
#' lesions with optional transcript annotations, chromosome context, genomic
#' coordinates, lesion-type legends, and gene-level GRIN2 statistics.
#'
#' The function supports both gene-centered plots and user-defined genomic
#' regions. Single-position lesions, or lesion types explicitly selected by the
#' user, are drawn at their exact genomic coordinates with user-controlled
#' marker size.
#'
#' @usage
#' lsn.transcripts.plot(
#'   grin.res,
#'   gene = NULL,
#'   chrom = NULL,
#'   plot.start = NULL,
#'   plot.end = NULL,
#'   transTrack = TRUE,
#'   show.ideogram = TRUE,
#'   show.axis = TRUE,
#'   show.transcript.labels = FALSE,
#'   transcript.mode = c("all", "selected"),
#'   selected.transcript = NULL,
#'   lesion.grp = NULL,
#'   lsn.clrs = NULL,
#'   spec.lsn.clr = NULL,
#'   point.lesion.groups = NULL,
#'   point.size.mm = 0.8,
#'   expand = 0.0005,
#'   hg38.transcripts = NULL,
#'   hg38.cytoband = NULL,
#'   genome = "hg38",
#'   transcript.height = 0.31,
#'   draw = TRUE
#' )
#'
#' @param grin.res A GRIN2 results object containing at least `lsn.data`.
#'   Gene-centered plots additionally require `gene.data`; gene-level recurrence
#'   and constellation statistics are displayed when `gene.hits` is available.
#'
#' @param gene Optional character scalar specifying the gene to plot. When
#'   supplied, genomic coordinates are obtained from `grin.res[["gene.data"]]`.
#'
#' @param chrom Chromosome for a user-defined regional plot. Used only when
#'   `gene = NULL`. Values with or without the `chr` prefix are accepted.
#'
#' @param plot.start,plot.end Numeric start and end coordinates for a
#'   user-defined regional plot. Both are required when `gene = NULL`.
#'
#' @param transTrack Logical; whether to display transcript annotations.
#'
#' @param show.ideogram Logical; whether to display a chromosome ideogram.
#'
#' @param show.axis Logical; whether to display a genomic-coordinate axis.
#'
#' @param show.transcript.labels Logical; whether transcript or gene labels are
#'   displayed within the transcript track.
#'
#' @param transcript.mode Character string specifying which transcript
#'   annotations to display. Use `"all"` for all overlapping transcripts or
#'   `"selected"` for transcript IDs supplied through `selected.transcript`.
#'
#' @param selected.transcript Character vector of transcript IDs retained when
#'   `transcript.mode = "selected"`. Version suffixes are ignored during
#'   matching, allowing identifiers such as `ENST00000342505` and
#'   `ENST00000342505.8` to match.
#'
#' @param lesion.grp Optional character vector specifying lesion types to plot.
#'
#' @param lsn.clrs Optional named character vector assigning colors to lesion
#'   types. Vector names must match the lesion-type values in
#'   `grin.res[["lsn.data"]]$lsn.type`. If an unnamed vector is supplied, it
#'   must contain exactly one color for each plotted lesion type. When
#'   `lsn.clrs = NULL`, colors are assigned automatically.
#'
#' @param spec.lsn.clr Optional color used when a single lesion type is
#'   selected through `lesion.grp`.
#'
#' @param point.lesion.groups Optional character vector specifying lesion types
#'   to display as point markers. When `NULL`, each lesion is classified from
#'   its genomic coordinates: lesions with identical start and end coordinates
#'   are displayed as circles, whereas lesions spanning multiple positions are
#'   displayed as intervals.
#'
#' @param point.size.mm Positive numeric scalar specifying the diameter, in
#'   millimeters, of lesions displayed as circles. The default is 0.8.
#'
#' @param expand Non-negative numeric proportion used to expand a gene-centered
#'   plotting interval upstream and downstream. It is ignored for user-defined
#'   regional plots, which use the supplied coordinates exactly.
#'
#' @param hg38.transcripts An `EnsDb` transcript-annotation object compatible
#'   with [ensembldb::getGeneRegionTrackForGviz()]. Required when
#'   `transTrack = TRUE`. A compatible GRCh38 transcript annotation can be
#'   obtained through AnnotationHub using the workflow shown in the examples.
#'
#' @param hg38.cytoband Cytoband annotation supplied to
#'   [Gviz::IdeogramTrack()]. Required when `show.ideogram = TRUE`. A GRCh38
#'   cytoband annotation is provided with GRIN2 as `hg38_cytoband` and can be
#'   loaded using `data(hg38_cytoband)`.
#'
#' @param genome Genome identifier passed to the Gviz annotation tracks.
#'
#' @param transcript.height Relative vertical space allocated to the transcript
#'   annotation track. Increasing this value provides additional space when
#'   many transcripts are displayed. This value is used only when
#'   `transTrack = TRUE`.
#'
#' @param draw Logical; if `TRUE`, draw the assembled plot on a new graphics
#'   page. The plot object is returned invisibly in either case.
#'
#' @details
#' In gene mode, the displayed interval is derived from the matching entry in
#' `gene.data` and can be expanded using `expand`. In regional mode, the exact
#' interval specified by `chrom`, `plot.start`, and `plot.end` is used.
#'
#' Transcript annotations can be shown in full or restricted to selected
#' transcript IDs. Selected IDs are matched after removal of version suffixes.
#' Unmatched IDs are reported in a warning when at least one requested
#' transcript is available; the function stops with an informative error when
#' none of the requested transcripts can be found.
#'
#' Point-like lesions are plotted at their exact genomic start coordinates and
#' use physical marker sizes measured in millimeters. Lesions are ordered first
#' by lesion type and then by interval size, with source order used to break
#' exact ties.
#'
#' When transcript annotations are not displayed, the vertical space allocated
#' to the ideogram and genomic-coordinate axis is reduced automatically so that
#' the lesion panel occupies the available plotting area more efficiently.
#'
#' For gene-centered plots, available lesion-specific subject counts and
#' transformed p- and q-values are displayed together with ordered
#' constellation statistics from `gene.hits`. Missing statistics are omitted
#' rather than preventing plot construction.
#'
#' @return Invisibly returns a `gtable`/`grob` object representing the assembled
#'   figure. When `draw = TRUE`, the figure is also drawn.
#'
#' @examples
#' data("grin.results")
#' data("hg38_cytoband")
#'
#' \donttest{
#' # Plot lesions across a user-defined genomic region without transcript
#' # annotations. This example uses only data included with GRIN2.
#' lsn.transcripts.plot(
#'   grin.res = grin.results,
#'   chrom = 9,
#'   plot.start = 19900000,
#'   plot.end = 25600000,
#'   transTrack = FALSE,
#'   show.ideogram = FALSE,
#'   lesion.grp = "loss",
#'   spec.lsn.clr = "blue"
#' )
#'
#' # Plot all lesion types across a larger genomic region without transcript
#' # annotations.
#' lsn.transcripts.plot(
#'   grin.res = grin.results,
#'   chrom = 9,
#'   plot.start = 1,
#'   plot.end = 141000000,
#'   transTrack = FALSE,
#'   show.ideogram = FALSE,
#'   point.size.mm = 1
#' )
#' }
#'
#' \dontrun{
#' # Transcript-track examples require an EnsDb annotation object retrieved
#' # through AnnotationHub and therefore require an internet connection.
#' ah <- AnnotationHub::AnnotationHub()
#'
#' # Ensembl release 110 transcript annotations for human GRCh38.
#' hg38.transcripts <- ah[["AH113665"]]
#'
#' # Define colors for the lesion types in the example dataset.
#' lesion.colors <- c(
#'   fusion = "black",
#'   gain = "red",
#'   loss = "blue",
#'   mutation = "olivedrab"
#' )
#'
#' # Plot all transcripts and genomic lesions affecting JAK3.
#' lsn.transcripts.plot(
#'   grin.res = grin.results,
#'   gene = "JAK3",
#'   transcript.mode = "all",
#'   lsn.clrs = lesion.colors,
#'   hg38.transcripts = hg38.transcripts,
#'   hg38.cytoband = hg38_cytoband
#' )
#'
#' # Display selected JAK1 transcripts and increase point-marker size.
#' lsn.transcripts.plot(
#'   grin.res = grin.results,
#'   gene = "JAK1",
#'   transcript.mode = "selected",
#'   selected.transcript = c(
#'     "ENST00000342505",
#'     "ENST00000672751"
#'   ),
#'   lsn.clrs = lesion.colors,
#'   point.size.mm = 0.7,
#'   transcript.height = 0.40,
#'   hg38.transcripts = hg38.transcripts,
#'   hg38.cytoband = hg38_cytoband
#' )
#'
#' # Plot transcript annotations and genomic lesions across a user-defined
#' # region while displaying gene names in the transcript track.
#' lsn.transcripts.plot(
#'   grin.res = grin.results,
#'   chrom = 9,
#'   plot.start = 21900000,
#'   plot.end = 21990000,
#'   transcript.mode = "all",
#'   show.transcript.labels = TRUE,
#'   lsn.clrs = lesion.colors,
#'   transcript.height = 0.46,
#'   hg38.transcripts = hg38.transcripts,
#'   hg38.cytoband = hg38_cytoband
#' )
#' }
#'
#' @seealso [grin.stats()], [Gviz::GeneRegionTrack()],
#'   [Gviz::GenomeAxisTrack()], [Gviz::IdeogramTrack()],
#'   [ensembldb::getGeneRegionTrackForGviz()]
#'
#' @export
#'
lsn.transcripts.plot=function(
    grin.res,

    # Gene-level or region-level plot
    gene=NULL,
    chrom=NULL,
    plot.start=NULL,
    plot.end=NULL,

    # Track options
    transTrack=TRUE,
    show.ideogram=TRUE,
    show.axis=TRUE,
    show.transcript.labels=FALSE,

    # Transcript-selection options
    transcript.mode=c("all","selected"),
    selected.transcript=NULL,

    # Lesion-selection and color options
    lesion.grp=NULL,
    lsn.clrs=NULL,
    spec.lsn.clr=NULL,

    # Point-like lesion options
    point.lesion.groups=NULL,
    point.size.mm=0.8,

    # Gene plotting-range expansion
    expand=0.0005,

    # Annotation objects
    hg38.transcripts=NULL,
    hg38.cytoband=NULL,
    genome="hg38",

    # Transcript-track height
    transcript.height=0.31,

    # Draw or return
    draw=TRUE)
{
  ####################################################################
  # 1. Required package checks
  ####################################################################

  required.packages=c("Gviz","grid","gridExtra","GenomeInfoDb","ensembldb",
                      "GenomicRanges","IRanges","S4Vectors")

  missing.packages=required.packages[
    !vapply(required.packages,requireNamespace,logical(1),quietly=TRUE)
  ]

  if (length(missing.packages)>0L)
    stop(paste0("The following required package(s) are not installed: ",
                paste(missing.packages,collapse=", "),
                ".\nInstall them using BiocManager::install()."),
         call.=FALSE)

  ####################################################################
  # 2. Internal graphical settings
  ####################################################################

  transcript.label="Transcripts"
  lesion.label="Genomic lesions"
  label.width=0.48

  row.background=c("white","grey92")
  interval.border=NA
  interval.height=0.82

  coordinate.cex=0.78
  title.cex=1.00
  legend.cex=0.85
  stats.cex=0.80

  # Align the lesion panel with the Gviz annotation panel
  lesion.track.inset.mm=3

  # Relative heights of fixed plot components
  title.height=0.045
  lesion.height=0.44
  legend.height=0.065
  lesion.stats.height=0.105
  constellation.stats.height=0.105

  ####################################################################
  # 3. General helper functions
  ####################################################################

  normalize.chromosome=function(x)
    sub("^chr","",as.character(x),ignore.case=TRUE)

  to.ucsc.chromosome=function(x)
    paste0("chr",normalize.chromosome(x))

  normalize.transcript.id=function(x)
    sub("\\.[0-9]+$","",as.character(x))

  safe.neglog10=function(x)
  {
    x=suppressWarnings(as.numeric(x))
    output=rep(NA_real_,length(x))
    positive=!is.na(x) & x>0
    zero=!is.na(x) & x==0

    output[positive]=-log10(x[positive])
    output[zero]=Inf
    output
  }

  format.statistic=function(x,digits=2L)
  {
    if (length(x)==0L || is.na(x))
      return("")

    if (is.infinite(x))
      return("Inf")

    format(round(x,digits),trim=TRUE,scientific=FALSE,nsmall=0)
  }

  detect.column=function(data,candidates,required=TRUE,description="column")
  {
    available=candidates[candidates %in% colnames(data)]

    if (length(available)>0L)
      return(available[1L])

    if (isTRUE(required))
      stop(paste0("Could not identify the ",description,
                  ". Expected one of: ",paste(candidates,collapse=", "),"."),call.=FALSE)

    NULL
  }

  detect.metadata.column=function(object,candidates,required=TRUE,
                                  description="metadata column")
  {
    metadata.names=colnames(S4Vectors::mcols(object))
    available=candidates[candidates %in% metadata.names]

    if (length(available)>0L)
      return(available[1L])

    if (isTRUE(required))
      stop(paste0("Could not identify the ",description,
                  " in the retrieved transcript annotation. Expected one of: ",
                  paste(candidates,collapse=", "),"."),call.=FALSE)

    NULL
  }

  obtain.lesion.colors=function(lesion.types,supplied.colors=NULL)
  {
    lesion.types=unique(as.character(lesion.types))

    if (length(lesion.types)==0L)
      return(character(0))

    # Use user-specified lesion colors when provided; otherwise use GRIN defaults
    if (is.null(supplied.colors))
      supplied.colors=default.grin.colors(lesion.types)

    # Assign lesion-type names when an unnamed color vector is provided
    if (is.null(names(supplied.colors)))
    {
      if (length(supplied.colors)!=length(lesion.types))
        stop("`lsn.clrs` must be a named color vector or contain exactly one color per plotted lesion type.",
             call.=FALSE)

      names(supplied.colors)=lesion.types
    }

    # Confirm that colors are available for all plotted lesion types
    missing.types=setdiff(lesion.types,names(supplied.colors))

    if (length(missing.types)>0L)
      stop(paste0("No color was provided for: ",
                  paste(missing.types,collapse=", "),"."),call.=FALSE)

    supplied.colors[lesion.types]
  }

  make.left.label.grob=function(label,fontsize=10,fontface="bold")
  {
    if (is.null(label) || !nzchar(label))
      return(grid::nullGrob())

    grid::textGrob(label=label,rot=90,
                   gp=grid::gpar(fontsize=fontsize,fontface=fontface))
  }

  add.fixed.label.column=function(plot.grob,label,width.inches)
  {
    gridExtra::arrangeGrob(
      make.left.label.grob(label=label),
      plot.grob,
      ncol=2,
      widths=grid::unit.c(grid::unit(width.inches,"in"),
                          grid::unit(1,"null")))
  }

  add.horizontal.inset=function(plot.grob,inset.mm=0)
  {
    if (length(inset.mm)!=1L || is.na(inset.mm) || inset.mm<0)
      stop("`inset.mm` must be one non-negative numeric value.",call.=FALSE)

    if (inset.mm==0)
      return(plot.grob)

    gridExtra::arrangeGrob(
      grid::nullGrob(),plot.grob,grid::nullGrob(),
      ncol=3,
      widths=grid::unit.c(grid::unit(inset.mm,"mm"),
                          grid::unit(1,"null"),
                          grid::unit(inset.mm,"mm")))
  }

  ####################################################################
  # 4. Transcript-selection helper
  ####################################################################

  filter.selected.transcript=function(transcript.object,selected.transcript)
  {
    if (is.null(selected.transcript) || length(selected.transcript)==0L)
      stop(paste0("`selected.transcript` must contain at least one non-empty ",
                  "transcript ID when `transcript.mode = \"selected\"`."),
           call.=FALSE)

    selected.transcript=trimws(as.character(selected.transcript))
    invalid.selected=is.na(selected.transcript) | !nzchar(selected.transcript)

    if (any(invalid.selected))
      stop(paste0("`selected.transcript` must contain only non-missing, ",
                  "non-empty transcript IDs."),call.=FALSE)

    transcript.metadata.column=detect.metadata.column(
      transcript.object,
      c("transcript","transcript_id","transcriptId","tx_id","txid","group"),
      required=TRUE,
      description="transcript-ID metadata column")

    requested.ids=unique(normalize.transcript.id(selected.transcript))
    available.ids=normalize.transcript.id(
      S4Vectors::mcols(transcript.object)[[transcript.metadata.column]]
    )

    matched.ids=intersect(requested.ids,unique(available.ids[!is.na(available.ids)]))
    missing.ids=setdiff(requested.ids,matched.ids)

    if (length(matched.ids)==0L)
    {
      available.unique=sort(unique(available.ids[!is.na(available.ids)]))

      stop(paste0(
        "None of the requested transcripts were found in the retrieved annotation: ",
        paste(requested.ids,collapse=", "),". ",
        if (length(available.unique)>0L)
        {
          paste0("Available transcript IDs include: ",
                 paste(utils::head(available.unique,10L),collapse=", "),
                 if (length(available.unique)>10L) ", ..." else "",
                 ".")
        } else {
          "No transcript IDs were available in the retrieved annotation."
        }),call.=FALSE)
    }

    if (length(missing.ids)>0L)
      warning(paste0("The following requested transcript ID(s) were not found ",
                     "and will be omitted: ",paste(missing.ids,collapse=", "),"."),
              call.=FALSE)

    keep=!is.na(available.ids) & available.ids %in% matched.ids
    filtered.object=transcript.object[keep]

    if (length(filtered.object)==0L)
      stop("No transcript annotations remained after transcript filtering.",
           call.=FALSE)

    filtered.object
  }

  ####################################################################
  # 5. Lesion-track constructor
  ####################################################################

  make.lesion.track.grob=function(lesion.data,from,to,colors,point.groups,
                                  point.size,interval.outline,
                                  interval.rectangle.height,backgrounds,
                                  coordinate.font.size)
  {
    number.lesions=nrow(lesion.data)

    if (number.lesions==0L)
      stop("No lesions are available for plotting.",call.=FALSE)

    lesion.data$plot.start=pmax(as.numeric(lesion.data$loc.start),from)
    lesion.data$plot.end=pmin(as.numeric(lesion.data$loc.end),to)
    children=list()

    # Alternating patient-row backgrounds
    for (i in seq_len(number.lesions))
    {
      y.center=number.lesions-i+0.5
      background.color=backgrounds[1L+((i-1L)%%length(backgrounds))]

      children[[length(children)+1L]]=grid::rectGrob(
        x=grid::unit((from+to)/2,"native"),
        y=grid::unit(y.center,"native"),
        width=grid::unit(to-from,"native"),
        height=grid::unit(1,"native"),
        gp=grid::gpar(fill=background.color,col=NA))
    }

    # Draw lesion intervals and point-like lesions
    for (i in seq_len(number.lesions))
    {
      lesion.type=as.character(lesion.data$lsn.type[i])
      lesion.color=unname(colors[lesion.type])
      y.center=number.lesions-i+0.5

      # When no lesion groups are specified, classify each lesion independently
      # according to its genomic span. Users may override this behavior by
      # supplying lesion-type names through point.lesion.groups.
      is.point.lesion=if (is.null(point.groups))
        lesion.data$loc.start[i]==lesion.data$loc.end[i]
      else
        lesion.type %in% point.groups

      if (isTRUE(is.point.lesion))
      {
        # Preserve the exact genomic coordinate for point-like lesions
        point.coordinate=as.numeric(lesion.data$loc.start[i])

        if (is.finite(point.coordinate) &&
            point.coordinate>=from &&
            point.coordinate<=to)
        {
          children[[length(children)+1L]]=grid::circleGrob(
            x=grid::unit(point.coordinate,"native"),
            y=grid::unit(y.center,"native"),
            r=grid::unit(point.size/2,"mm"),
            gp=grid::gpar(fill=lesion.color,col=lesion.color,lwd=0))
        }
      } else {
        interval.start=lesion.data$plot.start[i]
        interval.end=lesion.data$plot.end[i]

        if (is.finite(interval.start) &&
            is.finite(interval.end) &&
            interval.start<=to &&
            interval.end>=from)
        {
          interval.start=max(interval.start,from)
          interval.end=min(interval.end,to)
          interval.width=max(interval.end-interval.start,.Machine$double.eps)

          children[[length(children)+1L]]=grid::rectGrob(
            x=grid::unit((interval.start+interval.end)/2,"native"),
            y=grid::unit(y.center,"native"),
            width=grid::unit(interval.width,"native"),
            height=grid::unit(interval.rectangle.height,"native"),
            gp=grid::gpar(fill=lesion.color,col=interval.outline))
        }
      }
    }

    # Dashed boundaries at the selected range
    children[[length(children)+1L]]=grid::segmentsGrob(
      x0=grid::unit(c(from,to),"native"),
      x1=grid::unit(c(from,to),"native"),
      y0=grid::unit(c(0,0),"native"),
      y1=grid::unit(c(number.lesions,number.lesions),"native"),
      gp=grid::gpar(col="grey55",lty=2,lwd=0.7))

    # Start coordinate
    children[[length(children)+1L]]=grid::textGrob(
      label=format(from,scientific=FALSE,trim=TRUE),
      x=grid::unit(from,"native"),
      y=grid::unit(-0.18,"native"),
      just=c("left","top"),
      gp=grid::gpar(fontsize=coordinate.font.size))

    # End coordinate
    children[[length(children)+1L]]=grid::textGrob(
      label=format(to,scientific=FALSE,trim=TRUE),
      x=grid::unit(to,"native"),
      y=grid::unit(-0.18,"native"),
      just=c("right","top"),
      gp=grid::gpar(fontsize=coordinate.font.size))

    grid::grobTree(
      children=do.call(grid::gList,children),
      vp=grid::viewport(xscale=c(from,to),
                        yscale=c(-1.10,number.lesions+0.15),
                        clip="off"))
  }

  ####################################################################
  # 6. Legend constructor
  ####################################################################

  make.legend.grob=function(colors,fontsize)
  {
    lesion.types=names(colors)

    if (length(lesion.types)==0L)
      return(grid::nullGrob())

    legend.items=vector("list",length(lesion.types))
    legend.item.widths=vector("list",length(lesion.types))

    for (i in seq_along(lesion.types))
    {
      lesion.type=lesion.types[i]
      lesion.color=unname(colors[i])

      symbol.grob=grid::rectGrob(
        width=grid::unit(0.14,"in"),
        height=grid::unit(0.07,"in"),
        gp=grid::gpar(fill=lesion.color,col=NA))

      label.grob=grid::textGrob(
        label=lesion.type,x=0,just="left",
        gp=grid::gpar(fontsize=fontsize))

      label.width=grid::grobWidth(label.grob)

      legend.items[[i]]=gridExtra::arrangeGrob(
        symbol.grob,label.grob,ncol=2,
        widths=grid::unit.c(grid::unit(0.22,"in"),label.width))

      legend.item.widths[[i]]=grid::unit(0.22,"in")+
        label.width+grid::unit(0.12,"in")
    }

    legend.content=gridExtra::arrangeGrob(
      grobs=legend.items,nrow=1,
      widths=do.call(grid::unit.c,legend.item.widths))

    # Center the complete legend as one compact unit
    gridExtra::arrangeGrob(
      grid::nullGrob(),legend.content,grid::nullGrob(),
      ncol=3,
      widths=grid::unit.c(grid::unit(1,"null"),
                          grid::grobWidth(legend.content),
                          grid::unit(1,"null")))
  }

  ####################################################################
  # 7. Statistics-table constructors
  ####################################################################

  make.table.theme=function(base.size)
  {
    gridExtra::ttheme_minimal(
      base_size=base.size,
      padding=grid::unit(c(1.5,2.0),"mm"),
      core=list(fg_params=list(hjust=0.5,x=0.5)),
      colhead=list(fg_params=list(fontface="bold",hjust=0.5,x=0.5)))
  }

  make.gene.statistics.grobs=function(grin.object,requested.gene,
                                      lesion.types,base.size)
  {
    empty.output=list(lesion=grid::nullGrob(),
                      constellation=grid::nullGrob(),
                      lesion.available=FALSE,
                      constellation.available=FALSE)

    if (is.null(requested.gene) || is.null(grin.object[["gene.hits"]]))
      return(empty.output)

    gene.statistics=grin.object[["gene.hits"]]

    if (!"gene.name" %in% colnames(gene.statistics))
      return(empty.output)

    gene.statistics=gene.statistics[
      as.character(gene.statistics$gene.name)==requested.gene,,drop=FALSE
    ]

    if (nrow(gene.statistics)!=1L)
      return(empty.output)

    obtain.value=function(column.name)
    {
      if (!column.name %in% colnames(gene.statistics))
        return(NA_real_)

      suppressWarnings(as.numeric(gene.statistics[1,column.name,drop=TRUE]))
    }

    lesion.table=data.frame(
      Statistic=c("n","-log10(p)","-log10(q)"),
      stringsAsFactors=FALSE,
      check.names=FALSE)

    for (lesion.type in lesion.types)
    {
      n.value=obtain.value(paste0("nsubj.",lesion.type))
      p.value=obtain.value(paste0("p.nsubj.",lesion.type))
      q.value=obtain.value(paste0("q.nsubj.",lesion.type))

      lesion.table[[lesion.type]]=c(
        ifelse(is.na(n.value),"",
               format(n.value,trim=TRUE,scientific=FALSE)),
        format.statistic(safe.neglog10(p.value)),
        format.statistic(safe.neglog10(q.value)))
    }

    lesion.grob=gridExtra::tableGrob(
      lesion.table,rows=NULL,theme=make.table.theme(base.size))

    lesion.grob$widths=grid::unit.c(
      grid::unit(1.35,"in"),
      rep(grid::unit(0.80,"in"),ncol(lesion.table)-1L))

    p.columns=grep("^p[0-9]+\\.nsubj$",colnames(gene.statistics),value=TRUE)
    q.columns=grep("^q[0-9]+\\.nsubj$",colnames(gene.statistics),value=TRUE)

    extract.order=function(column.names)
      suppressWarnings(as.integer(
        sub("^[pq]([0-9]+)\\.nsubj$","\\1",column.names)
      ))

    p.orders=extract.order(p.columns)
    q.orders=extract.order(q.columns)
    common.orders=sort(intersect(p.orders,q.orders))

    if (length(common.orders)==0L)
      return(list(lesion=lesion.grob,
                  constellation=grid::nullGrob(),
                  lesion.available=TRUE,
                  constellation.available=FALSE))

    constellation.table=data.frame(
      Statistic=c("const.p","const.q"),
      stringsAsFactors=FALSE,
      check.names=FALSE)

    for (current.order in common.orders)
    {
      p.column=paste0("p",current.order,".nsubj")
      q.column=paste0("q",current.order,".nsubj")

      constellation.table[[paste0("const",current.order,".typ")]]=c(
        format.statistic(safe.neglog10(obtain.value(p.column))),
        format.statistic(safe.neglog10(obtain.value(q.column))))
    }

    constellation.grob=gridExtra::tableGrob(
      constellation.table,rows=NULL,theme=make.table.theme(base.size))

    constellation.grob$widths=grid::unit.c(
      grid::unit(1.35,"in"),
      rep(grid::unit(0.80,"in"),ncol(constellation.table)-1L))

    list(lesion=lesion.grob,
         constellation=constellation.grob,
         lesion.available=TRUE,
         constellation.available=TRUE)
  }

  ####################################################################
  # 8. Gviz-track constructor
  ####################################################################

  make.gviz.track.grob=function(chromosome,from,to,transcript.object,
                                cytoband.object,include.transcripts,
                                include.ideogram,include.axis,genome.name)
  {
    tracks=list()
    track.sizes=numeric(0)

    if (isTRUE(include.ideogram))
    {
      if (is.null(cytoband.object))
        stop(paste0("`hg38.cytoband` must be supplied when ",
                    "`show.ideogram = TRUE`."),call.=FALSE)

      ideogram.track=Gviz::IdeogramTrack(
        genome=genome.name,
        chromosome=chromosome,
        bands=cytoband.object,
        showId=FALSE)

      tracks=c(tracks,list(ideogram.track))
      track.sizes=c(track.sizes,1)
    }

    if (isTRUE(include.axis))
    {
      genome.axis.track=Gviz::GenomeAxisTrack(cex=0.6)
      tracks=c(tracks,list(genome.axis.track))
      track.sizes=c(track.sizes,1.25)
    }

    if (isTRUE(include.transcripts))
    {
      if (is.null(transcript.object) || length(transcript.object)==0L)
        stop("No transcript annotations are available for the selected region.",
             call.=FALSE)

      transcript.track=Gviz::GeneRegionTrack(
        transcript.object,
        chromosome=chromosome,
        genome=genome.name,
        name="",
        showId=show.transcript.labels,
        geneSymbol=show.transcript.labels,
        just.group="right",
        size=0.5,
        transcriptAnnotation=if (show.transcript.labels) "symbol" else NULL)

      tracks=c(tracks,list(transcript.track))
      track.sizes=c(track.sizes,5.5)
    }

    if (length(tracks)==0L)
      return(grid::nullGrob())

    grid::grid.grabExpr(
      Gviz::plotTracks(trackList=tracks,
                       chromosome=chromosome,
                       from=from,
                       to=to,
                       sizes=track.sizes,
                       showTitle=FALSE,
                       background.title="white",
                       col.title="black",
                       add=FALSE),
      wrap=TRUE)
  }

  ####################################################################
  # 9. Validate arguments
  ####################################################################

  transcript.mode=match.arg(transcript.mode)

  if (transcript.mode=="selected" && !isTRUE(transTrack))
    stop("`transcript.mode = \"selected\"` requires `transTrack = TRUE`.",
         call.=FALSE)

  gene.mode=!is.null(gene)

  if (gene.mode)
  {
    if (length(gene)!=1L)
      stop("Exactly one gene must be specified.",call.=FALSE)
  } else {
    if (is.null(chrom) || is.null(plot.start) || is.null(plot.end))
      stop(paste0("Either `gene` or all of `chrom`, `plot.start`, ",
                  "and `plot.end` must be supplied."),call.=FALSE)
  }

  if (!is.null(point.lesion.groups))
  {
    point.lesion.groups=as.character(point.lesion.groups)
    invalid.point.groups=is.na(point.lesion.groups) |
      !nzchar(trimws(point.lesion.groups))

    if (any(invalid.point.groups))
      stop(paste0("`point.lesion.groups` must be `NULL` or contain only ",
                  "non-missing, non-empty lesion-type names."),call.=FALSE)

    point.lesion.groups=unique(trimws(point.lesion.groups))
  }

  if (!is.numeric(point.size.mm) ||
      length(point.size.mm)!=1L ||
      is.na(point.size.mm) ||
      point.size.mm<=0)
    stop("`point.size.mm` must be one positive numeric value.",call.=FALSE)

  if (!is.numeric(expand) ||
      length(expand)!=1L ||
      is.na(expand) ||
      expand<0)
    stop("`expand` must be one non-negative numeric value.",call.=FALSE)

  if (!is.numeric(transcript.height) ||
      length(transcript.height)!=1L ||
      is.na(transcript.height) ||
      transcript.height<=0)
    stop("`transcript.height` must be one positive numeric value.",call.=FALSE)

  ####################################################################
  # 10. Determine plotting coordinates
  ####################################################################

  gene.data=grin.res[["gene.data"]]

  if (gene.mode)
  {
    if (is.null(gene.data))
      stop("`grin.res[['gene.data']]` was not found.",call.=FALSE)

    required.gene.columns=c("gene.name","chrom","loc.start","loc.end")
    missing.gene.columns=setdiff(required.gene.columns,colnames(gene.data))

    if (length(missing.gene.columns)>0L)
      stop(paste0("`gene.data` is missing: ",
                  paste(missing.gene.columns,collapse=", "),"."),
           call.=FALSE)

    gene.match=which(as.character(gene.data$gene.name)==gene)

    if (length(gene.match)==0L)
      stop(paste0(gene," was not found in `gene.data`."),call.=FALSE)

    if (length(gene.match)>1L)
      stop(paste0("Multiple entries for ",gene,
                  " were found in `gene.data`."),call.=FALSE)

    chromosome.raw=gene.data$chrom[gene.match]
    region.start=as.numeric(gene.data$loc.start[gene.match])
    region.end=as.numeric(gene.data$loc.end[gene.match])
    region.size=region.end-region.start+1

    from=floor(region.start-expand*region.size)
    to=ceiling(region.end+expand*region.size)
  } else {
    chromosome.raw=chrom

    # Region plots use exactly the user-specified range
    from=as.numeric(plot.start)
    to=as.numeric(plot.end)
  }

  if (!is.finite(from) || !is.finite(to) || from>=to)
    stop("Invalid plotting coordinates.",call.=FALSE)

  chromosome.normalized=normalize.chromosome(chromosome.raw)
  chromosome.ucsc=to.ucsc.chromosome(chromosome.raw)

  if (gene.mode)
  {
    panel.title=paste0(
      gene," | chr",chromosome.normalized,": ",
      format(from,big.mark=",",scientific=FALSE),"-",
      format(to,big.mark=",",scientific=FALSE))
  } else {
    panel.title=paste0(
      if (!is.null(lesion.grp))
        paste0(paste(lesion.grp,collapse=", ")," | ")
      else "",
      "chr",chromosome.normalized,": ",
      format(from,big.mark=",",scientific=FALSE),"-",
      format(to,big.mark=",",scientific=FALSE))
  }

  ####################################################################
  # 11. Prepare and order lesion data
  ####################################################################

  lesion.data=grin.res[["lsn.data"]]

  if (is.null(lesion.data))
    stop("`grin.res[['lsn.data']]` was not found.",call.=FALSE)

  required.lesion.columns=c("chrom","loc.start","loc.end","lsn.type")
  missing.lesion.columns=setdiff(required.lesion.columns,colnames(lesion.data))

  if (length(missing.lesion.columns)>0L)
    stop(paste0("`lsn.data` is missing: ",
                paste(missing.lesion.columns,collapse=", "),"."),
         call.=FALSE)

  # Save source ordering before filtering and sorting
  lesion.data$.original.order=seq_len(nrow(lesion.data))

  lesion.data=lesion.data[
    normalize.chromosome(lesion.data$chrom)==chromosome.normalized,
    ,drop=FALSE
  ]

  if (!is.null(lesion.grp))
    lesion.data=lesion.data[
      as.character(lesion.data$lsn.type) %in% lesion.grp,,drop=FALSE
    ]

  lesion.data$loc.start=suppressWarnings(as.numeric(lesion.data$loc.start))
  lesion.data$loc.end=suppressWarnings(as.numeric(lesion.data$loc.end))

  lesion.data=lesion.data[
    is.finite(lesion.data$loc.start) &
      is.finite(lesion.data$loc.end) &
      lesion.data$loc.start<=to &
      lesion.data$loc.end>=from,
    ,drop=FALSE
  ]

  if (nrow(lesion.data)==0L)
    stop("No lesions overlap the selected plotting region.",call.=FALSE)

  lesion.data$size=lesion.data$loc.end-lesion.data$loc.start+1

  # Primary order: lesion category
  # Within each category: lesion size
  # Original source order is used only to break exact ties
  lesion.data=lesion.data[
    order(as.character(lesion.data$lsn.type),
          lesion.data$size,
          lesion.data$.original.order),
    ,drop=FALSE
  ]

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

  if (!is.null(lesion.grp) &&
      length(lesion.grp)==1L &&
      !is.null(spec.lsn.clr))
  {
    lesion.colors=spec.lsn.clr
    names(lesion.colors)=lesion.grp
  } else {
    lesion.colors=obtain.lesion.colors(
      lesion.types=lesion.types,
      supplied.colors=lsn.clrs
    )
  }

  ####################################################################
  # 12. Retrieve and filter transcript annotations
  ####################################################################

  transcript.object=NULL

  if (isTRUE(transTrack))
  {
    if (is.null(hg38.transcripts))
      stop(paste0("`hg38.transcripts` must be supplied when ",
                  "`transTrack = TRUE`."),call.=FALSE)

    annotation.database=hg38.transcripts

    try(GenomeInfoDb::seqlevelsStyle(annotation.database) <- "UCSC",
        silent=TRUE)

    if (gene.mode)
    {
      transcript.object=ensembldb::getGeneRegionTrackForGviz(
        annotation.database,
        filter=~ genename==gene)
    } else {
      transcript.object=ensembldb::getGeneRegionTrackForGviz(
        annotation.database,
        chromosome=chromosome.ucsc,
        start=from,
        end=to)
    }

    if (is.null(transcript.object) || length(transcript.object)==0L)
      stop("No transcript annotations were found for the selected region.",
           call.=FALSE)

    # Keep transcript elements overlapping the exact plotting range
    transcript.range=GenomicRanges::GRanges(
      seqnames=chromosome.ucsc,
      ranges=IRanges::IRanges(start=from,end=to))

    overlap.hits=GenomicRanges::findOverlaps(
      transcript.object,
      transcript.range,
      ignore.strand=TRUE)

    query.indices=unique(S4Vectors::queryHits(overlap.hits))
    transcript.object=transcript.object[query.indices]

    if (length(transcript.object)==0L)
      stop(paste0("Transcript annotations were retrieved, but none overlap ",
                  "the exact selected plotting range."),call.=FALSE)

    if (transcript.mode=="selected")
      transcript.object=filter.selected.transcript(
        transcript.object=transcript.object,
        selected.transcript=selected.transcript
      )
  }

  ####################################################################
  # 13. Construct genomic annotation tracks
  ####################################################################

  upper.track.grob=NULL

  if (isTRUE(transTrack) || isTRUE(show.ideogram) || isTRUE(show.axis))
  {
    upper.track.grob=make.gviz.track.grob(
      chromosome=chromosome.ucsc,
      from=from,
      to=to,
      transcript.object=transcript.object,
      cytoband.object=hg38.cytoband,
      include.transcripts=transTrack,
      include.ideogram=show.ideogram,
      include.axis=show.axis,
      genome.name=genome)

    upper.track.grob=add.fixed.label.column(
      plot.grob=upper.track.grob,
      label=if (isTRUE(transTrack)) transcript.label else "",
      width.inches=label.width)
  }

  ####################################################################
  # 14. Construct lesion track
  ####################################################################

  lesion.track.grob=make.lesion.track.grob(
    lesion.data=lesion.data,
    from=from,
    to=to,
    colors=lesion.colors,
    point.groups=point.lesion.groups,
    point.size=point.size.mm,
    interval.outline=interval.border,
    interval.rectangle.height=interval.height,
    backgrounds=row.background,
    coordinate.font.size=10*coordinate.cex)

  # Reduce the physical width of the lesion panel without changing
  # genomic coordinates
  lesion.track.grob=add.horizontal.inset(
    plot.grob=lesion.track.grob,
    inset.mm=lesion.track.inset.mm)

  lesion.track.grob=add.fixed.label.column(
    plot.grob=lesion.track.grob,
    label=lesion.label,
    width.inches=label.width)

  ####################################################################
  # 15. Construct legend and statistics
  ####################################################################

  legend.grob=make.legend.grob(
    colors=lesion.colors,
    fontsize=10*legend.cex)

  legend.grob=add.fixed.label.column(
    plot.grob=legend.grob,
    label="",
    width.inches=label.width)

  statistics.grobs=make.gene.statistics.grobs(
    grin.object=grin.res,
    requested.gene=if (gene.mode) gene else NULL,
    lesion.types=names(lesion.colors),
    base.size=10*stats.cex)

  lesion.statistics.grob=add.fixed.label.column(
    plot.grob=statistics.grobs$lesion,
    label="",
    width.inches=label.width)

  constellation.statistics.grob=add.fixed.label.column(
    plot.grob=statistics.grobs$constellation,
    label="",
    width.inches=label.width)

  ####################################################################
  # 16. Assemble final layout
  ####################################################################

  title.grob=grid::textGrob(
    label=panel.title,
    gp=grid::gpar(fontsize=13*title.cex,fontface="bold"))

  layout.grobs=list(title.grob)
  layout.heights=c(title.height)

  if (!is.null(upper.track.grob))
  {
    # A full transcript panel requires substantially more vertical space than
    # an ideogram or genomic-coordinate axis displayed without transcripts.
    if (isTRUE(transTrack))
      upper.track.height=transcript.height
    else if (isTRUE(show.ideogram) && isTRUE(show.axis))
      upper.track.height=0.12
    else
      upper.track.height=0.07

    layout.grobs=c(layout.grobs,list(upper.track.grob))
    layout.heights=c(layout.heights,upper.track.height)
  }

  layout.grobs=c(layout.grobs,list(lesion.track.grob,legend.grob))
  layout.heights=c(layout.heights,lesion.height,legend.height)

  if (statistics.grobs$lesion.available)
  {
    layout.grobs=c(layout.grobs,list(lesion.statistics.grob))
    layout.heights=c(layout.heights,lesion.stats.height)
  }

  if (statistics.grobs$constellation.available)
  {
    layout.grobs=c(layout.grobs,list(constellation.statistics.grob))
    layout.heights=c(layout.heights,constellation.stats.height)
  }

  layout.heights=layout.heights/sum(layout.heights)

  final.grob=gridExtra::arrangeGrob(
    grobs=layout.grobs,
    ncol=1,
    heights=layout.heights)

  ####################################################################
  # 17. Draw and return
  ####################################################################

  if (isTRUE(draw))
  {
    grid::grid.newpage()
    grid::grid.draw(final.grob)
  }

  invisible(final.grob)
}