diff --git a/src/htrsim/input_estimation.R b/src/htrsim/input_estimation.R
new file mode 100644
index 0000000000000000000000000000000000000000..18cc20c3dd8852e3321113814cd29b524a1d94d5
--- /dev/null
+++ b/src/htrsim/input_estimation.R
@@ -0,0 +1,45 @@
+###################      Estimate alpha per gene        ########################
+
+estim.alpha <- function(dds){
+    #N.B: alpha = dispersion per gene
+    #dds  <- DESeq2::estimateDispersions(dds, quiet = F)
+    #dispersion estimation
+    dispersion_estimate <- dispersions(dds) 
+    
+    ## Shape and export
+    names(dispersion_estimate) <- tabl_cnts %>% rownames()
+    
+    ## drop NA in dispersion estimate (link to unexpress gene)
+    ### and convert to lovely dataframe
+    expressed_gene_dispersion <- dispersion_estimate[!is.na(dispersion_estimate)] %>%       
+      data.frame() %>% 
+      rownames_to_column() %>% 
+      rename(., "alpha" = ., gene_id = "rowname")
+    return(expressed_gene_dispersion)
+    
+    #disp_gene_express %>% dim
+    #write_tsv(disp_gene_express, 'results/2022-03-03/estimate_dispersion.tsv')
+}
+
+
+#################     Estimate mu distribution        #########################
+estim.mu <- function(dds){
+  mu_estimate <- dds@assays@data$mu
+  #dds@assays@data$mu %>% dim()
+  #mu_estimate %>% dim()
+  rownames(mu_estimate) = rownames(dds@assays@data$counts)
+  ## drop NA in dispersion estimate (link to unexpress gene)
+  ### and convert to lovely dataframe
+  mu_gene_express = mu_estimate %>% 
+                        na.omit()  %>%
+                        data.frame() 
+  
+  colnames(mu_gene_express) <- colnames(tabl_cnts)
+  mu_gene_express
+  mu_gene_express <- mu_gene_express %>% 
+                        mutate(gene_id = rownames(.)) %>% 
+                        select(gene_id, everything())
+  return(mu_gene_express)
+  
+  #write_tsv(mu_gene_express, 'results/2022-03-03/estimate_mu.tsv')
+}
\ No newline at end of file
diff --git a/src/htrsim/launch_deseq.R b/src/htrsim/launch_deseq.R
new file mode 100644
index 0000000000000000000000000000000000000000..7ad08f7158ba8cb9705602dce6939b3114734e5d
--- /dev/null
+++ b/src/htrsim/launch_deseq.R
@@ -0,0 +1,9 @@
+########### LAUNCH DESEQ #############
+
+run.deseq <- function(tabl_cnts, bioDesign ){
+  ## Design model
+  dds = DESeqDataSetFromMatrix( countData = round(tabl_cnts), colData = bioDesign , design = ~ mutant + env + mutant:env)
+  ## DESEQ standard analysis
+  dds <- DESeq(dds)
+  return(dds)
+}
\ No newline at end of file
diff --git a/src/htrsim/main.R b/src/htrsim/main.R
new file mode 100644
index 0000000000000000000000000000000000000000..505fe7594e4f2bc773f7496c84ce136d96d9f580
--- /dev/null
+++ b/src/htrsim/main.R
@@ -0,0 +1,26 @@
+######################### HTRSIM MAIN SRC ###################################
+
+htrsim <- function(countData, bioDesign, N_replicates){
+  
+  
+  source(file = "/home/rstudio/mydatalocal/counts_simulation/src/htrsim/launch_deseq.R")
+  dds <- run.deseq(countData, bioDesign)
+  
+  source(file="/home/rstudio/mydatalocal/counts_simulation/src/htrsim/input_estimation.R")
+  mu.input = estim.mu(dds)
+  alpha.input = estim.alpha(dds)
+  
+
+  source(file="/home/rstudio/mydatalocal/counts_simulation/src/htrsim/setup_cntsGenerator.R")
+  input = reshape_input2setup(mu.dtf = mu.input, alpha.dtf = alpha.input)
+  
+  setup.simulation <- setup_countGener(bioSample_id = input$bioSample_id,
+                                       n_rep = N_replicates,
+                                       alpha = input$alpha, 
+                                       gene_id = input$gene_id,
+                                       mu = input$mu)
+  ## Generate counts
+  source(file= "/home/rstudio/mydatalocal/counts_simulation/src/htrsim/counts_generator.R" )
+  htrs <- generate_counts(setup.simulation)
+  return(list(simul_cnts = htrs, mu.input = mu.input, alpha.input = alpha.input, setup = setup.simulation))
+}