diff --git a/src/simulators.R b/src/simulators.R deleted file mode 100644 index ed2fbcea96e00a869c26080bc679d436d5d7f09a..0000000000000000000000000000000000000000 --- a/src/simulators.R +++ /dev/null @@ -1,100 +0,0 @@ -############################## FUNCTIONS ################################### - -## count_generator(int, int, int) -> vec of length n_value -count_generator <- function(n_value, mu_theo, size_theo){ - rnbinom(n=n_value, mu = mu_theo, size = size_theo) -} - -## MATRIX_generator(int, int) -> matrice de dim(Ncol = N_cond*N_rep, Nrow = N_gene) -matrix_generator <- function(mu_theo, size_theo){ - n_value = N_gene*N_cond*N_rep #number of counts expected - mtx <- matrix(count_generator(n= n_value , mu = mu_theo, size = size_theo), ncol= N_cond*N_rep) - return(mtx) -} - - - - -mu_effect <- function(alpha, vec_of_mu){ - mu_observ <- c() - var_observ <- c() - statistical_power <- c() ## Init results of Differential expression analysis - res_DEA <- c() - for (mu in vec_of_mu){ - - # Print advancement message - cat(sprintf("Simulation for mu = %d\n", mu)) - - cnts <- matrix_generator(mu, alpha) - cond <- factor(rep(1:2, each=N_rep)) - dds <- DESeqDataSetFromMatrix(cnts, DataFrame(cond), ~ cond) - - # standard analysis - dds <- DESeq(dds, fitType='local') - res <- results(dds) - - #mu_observed - mu_observ <- c(mu_observ, mean(cnts)) - #var - var_observ <- c(var_observ, mean(rowVars(cnts))) - - - # results of DEA - cat(sprintf("Length table = %d\n", length(table(res$padj < 0.05)))) - if (dim(table(res$padj < 0.05)) == 1){ - cat(sprintf("NO DEG = %d\n", mu)) - cat(table(res$padj < 0.05)) - - res_DEA <- c(res_DEA, 0) ## case 1 : no DEG found by DESEQ2 - statistical_power <- c(statistical_power, NA) - } - else { - res_DEA <- c(res_DEA, table(res$padj < 0.05)[["TRUE"]]) ## case 2 : Nb DEG found by deseq2 - statistical_power <- c(statistical_power, min(abs(res$log2FoldChange[res$padj < 0.05]),na.rm=TRUE)) - } - } - return (data.frame(vec_of_mu, mu_observ, res_DEA, statistical_power, var_observ)) -} - - - -size_effect <- function(mu,vec_of_alpha){ - mu_observ <- c() - var_observ <- c() - statistical_power <- c() ## Init results of Differential expression analysis - res_DEA <- c() - for (alpha_params in vec_of_alpha){ - - # Print advancement message - cat(sprintf("Simulation for alpha = %f\n", alpha_params)) - - cnts <- matrix_generator(mu, alpha_params) - cond <- factor(rep(1:2, each=N_rep)) - dds <- DESeqDataSetFromMatrix(cnts, DataFrame(cond), ~ cond) - - # standard analysis - dds <- DESeq(dds, fitType='local') - res <- results(dds) - - #mu_observed - mu_observ <- c(mu_observ, mean(cnts)) - #var - var_observ <- c(var_observ, mean(rowVars(cnts))) - - - # results of DEA - cat(sprintf("Length table = %d\n", length(table(res$padj < 0.05)))) - if (dim(table(res$padj < 0.05)) == 1){ - cat(sprintf("NO DEG = %d\n", mu)) - cat(table(res$padj < 0.05)) - - res_DEA <- c(res_DEA, 0) ## case 1 : no DEG found by DESEQ2 - statistical_power <- c(statistical_power, NA) - } - else { - res_DEA <- c(res_DEA, table(res$padj < 0.05)[["TRUE"]]) ## case 2 : Nb DEG found by deseq2 - statistical_power <- c(statistical_power, min(abs(res$log2FoldChange[res$padj < 0.05]),na.rm=TRUE)) - } - } - return (data.frame(vec_of_alpha, mu_observ, res_DEA, statistical_power, var_observ)) -}