R/find.gene.lsn.overlaps.R

In GRIN2: Genomic Random Interval (GRIN)

#' Find Gene-Lesion Overlaps
#'
#' @description
#' Identifies overlaps between genes and genomic lesions using the output from
#' `prep.gene.lsn.data()`. The function detects instances in which a genomic
#' lesion spans or intersects the genomic coordinates of a gene.
#'
#' @usage
#' find.gene.lsn.overlaps(gl.data)
#'
#' @param gl.data A list returned by `prep.gene.lsn.data()` containing
#' processed and indexed gene and lesion data. When exon-level analysis was
#' requested, the list also contains gene- and chromosome-level exon target
#' sizes and the lesion types designated for exon-level analysis.
#'
#' @details
#' The function scans the combined gene and lesion position table generated by
#' `prep.gene.lsn.data()` and identifies lesions that overlap the genomic
#' coordinates of each gene. Each detected gene-lesion overlap is recorded for
#' subsequent counting and statistical analysis.
#'
#' Gene-lesion overlap detection is unchanged when `exon_level` is specified
#' and continues to use the complete genomic coordinates of each gene.
#' Exon-level gene and chromosome target sizes do not alter whether a lesion is
#' considered to overlap a gene. Instead, these objects are carried forward
#' for use in downstream GRIN probability calculations. Lesion types specified
#' in `exon_level` should therefore contain only exonic lesions, as described
#' in `prep.gene.lsn.data()`.
#'
#' All genes represented in `gene.data` are included in overlap detection,
#' regardless of whether they have a matching exon annotation. Genes without
#' a valid exon target size remain available for standard GRIN analyses but
#' will not receive probability estimates for lesion types specified in
#' `exon_level`.
#'
#' @return
#' A list containing the following components:
#' \describe{
#'   \item{lsn.data}{Processed lesion data.}
#'   \item{gene.data}{Processed gene annotation data.}
#'   \item{gene.lsn.data}{A `data.frame` ordered by chromosome and genomic
#'   position that contains both gene and lesion boundaries. The `cty` column
#'   identifies the position type: 1 = gene start, 2 = lesion start,
#'   3 = lesion end, and 4 = gene end.}
#'   \item{gene.lsn.hits}{A `data.frame` in which each row represents a gene
#'   overlapped by a genomic lesion. It contains the gene and lesion row
#'   indices, gene identifier and coordinates, patient or sample identifier,
#'   lesion coordinates, and lesion type.}
#'   \item{gene.index}{A `data.frame` indexing the rows corresponding to genes
#'   on each chromosome.}
#'   \item{lsn.index}{A `data.frame` indexing 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 have a value of `NA`. These genes remain
#'   available for standard GRIN analyses but are excluded from probability
#'   calculations for lesion types specified in `exon_level`. Returns `NULL`
#'   when exon-level analysis was not requested.}
#'   \item{exon.chrom.size}{A `data.frame` containing the genome-wide annotated
#'   exon target size for each chromosome. This object is carried forward for
#'   use in downstream exon-level probability calculations. Returns `NULL`
#'   when exon-level analysis was not requested.}
#'   \item{exon_level}{Character vector specifying the lesion types designated
#'   for exon-level analysis. Returns `NULL` when exon-level analysis was 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{prep.gene.lsn.data}},
#' \code{\link{count.hits}},
#' \code{\link{prob.hits}}
#'
#' @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"
#' )
#'
#' # Identify genes overlapped by genomic lesions
#' gene.lsn.overlap <- find.gene.lsn.overlaps(prep.gene.lsn)
#'
find.gene.lsn.overlaps=function(gl.data) # output of the prep.gene.lsn.data function
{
  gene.data=gl.data$gene.data
  lsn.data=gl.data$lsn.data
  lsn.index=gl.data$lsn.index
  gene.index=gl.data$gene.index
  gene.lsn.data=gl.data$gene.lsn.data
  gene.exon.size=gl.data$gene.exon.size
  exon.chrom.size=gl.data$exon.chrom.size
  exon_level=gl.data$exon_level

  m=nrow(gene.lsn.data)

  message(paste0("Scanning through combined lesion and gene data to find gene-lesion overlaps: ",date()))

  gene.row.mtch=NULL  # initialize vector for rows of gene data matched to rows of lesion data
  lsn.row.mtch=NULL   # initialize vector for rows of lesion data matched to rows of gene data
  current.genes=NULL  # initialize vector of genes overlapping the current point of the scan
  current.lsns=NULL   # initialize vector of lesions overlapping the current point of the scan

  for (i in seq_len(m)) # loop over rows of gene.lsn.data
  {
    # enter a gene
    if (gene.lsn.data$cty[i]==1)
    {
      # add this gene to the set of current genes
      current.genes=c(current.genes,gene.lsn.data$gene.row[i])

      # match this gene to the set of current lesions
      lsn.row.mtch=c(lsn.row.mtch,current.lsns)
      gene.row.mtch=c(gene.row.mtch,
                      rep(gene.lsn.data$gene.row[i],length(current.lsns)))
    }

    # exit a gene
    if (gene.lsn.data$cty[i]==4)
    {
      # drop this gene from the set of current genes
      current.genes=setdiff(current.genes,gene.lsn.data$gene.row[i])
    }

    # enter a lesion
    if (gene.lsn.data$cty[i]==2)
    {
      # match this lesion to the set of current genes
      lsn.row.mtch=c(lsn.row.mtch,
                     rep(gene.lsn.data$lsn.row[i],length(current.genes)))
      gene.row.mtch=c(gene.row.mtch,current.genes)

      # add this lesion to the set of current lesions
      current.lsns=c(current.lsns,gene.lsn.data$lsn.row[i])
    }

    # exit a lesion
    if (gene.lsn.data$cty[i]==3)
    {
      # drop this lesion from the set of current lesions
      current.lsns=setdiff(current.lsns,gene.lsn.data$lsn.row[i])
    }
  }

  message(paste0("Completed scan of combined gene-lesion data: ",date()))

  # Generate the gene-lesion hit data

  gene.lsn.hits=cbind.data.frame(
    gene.data[gene.row.mtch,
              c("gene.row","gene","chrom","loc.start","loc.end")],
    lsn.data[lsn.row.mtch,
             c("lsn.row","ID","chrom","loc.start","loc.end","lsn.type")])

  colnames(gene.lsn.hits)=c("gene.row","gene","gene.chrom","gene.loc.start","gene.loc.end",
                            "lsn.row","ID","lsn.chrom","lsn.loc.start","lsn.loc.end","lsn.type")

  res=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.lsn.hits=gene.lsn.hits,   # gene-lesion overlap records
           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

  return(res)
}