R/prob.hits.R

In GRIN2: Genomic Random Interval (GRIN)

#' Find Probability of Locus Hit
#'
#' @description
#' Computes the probability that each genomic locus, such as a gene or
#' regulatory region, is affected by one or more types of genomic lesions.
#' Statistical significance for lesion enrichment is estimated using a
#' convolution of independent but non-identically distributed Bernoulli
#' probabilities.
#'
#' @usage
#' prob.hits(hit.cnt, chr.size = NULL)
#'
#' @param hit.cnt A list returned by `count.hits()` containing the number of
#' subjects and lesions affecting each locus by lesion type, together with the
#' processed lesion and gene annotations. When exon-level analysis was
#' requested in `prep.gene.lsn.data()`, this object also contains the gene- and
#' chromosome-level exon target sizes and the lesion types designated for
#' exon-level analysis.
#'
#' @param chr.size A `data.frame` containing chromosome sizes for the
#' 22 autosomes and the X and Y chromosomes. It must contain the following
#' columns:
#' \describe{
#'   \item{chrom}{Chromosome identifier.}
#'   \item{size}{Chromosome length in base pairs.}
#' }
#'
#' @details
#' For each lesion type and chromosome, the function calculates the probability
#' that each lesion overlaps each affected locus according to the relative
#' sizes of the lesion, locus, and chromosome. These probabilities are combined
#' across lesions and subjects using the Poisson-binomial distribution to
#' estimate the probability of observing the recorded number of lesions or
#' affected subjects by chance.
#'
#' For lesion types not designated for exon-level analysis, locus size is
#' calculated using the complete genomic boundaries in `gene.data`, and
#' chromosome length is obtained from `chr.size`.
#'
#' For lesion types designated for exon-level analysis, the function instead
#' uses the total annotated exon size of each gene and the total annotated exon
#' target size of the corresponding chromosome. Lesion types specified in
#' `exon_level` must therefore contain only exonic lesions. For example, when
#' mutations are analyzed at the exon level, 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 exon-level specification and target sizes are carried through the GRIN
#' workflow within `hit.cnt`; users therefore specify the exon-level lesion
#' types only once when calling `prep.gene.lsn.data()`.
#'
#' Genes without a valid matching exon annotation remain in analyses of lesion
#' types that use standard genomic coordinates. For lesion types designated for
#' exon-level analysis, their lesion-specific p-values and q-values are returned
#' as `NA`. Ordered constellation statistics are also returned as `NA` for
#' genes that do not have a complete set of lesion-specific probabilities.
#'
#' The function calculates lesion-specific p-values for:
#' \describe{
#'   \item{Number of affected subjects}{The probability of observing the
#'   recorded number of unique subjects with at least one lesion affecting
#'   the locus.}
#'   \item{Number of lesion hits}{The probability of observing the recorded
#'   total number of lesions affecting the locus.}
#' }
#'
#' False discovery rate-adjusted q-values are calculated using the method of
#' Pounds and Cheng. When more than one lesion type is present, the function
#' also calculates ordered p- and q-values representing constellations of one
#' or more lesion types affecting the same locus.
#'
#' @return
#' A list containing the following components:
#' \describe{
#'   \item{gene.hits}{A `data.frame` containing locus annotations; the numbers
#'   of lesion hits and affected subjects by lesion type; lesion-specific
#'   p-values and q-values for hit counts (`p.nhit.*` and `q.nhit.*`) and
#'   affected-subject counts (`p.nsubj.*` and `q.nsubj.*`); and, when multiple
#'   lesion types are present, ordered constellation p-values and q-values
#'   (`p1.nhit`, `p2.nhit`, ..., `q1.nhit`, `q2.nhit`, ... and the
#'   corresponding `nsubj` columns).}
#'   \item{lsn.data}{Processed input lesion data.}
#'   \item{gene.data}{Processed input gene annotation data.}
#'   \item{gene.lsn.data}{A `data.frame` in which each row represents a locus
#'   overlapped by a genomic lesion.}
#'   \item{chr.size}{Chromosome size information used for standard
#'   chromosome-level calculations.}
#'   \item{gene.index}{A `data.frame` indexing loci by chromosome.}
#'   \item{lsn.index}{A `data.frame` indexing lesion groups by lesion type,
#'   chromosome, and subject.}
#'   \item{gene.exon.size}{A numeric vector containing the exon target size for
#'   each gene. Genes without a valid matching exon annotation have a value of
#'   `NA`; exon-level probability calculations are not performed for those
#'   genes. 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. These values provide the
#'   chromosome-level exonic search space used for lesion types designated for
#'   exon-level analysis. Returns `NULL` when exon-level analysis was not
#'   requested.}
#'   \item{exon_level}{A character vector specifying the lesion types evaluated
#'   using exon-level target sizes. Returns `NULL` when exon-level analysis was
#'   not requested.}
#' }
#'
#' @export
#'
#' @importFrom stats p.adjust
#'
#' @references
#' Pounds, S., et al. (2013). A genomic random interval model for statistical
#' analysis of genomic lesion data.
#'
#' Pounds, S., and Cheng, C. (2006). Robust estimation of the false discovery
#' rate.
#'
#' 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{count.hits}}
#'
#' @examples
#' data(lesion_data)
#' data(hg38_gene_annotation)
#' data(hg38_chrom_size)
#' data(example_exon_annotation)
#' data(hg38_exon_chrom_size)
#'
#' # 1) 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"
#' )
#'
#' # 2) Identify overlapping gene-lesion events
#' gene.lsn.overlap <- find.gene.lsn.overlaps(prep.gene.lsn)
#'
#' # 3) Count subjects and lesions affecting each gene
#' count.nsubj.nhits <- count.hits(gene.lsn.overlap)
#'
#' # 4) Compute p- and q-values for lesion enrichment
#' hits.prob <- prob.hits(count.nsubj.nhits, hg38_chrom_size)
#'
prob.hits=function(hit.cnt,    # output results of the count.hits function
                   chr.size=NULL) # chromosome size table for standard genomic calculations
{
  lsn.data=hit.cnt$lsn.data
  num.lsn=unique(lsn.data$lsn.type)

  gene.exon.size=hit.cnt$gene.exon.size
  exon.chrom.size=hit.cnt$exon.chrom.size
  exon_level=hit.cnt$exon_level

  # Convert chromosome exon sizes to a named numeric vector

  if (!is.null(exon.chrom.size))
  {
    exon.chrom.size.vec=exon.chrom.size$size
    names(exon.chrom.size.vec)=as.character(exon.chrom.size$chrom)
  } else {
    exon.chrom.size.vec=NULL
  }

  # Verify that exon-level target sizes are available when requested

  if (!is.null(exon_level))
  {
    if (is.null(gene.exon.size) || is.null(exon.chrom.size.vec))
      stop("Gene- and chromosome-level exon target sizes are required when ",
           "'exon_level' is specified. Run prep.gene.lsn.data() with ",
           "'exons.annotation', 'exon.chrom.size', and 'exon_level'.",
           call.=FALSE)
  }

  ###################
  # Order and index gene.lsn.data
  ###################

  ord=order(hit.cnt$gene.lsn.data$lsn.type,
            hit.cnt$gene.lsn.data$lsn.chrom,
            hit.cnt$gene.lsn.data$gene.row,
            hit.cnt$gene.lsn.data$ID)

  hit.cnt$gene.lsn.data=hit.cnt$gene.lsn.data[ord,]
  m=nrow(hit.cnt$gene.lsn.data)

  new.sect=which((hit.cnt$gene.lsn.data$gene.chrom[-1]!=hit.cnt$gene.lsn.data$gene.chrom[-m])|
                   (hit.cnt$gene.lsn.data$lsn.type[-1]!=hit.cnt$gene.lsn.data$lsn.type[-m])|
                   (hit.cnt$gene.lsn.data$gene.row[-1]!=hit.cnt$gene.lsn.data$gene.row[-m]))

  sect.start=c(1,new.sect+1)
  sect.end=c(new.sect,m)

  gene.lsn.index=cbind.data.frame(
    lsn.type=hit.cnt$gene.lsn.data$lsn.type[sect.start],
    chrom=hit.cnt$gene.lsn.data$gene.chrom[sect.start],
    gene.row=hit.cnt$gene.lsn.data$gene.row[sect.start],
    row.start=sect.start,
    row.end=sect.end,
    n.lsns=sect.end-sect.start+1)

  k=nrow(gene.lsn.index)

  new.chr=which(gene.lsn.index$chrom[-1]!=gene.lsn.index$chrom[-k])
  chr.start=c(1,new.chr+1)
  chr.end=c(new.chr,k)

  gene.lsn.chr.index=cbind.data.frame(
    lsn.type=gene.lsn.index$lsn.type[chr.start],
    chrom=gene.lsn.index$chrom[chr.start],
    row.start=chr.start,
    row.end=chr.end,
    n.rows=chr.end-chr.start+1)

  ###################
  # Create lesion chromosome index
  ###################

  nr.li=nrow(hit.cnt$lsn.index)

  new.chr=which((hit.cnt$lsn.index$lsn.type[-1]!=hit.cnt$lsn.index$lsn.type[-nr.li])|
                  (hit.cnt$lsn.index$chrom[-1]!=hit.cnt$lsn.index$chrom[-nr.li]))

  chr.start=c(1,new.chr+1)
  chr.end=c(new.chr,nr.li)

  lsn.chr.index=cbind.data.frame(
    lsn.type=hit.cnt$lsn.index$lsn.type[chr.start],
    chrom=hit.cnt$lsn.index$chrom[chr.start],
    row.start=chr.start,
    row.end=chr.end)

  b=nrow(gene.lsn.chr.index)
  g=nrow(hit.cnt$nhit.mtx)
  nlt=ncol(hit.cnt$nhit.mtx)

  p.nsubj=p.nhit=matrix(1,g,nlt)
  colnames(p.nsubj)=colnames(p.nhit)=colnames(hit.cnt$nhit.mtx)

  # Mark unavailable exon-level probabilities as missing

  if (!is.null(exon_level))
  {
    exon.columns=intersect(exon_level,colnames(p.nhit))
    missing.exon.genes=which(is.na(gene.exon.size) | gene.exon.size<=0)

    if (length(exon.columns)>0 && length(missing.exon.genes)>0)
    {
      p.nhit[missing.exon.genes,exon.columns]=NA_real_
      p.nsubj[missing.exon.genes,exon.columns]=NA_real_
    }
  }

  ###################
  # Compute lesion enrichment probabilities
  ###################

  for (i in seq_len(b))
  {
    # find rows for affected genes
    gli.start.row=gene.lsn.chr.index$row.start[i]
    gli.end.row=gene.lsn.chr.index$row.end[i]
    gld.start.row=gene.lsn.index$row.start[gli.start.row]
    gld.end.row=gene.lsn.index$row.end[gli.end.row]
    gld.rows=gld.start.row:gld.end.row

    gene.rows=unique(hit.cnt$gene.lsn.data$gene.row[gld.rows])
    n.genes=length(gene.rows)

    # find rows for lesions of this type on this chromosome
    lsn.chr.mtch=which((lsn.chr.index$lsn.type==gene.lsn.chr.index$lsn.type[i])&
                         (lsn.chr.index$chrom==gene.lsn.chr.index$chrom[i]))

    if (length(lsn.chr.mtch)>0)
    {
      lsn.index.start.row=lsn.chr.index$row.start[lsn.chr.mtch]
      lsn.index.end.row=lsn.chr.index$row.end[lsn.chr.mtch]
      lsn.start.row=hit.cnt$lsn.index$row.start[lsn.index.start.row]
      lsn.end.row=hit.cnt$lsn.index$row.end[lsn.index.end.row]
      lsn.rows=lsn.start.row:lsn.end.row
      n.lsns=length(lsn.rows)

      lsn.type=hit.cnt$lsn.data$lsn.type[lsn.start.row]
      chrom=gene.lsn.chr.index$chrom[i]

      message(paste0("Computing p-values for ",
                     n.genes," gene(s) on chromosome ",
                     chrom," affected by ",
                     n.lsns," ",lsn.type,
                     " (data block ",i," of ",b,"): ",date()))

      # obtain chromosome and gene target sizes
      if (!is.null(exon_level) && lsn.type %in% exon_level)
      {
        chrom.size=exon.chrom.size.vec[as.character(chrom)]
        gene.size=gene.exon.size[gene.rows]

        if (length(chrom.size)!=1 || is.na(chrom.size) || chrom.size<=0)
          stop("No valid exon target size was found for chromosome '",
               chrom,"' and lesion type '",lsn.type,"'.",call.=FALSE)

        valid.exon.genes=!is.na(gene.size) & gene.size>0
        gene.rows=gene.rows[valid.exon.genes]
        gene.size=gene.size[valid.exon.genes]
        n.genes=length(gene.rows)

        if (n.genes==0)
          next
      } else {
        chr.mtch=which(chrom==chr.size$chrom)
        chrom.size=chr.size$size[chr.mtch]

        gene.size=hit.cnt$gene.data$loc.end[gene.rows]-
          hit.cnt$gene.data$loc.start[gene.rows]+1

        if (length(chrom.size)!=1 || is.na(chrom.size) || chrom.size<=0)
          stop("No valid chromosome size was found for chromosome '",
               chrom,"'.",call.=FALSE)
      }

      # obtain lesion sizes
      lsn.size=hit.cnt$lsn.data$loc.end[lsn.rows]-
        hit.cnt$lsn.data$loc.start[lsn.rows]+1

      # compute gene-lesion overlap probabilities
      log.pr=log(rep(lsn.size,each=n.genes)+
                   rep(gene.size,times=n.lsns))-log(chrom.size)

      pr.gene.hit=matrix(exp(log.pr),n.genes,n.lsns)
      pr.gene.hit[pr.gene.hit>1]=1

      # compute subject-level probabilities
      lsn.subj.IDs=hit.cnt$lsn.data$ID[lsn.rows]
      pr.subj=row.prob.subj.hit(pr.gene.hit,lsn.subj.IDs)

      # compute hit- and subject-level p-values
      for (j in seq_len(n.genes))
      {
        nsubj=hit.cnt$nsubj.mtx[gene.rows[j],lsn.type]
        nhit=hit.cnt$nhit.mtx[gene.rows[j],lsn.type]

        p.nsubj[gene.rows[j],lsn.type]=rpbc(nsubj,pr.subj[j,])
        p.nhit[gene.rows[j],lsn.type]=rpbc(nhit,pr.gene.hit[j,])
      }
    }
  }

  rownames(p.nhit)=rownames(hit.cnt$nhit.mtx)
  rownames(p.nsubj)=rownames(hit.cnt$nsubj.mtx)

  colnames(hit.cnt$nhit.mtx)=paste0("nhit.",colnames(hit.cnt$nhit.mtx))
  colnames(hit.cnt$nsubj.mtx)=paste0("nsubj.",colnames(hit.cnt$nsubj.mtx))
  colnames(p.nhit)=paste0("p.",colnames(hit.cnt$nhit.mtx))
  colnames(p.nsubj)=paste0("p.",colnames(hit.cnt$nsubj.mtx))

  ###################
  # Compute q-values
  ###################

  message(paste0("Computing q-values: ",date()))

  q.nhit=matrix(NA_real_,nrow=nrow(p.nhit),ncol=ncol(p.nhit),
                dimnames=dimnames(p.nhit))
  q.nsubj=matrix(NA_real_,nrow=nrow(p.nsubj),ncol=ncol(p.nsubj),
                 dimnames=dimnames(p.nsubj))

  for (i in seq_len(ncol(q.nhit)))
  {
    valid.p=!is.na(p.nhit[,i])

    if (any(valid.p))
    {
      pi.hat=min(1,2*mean(p.nhit[valid.p,i]))
      q.nhit[valid.p,i]=pi.hat*stats::p.adjust(p.nhit[valid.p,i],method="fdr")
    }

    valid.p=!is.na(p.nsubj[,i])

    if (any(valid.p))
    {
      pi.hat=min(1,2*mean(p.nsubj[valid.p,i]))
      q.nsubj[valid.p,i]=pi.hat*stats::p.adjust(p.nsubj[valid.p,i],method="fdr")
    }
  }

  colnames(q.nhit)=paste0("q.",colnames(hit.cnt$nhit.mtx))
  colnames(q.nsubj)=paste0("q.",colnames(hit.cnt$nsubj.mtx))

  # Remove internal pointers and place row indices first

  gd.clms=setdiff(colnames(hit.cnt$gene.data),c("glp.row.start","glp.row.end"))
  lsn.clms=setdiff(colnames(hit.cnt$lsn.data),c("glp.row.start","glp.row.end"))

  gd.clms=c("gene.row",setdiff(gd.clms,"gene.row"))
  lsn.clms=c("lsn.row",setdiff(lsn.clms,"lsn.row"))

  ###################
  # Compute constellation statistics
  ###################

  if (length(num.lsn)>1)
  {
    message(paste0("Computing p-values for number of lesion types affecting genes: ",
                   date()))

    # compute ordered p-values only for genes with complete lesion-specific probabilities
    complete.nhit=stats::complete.cases(p.nhit)
    complete.nsubj=stats::complete.cases(p.nsubj)

    p.ord.nhit=matrix(NA_real_,nrow=nrow(p.nhit),ncol=ncol(p.nhit))
    p.ord.nsubj=matrix(NA_real_,nrow=nrow(p.nsubj),ncol=ncol(p.nsubj))

    if (any(complete.nhit))
      p.ord.nhit[complete.nhit,]=p.order(p.nhit[complete.nhit,,drop=FALSE])

    if (any(complete.nsubj))
      p.ord.nsubj[complete.nsubj,]=p.order(p.nsubj[complete.nsubj,,drop=FALSE])

    colnames(p.ord.nhit)=paste0("p",seq_len(ncol(p.nhit)),".nhit")
    colnames(p.ord.nsubj)=paste0("p",seq_len(ncol(p.nsubj)),".nsubj")

    # compute q-values of ordered p-values
    q.ord.nhit=matrix(NA_real_,nrow=nrow(p.ord.nhit),ncol=ncol(p.ord.nhit),
                      dimnames=dimnames(p.ord.nhit))
    q.ord.nsubj=matrix(NA_real_,nrow=nrow(p.ord.nsubj),ncol=ncol(p.ord.nsubj),
                       dimnames=dimnames(p.ord.nsubj))

    message(paste0("Computing q-values for number of lesion types affecting genes: ",
                   date()))

    for (i in seq_len(ncol(p.ord.nhit)))
    {
      valid.p=!is.na(p.ord.nhit[,i])

      if (any(valid.p))
      {
        pi.hat=min(1,2*mean(p.ord.nhit[valid.p,i]))
        q.ord.nhit[valid.p,i]=pi.hat*
          stats::p.adjust(p.ord.nhit[valid.p,i],method="fdr")
      }

      valid.p=!is.na(p.ord.nsubj[,i])

      if (any(valid.p))
      {
        pi.hat=min(1,2*mean(p.ord.nsubj[valid.p,i]))
        q.ord.nsubj[valid.p,i]=pi.hat*
          stats::p.adjust(p.ord.nsubj[valid.p,i],method="fdr")
      }
    }

    colnames(q.ord.nsubj)=paste0("q",seq_len(ncol(p.nsubj)),".nsubj")
    colnames(q.ord.nhit)=paste0("q",seq_len(ncol(p.nhit)),".nhit")

    gene.res=cbind.data.frame(hit.cnt$gene.data[,gd.clms],
                              hit.cnt$nsubj.mtx,
                              p.nsubj,
                              q.nsubj,
                              p.ord.nsubj,
                              q.ord.nsubj,
                              hit.cnt$nhit.mtx,
                              p.nhit,
                              q.nhit,
                              p.ord.nhit,
                              q.ord.nhit)
  }

  # if lesion data has only one lesion type, skip the constellation analysis

  else if (length(num.lsn)==1)
  {
    gene.res=cbind.data.frame(hit.cnt$gene.data[,gd.clms],
                              hit.cnt$nsubj.mtx,
                              p.nsubj,
                              q.nsubj,
                              hit.cnt$nhit.mtx,
                              p.nhit,
                              q.nhit)
  }

  ###################
  # Assemble results
  ###################

  res=list(
    gene.hits=gene.res,                    # GRIN statistical results
    lsn.data=hit.cnt$lsn.data[,lsn.clms], # processed lesion data
    gene.data=hit.cnt$gene.data[,gd.clms],# processed gene annotation data
    gene.lsn.data=hit.cnt$gene.lsn.data,  # gene-lesion overlap records
    chr.size=chr.size,                     # chromosome sizes for standard analysis
    gene.index=hit.cnt$gene.index,         # chromosome index for gene data
    lsn.index=hit.cnt$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)
}