From 23a127b9a6f94c58a2459eaf77d933bdc755796c Mon Sep 17 00:00:00 2001
From: aduvermy <arnaud.duvermy@ens-lyon.fr>
Date: Mon, 28 Mar 2022 15:28:10 +0000
Subject: [PATCH] add new src
---
src/htrsim/input_estimation.R | 45 +++++++++++++++++++++++++++++++++++
src/htrsim/launch_deseq.R | 9 +++++++
src/htrsim/main.R | 26 ++++++++++++++++++++
3 files changed, 80 insertions(+)
create mode 100644 src/htrsim/input_estimation.R
create mode 100644 src/htrsim/launch_deseq.R
create mode 100644 src/htrsim/main.R
diff --git a/src/htrsim/input_estimation.R b/src/htrsim/input_estimation.R
new file mode 100644
index 0000000..18cc20c
--- /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 0000000..7ad08f7
--- /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 0000000..505fe75
--- /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))
+}
--
GitLab