source("bandleFUN.R")

load("data/davies2018wt.rda")
load("data/davies2018ap4b1.rda")
load("data/davies2018ap4e1.rda")

davies2018wt <- filterNA(davies2018wt)
davies2018ap4b1 <- filterNA(davies2018ap4b1)
davies2018ap4e1 <- filterNA(davies2018ap4e1)

exprs(davies2018wt) <- log2(exprs(davies2018wt))
exprs(davies2018ap4b1) <- log2(exprs(davies2018ap4b1))
exprs(davies2018ap4e1) <- log2(exprs(davies2018ap4e1))

par(mfrow = c(2,3))
plot2D(davies2018wt[,1:5], dims = c(1,2))
plot2D(davies2018wt[,6:10], dims = c(1,2))
plot2D(davies2018ap4b1, dims = c(1,2))
plot2D(davies2018ap4e1[,1:5], dims = c(1,2))
getMarkerClasses(davies2018wt, fcol = "markers")

par(mfrow = c(1,1))

exprs(davies2018wt)

K <- length(getMarkerClasses(davies2018wt))
dir_prior = diag(rep(1, K)) + matrix(0.02, nrow = K, ncol = K)
prior_pred_dir(object = davies2018wt, dir_prior = dir_prior)

pc_prior <- matrix(NA, ncol = 3, K)
pc_prior[seq.int(1:K), ] <- matrix(rep(c(0.05, 60, 100), each = K), ncol = 3)

par(mfrow = c(4,4))
bandleres_daviese1 <- bandle(object_cond1 = c(davies2018wt[, 1:5], davies2018wt[, 6:10]),
                             object_cond2 = c(davies2018ap4e1[, 1:5], davies2018ap4e1[, 6:10]),
                             numIter = 20000, burnin = 10000, thin = 15, pg = FALSE, hypLearn = "LFBGS", 
                             maternCov = TRUE, nu = 2, pc_prior = pc_prior, numChains = 6, dir_prior = dir_prior)

save(bandleres_daviese1, file = "bandleres_daviese1.rda")
