R/row.prob.subj.hit.R

In GRIN2: Genomic Random Interval (GRIN)

## Internal helper function used by prob.hits().
##
## Computes the probability that each subject has at least one lesion affecting
## each gene. Rows of P represent genes, columns represent lesions, and IDs
## identifies the subject associated with each lesion column.
row.prob.subj.hit <- function(
    P,   # Matrix of lesion-hit probabilities; rows are genes and columns are lesions.
    IDs  # Subject identifier for each lesion; length must equal ncol(P).
) {
  if (length(IDs) != ncol(P)) {
    stop("length(IDs) must equal ncol(P).")
  }

  n.genes <- nrow(P)
  n.lesions <- ncol(P)

  # Arrange lesion columns so that lesions from the same subject are adjacent.
  ord <- order(IDs)
  IDs <- IDs[ord]
  P <- matrix(P[, ord], nrow = n.genes, ncol = n.lesions)

  # Identify the first and last lesion column for each subject.
  new.ID <- which(IDs[-1] != IDs[-n.lesions])
  ID.start <- c(1, new.ID + 1)
  ID.end <- c(new.ID, n.lesions)

  n.subjects <- length(ID.start)

  # Calculate gene-level hit probabilities for each subject.
  Pr <- matrix(NA_real_, nrow = n.genes, ncol = n.subjects)

  for (i in seq_len(n.subjects)) {
    pr.mtx <- matrix(
      P[, ID.start[i]:ID.end[i]],
      nrow = n.genes,
      ncol = ID.end[i] - ID.start[i] + 1
    )

    Pr[, i] <- 1 - exp(rowSums(log(1 - pr.mtx)))
  }

  return(Pr)
}