Création du Rmd

R/Method_on_realtree.Rmd

@@ -0,0 +1,111 @@
+---
+title: Méthode sur un vrai arbre
+---
+
+Ici nous appliquons les méthodes implémentées sur l'arbre de @chen2019.
+
+```{r import_modules, echo = FALSE}
+# Utiliser data.trans, ligne 883 voir RMD
+require(phylotools)
+require(phytools)
+require(phylolm)
+require(limma)
+require(edgeR)
+require(here)
+require(ggplot2)
+
+source("R/utils.R")
+```
+
+```{r import_donnees, echo = FALSE}
+### Data import
+cdata <- readRDS(here("data", "data_TER", "data", "chen2019_rodents_cpd.rds"))
+is.valid <- compcodeR:::check_phyloCompData(cdata)
+if (!(is.valid == TRUE)) stop("Not a valid phyloCompData object.")
+
+# Design
+design_formula <- as.formula(~condition)
+design_data <- compcodeR:::sample.annotations(cdata)[, "condition", drop = FALSE]
+design_data$condition <- factor(design_data$condition)
+design <- model.matrix(design_formula, design_data)
+
+# Normalisation
+nf <- edgeR::calcNormFactors(compcodeR:::count.matrix(cdata) / compcodeR:::length.matrix(cdata), method = "TMM")
+lib.size <- colSums(compcodeR:::count.matrix(cdata) / compcodeR:::length.matrix(cdata)) * nf
+data.norm <- sweep((compcodeR:::count.matrix(cdata) + 0.5) / compcodeR:::length.matrix(cdata), 2, lib.size + 1, "/")
+data.norm <- data.norm * 1e6
+
+# Transformation
+data.trans <- log2(data.norm)
+rownames(data.trans) <- rownames(compcodeR:::count.matrix(cdata))
+```
+
+```{r calcul_pvaleurs, echo = FALSE}
+### Pvalues computation
+
+#  computing pvalues vec for all genes
+pvalue_vec_vanilla <- sapply(seq(1, nrow(data.trans)), function(row_id) {
+    trait <- data.trans[row_id, ]
+    fit_phylo <- phylolm(trait ~ design_data$condition, phy = cdata@tree, measurement_error = TRUE)
+    compute_vanilla_pvalue(fit_phylo)
+})
+
+pvalue_vec_vanilla <- setNames(pvalue_vec_vanilla, rownames(data.trans))
+
+pvalue_vec_vanilla_adj <- p.adjust(pvalue_vec_vanilla, method = "BH")
+
+pvalue_vec_satterthwaite <- sapply(seq(1, nrow(data.trans)), function(row_id) {
+    trait <- data.trans[row_id, ]
+    fit_phylo <- phylolm(trait ~ design_data$condition, phy = cdata@tree, measurement_error = TRUE)
+    compute_satterthwaite_pvalue(fit_phylo, tree = cdata@tree)
+})
+
+pvalue_vec_satterthwaite <- setNames(pvalue_vec_satterthwaite, rownames(data.trans))
+
+pvalue_vec_satterthwaite_adj <- p.adjust(pvalue_vec_satterthwaite, method = "BH")
+
+
+pvalue_vec_lrt <- sapply(seq(1, nrow(data.trans)), function(row_id) {
+    trait <- data.trans[row_id, ]
+    fit_phylo <- phylolm(trait ~ design_data$condition, phy = cdata@tree, measurement_error = TRUE)
+    compute_lrt_pvalue(fit_phylo, tree = cdata@tree)
+})
+
+pvalue_vec_lrt <- setNames(pvalue_vec_lrt, rownames(data.trans))
+pvalue_vec_lrt_adj <- p.adjust(pvalue_vec_lrt, method = "BH")
+
+# REML
+pvalue_vec_satterthwaite.REML <- sapply(seq(1, nrow(data.trans)), function(row_id) {
+    trait <- data.trans[row_id, ]
+    fit_phylo <- phylolm(trait ~ design_data$condition, phy = cdata@tree, measurement_error = TRUE)
+    compute_satterthwaite_pvalue(fit_phylo, tree = cdata@tree, REML = TRUE)
+})
+
+pvalue_vec_satterthwaite.REML <- setNames(pvalue_vec_satterthwaite.REML, rownames(data.trans))
+
+pvalue_vec_satterthwaite_adj.REML <- p.adjust(pvalue_vec_satterthwaite.REML, method = "BH")
+```
+
+```{r dataframe_plot, echo = FALSE}
+## Préparation du dataframe
+pvalues_dataframe <- data.frame(
+    genes = rep(rownames(data.trans), 8),
+    pvalues = c(pvalue_vec_vanilla, pvalue_vec_vanilla_adj, 
+        pvalue_vec_satterthwaite, pvalue_vec_satterthwaite_adj, 
+        pvalue_vec_lrt, pvalue_vec_lrt_adj, pvalue_vec_satterthwaite.REML, 
+        pvalue_vec_satterthwaite_adj.REML),
+    test_method = rep(c("Vanilla", "VanillaAdj", "Satterthwaite", 
+    "SatterthwaiteAdj", "LRT", "LRTAdj", "SatterthwaiteREML", 
+    "SatterthwaiteAdjREML"), each = nrow(data.trans))
+)
+```
+
+```{r graphique_all_pvalues, echo = FALSE}
+## Graphiques
+ggplot(pvalues_dataframe) +
+    aes(x = genes, y = pvalues, fill = test_method) +
+    geom_bar(stat = "identity", position = "dodge") +
+    facet_wrap(~test_method) +
+    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
+```
+

R/realTreeMethod.R

@@ -108,9 +108,9 @@ ggplot(pvalues_dataframe) +
     geom_bar(stat = "identity", position = "dodge") +
     facet_wrap(~test_method)
 
-# TODO utiliser UpSetR pour diagramme de Venn
+# DONE utiliser UpSetR pour diagramme de Venn
 require(UpSetR)
-upset(pvalues_dataframe_wide)
+upset(pvalues_dataframe_wide, mainbar.y.label = "Nombre de gènes en commun")
 
 # TODO comparer avec le package evemodel, twothetatest
 # Comparer avec OU lrt