sim.counts {ssizeRNA}

RNA-seq Count Data Simulation from Negative-Binomial Distribution

Description

This function simulates count data from Negative-Binomial distribution for two-sample RNA-seq experiments with given mean, dispersion and fold change. A count data matrix is generated.

Usage

sim.counts(nGenes = 10000, pi0 = 0.8, m, mu, disp, fc, up = 0.5,
  replace = TRUE)

Arguments

nGenes

total number of genes, the default value is 10000.

pi0

proportion of non-differentially expressed genes, the default value is 0.8.

m

sample size per treatment group.

mu

a vector (or scalar) of mean counts in control group from which to simulate.

disp

a vector (or scalar) of dispersion parameter from which to simulate.

fc

a vector (or scalar, or a function that takes an integer n and generates a vector of length n) of fold change for differentially expressed (DE) genes.

up

proportion of up-regulated genes among all DE genes, the default value is 0.5.

replace

sample with or without replacement from given parameters. See Details for more information.

Details

If the total number of genes nGenes is larger than length of mu or disp, replace always equals TRUE.

Value

counts

RNA-seq count data matrix.

group

treatment group vector.

lambda0

mean counts in control group for each gene.

phi0

dispersion parameter for each gene.

de

differentially expressed genes indicator: 0 for non-differentially expressed genes, 1 for up-regulated genes, -1 for down-regulated genes.

delta

log2 fold change for each gene between treatment group and control group.

Author(s)

Ran Bi [e-mail], Peng Liu [e-mail]

Examples

m <- 3                    ## sample size per treatment group
mu <- 10                  ## mean counts in control group for all genes 
disp <- 0.1               ## dispersion for all genes
fc <- 2                   ## 2-fold change for DE genes

sim <- sim.counts(m = m, mu = mu, disp = disp, fc = fc)
sim$counts                ## count data matrix


## varying fold change
fc1 <- function(x){exp(rnorm(x, log(2), 0.5*log(2)))}
sim1 <- sim.counts(m = m, mu = mu, disp = disp, fc = fc1)