R/count.hits.R

In GRIN2: Genomic Random Interval (GRIN)

#' Count Gene-Lesion Hits
#'
#' @description
#' Computes the number of genomic lesions affecting each gene by lesion type
#' and the number of unique subjects whose lesions overlap each gene for each
#' lesion type.
#'
#' @usage
#' count.hits(ov.data)
#'
#' @param ov.data A list returned by `find.gene.lsn.overlaps()` containing
#' gene-lesion overlap data, processed gene and lesion annotations, and
#' supporting index objects. 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
#' This function summarizes the output of `find.gene.lsn.overlaps()` by
#' generating two matrices:
#'
#' \describe{
#'   \item{\code{nhit.mtx}}{The total number of overlapping lesions affecting
#'   each gene, categorized by lesion type. Multiple lesions of the same type
#'   in the same subject are counted separately.}
#'   \item{\code{nsubj.mtx}}{The number of unique subjects with at least one
#'   overlapping lesion affecting each gene, categorized by lesion type.
#'   Multiple lesions of the same type affecting the same gene in one subject
#'   are counted once.}
#' }
#'
#' For example, if three separate mutations from the same subject overlap
#' \emph{NOTCH1}, all three lesions are counted in `nhit.mtx`, whereas that
#' subject is counted once in `nsubj.mtx`.
#'
#' Exon-level target sizes do not alter the hit or affected-subject counts.
#' Gene-lesion overlaps are counted using the standard genomic coordinates of
#' each gene, including for lesion types specified in `exon_level`. The
#' gene-level exon target sizes, chromosome-level exon target sizes, and
#' exon-level lesion-type specification are retained in the returned list for
#' downstream GRIN probability calculations.
#'
#' All genes represented in `gene.data` remain included in `nhit.mtx` and
#' `nsubj.mtx`, regardless of whether they have a valid matching exon
#' annotation. Genes without valid exon target sizes remain available for
#' standard GRIN analyses but do not receive exon-level probability estimates
#' for lesion types specified in `exon_level`.
#'
#' @return
#' A list containing the following components:
#' \describe{
#'   \item{lsn.data}{Processed lesion data.}
#'   \item{lsn.index}{A `data.frame` indexing lesion groups defined by lesion
#'   type, chromosome, and subject.}
#'   \item{gene.data}{Processed gene annotation data.}
#'   \item{gene.index}{A `data.frame` indexing genes by chromosome.}
#'   \item{nhit.mtx}{A numeric matrix in which rows correspond to genes and
#'   columns correspond to lesion types. Each value is the number of lesions
#'   of the specified type affecting the gene.}
#'   \item{nsubj.mtx}{A numeric matrix with the same dimensions as `nhit.mtx`.
#'   Each value is the number of unique subjects with at least one lesion of
#'   the specified type affecting the gene.}
#'   \item{gene.lsn.data}{A `data.frame` in which each row represents a gene
#'   overlapped by a genomic lesion.}
#'   \item{glp.data}{The combined gene and lesion position table. The `cty`
#'   column identifies the boundary type: 1 = gene start, 2 = lesion start,
#'   3 = lesion end, and 4 = gene end.}
#'   \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
#'   included in the hit and affected-subject counts and in standard GRIN
#'   analyses, but exon-level probability calculations are not performed 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 retained 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{find.gene.lsn.overlaps}},
#' \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 overlapping gene-lesion events
#' gene.lsn.overlap <- find.gene.lsn.overlaps(prep.gene.lsn)
#'
#' # Count lesions and affected subjects for each gene and lesion type
#' count.nsubj.nhits <- count.hits(gene.lsn.overlap)
#'
count.hits=function(ov.data) # output results of find.gene.lsn.overlaps function
{
  lsn.data=ov.data$lsn.data
  lsn.index=ov.data$lsn.index
  gene.lsn.hits=ov.data$gene.lsn.hits
  gene.lsn.data=ov.data$gene.lsn.data
  gene.data=ov.data$gene.data
  gene.index=ov.data$gene.index
  gene.exon.size=ov.data$gene.exon.size
  exon.chrom.size=ov.data$exon.chrom.size
  exon_level=ov.data$exon_level

  g=nrow(gene.data)

  # Compute the number of hits matrix

  lsn.types=sort(unique(lsn.index[,"lsn.type"]))
  k=length(lsn.types)

  nhit.mtx=matrix(0,g,k)
  colnames(nhit.mtx)=lsn.types

  if (nrow(gene.lsn.hits)>0)
  {
    nhit.tbl=table(gene.lsn.hits$gene.row,
                   gene.lsn.hits$lsn.type)
    nhit.rows=as.numeric(rownames(nhit.tbl))

    for (i in seq_len(ncol(nhit.tbl)))
      nhit.mtx[nhit.rows,colnames(nhit.tbl)[i]]=nhit.tbl[,i]
  }

  # Compute the matrix of the number of subjects with a hit

  nsubj.mtx=matrix(0,g,k)
  colnames(nsubj.mtx)=lsn.types

  if (nrow(gene.lsn.hits)>0)
  {
    # identify unique gene-subject-lesion type combinations
    gene.subj.type=paste0(gene.lsn.hits$gene.row,"_",
                          gene.lsn.hits$ID,"_",
                          gene.lsn.hits$lsn.type)
    dup.gene.subj.type=duplicated(gene.subj.type)

    # retain one record per subject for each gene and lesion type
    subj.gene.hits=gene.lsn.hits[!dup.gene.subj.type,,drop=FALSE]

    nsubj.tbl=table(subj.gene.hits$gene.row,
                    subj.gene.hits$lsn.type)
    nsubj.rows=as.numeric(rownames(nsubj.tbl))

    for (i in seq_len(ncol(nsubj.tbl)))
      nsubj.mtx[nsubj.rows,colnames(nsubj.tbl)[i]]=nsubj.tbl[,i]
  }

  res=list(lsn.data=lsn.data,             # processed lesion data
           lsn.index=lsn.index,           # lesion type-chromosome-subject index
           gene.data=gene.data,           # processed gene annotation data
           gene.index=gene.index,         # chromosome index for gene data
           nhit.mtx=nhit.mtx,             # number of lesions affecting each gene by lesion type
           nsubj.mtx=nsubj.mtx,           # number of affected subjects by lesion type
           gene.lsn.data=gene.lsn.hits,   # gene-lesion overlap records
           glp.data=gene.lsn.data,        # combined gene and lesion position data
           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)
}