diff --git a/src/counts_matrix_generator.R b/src/counts_matrix_generator.R index 5e903c2efffcc297203a9e3f58ef6ebc28c33655..2ca4417c437d4dcda82288d22f4f7723a257c3ff 100644 --- a/src/counts_matrix_generator.R +++ b/src/counts_matrix_generator.R @@ -1,3 +1,10 @@ +############################# PCKGE REQUIRED ############################## +library(DESeq2) +library(ggplot2) +library(tydiverse) +### maybe others ### + + # fix seed set.seed(123) @@ -6,50 +13,17 @@ set.seed(123) source("mydatalocal/counts_simulation/src/simulators.R") -# visualization functions -mu_effect_visualization <- function(mu_effect_res){ - label_wrap <- c("mu observed", "N gene DE", "min(|logFC|)", "var observed") - names(label_wrap) <- c("mu_observ", "res_DEA", "statistical_power", "var_observ") - figure = mu_effect_res %>% ggplot(., aes(x=vec_of_mu, y = value, col=factor(N_rep))) + - geom_point() + facet_wrap(~variable, scales = "free_y", labeller = labeller(variable = label_wrap)) + labs(color = "N replicates") - return(figure) -} - -size_effect_visualization <- function(alpha_effect_res){ - label_wrap <- c("mu observed", "N gene DE", "min(|logFC|)", "var observed") - names(label_wrap) <- c("mu_observ", "res_DEA", "statistical_power", "var_observ") - figure = alpha_effect_res %>% ggplot(., aes(x=vec_of_alpha, y = value, col=factor(N_rep))) + - geom_point() + facet_wrap(~variable, scales = "free_y", labeller = labeller(variable = label_wrap)) + labs(color = "N replicates") - return(figure) -} +#visualization function +source("mydatalocal/counts_simulation/src/visualization_fun.R") -## main -# Params to specify design -min_rep = 2 #/!\ = 1 forbidden -max_rep = 10 +########################## INPUT PARAMS ##################################### N_cond = 2 N_gene = 6000 +n_rep_sim = seq(2, 5, by = 1) ### number of replicate to assessed -mu_simul = seq(100, 15000, by = 200) -mu_simul -mu_simul <- rep.int(1500, 8) -res_simul <- mu_effect(alpha = 2.9, mu_simul) -reshape_res_simul <- res_simul %>% reshape2::melt(.,id = c("vec_of_mu")) -mu_effect_visualization(reshape_res_simul) - - -alpha_simul = seq(0.01, 10, by = 0.1) -alpha_simul -res_simul2 <- size_effect(mu = 10000, alpha_simul) -res_simul2 -reshape_res_simul2 <- res_simul2 %>% reshape2::melt(.,id = c("vec_of_alpha")) -size_effect_visualization(reshape_res_simul2) - - -## replicate effect -n_rep_sim = seq(2, 5, by = 1) +############################ MU effect ####################################### mu_simul_dtf_res <- data.frame() for (N_rep in n_rep_sim){ @@ -62,16 +36,17 @@ for (N_rep in n_rep_sim){ mu_simul_dtf_res <- rbind(mu_simul_dtf_res, tmp_reshape_res_simul) } -## LOG transform +###### LOG transform ####### + +# -> SEE linearity of var observed & mu + #mu_simul_dtf_res$value[mu_simul_dtf_res$variable=="var_observ"]<-log(mu_simul_dtf_res$value[mu_simul_dtf_res$variable=="var_observ"]) #mu_simul_dtf_res$vec_of_mu <- log(mu_simul_dtf_res$vec_of_mu) -## Visualization +###### Visualization ###### figure_mu_effect <- mu_effect_visualization(mu_simul_dtf_res) figure_mu_effect -svg("mydatalocal/counts_simulation/img/fig_mu_effect.svg") -figure_mu_effect -dev.off() + ########################### ALPHA effect #################################### n_rep_sim = seq(2, 5, by = 1) @@ -86,10 +61,47 @@ for (N_rep in n_rep_sim){ tmp_reshape_res_simul <- res_simul %>% reshape2::melt(.,id = c("vec_of_alpha", "N_rep")) alpha_simul_dtf_res <- rbind(alpha_simul_dtf_res, tmp_reshape_res_simul) } + + +###### Visualization ###### alpha_simul_dtf_res figure_alpha_effect <- size_effect_visualization(alpha_simul_dtf_res) +figure_alpha_effect + + + +########################### EXPORT RESULTS ################################# + +svg("mydatalocal/counts_simulation/img/fig_mu_effect.svg") +figure_mu_effect +dev.off() svg("mydatalocal/counts_simulation/img/fig_size_effect.svg") figure_alpha_effect dev.off() + + +########################### Beta test ##################################### + +## main +# Params to specify design +min_rep = 2 #/!\ = 1 forbidden +max_rep = 10 + + + +mu_simul = seq(100, 15000, by = 200) +mu_simul +mu_simul <- rep.int(1500, 8) +res_simul <- mu_effect(alpha = 2.9, mu_simul) +reshape_res_simul <- res_simul %>% reshape2::melt(.,id = c("vec_of_mu")) +mu_effect_visualization(reshape_res_simul) + + +alpha_simul = seq(0.01, 10, by = 0.1) +alpha_simul +res_simul2 <- size_effect(mu = 10000, alpha_simul) +res_simul2 +reshape_res_simul2 <- res_simul2 %>% reshape2::melt(.,id = c("vec_of_alpha")) +size_effect_visualization(reshape_res_simul2)