From eb82537dfff3e9c9afd9030dede1d076a5e037a7 Mon Sep 17 00:00:00 2001
From: aduvermy <arnaud.duvermy@ens-lyon.fr>
Date: Wed, 21 Sep 2022 09:55:19 +0200
Subject: [PATCH] mv results v1
---
results/{ => v1}/2022-04-22_batch_effect.Rmd | 0
results/{ => v1}/2022-04-22_htrsim.Rmd | 856 +++++++++---------
.../{ => v1}/2022-04-22_lowCnts_effect.Rmd | 0
.../{ => v1}/2022-04-24_epsilon_effect.Rmd | 0
.../{ => v1}/2022-04-28_beta_SE_filter.Rmd | 0
results/{ => v1}/2022-05-22_alpha_effect.Rmd | 658 +++++++-------
results/{ => v1}/2022-06-24_subsampling.Rmd | 0
.../{ => v1}/2022-06-28_kallistoEffect.Rmd | 0
results/{ => v1}/2022_06_08_investigation.Rmd | 0
results/{ => v1}/SVA.R | 0
results/{ => v1}/SVA.Rmd | 0
results/{ => v1}/css/footer.css | 0
results/{ => v1}/css/template.css | 0
results/{ => v1}/lab-meeting_preparation.Rmd | 122 +--
14 files changed, 818 insertions(+), 818 deletions(-)
rename results/{ => v1}/2022-04-22_batch_effect.Rmd (100%)
rename results/{ => v1}/2022-04-22_htrsim.Rmd (96%)
rename results/{ => v1}/2022-04-22_lowCnts_effect.Rmd (100%)
rename results/{ => v1}/2022-04-24_epsilon_effect.Rmd (100%)
rename results/{ => v1}/2022-04-28_beta_SE_filter.Rmd (100%)
rename results/{ => v1}/2022-05-22_alpha_effect.Rmd (96%)
rename results/{ => v1}/2022-06-24_subsampling.Rmd (100%)
rename results/{ => v1}/2022-06-28_kallistoEffect.Rmd (100%)
rename results/{ => v1}/2022_06_08_investigation.Rmd (100%)
rename results/{ => v1}/SVA.R (100%)
rename results/{ => v1}/SVA.Rmd (100%)
rename results/{ => v1}/css/footer.css (100%)
rename results/{ => v1}/css/template.css (100%)
rename results/{ => v1}/lab-meeting_preparation.Rmd (96%)
diff --git a/results/2022-04-22_batch_effect.Rmd b/results/v1/2022-04-22_batch_effect.Rmd
similarity index 100%
rename from results/2022-04-22_batch_effect.Rmd
rename to results/v1/2022-04-22_batch_effect.Rmd
diff --git a/results/2022-04-22_htrsim.Rmd b/results/v1/2022-04-22_htrsim.Rmd
similarity index 96%
rename from results/2022-04-22_htrsim.Rmd
rename to results/v1/2022-04-22_htrsim.Rmd
index 62575c9..10be5f9 100644
--- a/results/2022-04-22_htrsim.Rmd
+++ b/results/v1/2022-04-22_htrsim.Rmd
@@ -1,428 +1,428 @@
----
-title: "HTRSIM"
-date: '2022-04-21'
-output:
- html_document:
-
-
-css:
- - css/template.css
- - css/footer.css
-
----
-
-
-## Introduction
-
-
-In living world, phenotypes are understanding as a mixture between a genotype effect, an environment effect and an interaction between G&E.
-$$Phenotype = Genotype + Environment + Genotype.Environment$$
-The quantification of each strengths (G,E; G&E) can be estimate by a coefficient $\beta$.
-Then, our expression becomes:
-$$Phenotype = \beta_{G} * Genotype + \beta_{E}*Environment + \beta_{G*E} * Genotype.Environment + \epsilon$$
-Notice that $\beta$ is specific of each component. Furthermore, we introduced above $\epsilon$. It's the residual of the model. $\epsilon$ can be seen as the difference between observed values and values predicted by the model.
-
-Genes expression can be also considered as a phenotype. <br>
-According to this, the quantification of $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$ for a given gene in a given condition may open the possibility to assess differences between the strengths in presence in different conditions.
-
-That's the purpose of Htrsim !
-
-## Htrsim
-
-##### Model
-
-In this aim, Htrsim is based on a model. <br>
-Because of is easy of use this model is managed by DESEQ2.
-Then, $K_{ij}$ for gene i, sample j are modeled using a Negative Binomial distribution with fitted mean $\mu_{ij}$ and a gene-specific dispersion parameter $\alpha_i$.
-
-$$
-K_{ij} \sim {\sf NB}(\mu_{ij} ; \sigma_i)
-$$
-$$
-\mu_{ij} = s_jq_{ij}
-$$
-$$
-log_2(q_{ij}) = x_j*\beta_i
-$$
-The fitted mean is composed of a sample-specific size factor $s_j$ and a parameter qij proportional to the expected true concentration of fragments for sample j.
-The coefficients $\beta_i$ give the log2 fold changes for gene i for each column of the model matrix X. The sample-specific size factors can be replaced by gene-specific normalization factors for each sample using normalizationFactors.
-
-
-
-According to the DESEQ2 GLM and our purpose, we can write:
-$$
-log_2(q_{ij}) = \beta_{G}*G + \beta_{E}*E + \beta_{G*E}*G.E + \beta_{0} + \epsilon_{ij}
-$$
-
-According to this generalized linear model, we wish to estimate $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$ for a given gene i, in a given condition j. Achieve this, would allow us to quantify each strengths (G, E, G&E) for a given gene i, in a given condition j.
-
-
-##### Required
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-library(htrsim)
-library(tidyverse)
-library(reshape2)
-```
-
-##### Worklow
-
-Using public libraries (from BioProject PRJNA675209b - chinese paper), and an usual RNA-seq pipeline, we build actual RNA-seq counts per genes for 3 genotypes and 2 environments.<br>
-<br>
-Using htrsim (in particular DESEQ2) and this count table, we are able to estimate $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$.
-
-
-a. Input
-
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## Import & reshape table counts
-fn = system.file("extdata/", "public_tablCnts.tsv", package = "htrsim")
-tabl_cnts <- read.table(file = fn, header = TRUE)
-rownames(tabl_cnts) <- tabl_cnts$gene_id
-tabl_cnts <- tabl_cnts %>% select(-gene_id)##suppr colonne GeneID
-tabl_cnts <- tabl_cnts %>% select(-gene_name) ##suppr colonne GeneName
-
-## import design of bioProject
-fn = system.file("extdata/", "public_bioDesign.csv", package = "htrsim")
-bioDesign <- read.table(file = fn, header = T, sep = ';')
-
-```
-
-b. Launch DESEQ2
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-dds = run.deseq(tabl_cnts = tabl_cnts, bioDesign = bioDesign)
-```
-
-
-DESEQ returns a dds object which contains many, many things ... <br>
-In particular it contains the $\beta$ coefficients. <br>
-<br>
-You can access to beta coefficients using:
-
-```{r}
-dds.mcols = S4Vectors::mcols(dds,use.names=TRUE)
-```
-
-c. $\mu_{ij}$
-
-Following our model, we can estimate $log_2(\mu_{ij]})$ from $\beta$ coefficients inferred by DESEQ2,
-
-$$
-log_2(q_{ij}) = \beta_{G}*G + \beta_{E}*E + \beta_{G*E}*G.E + \beta_{0} + \epsilon_{ij}
-$$
-
-Then, $\mu_{ij]}$ can be estimate
-
-$$
-\mu_{ij} = s_j * 2^{log_2(q_{ij})}
-$$
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## Model matrix per samples
-mm <- model.matrix(~genotype + env + genotype:env, bioDesign)
-
-## Input estimation
-estim_mu = estim.mu(dds, mm)
-mu.input = estim_mu$mu
-```
-
-d. $K_{ij}$
-
-As defined by our model, counts $K_{ij}$ for gene i, sample j are modeled using a Negative Binomial distribution with fitted mean $\mu_{ij}$ and a gene-specific dispersion parameter $\alpha_i$.
-
-$$
-K_{ij} \sim {\sf NB}(\mu_{ij} ; \alpha_i)
-$$
-The gene-specific dispersion parameter $\alpha_i$ is also stored in the dds object.<br>
-You can access to $\alpha_i$ using:
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-alpha.input = estim.alpha(dds)
-```
-
-Knowing $\alpha_i$ and $\mu_{ij]}$ for each gene and each condition given by the BioProject PRJNA675209b - chinese paper.
-We are now able to simulate $K_{ij}$ for each gene and each condition given by the BioProject PRJNA675209b.
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-# Setup simulation
-input = reshape_input2setup(mu.dtf = mu.input, alpha.dtf = alpha.input, average_rep = FALSE)
-
-#input$gene_id
-setup.simulation <- setup_countGener(bioSample_id = input$bioSample_id,
- n_rep = 1,
- alpha = input$alpha,
- gene_id = input$gene_id,
- mu = input$mu)
-
-#setup.simulation %>% dim()
-# Simulate counts
-htrs <- generate_counts(setup.simulation)
-
-```
-
-
-```{r, error=TRUE}
-sample = htrs %>% select(where(is.numeric)) %>% colnames()
-genotype = htrs %>%
- select(where(is.numeric)) %>%
- colnames() %>%
- map(., ~str_split(.,pattern = '_', simplify = T)[1]) %>% unlist()
-
-env = htrs %>%
- select(where(is.numeric)) %>%
- colnames() %>%
- map(., ~str_split(.,pattern = '_', simplify = T)[2]) %>% unlist()
-
-designSimu = cbind(sample, env, genotype) %>% data.frame()
-
-## RESHAPE HTRS
-htrs.reshape = htrs
-rownames(htrs.reshape) = htrs.reshape$gene_id
-htrs.reshape = htrs.reshape %>% select(-gene_id)
-########### LAUNCH DESEQ #############
-## Design model - specify reference
-designSimu$genotype <- factor(x = designSimu$genotype,levels = c('WT','Msn2D', 'Msn4D'))
-designSimu$env <- factor(x = designSimu$env,levels = c('control', 'KCl'))
-
-
-k_ij.simulation = htrs.reshape
-
-## DESEQ standard analysis
-dds_simu = DESeq2::DESeqDataSetFromMatrix( countData = k_ij.simulation ,
- colData = designSimu ,
- design = ~ genotype + env + genotype:env)
-
-max(k_ij.simulation)
-.Machine$integer.max
-
-```
-
-##### Evaluation
-
-```{r message=FALSE, warning=FALSE, include=TRUE}
-k_ij.simu = htrs %>% select(-gene_id) %>% flatten() %>% unlist()
-k_ij.actual = tabl_cnts %>% flatten() %>% unlist()
-
-df = cbind(k_ij.actual, k_ij.simu) %>% reshape2::melt(., value.name = "k_ij", variable.name = "origin")
-df$origin = df$Var2
-df = df %>% select(-Var2)
-max_k_ij.simu = df %>% filter(origin == "k_ij.simu") %>% select(k_ij) %>% max()
-max_k_ij.actual = df %>% filter(origin == "k_ij.actual") %>% select(k_ij) %>% max()
-
-ggplot(df, aes(x = k_ij, fill= origin )) + geom_density(bins = 100, alpha = 0.5) +
- geom_vline(xintercept = max_k_ij.actual, col = "#F8766D" ) +
- geom_vline(xintercept = (max_k_ij.simu), col= "#00BFC4" ) +
- scale_x_log10()
-```
-
-Comment: $K_{ij}$ simulated are abnormally huge !
-Comment: $K_{ij}$ simulated are slightly different from the actual K_{ij} !
-
-
-## Why so much differences
-
-b. $\epsilon$
-
-In our model, we define as follow:
-$$
-\epsilon_{ij} \sim {\sf N}(0 ; deviance_i)
-$$
-
-Let's see the distribution of $deviance_{i}$.
-
-
-```{r warning=FALSE}
-#estim_mu$beta.matrix
-
-deviance_i = estim_mu$deviance.sqrt[!is.na(estim_mu$deviance.sqrt)]^2
-#epsilon_ij <- mm[,1] %>% map(., ~rnorm(deviance_i.sqrt, mean = 0, sd = deviance_i.sqrt )) %>% data.frame() %>% flatten() %>% unlist()
-
-
-# Histogram logarithmic y axis
-ggplot(data.frame(deviance_i), aes(deviance_i)) +
- geom_histogram(bins = 100) #+ scale_x_log10()
-
-
-```
-
-
-The deviance is also inferred by DESEQ while computing its model.
-deviance is mostly inferred between 100 and 200.
-
-
-
-$$
-log_2(q_{ij}) = \beta_{G}*G + \beta_{E}*E + \beta_{G*E}*G.E + \beta_{0} + \epsilon_{ij}
-$$
-
-```{r warning=FALSE}
-#estim_mu$beta.matrix
-
-deviance_i.sqrt = estim_mu$deviance.sqrt[!is.na(estim_mu$deviance.sqrt)]
-epsilon_ij <- mm[,1] %>% map(., ~rnorm(deviance_i.sqrt, mean = 0, sd = deviance_i.sqrt )) %>% data.frame() %>% flatten() %>% unlist()
-#epsilon_ij <- mm[,1] %>% map(., ~rnorm(deviance_i.sqrt, mean = 0, sd = 0 )) %>% data.frame() %>% flatten() %>% unlist()
-
-# Histogram logarithmic y axis
-ggplot(data.frame(epsilon_ij), aes(epsilon_ij)) +
- geom_histogram(bins = 100) #+ scale_x_log10()
-
-
-```
-
-
-Comment: Some $\epsilon_{ij}$ are huge !
-Recall: $\epsilon$ can be seen as the difference between observed values and values predicted by the model.
-
-A large panel of $\epsilon$ mean that the model doesn't fit well with the observed data.
-
-It means that even if $\beta$ coefficients are well estimate. $log_2(q_{ij})$ will vary around them with a large panel of values (+/- 40)
-
-
-
-
-```{r warning=FALSE}
-beta.dtf = estim_mu$beta.matrix %>% data.frame()
-beta.dtf.long = beta.dtf %>% reshape2::melt(., value.name = "beta", variable.name = "origin")
-
-ggplot(beta.dtf.long, aes(x = beta )) + geom_density(bins = 100, alpha = 0.5, fill = 'grey') + facet_grid(~origin, scales = "free_x")
-
-```
-```{r warning=FALSE}
-beta.dtf = estim_mu$beta.matrix %>% data.frame()
-beta.dtf.long = beta.dtf %>% reshape2::melt(., value.name = "beta", variable.name = "origin")
-
-
-
-## Standard Error
-B0 = estim_mu$dds.mcols$SE_Intercept
-B1 = estim_mu$dds.mcols$SE_genotype_msn2D_vs_wt
-B2 <- estim_mu$dds.mcols$SE_genotype_msn4D_vs_wt
-B3 <- estim_mu$dds.mcols$SE_env_kcl_vs_control
-B4 <- estim_mu$dds.mcols$SE_genotypemsn2D.envkcl
-B5 <- estim_mu$dds.mcols$SE_genotypemsn4D.envkcl
-
-
-SE_B.dtf <- cbind(B0, B1, B2, B3, B4, B5) %>% data.frame()
-SE_B.dtf.long = SE_B.dtf %>% reshape2::melt(., value.name = "SE_beta", variable.name = "origin")
-
-ggplot(SE_B.dtf.long, aes(x = SE_beta, fill= origin )) + geom_density(bins = 100, alpha = 0.5) + facet_grid(~origin)
-```
-
-```{r warning=FALSE}
-bind_dtf<- cbind(SE_B.dtf.long, beta.dtf.long %>% select(-origin))
-ggplot(bind_dtf, aes(x = beta, y= SE_beta, fill= origin )) + geom_point(alpha = 0.1) + facet_grid(~origin)
-
-
-#new <- bind_dtf %>% mutate(annot = ifelse(origin == "B4 | B5" && SE_beta > 6 , TRUE, FALSE ))
-#new <- bind_dtf %>% tail
-#new %>% filter(beta == "B4")
-```
-
-
-```{r warning=FALSE}
-#dim(htrs)
-
-#new <- bind_dtf %>% mutate(annot = ifelse(((origin == "B4") | (origin == "B5")) & (SE_beta > 6) , TRUE, FALSE ))
-### WARNING
-#new %>% dcast(., annot ~ origin)
-
-
-SE_threshold = 6
-SE_B.dtf.annot = SE_B.dtf %>% mutate(annot = ifelse((B4 > SE_threshold) | (B5 > SE_threshold) , TRUE, FALSE ))
-SE_B.dtf.annot %>% group_by(annot) %>% tally()
-SE_B.dtf.annot.long = SE_B.dtf.annot %>% reshape2::melt(., value.name = "SE_beta", variable.name = "origin")
-
-
-bind_dtf.annot<- cbind(SE_B.dtf.annot.long, beta.dtf.long %>% select(-origin))
-bind_dtf.annot = bind_dtf.annot %>% filter(!is.na(annot))
-ggplot(bind_dtf.annot, aes(x = beta, y= SE_beta, col = annot )) + geom_point(alpha = 0.1, na.rm = T) + facet_grid(~origin)
-
-
-```
-
-
-## Beta0 vs SE & deviance
-
-
-
-```{r}
-
-B0 = beta.dtf$B0
-
-
-L = SE_B.dtf %>% colnames() %>% length()
-
-
-B0_vector = replicate(L, B0) %>% as.data.frame() %>% flatten() %>% unlist()
-deviance.sqrt_vec = replicate(L , estim_mu$deviance.sqrt )%>% as.data.frame() %>% flatten() %>% unlist()
-SE_B.dtf.annot.long$B0 = B0_vector
-SE_B.dtf.annot.long$deviance.sqrt = estim_mu$deviance.sqrt
-
-ggplot(SE_B.dtf.annot.long, aes(x = B0, y = SE_beta, col= annot )) + geom_point(alpha = 0.1, na.rm = T) + facet_grid(~origin)
-
-
-SE_B.dtf.annot.long_B0= SE_B.dtf.annot.long %>% filter(origin == "B0") %>% filter(!is.na(annot))
-
-ggplot( SE_B.dtf.annot.long_B0, aes(x = B0, y = 2^deviance.sqrt, col = annot)) + geom_point(alpha = 0.1, na.rm = TRUE) + scale_y_log10()
-
-
-```
-
-
-
-
-```{undefined eval=FALSE, include=FALSE}
-hist(dds.mcols$dispFit)
-hist(log_qij, )
-hist(log10(alpha.input$alpha))
-max(dds.mcols$dispersion, na.rm = T)
-hist(log(dds.mcols$dispersion))
-max(dds.mcols$deviance, na.rm = T)
-hist(dds.mcols$deviance)
-max(dds.mcols$dispFit, na.rm = T)
-DESeq2::design(dds)
-
-fitted.common.scale = t(t(dds@assays@data$mu)/dds$sizeFactor)
-
-hist(t(t(dds@assays@data$mu)/dds$sizeFactor))
-hist(residual_deseq)
-max(residual_deseq, na.rm = T)
-residual_deseq = (DESeq2::counts(dds, normalized=TRUE) - fitted.common.scale )
-
-
-w =DESeq2::nbinomWaldTest(dds)
-w@dispersionFunction()
-#S4Vectors::assays(dds)[["mu"]]
-
-
-vst = DESeq2::varianceStabilizingTransformation(dds)
-vst@assays@data$
-```
-
-
-
-https://support.bioconductor.org/p/123305/
-https://support.bioconductor.org/p/60567/
-http://bioconductor.org/packages/release/bioc/vignettes/RUVSeq/inst/doc/RUVSeq.pdf
-https://bioinformatics-core-shared-training.github.io/RNAseq-R/slides/LinearModels.pdf
-
-
-```{r eval=FALSE, include=FALSE}
-#install.packages("RUVSeq")
-#BiocManager::install("RUVSeq")
-library(RUVSeq)
-
-mm
-y <- DGEList(counts=counts(tabl_cnts), group=x)
-y <- calcNormFactors(y, method="upperquartile")
-y <- estimateGLMCommonDisp(y, design)
-y <- estimateGLMTagwiseDisp(y, design)
-fit <- glmFit(y, design)
-res <- residuals(fit, type="deviance")
-
-
-```
-
-
-
-
+---
+title: "HTRSIM"
+date: '2022-04-21'
+output:
+ html_document:
+
+
+css:
+ - css/template.css
+ - css/footer.css
+
+---
+
+
+## Introduction
+
+
+In living world, phenotypes are understanding as a mixture between a genotype effect, an environment effect and an interaction between G&E.
+$$Phenotype = Genotype + Environment + Genotype.Environment$$
+The quantification of each strengths (G,E; G&E) can be estimate by a coefficient $\beta$.
+Then, our expression becomes:
+$$Phenotype = \beta_{G} * Genotype + \beta_{E}*Environment + \beta_{G*E} * Genotype.Environment + \epsilon$$
+Notice that $\beta$ is specific of each component. Furthermore, we introduced above $\epsilon$. It's the residual of the model. $\epsilon$ can be seen as the difference between observed values and values predicted by the model.
+
+Genes expression can be also considered as a phenotype. <br>
+According to this, the quantification of $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$ for a given gene in a given condition may open the possibility to assess differences between the strengths in presence in different conditions.
+
+That's the purpose of Htrsim !
+
+## Htrsim
+
+##### Model
+
+In this aim, Htrsim is based on a model. <br>
+Because of is easy of use this model is managed by DESEQ2.
+Then, $K_{ij}$ for gene i, sample j are modeled using a Negative Binomial distribution with fitted mean $\mu_{ij}$ and a gene-specific dispersion parameter $\alpha_i$.
+
+$$
+K_{ij} \sim {\sf NB}(\mu_{ij} ; \sigma_i)
+$$
+$$
+\mu_{ij} = s_jq_{ij}
+$$
+$$
+log_2(q_{ij}) = x_j*\beta_i
+$$
+The fitted mean is composed of a sample-specific size factor $s_j$ and a parameter qij proportional to the expected true concentration of fragments for sample j.
+The coefficients $\beta_i$ give the log2 fold changes for gene i for each column of the model matrix X. The sample-specific size factors can be replaced by gene-specific normalization factors for each sample using normalizationFactors.
+
+
+
+According to the DESEQ2 GLM and our purpose, we can write:
+$$
+log_2(q_{ij}) = \beta_{G}*G + \beta_{E}*E + \beta_{G*E}*G.E + \beta_{0} + \epsilon_{ij}
+$$
+
+According to this generalized linear model, we wish to estimate $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$ for a given gene i, in a given condition j. Achieve this, would allow us to quantify each strengths (G, E, G&E) for a given gene i, in a given condition j.
+
+
+##### Required
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+library(htrsim)
+library(tidyverse)
+library(reshape2)
+```
+
+##### Worklow
+
+Using public libraries (from BioProject PRJNA675209b - chinese paper), and an usual RNA-seq pipeline, we build actual RNA-seq counts per genes for 3 genotypes and 2 environments.<br>
+<br>
+Using htrsim (in particular DESEQ2) and this count table, we are able to estimate $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$.
+
+
+a. Input
+
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## Import & reshape table counts
+fn = system.file("extdata/", "public_tablCnts.tsv", package = "htrsim")
+tabl_cnts <- read.table(file = fn, header = TRUE)
+rownames(tabl_cnts) <- tabl_cnts$gene_id
+tabl_cnts <- tabl_cnts %>% select(-gene_id)##suppr colonne GeneID
+tabl_cnts <- tabl_cnts %>% select(-gene_name) ##suppr colonne GeneName
+
+## import design of bioProject
+fn = system.file("extdata/", "public_bioDesign.csv", package = "htrsim")
+bioDesign <- read.table(file = fn, header = T, sep = ';')
+
+```
+
+b. Launch DESEQ2
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+dds = run.deseq(tabl_cnts = tabl_cnts, bioDesign = bioDesign)
+```
+
+
+DESEQ returns a dds object which contains many, many things ... <br>
+In particular it contains the $\beta$ coefficients. <br>
+<br>
+You can access to beta coefficients using:
+
+```{r}
+dds.mcols = S4Vectors::mcols(dds,use.names=TRUE)
+```
+
+c. $\mu_{ij}$
+
+Following our model, we can estimate $log_2(\mu_{ij]})$ from $\beta$ coefficients inferred by DESEQ2,
+
+$$
+log_2(q_{ij}) = \beta_{G}*G + \beta_{E}*E + \beta_{G*E}*G.E + \beta_{0} + \epsilon_{ij}
+$$
+
+Then, $\mu_{ij]}$ can be estimate
+
+$$
+\mu_{ij} = s_j * 2^{log_2(q_{ij})}
+$$
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## Model matrix per samples
+mm <- model.matrix(~genotype + env + genotype:env, bioDesign)
+
+## Input estimation
+estim_mu = estim.mu(dds, mm)
+mu.input = estim_mu$mu
+```
+
+d. $K_{ij}$
+
+As defined by our model, counts $K_{ij}$ for gene i, sample j are modeled using a Negative Binomial distribution with fitted mean $\mu_{ij}$ and a gene-specific dispersion parameter $\alpha_i$.
+
+$$
+K_{ij} \sim {\sf NB}(\mu_{ij} ; \alpha_i)
+$$
+The gene-specific dispersion parameter $\alpha_i$ is also stored in the dds object.<br>
+You can access to $\alpha_i$ using:
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+alpha.input = estim.alpha(dds)
+```
+
+Knowing $\alpha_i$ and $\mu_{ij]}$ for each gene and each condition given by the BioProject PRJNA675209b - chinese paper.
+We are now able to simulate $K_{ij}$ for each gene and each condition given by the BioProject PRJNA675209b.
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+# Setup simulation
+input = reshape_input2setup(mu.dtf = mu.input, alpha.dtf = alpha.input, average_rep = FALSE)
+
+#input$gene_id
+setup.simulation <- setup_countGener(bioSample_id = input$bioSample_id,
+ n_rep = 1,
+ alpha = input$alpha,
+ gene_id = input$gene_id,
+ mu = input$mu)
+
+#setup.simulation %>% dim()
+# Simulate counts
+htrs <- generate_counts(setup.simulation)
+
+```
+
+
+```{r, error=TRUE}
+sample = htrs %>% select(where(is.numeric)) %>% colnames()
+genotype = htrs %>%
+ select(where(is.numeric)) %>%
+ colnames() %>%
+ map(., ~str_split(.,pattern = '_', simplify = T)[1]) %>% unlist()
+
+env = htrs %>%
+ select(where(is.numeric)) %>%
+ colnames() %>%
+ map(., ~str_split(.,pattern = '_', simplify = T)[2]) %>% unlist()
+
+designSimu = cbind(sample, env, genotype) %>% data.frame()
+
+## RESHAPE HTRS
+htrs.reshape = htrs
+rownames(htrs.reshape) = htrs.reshape$gene_id
+htrs.reshape = htrs.reshape %>% select(-gene_id)
+########### LAUNCH DESEQ #############
+## Design model - specify reference
+designSimu$genotype <- factor(x = designSimu$genotype,levels = c('WT','Msn2D', 'Msn4D'))
+designSimu$env <- factor(x = designSimu$env,levels = c('control', 'KCl'))
+
+
+k_ij.simulation = htrs.reshape
+
+## DESEQ standard analysis
+dds_simu = DESeq2::DESeqDataSetFromMatrix( countData = k_ij.simulation ,
+ colData = designSimu ,
+ design = ~ genotype + env + genotype:env)
+
+max(k_ij.simulation)
+.Machine$integer.max
+
+```
+
+##### Evaluation
+
+```{r message=FALSE, warning=FALSE, include=TRUE}
+k_ij.simu = htrs %>% select(-gene_id) %>% flatten() %>% unlist()
+k_ij.actual = tabl_cnts %>% flatten() %>% unlist()
+
+df = cbind(k_ij.actual, k_ij.simu) %>% reshape2::melt(., value.name = "k_ij", variable.name = "origin")
+df$origin = df$Var2
+df = df %>% select(-Var2)
+max_k_ij.simu = df %>% filter(origin == "k_ij.simu") %>% select(k_ij) %>% max()
+max_k_ij.actual = df %>% filter(origin == "k_ij.actual") %>% select(k_ij) %>% max()
+
+ggplot(df, aes(x = k_ij, fill= origin )) + geom_density(bins = 100, alpha = 0.5) +
+ geom_vline(xintercept = max_k_ij.actual, col = "#F8766D" ) +
+ geom_vline(xintercept = (max_k_ij.simu), col= "#00BFC4" ) +
+ scale_x_log10()
+```
+
+Comment: $K_{ij}$ simulated are abnormally huge !
+Comment: $K_{ij}$ simulated are slightly different from the actual K_{ij} !
+
+
+## Why so much differences
+
+b. $\epsilon$
+
+In our model, we define as follow:
+$$
+\epsilon_{ij} \sim {\sf N}(0 ; deviance_i)
+$$
+
+Let's see the distribution of $deviance_{i}$.
+
+
+```{r warning=FALSE}
+#estim_mu$beta.matrix
+
+deviance_i = estim_mu$deviance.sqrt[!is.na(estim_mu$deviance.sqrt)]^2
+#epsilon_ij <- mm[,1] %>% map(., ~rnorm(deviance_i.sqrt, mean = 0, sd = deviance_i.sqrt )) %>% data.frame() %>% flatten() %>% unlist()
+
+
+# Histogram logarithmic y axis
+ggplot(data.frame(deviance_i), aes(deviance_i)) +
+ geom_histogram(bins = 100) #+ scale_x_log10()
+
+
+```
+
+
+The deviance is also inferred by DESEQ while computing its model.
+deviance is mostly inferred between 100 and 200.
+
+
+
+$$
+log_2(q_{ij}) = \beta_{G}*G + \beta_{E}*E + \beta_{G*E}*G.E + \beta_{0} + \epsilon_{ij}
+$$
+
+```{r warning=FALSE}
+#estim_mu$beta.matrix
+
+deviance_i.sqrt = estim_mu$deviance.sqrt[!is.na(estim_mu$deviance.sqrt)]
+epsilon_ij <- mm[,1] %>% map(., ~rnorm(deviance_i.sqrt, mean = 0, sd = deviance_i.sqrt )) %>% data.frame() %>% flatten() %>% unlist()
+#epsilon_ij <- mm[,1] %>% map(., ~rnorm(deviance_i.sqrt, mean = 0, sd = 0 )) %>% data.frame() %>% flatten() %>% unlist()
+
+# Histogram logarithmic y axis
+ggplot(data.frame(epsilon_ij), aes(epsilon_ij)) +
+ geom_histogram(bins = 100) #+ scale_x_log10()
+
+
+```
+
+
+Comment: Some $\epsilon_{ij}$ are huge !
+Recall: $\epsilon$ can be seen as the difference between observed values and values predicted by the model.
+
+A large panel of $\epsilon$ mean that the model doesn't fit well with the observed data.
+
+It means that even if $\beta$ coefficients are well estimate. $log_2(q_{ij})$ will vary around them with a large panel of values (+/- 40)
+
+
+
+
+```{r warning=FALSE}
+beta.dtf = estim_mu$beta.matrix %>% data.frame()
+beta.dtf.long = beta.dtf %>% reshape2::melt(., value.name = "beta", variable.name = "origin")
+
+ggplot(beta.dtf.long, aes(x = beta )) + geom_density(bins = 100, alpha = 0.5, fill = 'grey') + facet_grid(~origin, scales = "free_x")
+
+```
+```{r warning=FALSE}
+beta.dtf = estim_mu$beta.matrix %>% data.frame()
+beta.dtf.long = beta.dtf %>% reshape2::melt(., value.name = "beta", variable.name = "origin")
+
+
+
+## Standard Error
+B0 = estim_mu$dds.mcols$SE_Intercept
+B1 = estim_mu$dds.mcols$SE_genotype_msn2D_vs_wt
+B2 <- estim_mu$dds.mcols$SE_genotype_msn4D_vs_wt
+B3 <- estim_mu$dds.mcols$SE_env_kcl_vs_control
+B4 <- estim_mu$dds.mcols$SE_genotypemsn2D.envkcl
+B5 <- estim_mu$dds.mcols$SE_genotypemsn4D.envkcl
+
+
+SE_B.dtf <- cbind(B0, B1, B2, B3, B4, B5) %>% data.frame()
+SE_B.dtf.long = SE_B.dtf %>% reshape2::melt(., value.name = "SE_beta", variable.name = "origin")
+
+ggplot(SE_B.dtf.long, aes(x = SE_beta, fill= origin )) + geom_density(bins = 100, alpha = 0.5) + facet_grid(~origin)
+```
+
+```{r warning=FALSE}
+bind_dtf<- cbind(SE_B.dtf.long, beta.dtf.long %>% select(-origin))
+ggplot(bind_dtf, aes(x = beta, y= SE_beta, fill= origin )) + geom_point(alpha = 0.1) + facet_grid(~origin)
+
+
+#new <- bind_dtf %>% mutate(annot = ifelse(origin == "B4 | B5" && SE_beta > 6 , TRUE, FALSE ))
+#new <- bind_dtf %>% tail
+#new %>% filter(beta == "B4")
+```
+
+
+```{r warning=FALSE}
+#dim(htrs)
+
+#new <- bind_dtf %>% mutate(annot = ifelse(((origin == "B4") | (origin == "B5")) & (SE_beta > 6) , TRUE, FALSE ))
+### WARNING
+#new %>% dcast(., annot ~ origin)
+
+
+SE_threshold = 6
+SE_B.dtf.annot = SE_B.dtf %>% mutate(annot = ifelse((B4 > SE_threshold) | (B5 > SE_threshold) , TRUE, FALSE ))
+SE_B.dtf.annot %>% group_by(annot) %>% tally()
+SE_B.dtf.annot.long = SE_B.dtf.annot %>% reshape2::melt(., value.name = "SE_beta", variable.name = "origin")
+
+
+bind_dtf.annot<- cbind(SE_B.dtf.annot.long, beta.dtf.long %>% select(-origin))
+bind_dtf.annot = bind_dtf.annot %>% filter(!is.na(annot))
+ggplot(bind_dtf.annot, aes(x = beta, y= SE_beta, col = annot )) + geom_point(alpha = 0.1, na.rm = T) + facet_grid(~origin)
+
+
+```
+
+
+## Beta0 vs SE & deviance
+
+
+
+```{r}
+
+B0 = beta.dtf$B0
+
+
+L = SE_B.dtf %>% colnames() %>% length()
+
+
+B0_vector = replicate(L, B0) %>% as.data.frame() %>% flatten() %>% unlist()
+deviance.sqrt_vec = replicate(L , estim_mu$deviance.sqrt )%>% as.data.frame() %>% flatten() %>% unlist()
+SE_B.dtf.annot.long$B0 = B0_vector
+SE_B.dtf.annot.long$deviance.sqrt = estim_mu$deviance.sqrt
+
+ggplot(SE_B.dtf.annot.long, aes(x = B0, y = SE_beta, col= annot )) + geom_point(alpha = 0.1, na.rm = T) + facet_grid(~origin)
+
+
+SE_B.dtf.annot.long_B0= SE_B.dtf.annot.long %>% filter(origin == "B0") %>% filter(!is.na(annot))
+
+ggplot( SE_B.dtf.annot.long_B0, aes(x = B0, y = 2^deviance.sqrt, col = annot)) + geom_point(alpha = 0.1, na.rm = TRUE) + scale_y_log10()
+
+
+```
+
+
+
+
+```{undefined eval=FALSE, include=FALSE}
+hist(dds.mcols$dispFit)
+hist(log_qij, )
+hist(log10(alpha.input$alpha))
+max(dds.mcols$dispersion, na.rm = T)
+hist(log(dds.mcols$dispersion))
+max(dds.mcols$deviance, na.rm = T)
+hist(dds.mcols$deviance)
+max(dds.mcols$dispFit, na.rm = T)
+DESeq2::design(dds)
+
+fitted.common.scale = t(t(dds@assays@data$mu)/dds$sizeFactor)
+
+hist(t(t(dds@assays@data$mu)/dds$sizeFactor))
+hist(residual_deseq)
+max(residual_deseq, na.rm = T)
+residual_deseq = (DESeq2::counts(dds, normalized=TRUE) - fitted.common.scale )
+
+
+w =DESeq2::nbinomWaldTest(dds)
+w@dispersionFunction()
+#S4Vectors::assays(dds)[["mu"]]
+
+
+vst = DESeq2::varianceStabilizingTransformation(dds)
+vst@assays@data$
+```
+
+
+
+https://support.bioconductor.org/p/123305/
+https://support.bioconductor.org/p/60567/
+http://bioconductor.org/packages/release/bioc/vignettes/RUVSeq/inst/doc/RUVSeq.pdf
+https://bioinformatics-core-shared-training.github.io/RNAseq-R/slides/LinearModels.pdf
+
+
+```{r eval=FALSE, include=FALSE}
+#install.packages("RUVSeq")
+#BiocManager::install("RUVSeq")
+library(RUVSeq)
+
+mm
+y <- DGEList(counts=counts(tabl_cnts), group=x)
+y <- calcNormFactors(y, method="upperquartile")
+y <- estimateGLMCommonDisp(y, design)
+y <- estimateGLMTagwiseDisp(y, design)
+fit <- glmFit(y, design)
+res <- residuals(fit, type="deviance")
+
+
+```
+
+
+
+
diff --git a/results/2022-04-22_lowCnts_effect.Rmd b/results/v1/2022-04-22_lowCnts_effect.Rmd
similarity index 100%
rename from results/2022-04-22_lowCnts_effect.Rmd
rename to results/v1/2022-04-22_lowCnts_effect.Rmd
diff --git a/results/2022-04-24_epsilon_effect.Rmd b/results/v1/2022-04-24_epsilon_effect.Rmd
similarity index 100%
rename from results/2022-04-24_epsilon_effect.Rmd
rename to results/v1/2022-04-24_epsilon_effect.Rmd
diff --git a/results/2022-04-28_beta_SE_filter.Rmd b/results/v1/2022-04-28_beta_SE_filter.Rmd
similarity index 100%
rename from results/2022-04-28_beta_SE_filter.Rmd
rename to results/v1/2022-04-28_beta_SE_filter.Rmd
diff --git a/results/2022-05-22_alpha_effect.Rmd b/results/v1/2022-05-22_alpha_effect.Rmd
similarity index 96%
rename from results/2022-05-22_alpha_effect.Rmd
rename to results/v1/2022-05-22_alpha_effect.Rmd
index 2e78ece..d553f17 100644
--- a/results/2022-05-22_alpha_effect.Rmd
+++ b/results/v1/2022-05-22_alpha_effect.Rmd
@@ -1,329 +1,329 @@
----
-title: "Alpha effect"
-date: '2022-04-21'
-output:
- html_document:
-
-
-css:
- - css/template.css
-
----
-
-
-## Required
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-library(htrsim)
-library(tidyverse)
-library(reshape2)
-```
-
-
-
-## Worklow
-
-Using public libraries (from BioProject PRJNA675209b - chinese paper), and an usual RNA-seq pipeline, we build actual RNA-seq counts per genes for 3 genotypes and 2 environments.<br>
-<br>
-Using htrsim (in particular DESEQ2) and this count table, we are able to estimate $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$.
-
-
-##### a. Input
-
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## Import & reshape table counts
-fn = system.file("extdata/", "public_tablCnts.tsv", package = "htrsim")
-tabl_cnts <- read.table(file = fn, header = TRUE)
-rownames(tabl_cnts) <- tabl_cnts$gene_id
-tabl_cnts <- tabl_cnts %>% select(-gene_id)##suppr colonne GeneID
-tabl_cnts <- tabl_cnts %>% select(-gene_name) ##suppr colonne GeneName
-
-## import design of bioProject
-fn = system.file("extdata/", "public_bioDesign.csv", package = "htrsim")
-bioDesign <- read.table(file = fn, header = T, sep = ';')
-```
-
-##### Launch DESEQ
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-dds = run.deseq(tabl_cnts = tabl_cnts, bioDesign = bioDesign)
-```
-
-DESEQ returns a dds object which contains many, many things ... <br>
-In particular it contains the $\beta$ coefficients. <br>
-<br>
-You can access to beta coefficients using:
-
-```{r}
-dds.mcols = S4Vectors::mcols(dds,use.names=TRUE)
-```
-
-##### mu estimation
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## Model matrix per samples
-mm <- model.matrix(~genotype + env + genotype:env, bioDesign)
-
-## Input estimation
-estim_mu = estim.mu(dds, mm)
-mu.input = estim_mu$mu
-```
-
-d. $\alpha_{i}$
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-alpha.input = estim.alpha(dds)
-```
-
-
-```{r}
-hist(alpha.input$alpha)
-maximum_alpha = max(alpha.input$alpha, na.rm = T)
-maximum_alpha
-min(alpha.input$alpha, na.rm = T)
-
-```
-
-variance = n(1-alpha)/alpha^2.
-
-```{r}
-n = 1000
-alpha = 2.5
-variance = n*(1 - alpha)/alpha^2
-variance
-hist(rnbinom(n=1000, mu = 400, size = alpha))
-```
-
-```{r}
-n = 1000
-alpha = 0.01
-variance = n*(1 - alpha)/alpha^2
-variance
-hist(rnbinom(n=1000, mu = 400, size = alpha))
-```
-
-
-```{r}
-## filter on alphz
-#alpha.input.filtered = alpha.input %>% filter(gene_id %in% beta_filtered.short$gene_id)
-#alpha.input.filtered = alpha.input %>% filter(alpha > 10)
-
-
-### FIX ALPHA
-N_genes = length(alpha.input$gene_id)
-alpha.input$alpha = runif(n= N_genes, min = 0.1, max = maximum_alpha)
-```
-
-Knowing $\alpha_i$ and $\mu_{ij]}$ for each gene and each condition given by the BioProject PRJNA675209b - chinese paper.
-We are now able to simulate $K_{ij}$ for each gene and each condition given by the BioProject PRJNA675209b.
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## filter mu based on alpha filter
-#mu.input.filtered = mu.input %>% filter(gene_id %in% alpha.input.filtered$gene_id )
-
-# Setup simulation
-input = reshape_input2setup(mu.dtf = mu.input, alpha.dtf = alpha.input, average_rep = FALSE)
-
-#input$gene_id
-setup.simulation <- setup_countGener(bioSample_id = input$bioSample_id,
- n_rep = 1,
- alpha = input$alpha,
- gene_id = input$gene_id,
- mu = input$mu)
-
-#setup.simulation %>% dim()
-# Simulate counts
-htrs <- generate_counts(setup.simulation)
-
-```
-
-##### Evaluation
-
-```{r message=FALSE, warning=FALSE, include=TRUE}
-k_ij.simu = htrs %>% select(-gene_id) %>% flatten() %>% unlist()
-k_ij.actual = tabl_cnts %>% filter(rownames(.) %in% htrs$gene_id ) %>% flatten() %>% unlist()
-
-df = cbind(k_ij.actual, k_ij.simu) %>% reshape2::melt(., value.name = "k_ij", variable.name = "origin")
-df$origin = df$Var2
-df = df %>% select(-Var2)
-max_k_ij.simu = df %>% filter(origin == "k_ij.simu") %>% select(k_ij) %>% max()
-max_k_ij.actual = df %>% filter(origin == "k_ij.actual") %>% select(k_ij) %>% max()
-
-ggplot(df, aes(x = k_ij, fill= origin )) + geom_density(bins = 100, alpha = 0.5) +
- geom_vline(xintercept = max_k_ij.actual, col = "#F8766D" ) +
- geom_vline(xintercept = (max_k_ij.simu), col= "#00BFC4" ) +
- scale_x_log10()
-```
-
-
-
-## Without Epsilon
-
-##### mu estimation
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## Model matrix per samples
-mm <- model.matrix(~genotype + env + genotype:env, bioDesign)
-
-## Input estimation
-estim_mu = estim.mu(dds, mm, epsilon = FALSE)
-mu.input = estim_mu$mu
-```
-
-d. $\alpha_{i}$
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-alpha.input = estim.alpha(dds)
-```
-
-```{r}
-## filter on alphz
-#alpha.input.filtered = alpha.input %>% filter(gene_id %in% beta_filtered.short$gene_id)
-#alpha.input.filtered = alpha.input %>% filter(alpha > 10)
-
-
-### FIX ALPHA
-N_genes = length(alpha.input$gene_id)
-alpha.input$alpha = runif(n= N_genes, min = 0.1, max = maximum_alpha)
-#alpha.input$alpha = rnbinom(n = N_genes, mu = mean(alpha.input$alpha), size = 20)
-#hist(alpha.input$alpha)
-```
-
-Knowing $\alpha_i$ and $\mu_{ij]}$ for each gene and each condition given by the BioProject PRJNA675209b - chinese paper.
-We are now able to simulate $K_{ij}$ for each gene and each condition given by the BioProject PRJNA675209b.
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## filter mu based on alpha filter
-#mu.input.filtered = mu.input %>% filter(gene_id %in% alpha.input.filtered$gene_id )
-
-# Setup simulation
-input = reshape_input2setup(mu.dtf = mu.input, alpha.dtf = alpha.input, average_rep = FALSE)
-#input$gene_id %>% length()
-#input$gene_id
-setup.simulation <- setup_countGener(bioSample_id = input$bioSample_id,
- n_rep = 1,
- alpha = input$alpha,
- gene_id = input$gene_id,
- mu = input$mu)
-
-#setup.simulation %>% dim()
-# Simulate counts
-htrs <- generate_counts(setup.simulation)
-
-```
-
-##### Evaluation
-
-```{r message=FALSE, warning=FALSE, include=TRUE}
-k_ij.simu = htrs %>% select(-gene_id) %>% flatten() %>% unlist()
-k_ij.actual = tabl_cnts %>% filter(rownames(.) %in% htrs$gene_id ) %>% flatten() %>% unlist()
-
-df = cbind(k_ij.actual, k_ij.simu) %>% reshape2::melt(., value.name = "k_ij", variable.name = "origin")
-df$origin = df$Var2
-df = df %>% select(-Var2)
-max_k_ij.simu = df %>% filter(origin == "k_ij.simu") %>% select(k_ij) %>% max()
-max_k_ij.actual = df %>% filter(origin == "k_ij.actual") %>% select(k_ij) %>% max()
-
-ggplot(df, aes(x = k_ij, fill= origin )) + geom_density(bins = 100, alpha = 0.5) +
- geom_vline(xintercept = max_k_ij.actual, col = "#F8766D" ) +
- geom_vline(xintercept = (max_k_ij.simu), col= "#00BFC4" ) +
- scale_x_log10()
-```
-
-
-```{r}
-#colnames(htrs)
-##colnames(tabl_cnts)
-#tabl_cnts.geneSimulate <- tabl_cnts %>% filter(rownames(.) %in% htrs$gene_id )
-#actual = tabl_cnts.geneSimulate$WT_control_rep1
-#simu = htrs$WT_control_rep1_1
-ggplot() + geom_point(aes(x=k_ij.actual, y = k_ij.simu), alpha=0.5) + geom_abline(intercept = 0, slope = 1)
-#abline(a=0, b=1)
-```
-
-
-## Re-estimate beta
-
-##### Design
-
-```{r}
-sample = htrs %>% select(where(is.numeric)) %>% colnames()
-genotype = htrs %>%
- select(where(is.numeric)) %>%
- colnames() %>%
- map(., ~str_split(.,pattern = '_', simplify = T)[1]) %>% unlist()
-
-env = htrs %>%
- select(where(is.numeric)) %>%
- colnames() %>%
- map(., ~str_split(.,pattern = '_', simplify = T)[2]) %>% unlist()
-
-designSimu = cbind(sample, env, genotype) %>% data.frame()
-```
-
-
-##### Launch DESEQ
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-## RESHAPE HTRS
-htrs.reshape = htrs
-rownames(htrs.reshape) = htrs.reshape$gene_id
-htrs.reshape = htrs.reshape %>% select(-gene_id)
-########### LAUNCH DESEQ #############
-## Design model - specify reference
-designSimu$genotype <- factor(x = designSimu$genotype,levels = c('WT','Msn2D', 'Msn4D'))
-designSimu$env <- factor(x = designSimu$env,levels = c('control', 'KCl'))
-
-
-k_ij.simulate = htrs.reshape
-
-## DESEQ standard analysis
-dds_simu = DESeq2::DESeqDataSetFromMatrix( countData = k_ij.simulate , colData = designSimu , design = ~ genotype + env + genotype:env)
-
-max(htrs.reshape)
-.Machine$integer.max
-
-dds_simu <- DESeq2::DESeq(dds_simu)
-```
-
-DESEQ returns a dds object which contains many, many things ... <br>
-In particular it contains the $\beta$ coefficients. <br>
-<br>
-You can access to beta coefficients using:
-
-```{r}
-dds_simu.mcols = S4Vectors::mcols(dds_simu,use.names=TRUE)
-#rownames(tabl_cnts)
-```
-
-
-## Evaluation of beta inference
-
-
-```{r}
-## BETA INPUT
-beta_input = estim_mu$beta.matrix %>% as.data.frame()
-idx_nonNA = which(!is.na(beta_input$B0))
-beta_input = beta_input[idx_nonNA,]
-beta_input$gene_id = input$gene_id
-## BETA SIMU
-B0 <- dds_simu.mcols$Intercept
-B1 <- dds_simu.mcols$genotype_Msn2D_vs_WT
-B2 <- dds_simu.mcols$genotype_Msn4D_vs_WT
-B3 <- dds_simu.mcols$env_KCl_vs_control
-B4 <- dds_simu.mcols$genotypeMsn2D.envKCl
-B5 <- dds_simu.mcols$genotypeMsn4D.envKCl
-
-
-beta.dtf = cbind(B0, B1,B2,B3,B4,B5) %>% as.data.frame()
-beta.dtf$gene_id = input$gene_id
-beta.dtf$origin = "inference"
-beta_input$origin = "input"
-dtf.merged = rbind(beta.dtf, beta_input)
-
-dtf.merged.long.tmp = dtf.merged %>% reshape2::melt(., value.name = "value", variable.name= "beta")#, variable.name = "origin")()
-dtf.merged.long = dtf.merged.long.tmp %>% reshape2::dcast(., gene_id + beta ~ origin)
-
-ggplot(dtf.merged.long) + geom_point(aes(x=input, y = inference),alpha =0.1)+ geom_abline(intercept = 0, slope = 1) + facet_grid(~beta)
-```
+---
+title: "Alpha effect"
+date: '2022-04-21'
+output:
+ html_document:
+
+
+css:
+ - css/template.css
+
+---
+
+
+## Required
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+library(htrsim)
+library(tidyverse)
+library(reshape2)
+```
+
+
+
+## Worklow
+
+Using public libraries (from BioProject PRJNA675209b - chinese paper), and an usual RNA-seq pipeline, we build actual RNA-seq counts per genes for 3 genotypes and 2 environments.<br>
+<br>
+Using htrsim (in particular DESEQ2) and this count table, we are able to estimate $\beta_{G}$, $\beta_{E}$ and $\beta_{G*E}$.
+
+
+##### a. Input
+
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## Import & reshape table counts
+fn = system.file("extdata/", "public_tablCnts.tsv", package = "htrsim")
+tabl_cnts <- read.table(file = fn, header = TRUE)
+rownames(tabl_cnts) <- tabl_cnts$gene_id
+tabl_cnts <- tabl_cnts %>% select(-gene_id)##suppr colonne GeneID
+tabl_cnts <- tabl_cnts %>% select(-gene_name) ##suppr colonne GeneName
+
+## import design of bioProject
+fn = system.file("extdata/", "public_bioDesign.csv", package = "htrsim")
+bioDesign <- read.table(file = fn, header = T, sep = ';')
+```
+
+##### Launch DESEQ
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+dds = run.deseq(tabl_cnts = tabl_cnts, bioDesign = bioDesign)
+```
+
+DESEQ returns a dds object which contains many, many things ... <br>
+In particular it contains the $\beta$ coefficients. <br>
+<br>
+You can access to beta coefficients using:
+
+```{r}
+dds.mcols = S4Vectors::mcols(dds,use.names=TRUE)
+```
+
+##### mu estimation
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## Model matrix per samples
+mm <- model.matrix(~genotype + env + genotype:env, bioDesign)
+
+## Input estimation
+estim_mu = estim.mu(dds, mm)
+mu.input = estim_mu$mu
+```
+
+d. $\alpha_{i}$
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+alpha.input = estim.alpha(dds)
+```
+
+
+```{r}
+hist(alpha.input$alpha)
+maximum_alpha = max(alpha.input$alpha, na.rm = T)
+maximum_alpha
+min(alpha.input$alpha, na.rm = T)
+
+```
+
+variance = n(1-alpha)/alpha^2.
+
+```{r}
+n = 1000
+alpha = 2.5
+variance = n*(1 - alpha)/alpha^2
+variance
+hist(rnbinom(n=1000, mu = 400, size = alpha))
+```
+
+```{r}
+n = 1000
+alpha = 0.01
+variance = n*(1 - alpha)/alpha^2
+variance
+hist(rnbinom(n=1000, mu = 400, size = alpha))
+```
+
+
+```{r}
+## filter on alphz
+#alpha.input.filtered = alpha.input %>% filter(gene_id %in% beta_filtered.short$gene_id)
+#alpha.input.filtered = alpha.input %>% filter(alpha > 10)
+
+
+### FIX ALPHA
+N_genes = length(alpha.input$gene_id)
+alpha.input$alpha = runif(n= N_genes, min = 0.1, max = maximum_alpha)
+```
+
+Knowing $\alpha_i$ and $\mu_{ij]}$ for each gene and each condition given by the BioProject PRJNA675209b - chinese paper.
+We are now able to simulate $K_{ij}$ for each gene and each condition given by the BioProject PRJNA675209b.
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## filter mu based on alpha filter
+#mu.input.filtered = mu.input %>% filter(gene_id %in% alpha.input.filtered$gene_id )
+
+# Setup simulation
+input = reshape_input2setup(mu.dtf = mu.input, alpha.dtf = alpha.input, average_rep = FALSE)
+
+#input$gene_id
+setup.simulation <- setup_countGener(bioSample_id = input$bioSample_id,
+ n_rep = 1,
+ alpha = input$alpha,
+ gene_id = input$gene_id,
+ mu = input$mu)
+
+#setup.simulation %>% dim()
+# Simulate counts
+htrs <- generate_counts(setup.simulation)
+
+```
+
+##### Evaluation
+
+```{r message=FALSE, warning=FALSE, include=TRUE}
+k_ij.simu = htrs %>% select(-gene_id) %>% flatten() %>% unlist()
+k_ij.actual = tabl_cnts %>% filter(rownames(.) %in% htrs$gene_id ) %>% flatten() %>% unlist()
+
+df = cbind(k_ij.actual, k_ij.simu) %>% reshape2::melt(., value.name = "k_ij", variable.name = "origin")
+df$origin = df$Var2
+df = df %>% select(-Var2)
+max_k_ij.simu = df %>% filter(origin == "k_ij.simu") %>% select(k_ij) %>% max()
+max_k_ij.actual = df %>% filter(origin == "k_ij.actual") %>% select(k_ij) %>% max()
+
+ggplot(df, aes(x = k_ij, fill= origin )) + geom_density(bins = 100, alpha = 0.5) +
+ geom_vline(xintercept = max_k_ij.actual, col = "#F8766D" ) +
+ geom_vline(xintercept = (max_k_ij.simu), col= "#00BFC4" ) +
+ scale_x_log10()
+```
+
+
+
+## Without Epsilon
+
+##### mu estimation
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## Model matrix per samples
+mm <- model.matrix(~genotype + env + genotype:env, bioDesign)
+
+## Input estimation
+estim_mu = estim.mu(dds, mm, epsilon = FALSE)
+mu.input = estim_mu$mu
+```
+
+d. $\alpha_{i}$
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+alpha.input = estim.alpha(dds)
+```
+
+```{r}
+## filter on alphz
+#alpha.input.filtered = alpha.input %>% filter(gene_id %in% beta_filtered.short$gene_id)
+#alpha.input.filtered = alpha.input %>% filter(alpha > 10)
+
+
+### FIX ALPHA
+N_genes = length(alpha.input$gene_id)
+alpha.input$alpha = runif(n= N_genes, min = 0.1, max = maximum_alpha)
+#alpha.input$alpha = rnbinom(n = N_genes, mu = mean(alpha.input$alpha), size = 20)
+#hist(alpha.input$alpha)
+```
+
+Knowing $\alpha_i$ and $\mu_{ij]}$ for each gene and each condition given by the BioProject PRJNA675209b - chinese paper.
+We are now able to simulate $K_{ij}$ for each gene and each condition given by the BioProject PRJNA675209b.
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## filter mu based on alpha filter
+#mu.input.filtered = mu.input %>% filter(gene_id %in% alpha.input.filtered$gene_id )
+
+# Setup simulation
+input = reshape_input2setup(mu.dtf = mu.input, alpha.dtf = alpha.input, average_rep = FALSE)
+#input$gene_id %>% length()
+#input$gene_id
+setup.simulation <- setup_countGener(bioSample_id = input$bioSample_id,
+ n_rep = 1,
+ alpha = input$alpha,
+ gene_id = input$gene_id,
+ mu = input$mu)
+
+#setup.simulation %>% dim()
+# Simulate counts
+htrs <- generate_counts(setup.simulation)
+
+```
+
+##### Evaluation
+
+```{r message=FALSE, warning=FALSE, include=TRUE}
+k_ij.simu = htrs %>% select(-gene_id) %>% flatten() %>% unlist()
+k_ij.actual = tabl_cnts %>% filter(rownames(.) %in% htrs$gene_id ) %>% flatten() %>% unlist()
+
+df = cbind(k_ij.actual, k_ij.simu) %>% reshape2::melt(., value.name = "k_ij", variable.name = "origin")
+df$origin = df$Var2
+df = df %>% select(-Var2)
+max_k_ij.simu = df %>% filter(origin == "k_ij.simu") %>% select(k_ij) %>% max()
+max_k_ij.actual = df %>% filter(origin == "k_ij.actual") %>% select(k_ij) %>% max()
+
+ggplot(df, aes(x = k_ij, fill= origin )) + geom_density(bins = 100, alpha = 0.5) +
+ geom_vline(xintercept = max_k_ij.actual, col = "#F8766D" ) +
+ geom_vline(xintercept = (max_k_ij.simu), col= "#00BFC4" ) +
+ scale_x_log10()
+```
+
+
+```{r}
+#colnames(htrs)
+##colnames(tabl_cnts)
+#tabl_cnts.geneSimulate <- tabl_cnts %>% filter(rownames(.) %in% htrs$gene_id )
+#actual = tabl_cnts.geneSimulate$WT_control_rep1
+#simu = htrs$WT_control_rep1_1
+ggplot() + geom_point(aes(x=k_ij.actual, y = k_ij.simu), alpha=0.5) + geom_abline(intercept = 0, slope = 1)
+#abline(a=0, b=1)
+```
+
+
+## Re-estimate beta
+
+##### Design
+
+```{r}
+sample = htrs %>% select(where(is.numeric)) %>% colnames()
+genotype = htrs %>%
+ select(where(is.numeric)) %>%
+ colnames() %>%
+ map(., ~str_split(.,pattern = '_', simplify = T)[1]) %>% unlist()
+
+env = htrs %>%
+ select(where(is.numeric)) %>%
+ colnames() %>%
+ map(., ~str_split(.,pattern = '_', simplify = T)[2]) %>% unlist()
+
+designSimu = cbind(sample, env, genotype) %>% data.frame()
+```
+
+
+##### Launch DESEQ
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+## RESHAPE HTRS
+htrs.reshape = htrs
+rownames(htrs.reshape) = htrs.reshape$gene_id
+htrs.reshape = htrs.reshape %>% select(-gene_id)
+########### LAUNCH DESEQ #############
+## Design model - specify reference
+designSimu$genotype <- factor(x = designSimu$genotype,levels = c('WT','Msn2D', 'Msn4D'))
+designSimu$env <- factor(x = designSimu$env,levels = c('control', 'KCl'))
+
+
+k_ij.simulate = htrs.reshape
+
+## DESEQ standard analysis
+dds_simu = DESeq2::DESeqDataSetFromMatrix( countData = k_ij.simulate , colData = designSimu , design = ~ genotype + env + genotype:env)
+
+max(htrs.reshape)
+.Machine$integer.max
+
+dds_simu <- DESeq2::DESeq(dds_simu)
+```
+
+DESEQ returns a dds object which contains many, many things ... <br>
+In particular it contains the $\beta$ coefficients. <br>
+<br>
+You can access to beta coefficients using:
+
+```{r}
+dds_simu.mcols = S4Vectors::mcols(dds_simu,use.names=TRUE)
+#rownames(tabl_cnts)
+```
+
+
+## Evaluation of beta inference
+
+
+```{r}
+## BETA INPUT
+beta_input = estim_mu$beta.matrix %>% as.data.frame()
+idx_nonNA = which(!is.na(beta_input$B0))
+beta_input = beta_input[idx_nonNA,]
+beta_input$gene_id = input$gene_id
+## BETA SIMU
+B0 <- dds_simu.mcols$Intercept
+B1 <- dds_simu.mcols$genotype_Msn2D_vs_WT
+B2 <- dds_simu.mcols$genotype_Msn4D_vs_WT
+B3 <- dds_simu.mcols$env_KCl_vs_control
+B4 <- dds_simu.mcols$genotypeMsn2D.envKCl
+B5 <- dds_simu.mcols$genotypeMsn4D.envKCl
+
+
+beta.dtf = cbind(B0, B1,B2,B3,B4,B5) %>% as.data.frame()
+beta.dtf$gene_id = input$gene_id
+beta.dtf$origin = "inference"
+beta_input$origin = "input"
+dtf.merged = rbind(beta.dtf, beta_input)
+
+dtf.merged.long.tmp = dtf.merged %>% reshape2::melt(., value.name = "value", variable.name= "beta")#, variable.name = "origin")()
+dtf.merged.long = dtf.merged.long.tmp %>% reshape2::dcast(., gene_id + beta ~ origin)
+
+ggplot(dtf.merged.long) + geom_point(aes(x=input, y = inference),alpha =0.1)+ geom_abline(intercept = 0, slope = 1) + facet_grid(~beta)
+```
diff --git a/results/2022-06-24_subsampling.Rmd b/results/v1/2022-06-24_subsampling.Rmd
similarity index 100%
rename from results/2022-06-24_subsampling.Rmd
rename to results/v1/2022-06-24_subsampling.Rmd
diff --git a/results/2022-06-28_kallistoEffect.Rmd b/results/v1/2022-06-28_kallistoEffect.Rmd
similarity index 100%
rename from results/2022-06-28_kallistoEffect.Rmd
rename to results/v1/2022-06-28_kallistoEffect.Rmd
diff --git a/results/2022_06_08_investigation.Rmd b/results/v1/2022_06_08_investigation.Rmd
similarity index 100%
rename from results/2022_06_08_investigation.Rmd
rename to results/v1/2022_06_08_investigation.Rmd
diff --git a/results/SVA.R b/results/v1/SVA.R
similarity index 100%
rename from results/SVA.R
rename to results/v1/SVA.R
diff --git a/results/SVA.Rmd b/results/v1/SVA.Rmd
similarity index 100%
rename from results/SVA.Rmd
rename to results/v1/SVA.Rmd
diff --git a/results/css/footer.css b/results/v1/css/footer.css
similarity index 100%
rename from results/css/footer.css
rename to results/v1/css/footer.css
diff --git a/results/css/template.css b/results/v1/css/template.css
similarity index 100%
rename from results/css/template.css
rename to results/v1/css/template.css
diff --git a/results/lab-meeting_preparation.Rmd b/results/v1/lab-meeting_preparation.Rmd
similarity index 96%
rename from results/lab-meeting_preparation.Rmd
rename to results/v1/lab-meeting_preparation.Rmd
index 8d8498e..8bc7cae 100644
--- a/results/lab-meeting_preparation.Rmd
+++ b/results/v1/lab-meeting_preparation.Rmd
@@ -1,61 +1,61 @@
----
-title: "lab-meeting"
-output:
- html_document
----
-
-
-```{r}
-library(tidyverse)
-
-
-x=runif(100,0,100)
-epsilon.simu = rnorm(100,0,16)
-y = 150 + 4*x + epsilon.simu
-mysimu1 = cbind(x, y) %>% data.frame()
-mysimu1$epsilon = "eps ~ N(mean = 0, sd = 16)"
-
-ggplot(mysimu1) + geom_point(aes(x=x, y=y),col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100)
-
-ggplot(mysimu1) + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100)
-
-epsilon.simu= epsilon.simu %>% data.frame() %>% mutate(., epsilon = .) %>% select(-.)
-ggplot(epsilon.simu) + geom_density(aes(x=epsilon))
-
-```
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-
-x=runif(100,0,100)
-epsilon.simu = rnorm(1000,0,160)
-y = 150 + 4*x + epsilon.simu
-mysimu2 = cbind(x, y) %>% data.frame()
-mysimu2$epsilon = "eps ~ N(mean = 0, sd = 160)"
-ggplot(mysimu2) + geom_point(aes(x=x, y=y), col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100)
-epsilon.simu= epsilon.simu %>% data.frame() %>% mutate(., epsilon = .) %>% select(-.)
-ggplot(epsilon.simu) + geom_density(aes(x=epsilon))
-
-```
-
-```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
-dtf.merged <- rbind(mysimu1, mysimu2)
-
-ggplot(dtf.merged) + geom_point(aes(x=x, y=y), col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100) + facet_grid(~epsilon)
-
-```
-
-```{r}
-dtf.merged <- rbind(mysimu1, mysimu2)
-
-ggplot(dtf.merged) + geom_point(aes(x=x, y=y), col= "#00BFC4") + ylim(100, 500) + xlim(0, 100) + facet_grid(~epsilon)
-```
-
-```{r}
-x=runif(10,40,70)
-epsilon.simu = rnorm(10,0,100)
-y = 150 + 4*x + epsilon.simu
-mysimu1 = cbind(x, y) %>% data.frame()
-mysimu1$epsilon = "eps ~ N(mean = 0, sd = 16)"
-
-ggplot(mysimu1) + geom_point(aes(x=x, y=y),col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D")
-```
+---
+title: "lab-meeting"
+output:
+ html_document
+---
+
+
+```{r}
+library(tidyverse)
+
+
+x=runif(100,0,100)
+epsilon.simu = rnorm(100,0,16)
+y = 150 + 4*x + epsilon.simu
+mysimu1 = cbind(x, y) %>% data.frame()
+mysimu1$epsilon = "eps ~ N(mean = 0, sd = 16)"
+
+ggplot(mysimu1) + geom_point(aes(x=x, y=y),col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100)
+
+ggplot(mysimu1) + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100)
+
+epsilon.simu= epsilon.simu %>% data.frame() %>% mutate(., epsilon = .) %>% select(-.)
+ggplot(epsilon.simu) + geom_density(aes(x=epsilon))
+
+```
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+
+x=runif(100,0,100)
+epsilon.simu = rnorm(1000,0,160)
+y = 150 + 4*x + epsilon.simu
+mysimu2 = cbind(x, y) %>% data.frame()
+mysimu2$epsilon = "eps ~ N(mean = 0, sd = 160)"
+ggplot(mysimu2) + geom_point(aes(x=x, y=y), col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100)
+epsilon.simu= epsilon.simu %>% data.frame() %>% mutate(., epsilon = .) %>% select(-.)
+ggplot(epsilon.simu) + geom_density(aes(x=epsilon))
+
+```
+
+```{r message=FALSE, warning=FALSE, include=TRUE, results="hide"}
+dtf.merged <- rbind(mysimu1, mysimu2)
+
+ggplot(dtf.merged) + geom_point(aes(x=x, y=y), col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D") + ylim(100, 500) + xlim(0, 100) + facet_grid(~epsilon)
+
+```
+
+```{r}
+dtf.merged <- rbind(mysimu1, mysimu2)
+
+ggplot(dtf.merged) + geom_point(aes(x=x, y=y), col= "#00BFC4") + ylim(100, 500) + xlim(0, 100) + facet_grid(~epsilon)
+```
+
+```{r}
+x=runif(10,40,70)
+epsilon.simu = rnorm(10,0,100)
+y = 150 + 4*x + epsilon.simu
+mysimu1 = cbind(x, y) %>% data.frame()
+mysimu1$epsilon = "eps ~ N(mean = 0, sd = 16)"
+
+ggplot(mysimu1) + geom_point(aes(x=x, y=y),col= "#00BFC4") + geom_abline(intercept = 150, slope = 4, col= "#F8766D")
+```
--
GitLab