Delete simulators.R

src/simulators.R

-############################## 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)) 
-}