estimation_fun: Estimation for MAP data set:

In meca7653/MAPTest: MAP Test for detection DE genes

Description

Estimation for MAP data set:

Usage

estimation_fun(n_control = 10, n_treat = 10, n_rep = 3, x, Y1,
  k = 4, nn = 300, phi = NULL, type = NULL, tttt)

Arguments

n_control	Number of time points in control group
n_treat	Number of time points in treatment group
n_rep	Number of replicates in both group
x	Time structured design for the simulated data
Y1	RNA-seq data set: G by (n_control + n_treat) * n_rep matrix
k	Always 4
nn	Number of QMC nodes used in likelihood estimation
phi	the dispersion parameter
type	dispersion parameter estimation type
tttt	Time points used in the design

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)] = λ_{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. We used latented negative binomial model with EM algorithm to estimate the paramters of mixture model.

Value

A list of parameters that we need to use for DE analysis.

data_use List includes the data information

• Y1 RNA-seq data set: G by (n_control + n_treat) x n_rep matrix

• start Starting point of the EM algorithm

• k always 4

• n_basis Number of basis function have been used to construct the design matrix

• X1 Vector indicate control data (0) or treatment data (1)

• x Design matrix been used for estimation

• tttt Time points used in the design
  

result1 List includes the Parameters estimations

• mu1 Mean parameter for η_{g1}

• sigma1_2 Variance parameter for η_{g1}

• sigma2 Variance parameter for τ_{g}

• sigma2_2 Variance parameter for η_{g2}

• p_k Proportion for each mixture component

• aa Dispertion parameter

Examples

library(matlib)
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
#------------gaussian basis construction------------------
# method <- "gaussian"
# n_basis <- 2
# a = 1/4
# b = 2
# c = sqrt(a^2 + 2 * a * b)
# A = a + b + c
# B = b/A
# phi_cal = function(k, x){
#   Lambda_k =sqrt(2 * a / A) * B^k
#   1/(sqrt(sqrt(a/c) * 2^k * gamma(k+1))) *
#      exp(-(c-a) * x^2) * hermite(sqrt(2 * c) * x, k, prob = F)
#
# }
# z = do.call(cbind, lapply(c(1:n_basis), phi_cal, x = (tttt - mean(tttt))))
# x = matrix(0, length(tttt), n_basis)
# if(n_basis == 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])
#
# }else{
#   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])
#   x[,3] = z[,3] - Proj(z[,3], rep(1, length(tttt))) - Proj(z[,3], x[,1]) - Proj(z[,3], x[,2])
# }

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