R/prep.gene.lsn.data.R

In GRIN2: Genomic Random Interval (GRIN)

#' Prepare Gene and Lesion Data for GRIN Analysis
#'
#' @description
#' Prepares and indexes gene and lesion data for downstream GRIN
#' (Genomic Random Interval) analysis. The function merges and orders gene
#' and lesion coordinates to support efficient computation of overlaps between
#' genes and genomic lesions. When requested, it also prepares gene- and
#' chromosome-level exon target sizes for lesion types restricted to exonic
#' regions.
#'
#' @usage
#' prep.gene.lsn.data(lsn.data,
#'                    gene.data,
#'                    exons.annotation = NULL,
#'                    exon.chrom.size = NULL,
#'                    exon_level = NULL,
#'                    mess.freq = 10)
#'
#' @param lsn.data A `data.frame` containing lesion data in GRIN-compatible
#' format. The following five columns are required:
#' \describe{
#'   \item{ID}{Unique patient identifier.}
#'   \item{chrom}{Chromosome on which the lesion is located.}
#'   \item{loc.start}{Start position of the lesion in base pairs.}
#'   \item{loc.end}{End position of the lesion in base pairs.}
#'   \item{lsn.type}{Type of genomic lesion, such as mutation, breakpoint,
#'   gain, amplification, heterozygous deletion, or homozygous deletion.}
#' }
#'
#' @param gene.data A `data.frame` containing gene annotation data with the
#' following four required columns:
#' \describe{
#'   \item{gene}{Ensembl gene identifier.}
#'   \item{chrom}{Chromosome on which the gene is located.}
#'   \item{loc.start}{Start position of the gene in base pairs.}
#'   \item{loc.end}{End position of the gene in base pairs.}
#' }
#'
#' @param exons.annotation An optional `data.frame` containing exon annotation
#' data. This argument is required when `exon_level` is specified. The exon
#' annotation may contain the same gene set as `gene.data` or a subset of
#' those genes. Genes in `gene.data` without a valid matching exon annotation
#' remain in standard GRIN analyses, but exon-level probabilities are not
#' computed for those genes. Compatible exon annotation data for supported
#' genome assemblies can be retrieved using [get.ensembl.annotation()].
#' The following four columns are required:
#' \describe{
#'   \item{gene}{Ensembl gene identifier of the gene to which each annotated
#'   exon belongs. Multiple rows may therefore share the same gene identifier,
#'   with each row representing a different exon of that gene.}
#'   \item{chrom}{Chromosome on which the exon is located.}
#'   \item{loc.start}{Start position of the exon in base pairs.}
#'   \item{loc.end}{End position of the exon in base pairs.}
#' }
#'
#' @param exon.chrom.size An optional `data.frame` containing the total
#' annotated exon target size for each chromosome. This argument is required
#' when `exon_level` is specified and must contain:
#' \describe{
#'   \item{chrom}{Chromosome identifier.}
#'   \item{size}{Total annotated exon target size for the chromosome in base
#'   pairs.}
#' }
#' The chromosome-level exon sizes should be derived from a genome-wide exon
#' annotation using the same genome assembly and representative-transcript
#' selection procedure used to construct `exons.annotation`. A compatible
#' chromosome-level exon-size dataset is provided with GRIN2 and can be loaded
#' using `data(hg38_exon_chrom_size)`. Corresponding annotation data can also
#' be obtained using [get.ensembl.annotation()].
#'
#' @param exon_level An optional character vector specifying the lesion types
#' that should use exon-level target sizes in downstream GRIN probability
#' calculations. Lesions belonging to these types must be restricted to exonic
#' regions before analysis. For example, when mutations are specified for
#' exon-level analysis, exonic variants such as missense, nonsense, frameshift,
#' and synonymous mutations may be included, whereas intronic and other
#' non-coding mutations should be excluded before running the analysis.
#' The default is `NULL`.
#'
#' @param mess.freq Integer specifying the frequency of progress messages.
#' The default is 10.
#'
#' @details
#' The function first orders and indexes the lesion data by lesion type,
#' chromosome, and subject and orders the gene annotation by chromosome and
#' genomic position. Gene and lesion boundaries are then combined into a
#' unified genomic position table.
#'
#' The `cty` column in the combined table identifies the boundary represented
#' by each row:
#' \describe{
#'   \item{1}{Gene start.}
#'   \item{2}{Lesion start.}
#'   \item{3}{Lesion end.}
#'   \item{4}{Gene end.}
#' }
#'
#' The resulting table and index objects are used by
#' `find.gene.lsn.overlaps()` to identify gene-lesion overlaps.
#'
#' When `exon_level` is specified, exon lengths are calculated from
#' `exons.annotation` and summed to determine the exon target size for each
#' matched gene. Chromosome-level exon target sizes are not calculated by
#' this function and must instead be supplied through `exon.chrom.size`.
#' `exons.annotation` may contain all genes in `gene.data` or only a subset.
#' Genes without matching exon annotations remain available for lesion types
#' analyzed using standard genomic coordinates, but their gene-level exon
#' target sizes are recorded as missing and exon-level probabilities are not
#' computed for those genes.
#'
#' Gene-lesion overlaps continue to be determined using the complete genomic
#' coordinates of each gene, including when `exon_level` is specified.
#' Exon annotation is used only to calculate exon-based target sizes for
#' downstream probability calculations. Therefore, lesion types specified in
#' `exon_level` must contain only exonic lesions; non-exonic lesions must be
#' removed by the user before running the analysis.
#'
#' Lesion types not included in `exon_level` continue to use the standard
#' chromosome-level GRIN workflow.
#'
#' Exact duplicate exon records are removed before gene-level exon sizes are
#' calculated. Incomplete or malformed exon records are excluded with a
#' warning.
#'
#' @return
#' A list with the following components:
#' \describe{
#'   \item{lsn.data}{Processed lesion data, including indices identifying the
#'   range of rows in `gene.lsn.data` corresponding to each lesion.}
#'   \item{gene.data}{Processed gene annotation data, including indices
#'   identifying the range of rows in `gene.lsn.data` corresponding to each
#'   gene.}
#'   \item{gene.lsn.data}{Combined and ordered `data.frame` of gene and lesion
#'   positions. The `cty` column encodes the position type: 1 = gene start,
#'   2 = lesion start, 3 = lesion end, and 4 = gene end.}
#'   \item{gene.index}{Index `data.frame` indicating the ordered start and end
#'   rows for each chromosome in the gene data.}
#'   \item{lsn.index}{Index `data.frame` indicating the ordered start and end
#'   rows for lesion groups defined by lesion type, chromosome, and subject.}
#'   \item{gene.exon.size}{Numeric vector containing the total annotated exon
#'   target size for each gene, aligned by `gene.row`. Genes without a valid
#'   matching exon annotation receive `NA` and are excluded from probability
#'   calculations for lesion types specified in `exon_level`. Returned as
#'   `NULL` when `exon_level = NULL`.}
#'   \item{exon.chrom.size}{A `data.frame` containing chromosome identifiers
#'   in the `chrom` column and total genome-wide annotated exon target sizes in
#'   the `size` column. Returned as `NULL` when `exon_level = NULL`.}
#'   \item{exon_level}{The character vector of lesion types designated for
#'   exon-level analysis. Returns `NULL` when exon-level preprocessing is not
#'   requested.}
#' }
#'
#' @export
#'
#' @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{order.index.gene.data}},
#' \code{\link{order.index.lsn.data}},
#' \code{\link{find.gene.lsn.overlaps}}
#'
#' @examples
#' data(lesion_data)
#' data(hg38_gene_annotation)
#' data(example_exon_annotation)
#' data(hg38_exon_chrom_size)
#'
#' # Prepare gene and lesion data using the optional arguments
#' # for exon-level analysis
#' prep.gene.lsn <- prep.gene.lsn.data(
#'   lsn.data = lesion_data,
#'   gene.data = hg38_gene_annotation,
#'   exons.annotation = example_exon_annotation,
#'   exon.chrom.size = hg38_exon_chrom_size,
#'   exon_level = "mutation"
#' )
#'
prep.gene.lsn.data=function(lsn.data,               # lesion data in GRIN-compatible format
                            gene.data,              # gene annotation data
                            exons.annotation=NULL,  # optional exon annotation for exon-level analysis
                            exon.chrom.size=NULL,   # optional chromosome-level exon target sizes
                            exon_level=NULL,        # lesion types to analyze using exon-level target sizes
                            mess.freq=10)           # message frequency
{
  old_opt <- options(stringsAsFactors = FALSE)
  on.exit(options(old_opt), add = TRUE)

  # order lesion data by type, chromosome, and subject
  lsn.dset=order.index.lsn.data(lsn.data)
  lsn.data=lsn.dset$lsn.data
  lsn.index=lsn.dset$lsn.index

  # order and index gene locus data by chromosome and position
  gene.dset=order.index.gene.data(gene.data)
  gene.data=gene.dset$gene.data
  gene.index=gene.dset$gene.index

  # Extract some basic information
  g=nrow(gene.data) # number of genes
  l=nrow(lsn.data)  # number of lesions

  # Initialize exon-level gene target sizes
  gene.exon.size=NULL

  # Prepare exon-level target sizes for selected lesion types
  if (!is.null(exon_level))
  {
    # Validate exon annotation inputs
    if (is.null(exons.annotation))
      stop("'exons.annotation' must be provided when 'exon_level' is specified.",
           call.=FALSE)

    if (is.null(exon.chrom.size))
      stop("'exon.chrom.size' must be provided when 'exon_level' is specified.",
           call.=FALSE)

    required.exon.columns=c("gene","chrom","loc.start","loc.end")
    missing.exon.columns=setdiff(required.exon.columns,colnames(exons.annotation))

    if (length(missing.exon.columns)>0)
      stop("'exons.annotation' must contain the following columns: ",
           paste(required.exon.columns,collapse=", "),".",call.=FALSE)

    required.exon.chrom.columns=c("chrom","size")
    missing.exon.chrom.columns=setdiff(required.exon.chrom.columns,
                                       colnames(exon.chrom.size))

    if (length(missing.exon.chrom.columns)>0)
      stop("'exon.chrom.size' must contain the following columns: ",
           paste(required.exon.chrom.columns,collapse=", "),".",call.=FALSE)

    warning(
      "For exon-level analysis, lesion types specified in 'exon_level' should ",
      "contain only exonic lesions. Intronic and other non-coding lesions should ",
      "be excluded before running the analysis.",
      call.=FALSE
    )

    message(paste0("Preparing exon-level gene and chromosome target sizes: ",date()))

    # Use consistent temporary formats for annotation matching
    ex=exons.annotation[,required.exon.columns,drop=FALSE]
    ex$gene=as.character(ex$gene)
    ex$chrom=as.character(ex$chrom)
    ex$loc.start=suppressWarnings(as.numeric(as.character(ex$loc.start)))
    ex$loc.end=suppressWarnings(as.numeric(as.character(ex$loc.end)))

    # Remove incomplete or malformed exon records
    valid.exons=!is.na(ex$gene) &
      nzchar(trimws(ex$gene)) &
      !is.na(ex$chrom) &
      nzchar(trimws(ex$chrom)) &
      !is.na(ex$loc.start) &
      !is.na(ex$loc.end) &
      ex$loc.start>=1 &
      ex$loc.end>=ex$loc.start

    n.invalid.exons=sum(!valid.exons)

    if (n.invalid.exons>0)
      warning(n.invalid.exons,
              " incomplete or malformed exon annotation row(s) were removed.",
              call.=FALSE)

    ex=ex[valid.exons,,drop=FALSE]

    if (nrow(ex)==0)
      stop("No valid exon records remained in 'exons.annotation'.",
           call.=FALSE)

    # Remove exact duplicate exon records and calculate exon lengths
    ex=unique(ex)
    ex$exon.len=ex$loc.end-ex$loc.start+1

    # Calculate the total exon target size for each gene
    gene.exon.summary=aggregate(exon.len~gene+chrom,data=ex,FUN=sum)
    colnames(gene.exon.summary)[
      colnames(gene.exon.summary)=="exon.len"
    ]="exon.size"

    # Match exon annotations to the ordered gene annotation
    exon.gene.key=paste(gene.exon.summary$gene,
                        gene.exon.summary$chrom,sep="\r")
    gene.data.key=paste(as.character(gene.data$gene),
                        as.character(gene.data$chrom),sep="\r")
    gene.match=match(exon.gene.key,gene.data.key)
    unmatched.exons=is.na(gene.match)

    if (any(unmatched.exons))
      warning(sum(unmatched.exons),
              " gene-level exon annotation record(s) could not be matched ",
              "to 'gene.data' and were excluded.",call.=FALSE)

    # Build an exon-size vector aligned with the ordered gene data
    gene.exon.size=rep(NA_real_,g)
    matched.exons=!unmatched.exons

    if (any(matched.exons))
    {
      gene.rows=gene.data$gene.row[gene.match[matched.exons]]
      gene.exon.size[gene.rows]=gene.exon.summary$exon.size[matched.exons]
    }

    # Identify genes without a valid exon target size
    missing.gene.exons=is.na(gene.exon.size) | gene.exon.size<=0

    if (any(missing.gene.exons))
      warning(sum(missing.gene.exons),
              " gene(s) in 'gene.data' did not have a valid matching exon ",
              "annotation. These genes will remain in the standard GRIN analyses, ",
              "but probability calculations will not be performed for lesion types ",
              "specified in 'exon_level'.",call.=FALSE)

    # Prepare and validate chromosome-level exon target sizes
    exon.chrom.size=exon.chrom.size[,required.exon.chrom.columns,drop=FALSE]
    exon.chrom.size$chrom=as.character(exon.chrom.size$chrom)
    exon.chrom.size$size=suppressWarnings(
      as.numeric(as.character(exon.chrom.size$size))
    )

    valid.exon.chrom.rows=!is.na(exon.chrom.size$chrom) &
      nzchar(trimws(exon.chrom.size$chrom)) &
      !is.na(exon.chrom.size$size) &
      exon.chrom.size$size>0

    if (!all(valid.exon.chrom.rows))
      stop("'exon.chrom.size' contains missing or invalid chromosome sizes.",
           call.=FALSE)

    if (anyDuplicated(exon.chrom.size$chrom)>0)
      stop("'exon.chrom.size' must contain only one row per chromosome.",
           call.=FALSE)

    # Confirm exon target sizes are available for chromosomes used in exon-level analysis
    exon.lsn.chrom=unique(as.character(
      lsn.data$chrom[lsn.data$lsn.type %in% exon_level]
    ))
    missing.exon.chrom=setdiff(exon.lsn.chrom,exon.chrom.size$chrom)

    if (length(missing.exon.chrom)>0)
      stop("'exon.chrom.size' does not contain exon target sizes for ",
           "the following chromosome(s) represented by lesion types in ",
           "'exon_level': ",paste(missing.exon.chrom,collapse=", "),".",
           call.=FALSE)
  }

  # Create gene position data
  message(paste0("Formatting gene position data for counting: ",date()))

  gene.pos.data=rbind.data.frame(
    cbind.data.frame(ID="", # gene start data
                     lsn.type="",
                     lsn.row=NA,
                     gene=gene.data[,"gene"],
                     gene.row=gene.data[,"gene.row"],
                     chrom=gene.data[,"chrom"],
                     pos=gene.data[,"loc.start"],
                     cty=1),
    cbind.data.frame(ID="", # gene end data
                     lsn.type="",
                     lsn.row=NA,
                     gene=gene.data[,"gene"],
                     gene.row=gene.data[,"gene.row"],
                     chrom=gene.data[,"chrom"],
                     pos=gene.data[,"loc.end"],
                     cty=4))

  # order gene position data
  ord=order(gene.pos.data[,"chrom"],
            gene.pos.data[,"pos"],
            gene.pos.data[,"cty"])
  gene.pos.data=gene.pos.data[ord,]

  # Create lesion position data with one row for each edge of each lesion
  message(paste0("Formatting lesion position data for counting: ",date()))

  lsn.pos.data=rbind.data.frame(
    cbind.data.frame(ID=lsn.data[,"ID"],
                     lsn.type=lsn.data[,"lsn.type"],
                     lsn.row=lsn.data[,"lsn.row"],
                     gene="",
                     gene.row=NA,
                     chrom=lsn.data[,"chrom"],
                     pos=lsn.data[,"loc.start"],
                     cty=2),
    cbind.data.frame(ID=lsn.data[,"ID"],
                     lsn.type=lsn.data[,"lsn.type"],
                     lsn.row=lsn.data[,"lsn.row"],
                     gene="",
                     gene.row=NA,
                     chrom=lsn.data[,"chrom"],
                     pos=lsn.data[,"loc.end"],
                     cty=3))

  # order lesion position data
  ord=order(lsn.pos.data[,"chrom"],
            lsn.pos.data[,"pos"],
            lsn.pos.data[,"cty"])
  lsn.pos.data=lsn.pos.data[ord,]

  # Combine gene & lesion data
  message(paste0("Combining formatted gene and lesion position data: ",date()))
  gene.lsn.data=rbind.data.frame(gene.pos.data,lsn.pos.data)

  # Order and index gene & lesion data
  ord=order(gene.lsn.data[,"chrom"],
            gene.lsn.data[,"pos"],
            gene.lsn.data[,"cty"])
  gene.lsn.data=gene.lsn.data[ord,]
  m=nrow(gene.lsn.data)
  gene.lsn.data[,"glp.row"]=1:m

  # compute vector to order gene.lsn.data by lsn.row and gene.row
  ord=order(gene.lsn.data[,"lsn.row"],
            gene.lsn.data[,"gene.row"],
            gene.lsn.data[,"cty"])

  # use that vector to add gene.lsn.data row.start and row.end indices to lsn.data
  lsn.pos=gene.lsn.data[ord[1:(2*l)],]
  lsn.data[,"glp.row.start"]=lsn.pos[2*(1:l)-1,"glp.row"]
  lsn.data[,"glp.row.end"]=lsn.pos[2*(1:l),"glp.row"]

  # use that vector to add gene.lsn.data row.start and row.end indices to gene.data
  gene.pos=gene.lsn.data[ord[-(1:(2*l))],]
  gene.data[,"glp.row.start"]=gene.pos[2*(1:g)-1,"glp.row"]
  gene.data[,"glp.row.end"]=gene.pos[2*(1:g),"glp.row"]

  # Double-check table pointers from lsn.data and gene.data to gene.lsn.data
  message(paste0("Verifying structure of combined gene and lesion data: ",date()))

  glp.gene.start=gene.lsn.data[gene.data$glp.row.start,c("gene","chrom","pos")]
  colnames(glp.gene.start)=c("gene","chrom","loc.start")
  ok.glp.gene.start=all(
    glp.gene.start==gene.data[,c("gene","chrom","loc.start")]
  )

  glp.gene.end=gene.lsn.data[gene.data$glp.row.end,c("gene","chrom","pos")]
  colnames(glp.gene.end)=c("gene","chrom","loc.end")
  ok.glp.gene.end=all(
    glp.gene.end==gene.data[,c("gene","chrom","loc.end")]
  )

  glp.lsn.start=gene.lsn.data[
    lsn.data$glp.row.start,c("ID","chrom","pos","lsn.type")
  ]
  colnames(glp.lsn.start)=c("ID","chrom","loc.start","lsn.type")
  ok.glp.lsn.start=all(
    glp.lsn.start==lsn.data[,c("ID","chrom","loc.start","lsn.type")]
  )

  glp.lsn.end=gene.lsn.data[
    lsn.data$glp.row.end,c("ID","chrom","pos","lsn.type")
  ]
  colnames(glp.lsn.end)=c("ID","chrom","loc.end","lsn.type")
  ok.glp.lsn.end=all(
    glp.lsn.end==lsn.data[,c("ID","chrom","loc.end","lsn.type")]
  )

  # Double-check table pointers from gene.lsn.data to gene.data and lsn.data
  glp.gene.start=gene.lsn.data[
    gene.lsn.data$cty==1,c("gene.row","gene","chrom","pos")
  ]
  ok.gene.start=all(
    glp.gene.start[,c("gene","chrom","pos")]==
      gene.data[glp.gene.start$gene.row,c("gene","chrom","loc.start")]
  )

  glp.gene.end=gene.lsn.data[
    gene.lsn.data$cty==4,c("gene.row","gene","chrom","pos")
  ]
  ok.gene.end=all(
    glp.gene.end[,c("gene","chrom","pos")]==
      gene.data[glp.gene.end$gene.row,c("gene","chrom","loc.end")]
  )

  glp.lsn.start=gene.lsn.data[
    gene.lsn.data$cty==2,c("lsn.row","ID","chrom","pos","lsn.type")
  ]
  ok.lsn.start=all(
    glp.lsn.start[,c("ID","chrom","pos","lsn.type")]==
      lsn.data[glp.lsn.start$lsn.row,c("ID","chrom","loc.start","lsn.type")]
  )

  glp.lsn.end=gene.lsn.data[
    gene.lsn.data$cty==3,c("lsn.row","ID","chrom","pos","lsn.type")
  ]
  ok.lsn.end=all(
    glp.lsn.end[,c("ID","chrom","pos","lsn.type")]==
      lsn.data[glp.lsn.end$lsn.row,c("ID","chrom","loc.end","lsn.type")]
  )

  all.ok=all(c(ok.glp.gene.start,ok.glp.gene.end,
               ok.glp.lsn.start,ok.glp.lsn.end,
               ok.gene.start,ok.gene.end,
               ok.lsn.start,ok.lsn.end))

  if (!all.ok)
    stop("Error in constructing and indexing combined lesion and gene data.")

  message(paste0(
    "Verified correct construction and indexing of combined lesion and gene data: ",
    date()
  ))

  return(list(lsn.data=lsn.data,             # processed lesion data
              gene.data=gene.data,           # processed gene annotation data
              gene.lsn.data=gene.lsn.data,   # combined gene and lesion position data
              gene.index=gene.index,         # chromosome index for gene data
              lsn.index=lsn.index,           # lesion type-chromosome-subject index
              gene.exon.size=gene.exon.size, # gene-level exon target sizes
              exon.chrom.size=exon.chrom.size, # chromosome-level exon target sizes
              exon_level=exon_level))        # lesion types using exon-level analysis
}