Summary_MAP: Summary for MAP Test

In meca7653/MAPTest: MAP Test for detection DE genes

Description

Summary for MAP Test

Usage

Summary_MAP(test_result, alpha = 0.05)

Arguments

test_result	MAP_test result: A list of results for DE analysis
alpha	Nominal level of FDR

Details

The vector of read counts for gene g, treatment group i, replicate j, at time point t,Y_{gij}(t), follows a Negative Binomial distribution parameterized mean λ_{gi} and φ_g, where E[Y_{gij}(t)] = lambda_{gi}(t). λ_{gi}(t) is further modeled as λ_{gi}(t) = S_{ij} exp[η_{g1}I_{i = 2} + B'(t)η_{g2}I_{i = 2} + B'(t)τ_{g}] We have B'(t) are design matrix, which is constructed by 2 orthorgonal polynomial bases.

• t = 1,..., n_treat (or n_control if control group);

• j = 1,..., n_rep;

• g = 1,...,G; and

• [η_{g1}, η_{g2}, τ_{g}] ~ 4-component gausssian mixture model. After the parameter estimation we run the MAP test to detect the DE genes. At given nominal level, a list of differential genes are returned at a given hypothesis.

Value

A list of results for DE analysis.

• Type Type of hypothesis

• Reject_index DE genes list

• FDR_hat Estimated FDR

Examples

library(matlib)
set.seed(1)
n_basis = 2
n_control = 10
n_treat   = 10
n_rep = 3
tt_treat  = c(1:n_treat)/n_treat
nt = length(tt_treat)
ind_t = sort(sample(c(1:nt), n_control))
tt = tt_treat[ind_t]
tttt = c(rep(tt, n_rep), rep(tt_treat, n_rep))
z = x = matrix(0, length(tttt), n_basis)
z[,1] = 1.224745*tttt
z[,2] = -0.7905694 + 2.371708*tttt^2
x[,1] = z[,1] - Proj(z[,1], rep(1, length(tttt)))
x[,2] = z[,2] - Proj(z[,2], rep(1, length(tttt))) - Proj(z[,2], x[,1])
Y1 = data_generation(G = 100,
                     n_control = n_control,
                     n_treat   = n_treat,
                     n_rep     = n_rep,
                     k_real = 4,
                     sigma2_r = rep(1, 2),
                     sigma1_2_r = 1,
                     sigma2_2_r = c(3,2),
                     mu1_r = 4,
                     phi_g_r = rep(1,100),
                     p_k_real = c(0.7, 0.1, 0.1, 0.1),
                     x = x)
aaa <- proc.time()
est_result <- estimation_fun(n_control = n_control,
                             n_treat = n_treat,
                             n_rep = n_rep,
                             x = x,
                             Y1 = Y1,
                             nn = 300,
                             k = 4,
                             phi = NULL,
                             type = 2,
                             tttt = tttt)
aaa1 <- proc.time()
aaa1 - aaa

G <- 100
k_real <- 4
p_k_real <- c(0.7, 0.1, 0.1, 0.1)
dd = rep(c(0:(k_real-1)), p_k_real * G)
result = MAP_test(est_result = est_result, Type = c(1:6), dd = dd, nn = 300)
ss <- Summary_MAP(result)

meca7653/MAPTest documentation built on June 20, 2019, 2 p.m.