Simulate molecular count data

This function generates a matrix of negative-binomial count data using the nonlinear GxT parametrization consistent with the models fit by map_qtl when fn = "nonlinear". Each feature is associated with its own SNP. When geno is NULL, a separate genotype vector is simulated per feature.

make_count_data(
  anno,
  n.feat,
  n.sample,
  geno = NULL,
  sd = c(0, 0, 0),
  coef = rep(NA, 3),
  lb.maf = 0.05,
  ub.maf = 0.5,
  filter.geno = FALSE,
  intercept.mean = 4,
  intercept.sd = 2,
  disp.mean.rel = function(x) 4/x + 0.1,
  size.factors = rep(1, n.sample)
)

Arguments

anno

A data frame containing sample IDs (character strings or integers), subject IDs (character strings or integers), and treatment conditions (0 or 1). The columns must be named as "sample", "subject", and "condition".

n.feat

An integer specifying the number of features (e.g., genes and ATAC-seq peaks).

n.sample

An integer specifying the number of samples.

geno

An n.feat x n.sub matrix of genotype values (0, 1, 2). Column names must match unique(anno$subject). Row i gives the subject-level genotypes for feature i. When NULL (the default), a separate genotype vector is simulated for each feature using lb.maf and ub.maf.

sd

A numeric vector of length three specifying the standard deviations of the genotype, treatment, and interaction effects, respectively. Each effect is drawn independently per feature from Normal(0, sd). Default c(0, 0, 0) gives no effects. For each effect, exactly one of the corresponding elements of sd and coef must be non-NA.

coef

A numeric vector of length three specifying fixed values for the genotype, treatment, and interaction effects, respectively. For each effect, exactly one of the corresponding elements of sd and coef must be non-NA.

lb.maf

A scalar specifying the lower bound of the MAF used when simulating genotypes (geno = NULL).

ub.maf

A scalar specifying the upper bound of the MAF used when simulating genotypes (geno = NULL).

filter.geno

A logical scalar. When TRUE and geno = NULL, genotype simulation for each feature is repeated until all three levels (0, 1, 2) are present.

intercept.mean

A scalar specifying the mean of the per-feature intercept b0 drawn from Normal(intercept.mean, intercept.sd).

intercept.sd

A scalar specifying the standard deviation of the per-feature intercept.

disp.mean.rel

A function that maps exp(b0) (the baseline mean count) to the negative-binomial dispersion (size = 1 / dispersion in rnbinom).

size.factors

A numeric vector of length n.sample specifying multiplicative size factors applied to the count-scale mean.

Value

A matrix of non-negative integers with n.feat rows and n.sample columns. Row names are "f1", "f2", ... and column names are the sample IDs in anno$sample.

Examples

anno <- data.frame(
    sample = paste0("s", 1:20),
    subject = rep(paste0("d", 1:10), each = 2),
    condition = rep(c(0, 1), 10)
)
counts <- make_count_data(
    anno = anno, n.feat = 5, n.sample = 20,
    sd = c(0, 0, 0)
)