diff --git a/src/tutorial_htrsim.Rmd b/src/tutorial_htrsim.Rmd index 23fbe4a325adcc8d6ed5c8552decba10100a8415..8ddc0c5278709881116c40441fa15394635eea21 100644 --- a/src/tutorial_htrsim.Rmd +++ b/src/tutorial_htrsim.Rmd @@ -9,12 +9,13 @@ This is a tutorial for *htrsim* utilization ```{r required} library(data.table) library(tidyverse) -#library("DESeq2") +library(DESeq2) ``` ```{r setworkdir} -setwd("~/mydatalocal/counts_simulation/src/") +# on berthollet +setwd("/home2/aduvermy/counts_simulation/src/") ``` @@ -108,7 +109,7 @@ Next, plot your distribution of simulated counts per gene.</br> ```{r plot hist, echo = TRUE} data2plot <- htrs %>% reshape2::melt( ., id=c("name_gene"), variable.name = "Run") -data2plot <- data2plot %>% group_by(name_gene) %>% mutate(mean_obs = mean(value)) %>% ungroup() +data2plot <- data2plot %>% dplyr::group_by(name_gene) %>% dplyr::mutate(mean_obs = mean(value)) %>% dplyr::ungroup() figure = data2plot %>% ggplot(., aes(x=value)) + geom_histogram(fill= "grey", binwidth = 30) + @@ -327,11 +328,11 @@ We will mimic an existing design. ```{r message=FALSE, warning=FALSE} ## import mu(ij) params for each gene & each sample -mu_params = readr::read_tsv(file="../../rna-seq_public_library_investigations/results/2022-03-03/estimate_mu.tsv") +mu_params = readr::read_tsv(file="../../rna-seq_public_library_investigations/results/2022-03-03/estimate_mu.tsv", show_col_type=F) mu_params <- mu_params %>% select(., contains("rep1")) ## import alpha(i) for each genes -alpha_params = readr::read_tsv(file="../../rna-seq_public_library_investigations/results/2022-03-03/estimate_dispersion.tsv") +alpha_params = readr::read_tsv(file="../../rna-seq_public_library_investigations/results/2022-03-03/estimate_dispersion.tsv", show_col_type=F) ## Defining sample names @@ -366,6 +367,7 @@ genes <- list(name = nameGene.set , alpha = alphaGene.set) %>% as.data.frame() source(file = "htrsim/setup_cntsGenerator.R") ## with mu dataframe setup.simulation <- setup_countGener(sample_names= samples$name , n_rep= samples$n_rep , gene_dispersion = genes$alpha, gene_names = genes$name, mu = mu_params) +setup.simulation %>% head source(file= "htrsim/counts_generator.R" ) htrs <- generate_counts(setup.simulation) ```