R/mutation.proba.r

In IOstrovnaya/Clonality: Clonality testing

# mutation.proba: Probability of being clonal

#' Probability of being clonal
#'
#' @description This function uses the results from mutation.rem to estimate the diagnostic probability of clonal relatedness for new cases. 
#' It is obtained from Bayes theorem using the prior probability of clonal relatedness (pi) and the contributions to the 
#' likelihood based on the mutations observed for the case. We recommand to use this function to estimate probabilities of clonality for 
#' new subjects, ie who are not used for the model estimation. To obtain estimate for the subjects on which the model estimation is based, 
#' the option "proba=TRUE" can be used in the mutation.rem function.
#' 
#' @usage mutation.proba(para, likmat, xigrid = c(0, seq(0.0005, 0.9995, by=0.001)))
#' @param para Value of the model parameters, in the form c(mu, sigma, pi). 
#' @param likmat Grid of conditional probabilities for each tumor pair (rows) and each value of xi (columns). 
#' This matrix is generated by the auxiliary function grid.lik, and returned as a parameter by the principal function
#' mutation.rem.
#' @param xigrid Grid of the values of xi, corresponding to its domain of definition. The default is c(0, seq(0.0005, 0.9995, by=0.001)).
#' @return Returns the vectors of probability of clonality for each pairs of tumors contained in the matrix likmat 
#' (the number of pairs is the number of rows of the matrix).
#' @author Audrey Mauguen \email{mauguena@mskcc.org} and Venkatraman E. Seshan.
#' @references Mauguen A, Seshan VE, Ostrovnaya I, Begg CB. Estimating the Probability of Clonal Relatedness of 
#' Pairs of Tumors in Cancer Patients. Submitted.
#' @examples
#' #___ Analysis of LCIS data
#' data(lcis)
#' 
#' #__ Parameters estimation
#' mod <- mutation.rem(lcis)
#' mod
#' 
#' #__ Probability of being clonal for a new subject
#' # generate a case with 30 mutations
#' # probabilities of each observed mutation
#' pi <- runif(30,0.001,0.13)
#' # mutation 1=shared or 2=private
#' newpair <- cbind(pi,rbinom(30,1,1-pi^2)+1)
#' # generate the matrix of likelihood values
#' new.likmat <- grid.lik(xigrid=c(0, seq(0.0005, 0.9995, by=0.001)), as.matrix(newpair[,c(-1)]), newpair[,1])
#' # probability of being clonal using the model previoulsy estimated
#' proba <- mutation.proba(c(mod$mu, mod$sigma, mod$pi), t(as.matrix(new.likmat)) )
#' @export

mutation.proba <- function(para, likmat, xigrid = c(0, seq(0.0005, 0.9995, by=0.001))){
  pclon <- para[3] * rowSums(likmat[,-1] %*% xidens(para[1],para[2], xigrid)[-1] ) * 0.001
  pnclon <- (1-para[3]) * likmat[,1]
  out <- pclon / (pclon+pnclon)
  class(out) <- "mutation.proba"
  out
}