diff --git a/Rnw/donnees-reelles.Rnw b/Rnw/donnees-reelles.Rnw index 5457e8f..87d3683 100644 --- a/Rnw/donnees-reelles.Rnw +++ b/Rnw/donnees-reelles.Rnw @@ -26,23 +26,24 @@ if (!all(necessary_packages %in% installed.packages())) { install.packages("remotes") remotes::install_gitlab("sandve-lab/evemodel") remotes::install_github("pbastide/compcodeR", ref = "phylolimma") -} else { - require(phylotools) - require(phytools) - require(phylolm) - require(limma) - require(edgeR) - require(here) - require(ggplot2) - require(dplyr) - require(tidyr) - require(UpSetR) - require(evemodel) - require(compcodeR) - require(mvSLOUCH) - require(patchwork) } +require(phylotools) +require(phytools) +require(phylolm) +require(limma) +require(edgeR) +require(here) +require(ggplot2) +require(dplyr) +require(tidyr) +require(UpSetR) +require(evemodel) +require(compcodeR) +require(mvSLOUCH) +require(patchwork) + + source(here("R","utils.R")) @ @@ -120,7 +121,7 @@ if (!file.exists(here("data",pvalues_data)) | !file.exists(here("data",pvalues_a 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, REML = TRUE, measurement_error = TRUE) - compute_satterthwaite_pvalue(fit_phylo, tree = cdata@tree, REML = TRUE) + compute_satterthwaite_pvalue(fit_phylo, tree = cdata@tree) }) pvalue_vec_satterthwaite.REML <- setNames(pvalue_vec_satterthwaite.REML, rownames(data.trans)) @@ -136,9 +137,9 @@ if (!file.exists(here("data",pvalues_data)) | !file.exists(here("data",pvalues_a pvalue_vec_satterthwaite_adj, pvalue_vec_lrt_adj, pvalue_vec_satterthwaite_adj.REML), - test_method = rep(c( "VanillaAdj", "VanillaAdjREML", - "SatterthwaiteAdj", "LRTAdj", - "SatterthwaiteAdjREML"), each = nrow(data.trans)) + test_method = rep(c( "ANOVA Phylo Ajustée", "ANOVA Phylo REML Ajustée", + "ANOVA Phylo Satterthwaite Ajustée", "LRT Ajusté", + "ANOVA Phylo Satterthwaite REML Ajustée"), each = nrow(data.trans)) ) pvalues_dataframe <- data.frame( @@ -148,9 +149,9 @@ if (!file.exists(here("data",pvalues_data)) | !file.exists(here("data",pvalues_a pvalue_vec_satterthwaite, pvalue_vec_lrt, pvalue_vec_satterthwaite.REML), - test_method = rep(c( "Vanilla", "VanillaREML", - "Satterthwaite", "LRT", - "SatterthwaiteREML"), each = nrow(data.trans)) + test_method = rep(c( "ANOVA Phylo", "ANOVA Phylo REML", + "ANOVA Phylo Satterthwaite", "LRT", + "ANOVA Phylo Satterthwaite REML"), each = nrow(data.trans)) ) pvalues_dataframe$test_method <- as.factor(pvalues_dataframe$test_method) @@ -167,11 +168,12 @@ if (!file.exists(here("data",pvalues_data)) | !file.exists(here("data",pvalues_a } @ -Nous appliquons les différentes méthodes que nous avons implémentés dans le +Nous appliquons les différentes méthodes que nous avons implémentées dans le code. -Ci-dessous la figure~\ref{fig:pval-methods} présente les p-values des -différentes méthodes. Il est important de noter que ce graphique présente les +Ci-dessous la figure~\ref{fig:pval-methods} présente les p-values ordonnées des +différentes méthodes. Il s'agit d'une visualisation classique pour les données +RNA-seq. Il est important de noter que ce graphique présente les p-values \emph{non ajustées}. \begin{figure}[H] @@ -193,7 +195,7 @@ p-values \emph{non ajustées}. axis_titles = "collect", tag_level = "new") + plot_annotation(title = "Selected genes by tested methods") @ - \caption{\emph{p-values} pour les différents tests} + \caption{\emph{p-values} ordonnées pour les différents tests} \label{fig:pval-methods} \end{figure} @@ -202,6 +204,7 @@ pour les test multiples, nommément la correction de \cite{benjaminiControllingF << wide_data, echo = FALSE>>= pvalues_adj_dataframe_wide <- pvalues_adj_dataframe %>% + filter(test_method != "ANOVA Phylo Satterthwaite Ajustée") %>% pivot_wider(id_cols = gene, names_from = test_method, values_from = selected) %>% @@ -211,26 +214,15 @@ pvalues_adj_dataframe_wide <- pvalues_adj_dataframe %>% Une fois ces corrections appliquées, nous allons comparer les gènes sélectionnés, c'est-à-dire différentiellement exprimés. - -Ces résultats sont présentés dans le diagramme de Venn (figure~\ref{fig:venn-all-methods}). - - On peut voir que la méthode de Satterthwaite sans REML a sélectionné énormément de gènes, -$\Sexpr{sum(pvalues_adj_dataframe[pvalues_adj_dataframe$test_method == "SatterthwaiteAdj",]$selected)}$ -comme étant différentiellement exprimés. -% TODO Comprendre la sur-sélection de Satterthwaite et pas de SatterthwaiteREML +$\Sexpr{sum(pvalues_adj_dataframe[pvalues_adj_dataframe$test_method == "ANOVA Phylo Satterthwaite Ajustée",]$selected)}$ +comme étant différentiellement exprimés. -\begin{figure}[H] - << upset_selection, echo = FALSE>>= - upset(pvalues_adj_dataframe_wide, - nsets = 5, - mainbar.y.label = "Nombre de gènes en commun", - sets.x.label = "Nombre de gènes sélectionnés") - @ - \caption{Diagramme de Venn comparant les gènes sélectionnés selon les méthodes} - \label{fig:venn-all-methods} -\end{figure} +Ce résultat n'étant pas biologiquement crédible, nous préférons ne pas +l'afficher dans le diagramme de Venn, figure~\ref{fig:venn-all-methods-eve}. +% DONE Préciser que Satterthwaite sur-sélectionne et le retirer des graphiques * +% pour cette raison \subsection*{EVEmodel} @@ -293,6 +285,10 @@ Dans l'article \cite{rohlfsPhylogeneticANOVAExpression2015}, les auteurs introduisent une méthode de détection des gènes différentiellement exprimés. Cette méthode est à l'heure actuelle très utilisée. +Son principe de fonctionnement suppose que les traits évoluent selon un +processus d'Ornstein-Uhlenbeck et le test réalisé est un \emph{Likelihood +Ratio test}. + \emph{Remarque :} La méthode a produit des \texttt{NA} pour certains gènes, d'après le message d'erreur, une optimisation n'a pas convergé. Ces gènes sont présentés dans le tableau~\ref{tab:na-evemodel}. @@ -307,6 +303,7 @@ méthodes. pvalueseve_dataframe <- rbind(pvalues_adj_dataframe, evemodel_dataframe) pvalueseve_dataframe_wide <- pvalueseve_dataframe %>% + filter(test_method != "ANOVA Phylo Satterthwaite Ajustée") %>% pivot_wider(id_cols = gene, names_from = test_method, values_from = selected, values_fill = 0) %>% diff --git a/Rnw/simulations-methodes.Rnw b/Rnw/simulations-methodes.Rnw index 102f11e..e8a61cd 100644 --- a/Rnw/simulations-methodes.Rnw +++ b/Rnw/simulations-methodes.Rnw @@ -31,13 +31,14 @@ l'\emph{erreur de première espèce} et la \emph{puissance} obtenue. \end{itemize} -Pour les simulations qui ne se font pas selon la phylogénie, nous nous attendons -à ce que l'ANOVA classique obtienne de bons résultats puisque c'est une -situation correspondant à l'application du modèle. +En sélectionnant des espèces de manière aléatoire, nous cassons la structure +induite par la phylogénie. Nous nous attendons donc à ce que l'ANOVA réalise de +meilleurs résultats en ne prenant pas en compte l'information phylogénétique. -Pour les simulations qui se font selon l'information de l'arbre phylogénétique, -nous nous attendons à ce que l'ANOVA phylogénétique parvienne à mieux prendre en -compte l'information apportée par la phylogénie et à démêler son effet.\\ +Pour les simulations avec des groupes respectant la structure de l'arbre +phylogénétique, nous nous attendons à ce que l'ANOVA phylogénétique +parvienne à mieux prendre en compte l'information apportée par la phylogénie et +à démêler son effet.\\ <<'modules-simulations', include = FALSE, eval=TRUE>>= necessary_simu <- c("ape", "remotes", "phylolm", "phylolimma", "phytools", @@ -67,7 +68,7 @@ nb_species <- 43 plot_group_on_tree <- function(tree, groups, ...) { plot(tree, ...) - tiplabels(bg = groups, pch = 21) + tiplabels(bg = groups, pch = 21, cex = 1.5) } tree <- read.tree(here("R","chen2019.tree")) # Normalising tree edge length @@ -113,7 +114,7 @@ Enfin pour que notre analyse soit reproductible nous fixons la graine à \Sexpr{seed}. \begin{figure}[H] - \begin{subfigure}[H]{0.45\textwidth} + \begin{subfigure}[H]{0.49\textwidth} \centering <<'plot-groupes-mus'>>= plot_group_on_tree(tree, group = group_mus_rat_vs_other) @@ -121,12 +122,12 @@ Enfin pour que notre analyse soit reproductible nous fixons la graine à \caption{Groupes \emph{Mus} et rats contre les autres} \label{fig:simu-groupes-mus} \end{subfigure} - \begin{subfigure}[H]{0.45\textwidth} + \begin{subfigure}[H]{0.49\textwidth} \centering <<'plot-groupes-random'>>= plot_group_on_tree(tree, group = random_groups) @ - \caption{Groupes respectant les proportions} + \caption{Groupes sélectionnés sans respect de la phylogénie.} \label{fig:simu-groupes-prop} \end{subfigure} \caption{Arbre et groupes pour les simulations} @@ -179,7 +180,7 @@ plot_list <- lapply(seq_len(length(heri)), function(idx) if (!file.exists(filename)) { sim <- N_simulation_typeI_power(N, groups_list = list(RatMus = group_mus_rat_vs_other, - Proportions = random_groups), + Sélectionnés = random_groups), tree = tree, base_values = c(0, snr * total_variance), sigma2_phylo = her * total_variance, @@ -204,6 +205,14 @@ Sur toutes les sous-figures de la figure~\ref{fig:simus-results}, les étiquette A présentent les erreurs de type I commises par les méthodes et les étiquettes B présentent les puissances des mêmes méthodes. +L'erreur de type I est particulièrement importante à contrôler, en effet elle +indique le nombre de faux positifs et l'on veut pouvoir en déterminer le seuil +$\alpha$ avec comme seuil classique $0.05$. + +D'après nos simulations, les méthodes utilisant le REML contrôle toujours mieux +l'erreur de première espèce que les méthodes utilisant le maximum de +vraisemblance. + TODO Ajouter les commentaires sur les simulations \paragraph*{REML vs Maximum Likelihood (ML)} En comparant les méthodes selon diff --git a/data/chen2019pvalues.Rds b/data/chen2019pvalues.Rds index 67099aa..66ab9b7 100644 Binary files a/data/chen2019pvalues.Rds and b/data/chen2019pvalues.Rds differ diff --git a/data/chen2019pvaluesadj.Rds b/data/chen2019pvaluesadj.Rds index b5becf1..07d0fb4 100644 Binary files a/data/chen2019pvaluesadj.Rds and b/data/chen2019pvaluesadj.Rds differ diff --git a/data/simus/real_her_0.3_seed_1234.Rds b/data/simus/real_her_0.3_seed_1234.Rds index 7ebd303..a32eb8f 100644 Binary files a/data/simus/real_her_0.3_seed_1234.Rds and b/data/simus/real_her_0.3_seed_1234.Rds differ diff --git a/data/simus/real_her_0.5_seed_1234.Rds b/data/simus/real_her_0.5_seed_1234.Rds index 0c9fa7b..332dd40 100644 Binary files a/data/simus/real_her_0.5_seed_1234.Rds and b/data/simus/real_her_0.5_seed_1234.Rds differ diff --git a/data/simus/real_her_0.7_seed_1234.Rds b/data/simus/real_her_0.7_seed_1234.Rds index f96758a..8dabd7a 100644 Binary files a/data/simus/real_her_0.7_seed_1234.Rds and b/data/simus/real_her_0.7_seed_1234.Rds differ diff --git a/data/simus/real_her_0.9_seed_1234.Rds b/data/simus/real_her_0.9_seed_1234.Rds index 410e31f..f20b039 100644 Binary files a/data/simus/real_her_0.9_seed_1234.Rds and b/data/simus/real_her_0.9_seed_1234.Rds differ diff --git a/rapport.Rnw b/rapport.Rnw index d00c642..5353de0 100644 --- a/rapport.Rnw +++ b/rapport.Rnw @@ -193,8 +193,6 @@ ggplot(df) + % Besoin de le dire qu'on fait une régression linéaire matrice structurée, % figure avec le Brownien sur l'arbre à reprendre dans le chapitre de livre -TODO Etre assez concis sur l'histoire de la projection et le modèle et les différences avec l'ANOVA. - \subsection{ANOVA phylogénétique avec erreur de mesure} Dans la section précedente, on a supposé que la seule source de variabilité provenait du mouvement brownien sur l'arbre. On rajoute dans cette section une autre variabilité specifiée par $\sigma^2_{err}$ qui à partir de la formule précédente \eqref{eq:ANOVAphylo}, nous donne: diff --git a/rapport.pdf b/rapport.pdf index c1bf3d4..3236527 100644 Binary files a/rapport.pdf and b/rapport.pdf differ