simulate_read_count: Simulate read count (single SNP model)

In liangyy/mixqtl: Combine total and allele-specific counts for cis-QTL analysis

Description

Given gene configuration, genotype (single SNP), genetic effect size (log aFC), and other parameters, simulate read count.

Usage

simulate_read_count(gene, genotype, beta, L_read, y_dist)

Arguments

gene	gene instance generated from create_gene
genotype	genotype generated from create_genotype (the 'genotype' object in the list)
beta	log(aFC)
L_read	length of read
y_dist	distribution of Y^h | library size, relative abundance. It specifies the shape of the distribution whereas the mean is given by library size and relative abundance.

Value

read counts, where observed count include y1, y2 (AS count of each haplotypes) and ystar (total count - y1 - y2) along with library size Ti_lib. Also, it includes unobserved count as 'hidden' (y1star and y2star).

About sampling Y^h

To specify read count of each haplotype (Y^h) given library size and relative abundance. It is set in y_dist. Currently it supports three distribution types: 'poisson', 'lognormal', and 'negbinom'. Specifically, y_dist is a list with distribution name in 'type'. For y_dist$type = 'poisson', no other parameter is required. For y_dist$type = 'lognormal', 'sigma' needs to specify, and Y^h = round(library_size * relative_abundance * rlnorm(1, 0, y_dist$sigma)) For y_dist$type = 'negbinom', 'prob' and 'size_factor' need to specify, and Y^h = rnbinom(1, size = y_dist$size_factor * library_size * relative_abundance, prob = y_dist$prob)

Examples

gene = create_gene()
genotype = create_genotype(
  maf = c(0.05, 0.3),
  nsample = 300,
  nreplicate = 1000
)
simulate_read_count(
  gene = gene,
  genotype = genotype[[1]]$genotype,
  beta = 0.2,
  L_read = 75,
  y_dist = list(
    type = 'negbinom',
    prob = 2/3,
    size_factor = 2
  )
)

liangyy/mixqtl documentation built on Sept. 17, 2020, 11:36 a.m.