R/02-prior.R

In DeepCNV: Exploiting Normal Contamination to Infer Copy Number from Deep Sequencing

# Copyright (C) Kevin R. Coombes 2014-2017

###################################
# Really need two alternate forms of the prior on the frsaction 'nu'
# of normal cells.  The default is to take a (theoretical) beta-prior
# with paramerters 'alpha' and 'beta'.  At some point, we have to
# "realize" this prior by instantiating it along a grid of x-values
# in the interval [0,1] and computing the corresponding frequencies
# for the PDF.  We might even just want to take the second representation
# when the 'prior' has already been updated and is no longer a simple
# beta distribution.
#
# Made our lives simple: we just pass in a "resolution" along with the
# constructor and instantiate immediately.  Note that we may need to
# add another method or change this design to handle the multi-gene case.
setClass("CNVPrior",
         representation=list(
           pGain = "numeric",
           pLoss = "numeric",
           alpha = "numeric",
           beta = "numeric",
           grid = "numeric",
           pdf = "numeric"
           ),
         prototype=list(
           pGain=0.2, pLoss=0.2,
           alpha=1, beta=1, grid=numeric(), pdf=numeric()))

# how do we validate the structure of the discrete prior?
setValidity("CNVPrior", function(object) {
  (length(object@alpha) == 1 &
     length(object@beta) == 1 &
     length(object@grid) == length(object@pdf) &
     length(object@pGain) == 1 & length(object@pLoss) == 1 &
     object@pGain >= 0 & object@pLoss >= 0 &
    (object@pGain + object@pLoss <= 1)
  )
})

setMethod("plot", c("CNVPrior", "missing"), function(x, lwd=2, ...) {
  if (length(x@grid) > 1) {
    xx <- x@grid
    yy <- exp(x@pdf)
  } else {
    xx <- seq(0, 1, length=301)
    yy <- dbeta(xx, x@alpha, x@beta)
  }
  plot(xx, yy, type="l", lwd=lwd,
       main="Prior Distribution on Normal Fraction",
       xlab="Normal Fraction", ylab="Frequency", ...)
  invisible(x)
})

setMethod("summary", c("CNVPrior"), function(object, ...) {
  cat("Discrete Prior: \n")
  cat("Normal:", 1-object@pGain-object@pLoss,
      "Gained:", object@pGain,
      "Deleted:", object@pLoss, "\n")
})

setCNVPrior <- function(alpha = 1, beta = 1,
                        pAbnormal=2/3, pGain = NULL, pLoss = pGain,
                        resn=300) {
  if (is.null(pGain)) {
    pGain <- pAbnormal/2
    pLoss <- pAbnormal/2
  }
  # avoid infinities at boundary by removing the 0 and 1 entries
  grid <- seq(0, 1, length=2+resn)[c(-1, -(2+resn))]
  nuprior <- dbeta(grid, alpha, beta, log=TRUE)
  new("CNVPrior",
      pGain=pGain, pLoss=pLoss,
      alpha = alpha, beta = beta, grid=grid, pdf=nuprior)
}

cnPrior <- function(object) {
  v <- c(object@pLoss, object@pGain, 1 - object@pLoss - object@pGain)
  names(v) <- cnSet
  v
}