## ----setup, include = FALSE, warning = FALSE---------------------------------- knitr::opts_chunk$set(comment = FALSE, message = FALSE) ## ----------------------------------------------------------------------------- library(cellmig) library(ggplot2) library(ggforce) library(patchwork) library(rstan) ggplot2::theme_set(new = theme_bw(base_size = 10)) ## ----------------------------------------------------------------------------- # set seed for reproducible results set.seed(seed = 1253) # simulate data y_p <- gen_partial(control = list(N_biorep = 8, N_techrep = 5, N_cell = 50, delta = c(-0.4, -0.2, -0.1, 0, 0.1, 0.2, 0.4), sigma_bio = 0.1, sigma_tech = 0.05, offset = 4, prior_alpha_p_M = -0.5, prior_alpha_p_SD = 0.1, prior_kappa_mu_M = 1.5, prior_kappa_mu_SD = 0.1, prior_kappa_sigma_M = 0, prior_kappa_sigma_SD = 0.1)) ## ----------------------------------------------------------------------------- # see data structure str(y_p$y) ## ----fig.width=7, fig.height=3------------------------------------------------ ggplot(data = y_p$y)+ geom_sina(aes(x = paste0("t=", group), col = paste0("p=", plate), y = v, group = well), size = 0.5)+ theme_bw()+ theme(legend.position = "top", strip.text.x = element_text(margin = margin(0.03,0,0.03,0, "cm")))+ xlab(label = "dose")+ ylab(label = "migration speed")+ guides(color = guide_legend(override.aes = list(size = 3)))+ scale_color_grey(name = "plate")+ scale_y_log10() ## ----fig.width=7, fig.height=3------------------------------------------------ ggplot(data = aggregate(v~well+group+plate, data = y_p$y, FUN = mean))+ geom_sina(aes(x = paste0("t=", group), col = paste0("p=", plate), y = v, group = well), size = 1)+ theme_bw()+ theme(legend.position = "top", strip.text.x = element_text(margin = margin(0.03,0,0.03,0, "cm")))+ xlab(label = "dose")+ ylab(label = "migration speed")+ guides(color = guide_legend(override.aes = list(size = 3)))+ scale_y_log10()+ scale_color_grey(name = "plate") ## ----fig.width=7, fig.height=3.5---------------------------------------------- sim_data <- y_p$y # format simulated data to be used as input in cellmig sim_data$well <- as.character(sim_data$well) sim_data$compound <- as.character(sim_data$compound) sim_data$plate <- as.character(sim_data$plate) sim_data$offset <- 0 sim_data$offset[sim_data$group=="4"] <- 1 ## ----fig.width=7, fig.height=3.5---------------------------------------------- osd <- cellmig(x = sim_data, control = list(mcmc_warmup = 250, mcmc_steps = 800, mcmc_chains = 2, mcmc_cores = 2, mcmc_algorithm = "NUTS", adapt_delta = 0.8, max_treedepth = 10)) ## ----fig.width=4, fig.height=2------------------------------------------------ ggplot(data = osd$posteriors$alpha_p)+ geom_point(aes(y = plate_id, x = exp(mean)))+ geom_errorbarh(aes(y = plate_id, x = exp(mean), xmin = exp(X2.5.), xmax = exp(X97.5.)), height = 0.1)+ geom_point(data = data.frame(v = exp(y_p$par$alpha_p), plate = osd$posteriors$alpha_p$plate_id), aes(y = plate, x = v), col = "red")+ scale_y_continuous(breaks = osd$posteriors$alpha_p$plate_id)+ xlab(label = expression(alpha[p]*"'")) ## ----fig.width=4, fig.height=2.5---------------------------------------------- ggplot(data = osd$posteriors$delta_t)+ geom_point(aes(y = group_id, x = mean))+ geom_errorbarh(aes(y = group_id, x = mean, xmin = X2.5., xmax = X97.5.), height = 0.1)+ geom_point(data = data.frame(delta = c(-0.4, -0.2, -0.1, 0.1, 0.2, 0.4), group_id = 1:6), aes(y = group_id, x = delta), col = "red")+ xlab(label = expression(delta[t]))+ ylab(label = "treatment")+ scale_y_continuous(breaks = 1:8) ## ----fig.width=4, fig.height=2------------------------------------------------ plot(osd$fit, par = c("sigma_bio", "sigma_tech")) ## ----------------------------------------------------------------------------- control <- list(N_biorep = 4, N_techrep = 3, N_cell = 30, N_group = 8, prior_alpha_p_M = -0.5, prior_alpha_p_SD = 1.0, prior_kappa_mu_M = 1.5, prior_kappa_mu_SD = 1.0, prior_kappa_sigma_M = 0, prior_kappa_sigma_SD = 1.0, prior_sigma_bio_M = 0.0, prior_sigma_bio_SD = 1.0, prior_sigma_tech_M = 0.0, prior_sigma_tech_SD = 1.0, prior_delta_t_M = 0.0, prior_delta_t_SD = 1.0) ## ----------------------------------------------------------------------------- # generate data y_f <- gen_full(control = control) ## ----------------------------------------------------------------------------- str(y_f) ## ----fig.width=4, fig.height=4------------------------------------------------ w <- data.frame(v_f = y_f$y$v[1:2000], v_p = y_p$y$v[1:2000]) ggplot(data = w)+ geom_point(aes(x = v_f, y = v_p), size = 0.5)+ geom_density_2d(aes(x = v_f, y = v_p), col = "orange")+ scale_x_log10(name = "Speed from fully generative model", limits = c(0.01, 10^4))+ scale_y_log10(name = "Speed from partially generative model", limits = c(0.01, 10^4))+ annotation_logticks(base = 10, sides = "lb")+ theme_bw(base_size = 10) ## ----eval = FALSE, message=FALSE, warning=FALSE------------------------------- # # Constants # N_bioreps <- c(3, 6, 9) # N_sim <- 10 # true_deltas <- c(-0.3, -0.15, 0, 0.2, 0.4) # offset <- 3 # # # power analysis # deltas <- vector(mode = "list", length = length(N_bioreps)*N_sim) # i <- 1 # for(N_biorep in N_bioreps) { # for(b in 1:N_sim) { # # simulate data # y_p <- gen_partial(control = list(N_biorep = N_biorep, # N_techrep = 3, # N_cell = 40, # delta = true_deltas, # sigma_bio = 0.1, # sigma_tech = 0.05, # offset = offset, # prior_alpha_p_M = -0.5, # prior_alpha_p_SD = 0.1, # prior_kappa_mu_M = 1.5, # prior_kappa_mu_SD = 0.1, # prior_kappa_sigma_M = 0, # prior_kappa_sigma_SD = 0.1)) # # # format simulated data to be used as input in cellmig # sim_data <- y_p$y # sim_data$well <- as.character(sim_data$well) # sim_data$compound <- as.character(sim_data$compound) # sim_data$plate <- as.character(sim_data$plate) # sim_data$offset <- 0 # sim_data$offset[sim_data$group==offset] <- 1 # # # run cellmig # o <- cellmig(x = sim_data, # control = list(mcmc_warmup = 300, # mcmc_steps = 1000, # mcmc_chains = 1, # mcmc_cores = 1, # mcmc_algorithm = "NUTS", # adapt_delta = 0.8, # max_treedepth = 10)) # # delta <- o$posteriors$delta_t # delta$b <- b # delta$N_biorep <- N_biorep # delta$true_deltas <- true_deltas[-offset] # delta$TP <- (delta$X2.5. <= delta$true_deltas & # delta$X97.5. >= delta$true_deltas) & # !(delta$X2.5. <= 0 & delta$X97.5. >= 0) # deltas[[i]] <- delta # i <- i + 1 # } # } # # rm(i, delta, o, sim_data, y_p, b, N_biorep, offset, true_deltas) # # deltas <- do.call(rbind, deltas) ## ----eval = FALSE, fig.height=3, fig.width = 6-------------------------------- # ggplot(data = aggregate(TP~N_biorep+true_deltas, data = deltas, FUN = sum))+ # geom_point(aes(x = N_biorep, y = TP, col = abs(true_deltas), # group = as.factor(true_deltas)), size = 2, alpha = 0.5)+ # geom_line(aes(x = N_biorep, y = TP, col = abs(true_deltas), # group = as.factor(true_deltas)), alpha = 0.5)+ # ylab(label = "Number of true positives (TPs)")+ # scale_color_distiller(name = expression("|"*delta[t]*"|"),palette="Spectral")+ # theme_bw(base_size = 10) ## ----------------------------------------------------------------------------- sessionInfo()