HeritSeq (version 1.0.2)

getReadMatrix.CP: Simulate a read matrix from compound Poisson mixed effect models (CPMM).

Description

Simulate a (possibly unbalanced) read matrix from CPMM. For a compound Poisson (CP) random variable \(Y_{gsr}\) with mean \(\mu_{gs}\), its variance can be expressed as \(\phi_g\mu_{gs}^{p_g}\), for some \(1<p_g<2\). Under the CPMM, with a \(\log\)-link, the regression on the mean has the form: \(\log(\mu_{gs}) = \alpha_g+ b_{gs}, \;\;b_{gs}\sim N(0, \sigma^2_g).\)

Usage

getReadMatrix.CP(vec.num.rep, alphas, sigma2s, ps, phis)

Arguments

vec.num.rep

A vector of replicate numbers for each strain.

alphas

Intercept vector \(\alpha_g\)'s, \(1 \times \texttt{num.features}\).

sigma2s

Random effect variance vector \(\sigma^2_g\)'s, \(1 \times \texttt{num.features}\).

ps

Tweedie parameter in CP models, \(p_g\)'s, a \(1 \times \texttt{num.features}\) vector.

phis

Dispersion parameter in CP models, \(\phi_g\)'s, a \(1 \times \texttt{num.features}\) vector.

Value

A \(G \times N\) matrix with CP reads. \(N\) is the total number of samples; \(G\) is the number of features. Column names are sample names of the form "Ss_r", where S stands for sample, s is the strain number, r is the replicate number within the strain. Row names are the feature names of the form "Gene g", where g is the feature index.

Examples

# NOT RUN {
## Generate a sequencing dataset with 5 features and 6 strains. 
## Assign parameter values.
rep.num <- c(3, 5, 2, 3, 4, 2)
a0s <- c(-1, 1, 2, 5, 10)
sig2s <- c(10, 0.2, 0.1, 0.03, 0.01)
ps <- rep(1.5, 5)
phis <- c(1.5, 1, 0.5, 0.1, 0.1)

set.seed(1234)
## Generate reads:
cpData <- getReadMatrix.CP(rep.num, a0s, sig2s, ps, phis)
## Generate strain names:
str <- sapply(1:length(rep.num), function(x){
  str.x <- paste0("S", x)
  return(rep(str.x, rep.num[x]))
})
str <- do.call(c, str)
# }