R/pvalueBetaReg.R

In HengPan2007/MethSig: DNA Methylation Drivers Inference

#' @title Infer sample-specific hypermethylation
#' @description Expected DHcR is estimated through beta regression model and expected DHcR
#'     is tested against observed DHcR to evalute if there is promoter hypermethylation.
#' @param formula An object of class "formula" (or one that can be coerced to that class):
#'     a symbolic description of the model to be fitted. The name of y variable should be DHcR_Tumor_Beta.
#' @param data A data frame object generated by \emph{makeInputMatrix} containing following columns:
#'     Id, Hugo, Covariates included in \emph{matCV}, PDR_Tumor, Depth_Tumor, CpGs_Tumor and DHcR_Tumor.
#'     An example can be loaded by using \emph{invisible(inputMat)}.
#' @return A data frame summarizing information of sample-specific hypermethylation inference. Addtional columns
#'     are added to the input matrix.
#'     \itemize{
#'       \item DHcR_Tumor_Beta: Transformed data of DHcR_Tumor (from [0, 1] to (0, 1))
#'       \item Beta_Response: Response value estimated by beta regression
#'       \item Beta_Precision: Precision value estimated by beta regression
#'       \item Pvalue: P-value of sample-specific hypermethylation inference
#'     }
#' @examples
#' pvalueBetaReg(formula = 'DHcR_Tumor_Beta~DHcR_Normal+PDR_Normal+GEXP_Normal+Reptime+PDR_Tumor+Depth_Tumor+CpGs_Tumor',
#'               data = invisible(inputMat))
#'
pvalueBetaReg <- function(formula, data) {
  # data transformation from [0,1] to (0,1)
  scale.factor <- round(nrow(data)/length(unique(data$Id)))
  data$DHcR_Tumor_Beta <- (data$DHcR_Tumor*(scale.factor-1)+0.5)/scale.factor

  # fit beta regression model
  fit.beta <- betareg::betareg(formula, data=data)

  # estimate expected hypermethylation of tumor samples
  data$Beta_Response <- predict(fit.beta, data=data, type='response') # mu
  data$Beta_Precision <- predict(fit.beta, data=data, type='precision') # phi
  shape1 <- data$Beta_Response * data$Beta_Precision
  shape2 <- data$Beta_Precision - shape1
  ans <- data.frame(beta=data$DHcR_Tumor_Beta, shape1=shape1, shape2=shape2)

  # evaluate if observed DHcR is significantly higher than expected DHcR
  data$Pvalue <- apply(ans, 1, function(x) pbeta(as.numeric(x[1]), as.numeric(x[2]), as.numeric(x[3]), lower.tail=F))
  return(data)
}