R/langmuirMethArrays.R

In stephaniehicks/quantroSim: Supporting Data Simulation R-package for quantro

 #' @title Langmuir parameters for DNA Methylation Arrays
#'
#' @description These are the parameters used to simulate the DNA Methylation
#' arrays. The parameters are used in the Langmuir adsorption model.  
#'
#' @param nProbes number of probes
#' @param nSamps number of samples in each group
#' @param nGroups number of groups
#' @param mua hyperparameter: mean of log normal distribution for a 
#' @param siga hyperparameter: variance of log normal distribution for a 
#' @param mub hyperparameter: mean of log normal distribution for b
#' @param sigb hyperparameter: variance of log normal distribution for b
#' @param muOpt hyperparameter: mean of log normal distribution for optical 
#' noise
#' @param sigOpt hyperparameter: variance of log normal distribution for 
#' optical noise
#' @param muBG hyperparameter: mean of log normal distribution for d
#' @param sigBG hyperparameter: variance of log normal distribution for d
#' @param muERR hyperparameter: mean of log normal distribution for 
#' measurement error
#' @param sigERR hyperparameter: variance of log normal distribution for 
#' measurement error
#'
#' @return A list with elements
#' \item{gID}{Group number ID}
#' \item{a}{Simulated a values for the Langmuir Adsorption model}
#' \item{b}{Simulated b values for the Langmuir Adsorption model}
#' \item{d}{Simulated d values for the Langmuir Adsorption model}
#' \item{opt}{Simulated optical noise}
#' \item{errMeth}{Simulated measurement error for methylation values}
#' \item{errUnmeth}{Simulated measurement error for unmethylation values}
#' 
#' @import MASS
#' 
#' @author Stephanie Hicks
#' @export
#' @examples
#' langmuirMethArrays(nProbes = 10000, nSamps = 4, nGroups = 2)
langmuirMethArrays <- function(nProbes, nSamps, nGroups,
            mua = NULL, siga = NULL, mub = NULL, sigb = NULL, 
            muOpt = NULL, sigOpt = NULL, muBG = NULL, sigBG = NULL, 
            muERR = NULL, sigERR = NULL)
{
	totSamps <- nSamps * nGroups

	# Default parameters
	if(is.null(mua)){ mua = rep(16, totSamps) }
  if(is.null(siga)){ siga = 0.1 * diag(totSamps) }
  
	if(is.null(mub)){ mub = rep(22, totSamps) }
  if(is.null(sigb)){ sigb = 0.1 * diag(totSamps) }
  
	if(is.null(muOpt)){ muOpt = rep(5, totSamps) }
	if(is.null(sigOpt)){ sigOpt = 1 * diag(totSamps) }
  
  if(is.null(muBG)){ muBG = rep(5, totSamps)}
	if(is.null(sigBG)){ sigBG = 1 * diag(totSamps) }	
	
  if(is.null(muERR)){ muERR = rep(0, totSamps) }
	if(is.null(sigERR)){ sigERR = 1 * diag(totSamps) }	
	
	# Sample parameters: a, b, d, err
  a <- sweep(2^(mvrnorm(nProbes, mu = rep(0, totSamps), 
                          Sigma = 0.1 * diag(totSamps))), 2,
             2^(mvrnorm(n = 1, mu = mua, Sigma = siga)), FUN = "*") 
  b <- sweep(2^(mvrnorm(nProbes, mu = rep(0, totSamps), 
                        Sigma = 0.1 * diag(totSamps))), 2,
             2^(mvrnorm(n = 1, mu = mub, Sigma = sigb)), FUN = "*")
  opt <- sweep(2^(mvrnorm(n = nProbes, mu = rep(0, totSamps), 
                          Sigma = 0.1 * diag(totSamps))), 2, 
               2^(mvrnorm(n = 1, mu = muOpt, Sigma = sigOpt)), FUN="*")
	d <- 2^(mvrnorm(n = nProbes, mu = muBG, Sigma = sigBG))
	errMeth <- 2^(mvrnorm(n = nProbes, mu = muERR, Sigma = sigERR))
	errUnmeth <- 2^(mvrnorm(n = nProbes, mu = muERR, Sigma = sigERR))

	list("gID" = rep(1:nGroups, each = nSamps), 
	     "a" = a, "b" = b, "opt" = opt, "d" = d, 
         "errMeth" = errMeth, "errUnmeth" = errUnmeth)			 
}