mirror of
https://github.com/Polarolouis/anova-phylogenetique-projet-msv.git
synced 2026-06-17 10:15:25 +02:00
1173 lines
52 KiB
Text
1173 lines
52 KiB
Text
\documentclass[10pt]{beamer}
|
||
|
||
\usetheme[progressbar=frametitle]{metropolis}
|
||
\usecolortheme{aggie}
|
||
|
||
\usepackage{tikz}
|
||
\usepackage[compatibility=false]{caption}
|
||
\usepackage{subcaption}
|
||
\usepackage{graphicx}
|
||
\graphicspath{{img/}}
|
||
\usepackage{appendixnumberbeamer}
|
||
|
||
\setbeameroption{show notes on second screen}
|
||
|
||
\usepackage{booktabs}
|
||
\usepackage[scale=2]{ccicons}
|
||
\usepackage[style=authoryear-comp,backend=biber]{biblatex}
|
||
%== use and define color ==%
|
||
\AtEveryCite{\color{blue}}
|
||
\addbibresource{references.bib}
|
||
|
||
\usepackage{pgfplots}
|
||
\usepackage{amsmath}
|
||
\usepackage{bm}
|
||
\usepackage{listings}
|
||
\lstset{frame=tb,
|
||
language=R,
|
||
keywordstyle=\color{blue},
|
||
alsoletter={.}
|
||
}
|
||
\usepgfplotslibrary{dateplot}
|
||
|
||
% Customize footline to include current and total slide numbers
|
||
\makeatletter
|
||
|
||
\setbeamertemplate{footline}{
|
||
\begin{tikzpicture}[remember picture,overlay]
|
||
% Utiliser un compteur pour sélectionner différentes images avec modulo
|
||
\pgfmathtruncatemacro{\imageindex}{mod(\thepage-1,10)+1}
|
||
\node[anchor=south west, inner sep=0pt] at (current page.south west) {
|
||
% Insérer l'image à gauche
|
||
\includegraphics[height=0.75cm]{img/footer/image\imageindex}
|
||
};
|
||
\node[anchor=south east, inner sep=0pt] at (current page.south east) {
|
||
% Vous pouvez également ajouter d'autres éléments à droite du pied de page si nécessaire
|
||
\ifnum\c@framenumber>1%
|
||
\hfill%
|
||
{\fontsize{10}{12}\selectfont % Adjust the font size here
|
||
\color{gray}%
|
||
\raisebox{0pt}[\height-10pt][\depth]{\insertframenumber/\inserttotalframenumber}
|
||
}%
|
||
\fi%
|
||
};
|
||
\end{tikzpicture}
|
||
|
||
}
|
||
|
||
\makeatother
|
||
|
||
\usepackage{xspace}
|
||
\newcommand{\themename}{\textbf{\textsc{metropolis}}\xspace}
|
||
\newcommand{\dd}{\mathrm{d}}
|
||
\newcommand{\Var}{\mathbb{V}ar}
|
||
\newcommand{\normal}{\mathcal{N}}
|
||
|
||
\title{Projet: ANOVA Phylogénétique}
|
||
\subtitle{Présentation du Mardi 26 Mars 2024}
|
||
\date{}
|
||
\author{Alizée Geffroy, Louis Lacoste, encadrés par Mélina Gallopin et Paul Bastide}
|
||
\institute{M2 MathSV Université Paris-Saclay}
|
||
% \titlegraphic{\hfill\includegraphics[height=1.5cm]{logo.pdf}}
|
||
|
||
\begin{document}
|
||
|
||
\maketitle
|
||
|
||
<<'init', include=FALSE>>=
|
||
knitr::opts_knit$set(fig.align = "center", warnings = FALSE, echo = FALSE, format = "latex")
|
||
require("knitr", quietly = TRUE)
|
||
options(knitr.table.format = "latex")
|
||
@
|
||
|
||
<<'libraries', include=FALSE>>=
|
||
# "phytools", "phylotools"
|
||
necessary_packages <- c("ape", "here")
|
||
|
||
if (any(!(necessary_packages %in% installed.packages()))) {
|
||
install.packages(necessary_packages)
|
||
}
|
||
|
||
require(ape)
|
||
require(here)
|
||
|
||
|
||
source(here("R","utils.R"))
|
||
@
|
||
\begin{frame}{Sommaire}
|
||
\setbeamertemplate{section in toc}[sections numbered]
|
||
\tableofcontents[hideallsubsections]
|
||
\end{frame}
|
||
|
||
% \begin{frame}[allowframebreaks]{Idée structure}
|
||
% TODO Supprimer cette slide temporaire
|
||
% \begin{itemize}
|
||
% \item \textbf{Intro/Contexte} : biologique avec l'exemple de Chen (mettre l'arbre) + figure de l'article ? -> trouver les gènes différentiellement exprimés
|
||
% \item Il existe déjà des méthodes statistiques pour cette problématique (EVEmodel ? State of the Art)
|
||
% \item Transition avec le pourquoi du projet, trouver d'autres méthodes statistiques, adaptées de méthodes classiques qui pourraient bien marcher
|
||
% \item \textbf{Méthode pas par nous} : 1 slide par tiret
|
||
% \begin{itemize}
|
||
% \item Reprendre la forme matricielle de l'ANOVA phylo (mettre en rouge les diffs)
|
||
% \item Présenter le MB qui évolue sur l'arbre + lien matrice K
|
||
% \item Mettre la statistique de test (mettre en rouge la projection (donc diffs))
|
||
% \end{itemize}
|
||
% \item Transition vers notre travail
|
||
% \begin{itemize}
|
||
% \item Mettre la formule avec erreur de mesure avec justification de l'ajout de l'erreur de mesure, formule transfo $V_{\lambda}$, pointer la limite qui est l'erreur dûe à l'estimation du $\lambda$
|
||
% \end{itemize}
|
||
% \item \textbf{Méthode par nous} :
|
||
% \begin{itemize}
|
||
% \item Satterthwaite : préciser que c'est nos calculs à partir de résultats sur modèle mixte (faire slide en appendice) + stat approximée + df formule une méthode possible parmi tant d'autres: Kenward Roger classique
|
||
% \end{itemize}
|
||
% \item \textbf{Simulations} :
|
||
% \begin{itemize}
|
||
% \item les 2 arbres avec les groupes
|
||
% \item Modalités de simulations, bien préciser que l'idée de simuler c'est pour voir erreur de type I et puissance
|
||
% \item Les résultats de simulations: pour les résultats Mettre ANOVA , ANOVA phylo Satterthwaite LRT
|
||
% \end{itemize}
|
||
% \item \textbf{Applications aux données réelles} :
|
||
% \begin{itemize}
|
||
% \item Rappel du type de données, RNA-seq sur pleins de gènes (éventuellement un extrait du tableau ?)
|
||
% \item Mentionner toutes les méthodes rapidement et présenter l'UpSet diagramme avec son analyse et la remarque sur Satterthwaite ML qui sur-sélectionne
|
||
% \end{itemize}
|
||
% \item \textbf{Conclusions/Ouvertures}:
|
||
% \begin{itemize}
|
||
% \item \textbf{Conclusions} :
|
||
% \begin{itemize}
|
||
% \item Récap du projet sur son contenu scientifique
|
||
% \end{itemize}
|
||
% \item \textbf{Ouvertures} :
|
||
% \begin{itemize}
|
||
% \item Utiliser un autre processus stochastique Ornstein-Uhlenbeck
|
||
% \item Comprendre pourquoi Satterthwaite a sur-sélectionné dans l'application: mauvaise implémentation ? évaluer l'impact de l'approx
|
||
% \item Prendre un autre arbre ou ré-échantillonner les groupes dans les simus
|
||
% \item Agrandir le cadre de simulations
|
||
% \item Appliquer les méthodes à d'autres données
|
||
% \item modèle qui fait gène par gène: imaginer en prenant tous les gènes : Limma
|
||
% \end{itemize}
|
||
% \end{itemize}
|
||
% \end{itemize}
|
||
% \end{frame}
|
||
|
||
\section[Introduction]{Introduction}
|
||
|
||
\begin{frame}[fragile]{Contexte biologique}
|
||
\begin{figure}[H]
|
||
\centering
|
||
<<'plot-arbre-chen', out.width = "50%", echo = FALSE>>=
|
||
tree <- read.tree(here("R","chen2019.tree"))
|
||
# Normalising tree edge length
|
||
taille_tree <- diag(vcv(tree))[1]
|
||
tree$edge.length <- tree$edge.length / taille_tree
|
||
|
||
phytools::plotTree(tree, ftype="i")
|
||
@
|
||
\caption{Arbre phylogénétique de \cite{chenQuantitativeFrameworkCharacterizing2019}}
|
||
\label{fig:arbre-chen2019}
|
||
\end{figure}
|
||
|
||
\note{
|
||
\begin{itemize}
|
||
\item Avec l'avènement des données "-omiques" de masses, on a accès à des données quantitatives très nombreuses pour plusieurs espèces.
|
||
\item Le but derrière est de trouver les voies métaboliques et les gènes impliqués dans leur fonctionnement.
|
||
\item Problème : cela coûte cher de regarder tous les gènes.
|
||
\item On veut alors trouver les gènes qui s'expriment différentiellement en utilisant des méthodes statistiques en contrôlant l'erreur de type I, c'est-à-dire les faux positifs.
|
||
\end{itemize}
|
||
|
||
article de Chen:
|
||
\begin{itemize}
|
||
\item En compilant plusieurs jeux de données pour plusieurs espèces avec
|
||
parfois plusieurs individus \emph{cf} l'arbre les auteurs regardent par exemple
|
||
entre les espèces les gènes qui sont différentiellement exprimés dans le foie.
|
||
\item Pour voir s'il y a une différence entre 2 groupes d'espèces (\textit{RatMus vs Autres}). $\mu_1$, $\mu_2$ etc.
|
||
\end{itemize}
|
||
}
|
||
\end{frame}
|
||
\begin{frame}[fragile]{Mouvement brownien}
|
||
\begin{figure}[H]
|
||
\centering
|
||
<<'plot-MB', warnings = FALSE, message = FALSE, out.width = "75%", fig.height = 3.5, echo = FALSE>>=
|
||
source(here("simulations","mouvement_brownien.R"))
|
||
set.seed(12)
|
||
df <- generate_phylo_tree(5,100, 0.1)
|
||
df$id_branche <- as.factor(df$id_branche)
|
||
|
||
ggplot(df) +
|
||
aes(x = time, y = traj, color = as.factor(id_branche))+
|
||
geom_line() +
|
||
# geom_vline(data = df[df$is_spec_time,], aes(xintercept = time, color = .data$id_branche), linetype = "dashed")+
|
||
geom_point(data = df[df$is_spec_time,], aes(x = time, y = .data$traj), color = df[df$is_spec_time,]$id_branche, size = 4)+
|
||
labs(color = "Espèce", x = "Temps", y = "Trait") +
|
||
coord_fixed(ratio=0.75) +
|
||
theme_minimal()
|
||
@
|
||
\caption{Exemple d'un arbre phylogénétique dont le trait est généré selon un Mouvement Brownien}
|
||
\end{figure}
|
||
|
||
\note{
|
||
Pour un arbre phylo donné ça nous renseigne sur les instants de spéciation,
|
||
donc moment de divergence entre 2 espèces représenté ici par les ronds. \\
|
||
IMPORTANT : l'arbre phylogénétique est supposé connu, calibré\footnote{Il s'agit de pouvoir calibrer l'horloge moléculaire.} en temps et
|
||
on n'y touche pas, nous.
|
||
\begin{itemize}
|
||
\item Ici représenté l'évolution d'un trait cad d'une valeur quantitaive qu'on considère : ex comptage du nombre d'ARN exprimé pour un gène donné.
|
||
\item La valeur du trait peut diverger pour chaque espèce à partir du moment de spéciation.
|
||
\item Processus stochastique utilisé comme support de modélisation le mouvement brownien.
|
||
\item Finalement, on a juste accès aux observations aux feuilles, jamais à ce qu'il se passe jusqu'alors
|
||
\end{itemize}
|
||
}
|
||
|
||
\end{frame}
|
||
|
||
\begin{frame}
|
||
Pour un trait Y mesuré chez des espèces i et j, $Cov(Y_i, Y_j)= \sigma_{phylo}^2t_{i,j}$ où $t_{i,j}$ est le temps d'évolution commune.
|
||
\begin{center}
|
||
\includegraphics[width=0.7\textwidth]{matrix_K.png}\footcite{bastideContinuousTraitEvolution2022}
|
||
\end{center}
|
||
\note{
|
||
\begin{itemize}
|
||
\item attention pas l'arbre correspondant au trait simulé dans la slide d'avant
|
||
\item Une fois qu'on a les observations pour chaque, à partir du MB, on a cette covariance là
|
||
\item expliciter avec l'exemple
|
||
\item La matrice V ici porte alors l'information phylogénétique de l'arbre
|
||
\end{itemize}
|
||
}
|
||
\end{frame}
|
||
% \begin{frame}{partie maths}
|
||
% Le modèle de mouvement brownien va alors induire que les feuilles des arbres (nos observations) auront une distribution gausienne que l'on écrira sous la forme suivante
|
||
% % Passer au brownien sur l'arbre donc la dérive mélange tout, le problème est
|
||
% % dur et présenter l'ANOVA phylogénétique
|
||
% Le modèle le plus couramment utilisé pour ce type de données: Expression Variance and Evolution (EVE) modèle décrit dans \cite{rohlfsPhylogeneticANOVAExpression2015}
|
||
|
||
% Il suppose que les traits évoluent selon un
|
||
% processus d'Ornstein-Uhlenbeck et le test réalisé est un \emph{Likelihood Ratio test} (LRT).
|
||
|
||
% Notre projet s'inscrit dans un plus grand projet qui vise à trouver d'autres méthodes statistiques adaptées de méthodes classiques qui feraient aussi bien ou mieux
|
||
% \end{frame}
|
||
|
||
\section{État de l'art}
|
||
|
||
\note{
|
||
\begin{itemize}
|
||
\item On dispose des observations aux feuilles et de l'arbre, ou du moins de l'informtion phylo
|
||
\item Une méthode classique c'est l'ANOVA mais ça ne s'applique pas à nos données car elles ne sont pas indépendantes
|
||
\item Alors dans ce projet on va étudier et utiliser une méthode d'ANOVA phylogénétique que l'on va présenter
|
||
\end{itemize}
|
||
}
|
||
|
||
\begin{frame}{ANOVA phylogénétique}
|
||
\begin{equation}
|
||
Y = X\beta + u \text{, } u \sim \mathcal{N}_n(0, \textcolor{red}{\sigma^2_{phy}K})
|
||
\label{eq:ANOVAphylo}
|
||
\end{equation}
|
||
\[
|
||
\text{où } \mathbf{Y} = \begin{bmatrix} Y_{11} \\\vdots\\ Y_{1n_1} \\
|
||
Y_{21}\\ \vdots \\ Y_{2n_2} \end{bmatrix}\text{, }
|
||
\mathbf{X} = \begin{bmatrix} \mathbf{1} & \mathbf{1_{n_1}} \end{bmatrix}=\begin{bmatrix} 1 & 1 \\ \vdots & \vdots\\1 & 1 \\ 1 & 0\\ \vdots & \vdots\\1 & 0 \end{bmatrix} \text{, } \mathbf{\beta} = \begin{bmatrix} \beta_1 \\ \beta_2 \end{bmatrix} \text{, } n=n_1+n_2
|
||
\]
|
||
\note{
|
||
Rappel de l'ANOVA
|
||
\begin{itemize}
|
||
\item Pour rappel ici formule de 'ANOVA classique matricielle -> écrire au tableau
|
||
\item Au tableau pas oublier de dire $\beta1 = \mu1, \beta2 = \mu1 - \mu2$
|
||
\item ou $\mu_1$ et $\mu_2$ sont les moyennes du trait pour les groupes 1 et 2 que l'on considèrera
|
||
\end{itemize}
|
||
ANOVA phylogénétique
|
||
\begin{itemize}
|
||
|
||
\item L'anova phylo consiste à inclure l'information de l'arbre
|
||
\item Remarque, ici $K$ correspond à la matrice V présentée quand on parlait de l'arbre, la covariance des $Y_i, Y_j$
|
||
\end{itemize}
|
||
}
|
||
|
||
\end{frame}
|
||
|
||
|
||
\begin{frame}{Statistique de test}
|
||
Pour $\ell=\begin{bmatrix}0 \\1 \end{bmatrix}$ :
|
||
\[ H_0 : \beta_2 =0 \Leftrightarrow \ell^T\beta = \begin{bmatrix}0 \\0\end{bmatrix} \text{, les 2 groupes ont la même moyenne}\]
|
||
\[ H_1 : \beta_2\neq 0 \text{, les 2 groupes ont des moyennes différentes}\]
|
||
On a alors la statistique de test suivante : %venant de \cite{bastideContinuousTraitEvolution2022}
|
||
\begin{equation}
|
||
F_{ANOVAphylo}=\frac{||\hat{Y} - \bar{Y}||^2_{\textcolor{red}{K^{-1}}}(n-2)}{||Y - \hat{Y}||^2_{\textcolor{red}{K^{-1}}}} \underset{\mathcal{H}_0}{\sim}\mathcal{F}\text{isher} (1, n-2)
|
||
\end{equation}
|
||
\note{
|
||
\begin{itemize}
|
||
\item H0 : les 2 groupes ont la même moyenne, c'est à dire pour notre exemple que le gène n'est pas diff exprimés le gène a le même niveau d'expression
|
||
\item On peut noter ici par rapport à la stat de test pour l'anova classique, que cela revient à projeter l'écart à la moyenne et les erreurs sur l'inverse de K la matrice des temps de coévolution.
|
||
\item Pour la démo voir les slides en appendices. (Mettre slide 39 et 46 \cite{bastideContinuousTraitEvolution2022})
|
||
\item Donc c'est la projection orthogonale par rapport au produit scalaire $\langle u, v\rangle_{V^{-1}} = u^{T}V^{-1}v$.
|
||
\item Pourquoi les dl de la loi de Fisher : $ 1 = K - 1$ ici $K = 2$ et $n-2 = n - K$.
|
||
\end{itemize}
|
||
}
|
||
|
||
\end{frame}
|
||
|
||
\begin{frame}{ANOVA phylogénétique et erreur de mesure}
|
||
On ajoute une erreur de mesure qui correspond mieux à la réalité des données: erreur intraspécifique
|
||
\begin{equation}
|
||
Y = X\beta + u + \epsilon \text{, } \quad u \sim \mathcal{N}_n(0, \sigma^2_{phy}K) \text{,} \quad \epsilon \sim \mathcal{N}_n(0, \sigma^2_{err}I_n)
|
||
\label{eq:eq2err}
|
||
\end{equation}
|
||
En posant $\lambda = \frac{\sigma^2_{phy}}{\sigma^2_{err}}$ et $E=u+\epsilon$, on peut obtenir une nouvelle forme pour $Y$
|
||
\begin{align}
|
||
\label{eq:V_lambda}
|
||
&Y = X\beta + E \text{, où } Var(E)=V(\theta)=\sigma^2_{phy}(K-\lambda I_n)=\sigma^2_{phy}V_\lambda \\
|
||
&E \sim \mathcal{N}_n(0, \sigma^2_{phy}V_\lambda) \notag
|
||
\end{align}
|
||
Problème: $\lambda$ n'est en général pas connu et il faut l'estimer. Dans ce cas,
|
||
le test n'est pas exact et $F$ ne suit plus la même loi de Fisher.
|
||
|
||
\note{
|
||
\begin{itemize}
|
||
\item Erreur intra-spécifique : variabilité entre les observations
|
||
\item Donc on ne sait pas comment l'estimation de $\lambda$ fait évoluer les degrés de libertés et l'idée est donc ici d'utiliser l'approximation de Satterthwaite pour estimer les degrés de liberté.
|
||
\item Remarque : il est toujours possible de réaliser le test avec cette statistique mais l'on s'attend à ce que le test puisse se tromper.
|
||
\end{itemize}
|
||
}
|
||
|
||
\end{frame}
|
||
|
||
\section{Calculs}
|
||
\note{
|
||
\begin{itemize}
|
||
\item Jusqu'ici nous avons étudier le modèle d'ANOVA phylo,
|
||
ça a été un apprentissage. A partir d'ici ce sont nos calculs avec
|
||
pour but leur implémentation.
|
||
\end{itemize}
|
||
}
|
||
|
||
\begin{frame}{Calcul avec approximation de Satterthwaite}
|
||
Méthode pour approximer les véritables degrés de liberté quand $\lambda$ inconnu\\
|
||
Pour cela on peut voir le modèle comme un modèle mixte\footcite{kuznetsovaLmerTestPackageTests2017} :
|
||
\begin{align}
|
||
&F_{approx}=\frac{||\hat{Y} - \bar{Y}||^2_{V_\lambda^{-1}}df_{approx}}{||Y - \hat{Y}||^2_{V_\lambda^{-1}}} \underset{\mathcal{H}_0}{\sim}\mathcal{F}\text{isher} (1, df_{approx})\\
|
||
\text{Avec } &df_{approx} = \frac{2(f(\hat{\theta}))^2}{[\nabla f(\hat{\theta})]^T A[\nabla f(\hat{\theta})]} \\
|
||
\text{où } f(\theta) &= \ell^TC(\theta)\ell \text{ et A matrice de variance-covariance de } \hat{\theta}=(\hat{\sigma}^2_{phy}, \hat{\sigma}^2_{err}) \notag\\
|
||
C(\theta) &= (Cov(\beta_i , \beta_j))_{i,j}\notag \\
|
||
&= (X^TV(\theta)^{-1}X)^{-1} = (X^T(\sigma^2_{phy}K + \sigma^2_{err}I_n)^{-1}X)^{-1} \notag
|
||
\end{align}
|
||
|
||
\note{
|
||
\begin{itemize}
|
||
\item On obtient alors des degrés de liberté approximés, qui nous permettent d'obtenir une stat de test elle aussi approximée.
|
||
\item A partir de la doc du package lmerTest, en considérant le contexte de modèle mixte on a une formule approximée des degrés de liberté et donc nous avons calculé des formules explicites du gradient de f et de A. Et voilà.
|
||
\item Et nous les avons implémentées pour pouvoir réaliser les tests à partir de cette nouvelle statistique.
|
||
\end{itemize}
|
||
}
|
||
|
||
\end{frame}
|
||
|
||
|
||
\section{Simulations}
|
||
\note{
|
||
Après avoir implémenté, nous avons donc réalisé des simulations afin de regarder comment se comporte les méthodes :
|
||
ANOVA classique, ANOVA phylogénétique et ANOVA phylogénétique avec approximation de Satterthwaite.
|
||
}
|
||
<<'modules-simulations', include = FALSE, eval=TRUE>>=
|
||
necessary_simu <- c("ape", "remotes", "phylolm", "phylolimma", "phytools",
|
||
"latex2exp", "here")
|
||
|
||
if (any(!(necessary_simu %in% installed.packages()))){
|
||
install.packages(necessary_simu)
|
||
remotes::install_github("lamho86/phylolm", quiet = TRUE)
|
||
remotes::install_github("pbastide/phylolimma", quiet = TRUE)
|
||
}
|
||
|
||
library("ape")
|
||
library("phylolm")
|
||
library("phytools")
|
||
library("here")
|
||
library("tidyverse")
|
||
library("ggplot2")
|
||
library("patchwork")
|
||
library("latex2exp")
|
||
source(here("R","utils.R"))
|
||
@
|
||
|
||
<<'Import-arbre', include = FALSE>>=
|
||
K <- 2
|
||
|
||
nb_species <- 43
|
||
|
||
plot_group_on_tree <- function(tree, groups, ...) {
|
||
plot(tree, ...)
|
||
tiplabels(bg = groups, pch = 21, cex = 1.5)
|
||
}
|
||
tree <- read.tree(here("R","chen2019.tree"))
|
||
# Normalising tree edge length
|
||
taille_tree <- diag(vcv(tree))[1]
|
||
tree$edge.length <- tree$edge.length / taille_tree
|
||
@
|
||
|
||
<<'simus-groupes', include = FALSE>>=
|
||
seed <- 1234
|
||
set.seed(seed)
|
||
# Mus et Rat vs le reste
|
||
group_mus_rat_vs_other <- sapply(1:nb_species, function(tip) {
|
||
if (tip %in% getDescendants(tree = tree, 55)) {
|
||
return(1)
|
||
}
|
||
return(2)
|
||
})
|
||
|
||
rng_species <- c("chimp", "bonobo", "human", "orangutan", "marmoset",
|
||
"musMusculus", "rat", "dog", "ferret", "cow", "opossum")
|
||
random_groups <- rowSums(sapply(rng_species,
|
||
function(spec) grepl(spec, tree$tip.label)))
|
||
random_groups[random_groups == 0] <- 2
|
||
@
|
||
|
||
\begin{frame}[fragile]{Modalités de simulations}
|
||
\begin{figure}[H]
|
||
\begin{subfigure}[H]{0.49\textwidth}
|
||
\centering
|
||
<<'plot-groupes-mus', echo = FALSE>>=
|
||
plot_group_on_tree(tree, group = group_mus_rat_vs_other)
|
||
@
|
||
\caption{Groupes \emph{Mus} et rats contre les autres}
|
||
\label{fig:simu-groupes-mus}
|
||
\end{subfigure}
|
||
\hfill
|
||
\begin{subfigure}[H]{0.49\textwidth}
|
||
\centering
|
||
<<'plot-groupes-random', echo = FALSE>>=
|
||
plot_group_on_tree(tree, group = random_groups)
|
||
@
|
||
\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}
|
||
\label{fig:arbres-groupes}
|
||
\end{figure}
|
||
|
||
\note{
|
||
\begin{itemize}
|
||
\item Afin d'avoir une idée des performances des méthodes, nous avons choisis de les
|
||
comparer dans un contexte proche des cas d'application réels.
|
||
\item Nous reprenons l'arbre de \cite{chenQuantitativeFrameworkCharacterizing2019}
|
||
\item Et nous allons tester deux situations
|
||
\item \begin{enumerate}
|
||
\item RatMus contre les autres espèces Figure~\ref{fig:simu-groupes-mus}, donc l'information phylogénétique joue un rôle.
|
||
\item ET les espèces réparties en deux groupes sans lien avec la classification phylogénique
|
||
\end{enumerate}
|
||
\item Ce à quoi on s'attend \textbf{LE DESSIN}: un trait évoluant selon un processus stochastique, aboutissant à 2 valeurs différentes, l'ANOVA dit qu'il y a un saut, mais l'ANOVA phylogénétique connaissant la matrice $K$ dit que c'est normal c'est la dérive, on a divergé il y a longtemps.
|
||
\end{itemize}
|
||
}
|
||
\end{frame}
|
||
\begin{frame}[fragile]{Modalités de simulations}
|
||
<<'param-simulation', echo = FALSE>>=
|
||
# Generate data for rat&mus vs the rest
|
||
N <- 500
|
||
risk_threshold <- 0.05
|
||
|
||
## Standardized parameters
|
||
total_variance <- 1.0 # sigma2_phylo + sigma2_error, fixed [as tree_height = 1]
|
||
heri <- c(0.3, 0.5, 0.7, 0.9) # heritability her = sigma2_phylo / total_variance. 0 means only noise. 1 means only phylo.
|
||
snr <- 1 # signal to noise ratio snr = size_effect / total_variance
|
||
@
|
||
|
||
Nous re-paramétrisons :
|
||
|
||
$$v_{tot} = \sigma^2_{phylo} + \sigma^2_{err} = \Sexpr{total_variance}$$
|
||
|
||
Et alors les paramètres du modèle s'expriment :
|
||
|
||
\begin{align*}
|
||
\sigma^2_{phylo} &= h\times v_{tot},\\
|
||
\sigma^2_{err} &= (1-h) \times v_{tot} \\
|
||
\end{align*}
|
||
|
||
|
||
Ainsi, $h = 0$ signifie qu'il y a seulement du bruit, et $h = 1$
|
||
seulement de l'information phylogénétique.
|
||
|
||
Et nous avons réalisé des simulations pour $h \in \{\Sexpr{heri}\}$.
|
||
|
||
\note{
|
||
\begin{itemize}
|
||
\item A l'implémentation on ré-écrit le modèle pour n'avoir qu'un seul
|
||
paramètre à faire varier, $h$, l'héritabilité.
|
||
Qui se base sur le fait que la variance totale, $v_tot$ \textbf{la formule}.
|
||
\end{itemize}
|
||
}
|
||
|
||
% Choisir les figures que l'on veut mettre en valeur, les méthodes
|
||
|
||
% Montrer les erreurs de type I
|
||
% Pas présenter les ML
|
||
% Enlever ANOVA phylo, Satterthwaite
|
||
|
||
% Comparer ANOVA contre ANOVA phylo REML -> meilleur controle des faux positifs
|
||
|
||
% Parler de Satterthwaite, ANOVA et ANOVA phylo
|
||
|
||
% TODO faire retour sur
|
||
% La prise en compte de la structure de covariance même dans le groupe
|
||
% aléatoire aide car 2 espèce proches sont censées être proche et donc
|
||
% si différentes on peut se dire qu'ils sont bien différent, l'ANOVA phylo
|
||
% a un problème plus facile.
|
||
|
||
|
||
% Couper en deux : ANOVA vs ANOVA phylo
|
||
% ANOVA phylo vs Sattertwhaite
|
||
% Satterthwaite vs LRT
|
||
\end{frame}
|
||
\begin{frame}[fragile]{Résultats: Erreur de type I}
|
||
\begin{figure}[H]
|
||
\begin{subfigure}[H]{0.49\textwidth}
|
||
\centering
|
||
<<'typeI-data-plot-h01', echo = FALSE>>=
|
||
her <- 0.3
|
||
filename <- here("data", "simus",
|
||
paste0("real_her_", her, "_seed_", seed, ".Rds"))
|
||
sim <- readRDS(filename)
|
||
|
||
df_sim_plot <- compute_power_typeI(df = sim)
|
||
|
||
df_sim_plot_filtered <- df_sim_plot %>%
|
||
filter(tested_method %in% c("ANOVA", "ANOVA Phylo REML", "ANOVA Phylo Satterthwaite REML"))
|
||
|
||
ggplot(df_sim_plot_filtered) +
|
||
aes(x = group_type, y = errortypeI, fill = group_type) +
|
||
geom_bar(stat = "identity") +
|
||
scale_y_continuous(limits = c(0, 1)) +
|
||
ylab("Erreur type I") +
|
||
xlab("Type de groupe") +
|
||
# labs(fill = "Type de groupe") +
|
||
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
|
||
facet_wrap(vars(tested_method), ncol = 1) +
|
||
geom_text(aes(label = round(errortypeI, digits = 2)),
|
||
vjust = -0.5, position = position_dodge(width = 0.9)
|
||
) +
|
||
geom_hline(yintercept = 0.05) +
|
||
guides(fill="none") +
|
||
ggtitle("Erreur Type I") +
|
||
theme_minimal()
|
||
@
|
||
\caption{Erreur de type I pour les simulations avec $h=0.3$}
|
||
\label{fig:simu-errI-h03}
|
||
\end{subfigure}
|
||
\hfill
|
||
\begin{subfigure}[H]{0.49\textwidth}
|
||
\centering
|
||
<<'typeI-data-plot-h09', echo = FALSE>>=
|
||
her <- 0.9
|
||
filename <- here("data", "simus",
|
||
paste0("real_her_", her, "_seed_", seed, ".Rds"))
|
||
sim <- readRDS(filename)
|
||
|
||
df_sim_plot <- compute_power_typeI(df = sim)
|
||
|
||
df_sim_plot_filtered <- df_sim_plot %>%
|
||
filter(tested_method %in% c("ANOVA", "ANOVA Phylo REML", "ANOVA Phylo Satterthwaite REML"))
|
||
|
||
ggplot(df_sim_plot_filtered) +
|
||
aes(x = group_type, y = errortypeI, fill = group_type) +
|
||
geom_bar(stat = "identity") +
|
||
scale_y_continuous(limits = c(0, 1)) +
|
||
ylab("Erreur type I") +
|
||
xlab("Type de groupe") +
|
||
# labs(fill = "Type de groupe") +
|
||
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
|
||
facet_wrap(vars(tested_method), ncol = 1) +
|
||
geom_text(aes(label = round(errortypeI, digits = 2)),
|
||
vjust = -0.5, position = position_dodge(width = 0.9)
|
||
) +
|
||
geom_hline(yintercept = 0.05) +
|
||
guides(fill="none") +
|
||
ggtitle("Erreur Type I") +
|
||
theme_minimal()
|
||
@
|
||
\caption{Erreur de type I pour les simulations avec $h=0.9$}
|
||
\label{fig:simu-errI-h09}
|
||
\end{subfigure}
|
||
\caption{Erreurs de type I pour $h \in \{0.3 ; 0.9\}$}
|
||
\label{fig:simu-errI}
|
||
\end{figure}
|
||
|
||
\note{
|
||
\tiny
|
||
\begin{itemize}
|
||
\item Pour comparer les méthodes nous nous intéressons à deux métriques,
|
||
l'erreur de type I et la puissance.
|
||
\item L'erreur de type I est particulièrement importante à contrôler car
|
||
comme nous en avons parlé plus haut, des faux-positifs impliqueraient donc
|
||
des expériences inutiles et particulièrement coûteuses.
|
||
\item A noter ici que nous ne regardons que les méthodes qui minimisent le
|
||
critère REML (Restricted Maximum Likelihood) car ce sont elles qui
|
||
fournissent une estimation de la variance non biaisée. Ce critère améliore
|
||
sensiblement l'estimation de la variance.
|
||
\item Remarque: l'ANOVA telle qu'implémentée dans R utilise directement le
|
||
REML.
|
||
\end{itemize}
|
||
Analyse :
|
||
\begin{itemize}
|
||
\item erreur type I pour groupe phylogénétique: h 0.3 l'ANOVA classiques trompe bcp, l'ANOVA phylo se trompe,seule ANOVA phylo Satter est sous la barre des 0.5 = controle bien les faux positifs
|
||
\item h=0.9 : ANOVA classqiue se trompe toujours et plus, ANOVA phylo est au seuil 5 \% pas étonnant il y a plus d'info phylogénétique, avec Satterthwaite on a aucune erreur
|
||
\item au global l'ANOVA phylo avec Satterthaite controle dans les 2 cas l'erreur de type I, et comme on s'y attend l'ANOVA phylo fait mieux que l'ANOVA classique
|
||
\hfill
|
||
\item groupe pas phylo: h=0.3 l'ANOVA se trompe legerement, elle depasse le seuil, les autres sont en dessous à 0.03
|
||
\item pour h =0.9 l'ANOVA se trompe plus, elle depasse le seuil, les autres sont en dessous
|
||
\item touk, avec faible héritabilité on est dasn un résultat proche de l'attendu : l'ANOVA se trompe à peine, avec forte héritabilité l'erreur est plus marquée ce qui est étonnant au vu des groupes selectionnes
|
||
\item Tout d'abord nos données ne respectent les hypothèses de l'ANOVA.
|
||
On suspecte que la manière dont on a constitué les groupes n'a pas
|
||
suffisamment cassé la phylogénie.
|
||
\end{itemize}
|
||
}
|
||
|
||
\end{frame}
|
||
\begin{frame}{Résultats: puissance}
|
||
\begin{figure}[H]
|
||
\begin{subfigure}[H]{0.49\textwidth}
|
||
\centering
|
||
<<'puissance-data-plot-h01', echo = FALSE>>=
|
||
her <- 0.3
|
||
filename <- here("data", "simus",
|
||
paste0("real_her_", her, "_seed_", seed, ".Rds"))
|
||
sim <- readRDS(filename)
|
||
|
||
df_sim_plot <- compute_power_typeI(df = sim)
|
||
|
||
df_sim_plot_filtered <- df_sim_plot %>%
|
||
filter(tested_method %in% c("ANOVA", "ANOVA Phylo REML", "ANOVA Phylo Satterthwaite REML"))
|
||
|
||
ggplot(df_sim_plot_filtered) +
|
||
aes(x = group_type, y = power, fill = group_type) +
|
||
geom_bar(stat = "identity") +
|
||
scale_y_continuous(limits = c(0, 1)) +
|
||
ylab("Puissance") +
|
||
xlab("Type de groupe") +
|
||
# labs(fill = "Type de groupe") +
|
||
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
|
||
facet_wrap(vars(tested_method), ncol = 1) +
|
||
geom_text(aes(label = round(power, digits = 2)),
|
||
vjust = -0.5, position = position_dodge(width = 0.9)
|
||
) +
|
||
guides(fill="none") +
|
||
ggtitle("Puissance") +
|
||
theme_minimal()
|
||
@
|
||
\caption{Puissances pour les simulations avec $h=0.3$}
|
||
\label{fig:simu-power-h03}
|
||
\end{subfigure}
|
||
\hfill
|
||
\begin{subfigure}[H]{0.49\textwidth}
|
||
\centering
|
||
<<'power-data-plot-h09', echo = FALSE>>=
|
||
her <- 0.9
|
||
filename <- here("data", "simus",
|
||
paste0("real_her_", her, "_seed_", seed, ".Rds"))
|
||
sim <- readRDS(filename)
|
||
|
||
df_sim_plot <- compute_power_typeI(df = sim)
|
||
|
||
df_sim_plot_filtered <- df_sim_plot %>%
|
||
filter(tested_method %in% c("ANOVA", "ANOVA Phylo REML", "ANOVA Phylo Satterthwaite REML"))
|
||
|
||
ggplot(df_sim_plot_filtered) +
|
||
aes(x = group_type, y = power, fill = group_type) +
|
||
geom_bar(stat = "identity") +
|
||
scale_y_continuous(limits = c(0, 1)) +
|
||
ylab("Puissance") +
|
||
xlab("Type de groupe") +
|
||
# labs(fill = "Type de groupe") +
|
||
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
|
||
facet_wrap(vars(tested_method), ncol = 1) +
|
||
geom_text(aes(label = round(power, digits = 2)),
|
||
vjust = -0.5, position = position_dodge(width = 0.9)
|
||
) +
|
||
guides(fill="none") +
|
||
ggtitle("Puissance") +
|
||
theme_minimal()
|
||
@
|
||
\caption{Puissances pour les simulations avec $h=0.9$}
|
||
\label{fig:simu-power-h09}
|
||
\end{subfigure}
|
||
\caption{Puissances pour $h \in \{0.3 ; 0.9\}$}
|
||
\label{fig:simu-power}
|
||
\end{figure}
|
||
|
||
\note{
|
||
\begin{itemize}
|
||
\item puissance statistique: quantifie les vraies positifs il faut qu'elle soit grande
|
||
\item Situation groupe avec info phylogénétique: Revers de la médaille mauvaises puissances voire très mauvaises: CATASTROPHIQUE
|
||
\item h = 0.3 + groupes selectionnes: on a des puissance correctes ce qui nous rassure sur l'implémentation des Méthodes
|
||
\item h= 0.9 + groupes selectionnes : ANOVA et ANOVA phylo pas très bonne puissances, c'est inquiétant : toujours meme hypothèse il reste ude l'info phylo
|
||
\item pour SAtterthwaite ca parait bien c'est la meme chose
|
||
\end{itemize}
|
||
}
|
||
|
||
\end{frame}
|
||
|
||
|
||
\section{Application aux données réelles}
|
||
\note{
|
||
\begin{itemize}
|
||
\item On a des données pour environ 5400 genes pour 17 especes toujours selon le meme arbre
|
||
\item On passe a un test multiple sur tous les genes : on fait un test par gene et on corrige par la technique de Benjamini et Hochberg
|
||
\item Nous allons appliquer les différentes méthodes aux données compilées par
|
||
\cite{chenQuantitativeFrameworkCharacterizing2019}. Il s'agit de données de
|
||
RNA-seq chez 17 espèces et de l'arbre phylogénétique présenté
|
||
figure~\ref{fig:arbre-chen2019}.
|
||
\end{itemize}
|
||
}
|
||
|
||
\begin{frame}[fragile]{UpSet diagram}
|
||
<< 'knitr_options', echo = FALSE>>=
|
||
knitr::opts_knit$set(fig.pos = "HT", fig.width = 6, fig.height = 6,
|
||
fig.align = "center", echo = FALSE, format = "latex")
|
||
|
||
@
|
||
|
||
<< import_modules, echo = FALSE, include=FALSE>>=
|
||
necessary_packages <- c("phylotools", "phytools", "phylolm", "limma", "edgeR",
|
||
"here", "ggplot2", "patchwork", "dplyr", "tidyr", "evemodel",
|
||
"compcodeR", "mvSLOUCH", "ComplexUpset", "see")
|
||
|
||
if (!all(necessary_packages %in% installed.packages())) {
|
||
install.packages(necessary_packages)
|
||
install.packages("remotes")
|
||
remotes::install_gitlab("sandve-lab/evemodel")
|
||
remotes::install_github("pbastide/compcodeR", ref = "phylolimma")
|
||
}
|
||
|
||
require(phylotools)
|
||
require(phytools)
|
||
require(phylolm)
|
||
require(limma)
|
||
require(edgeR)
|
||
require(here)
|
||
require(ggplot2)
|
||
require(dplyr)
|
||
require(tidyr)
|
||
require(ComplexUpset)
|
||
require(see)
|
||
require(evemodel)
|
||
require(compcodeR)
|
||
require(mvSLOUCH)
|
||
require(patchwork)
|
||
|
||
|
||
source(here("R","utils.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))
|
||
@
|
||
<< calcul_pvaleurs, echo = FALSE, warning = FALSE>>=
|
||
### Pvalues computation
|
||
pvalues_data = "chen2019pvalues.Rds"
|
||
pvalues_adj_data = "chen2019pvaluesadj.Rds"
|
||
if (!file.exists(here("data",pvalues_data)) | !file.exists(here("data",pvalues_adj_data))){
|
||
# 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_vanilla.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_vanilla_pvalue(fit_phylo)
|
||
})
|
||
|
||
pvalue_vec_vanilla.REML <- setNames(pvalue_vec_vanilla.REML, rownames(data.trans))
|
||
|
||
pvalue_vec_vanilla_adj.REML <- p.adjust(pvalue_vec_vanilla.REML, 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, REML = TRUE, measurement_error = TRUE)
|
||
compute_satterthwaite_pvalue(fit_phylo, tree = cdata@tree)
|
||
})
|
||
|
||
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")
|
||
|
||
|
||
## Préparation du dataframe
|
||
pvalues_adj_dataframe <- data.frame(
|
||
gene = rep(rownames(data.trans), 5),
|
||
pvalue = c(pvalue_vec_vanilla_adj,
|
||
pvalue_vec_vanilla_adj.REML,
|
||
pvalue_vec_satterthwaite_adj,
|
||
pvalue_vec_lrt_adj,
|
||
pvalue_vec_satterthwaite_adj.REML),
|
||
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(
|
||
gene = rep(rownames(data.trans), 5),
|
||
pvalue = c(pvalue_vec_vanilla,
|
||
pvalue_vec_vanilla.REML,
|
||
pvalue_vec_satterthwaite,
|
||
pvalue_vec_lrt,
|
||
pvalue_vec_satterthwaite.REML),
|
||
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)
|
||
pvalues_dataframe <- pvalues_dataframe %>% mutate(selected = ifelse(pvalue < 0.05, 1, 0))
|
||
|
||
pvalues_adj_dataframe$test_method <- as.factor(pvalues_adj_dataframe$test_method)
|
||
pvalues_adj_dataframe <- pvalues_adj_dataframe %>% mutate(selected = ifelse(pvalue < 0.05, 1, 0))
|
||
|
||
save(pvalues_dataframe, file = here("data", pvalues_data))
|
||
save(pvalues_adj_dataframe, file = here("data", pvalues_adj_data))
|
||
} else {
|
||
load(here("data", pvalues_data))
|
||
load(here("data", pvalues_adj_data))
|
||
}
|
||
@
|
||
<<'evemodel', echo = FALSE, warning = FALSE, include = FALSE>>=
|
||
eve_data = "evechen2019pvalues.Rds"
|
||
if (!file.exists(here("data", eve_data))){
|
||
|
||
# TODO comparer avec le package evemodel, twothetatest
|
||
# Comparer avec OU lrt
|
||
# Arbre sans les replicats et les genes data
|
||
# remotes::install_gitlab("sandve-lab/evemodel")
|
||
# TODO Utiliser les infos de la ligne 83 du Rmd
|
||
|
||
cdataEVE <- readRDS(here("data", "data_TER", "data", "chen2019_rodents_cpd.rds"))
|
||
is.valid <- compcodeR:::check_phyloCompData(cdataEVE)
|
||
if (!(is.valid == TRUE)) stop('Not a valid phyloCompData object.')
|
||
tree_rep <- compcodeR:::getTree(cdataEVE)
|
||
tree_norep <- compcodeR:::getTreeEVE(cdataEVE)
|
||
theta_2_vec <- compcodeR:::getIsTheta2edge(cdataEVE, tree_norep)
|
||
#col_species <- tree_norep$tip.label[compcodeR:::sample.annotations(cdataEVE)$id.species]
|
||
col_species <- tree_norep$tip.label[cumsum(!duplicated(compcodeR:::sample.annotations(cdataEVE)$id.species))]
|
||
|
||
# Normalisation
|
||
nfEVE <- edgeR::calcNormFactors(compcodeR:::count.matrix(cdataEVE) / compcodeR:::length.matrix(cdataEVE), method = 'TMM')
|
||
lib.sizeEVE <- colSums(compcodeR:::count.matrix(cdataEVE) / compcodeR:::length.matrix(cdataEVE)) * nfEVE
|
||
data.normEVE <- sweep((compcodeR:::count.matrix(cdataEVE) + 0.5) / compcodeR:::length.matrix(cdataEVE), 2, lib.sizeEVE + 1, '/')
|
||
data.normEVE <- data.normEVE * 1e6
|
||
|
||
# Transformation
|
||
data.transEVE <- log2(data.normEVE)
|
||
rownames(data.transEVE) <- rownames(compcodeR:::count.matrix(cdataEVE))
|
||
|
||
# Analysis with EVE
|
||
evemodel.results_list <- evemodel::twoThetaTest(tree = tree_norep, gene.data = data.transEVE, isTheta2edge = theta_2_vec, colSpecies = col_species, upperBound = c(theta = Inf, sigma2 = Inf, alpha = log(2)/0.001/1))
|
||
|
||
result.table <- data.frame(pvalue = pchisq(evemodel.results_list$LRT, df = 1, lower.tail = FALSE), logFC = compcodeR:::getlogFCEVE(evemodel.results_list$twoThetaRes, theta_2_vec, tree_norep))
|
||
result.table$score <- 1 - result.table$pvalue
|
||
result.table$adjpvalue <- p.adjust(result.table$pvalue, 'BH')
|
||
|
||
rownames(result.table) <- rownames(compcodeR:::count.matrix(cdataEVE))
|
||
|
||
evemodel_dataframe <- data.frame(gene = rep(rownames(result.table)),
|
||
pvalue = c(result.table$adjpvalue),
|
||
test_method = rep("EVEAdj", each = nrow(result.table)))
|
||
save(evemodel_dataframe, file = here("data", eve_data))
|
||
} else {
|
||
load(file = here("data", eve_data))
|
||
}
|
||
|
||
evegenesNA <- (evemodel_dataframe%>% filter(is.na(pvalue)))$gene
|
||
|
||
evemodel_dataframe <- evemodel_dataframe %>%
|
||
filter(!is.na(pvalue)) %>%
|
||
mutate(selected = ifelse(pvalue < 0.05, 1, 0))
|
||
|
||
evemodel_dataframe$test_method <- as.factor(evemodel_dataframe$test_method)
|
||
@
|
||
<< pvalue_eve_upset, echo = FALSE>>=
|
||
pvalueseve_dataframe <- rbind(pvalues_adj_dataframe, evemodel_dataframe)
|
||
|
||
pvalueseve_dataframe_wide <- pvalueseve_dataframe %>%
|
||
filter(!(test_method %in% c("ANOVA Phylo Satterthwaite Ajustée", "ANOVA Phylo Ajustée"))) %>%
|
||
pivot_wider(id_cols = gene,
|
||
names_from = test_method,
|
||
values_from = selected, values_fill = 0) %>%
|
||
data.frame()
|
||
@
|
||
|
||
\begin{figure}[H]
|
||
\centering
|
||
<< full_plot, out.width = "70%",echo = FALSE>>=
|
||
sets <- colnames(pvalueseve_dataframe_wide)[-1]
|
||
sets_colors <- okabeito_colors()[-2][1:length(sets)]
|
||
highlight_intersections <- lapply(seq_along(sets), function(i) {
|
||
upset_query(set = sets[i], fill = sets_colors[i], color = sets_colors[i])
|
||
})
|
||
|
||
names(sets_colors) <- sets
|
||
|
||
plot <- upset(pvalueseve_dataframe_wide, sets,
|
||
name = "Méthode",
|
||
width_ratio=0.1,
|
||
base_annotations=list(
|
||
# 'Intersection size'=intersection_size(),
|
||
"Taille d'intersection" = intersection_size() +
|
||
scale_fill_venn_mix(
|
||
data=pvalueseve_dataframe_wide,
|
||
guide='none',
|
||
colors=sets_colors
|
||
)
|
||
),
|
||
queries = highlight_intersections,
|
||
set_sizes=(
|
||
upset_set_size() +
|
||
theme(axis.text.x=element_text(angle=90))
|
||
))
|
||
plot
|
||
|
||
# upset(pvalueseve_dataframe_wide,
|
||
# sets = sets,
|
||
# mainbar.y.label = "Nombre de gènes en commun",
|
||
# sets.x.label = "Nombre de gènes sélectionnés",
|
||
# sets.bar.color = sets_colors)
|
||
@
|
||
\caption{\emph{UpSet diagram} de toutes les méthodes en incluant la méthode EVE}
|
||
\label{fig:venn-all-methods-eve}
|
||
\end{figure}
|
||
\note{
|
||
\begin{itemize}
|
||
\item ANOVA phylo Satterthwait REML surselectionne des genes
|
||
\item mais il y a des recoupements
|
||
\item Lrt sous selectionne
|
||
\item la méthode EVE (Expression Variance and Evolution) c'est l'état de
|
||
l'art basé sur lrt
|
||
\end{itemize}
|
||
}
|
||
\end{frame}
|
||
\section{Conclusions et ouvertures}
|
||
\begin{frame}{Conclusions}
|
||
\begin{itemize}
|
||
\item La méthode d'ANOVA phylogénétique avec Satterthwaite parait
|
||
intéressante, notamment pour le contrôle de l'erreur de type I. Mais il
|
||
faudra creuser pour essayer de comprendre la dégradation de la puissance.
|
||
|
||
\item Au début nous utilisions une approximation de l'hessienne,
|
||
remplacée par la forme analytique une fois celle-ci obtenue.
|
||
\end{itemize}
|
||
\note{
|
||
\begin{itemize}
|
||
\item Approximation Hessienne : calculée par approximation numérique, méthode de
|
||
Richardson.
|
||
\end{itemize}
|
||
}
|
||
\end{frame}
|
||
\begin{frame}{Ouvertures}
|
||
|
||
\begin{itemize}
|
||
\item Changer de processus stochastique ? Le processus d'Ornstein-Uhlenbeck.
|
||
\item Comprendre pourquoi avec l'approximation de Satterthwaite sur les
|
||
données réelles il y a eu une sur-sélection.
|
||
\item Changer les conditions de simulations : prendre un autre arbre,
|
||
autres données, ou ré-échantillonner les groupes.
|
||
\item Ces méthodes test gène par gène puis corrige pour faire un test
|
||
multiple. On pourrait développer des méthodes qui font sur tous les gènes
|
||
en même temps (adaptation phylogénétique de la méthode LIMMA).
|
||
\end{itemize}
|
||
\end{frame}
|
||
|
||
\begin{frame}{Remerciements}
|
||
|
||
Merci pour votre attention.
|
||
|
||
|
||
Merci à nos encadrants pour leur accompagnement, leur disponibilité et
|
||
leur gentillesse.
|
||
\end{frame}
|
||
|
||
\section{Références et appendices}
|
||
\begin{frame}[allowframebreaks]{References}
|
||
|
||
\printbibliography
|
||
% \bibliographystyle{abbrv}
|
||
|
||
\end{frame}
|
||
|
||
\appendix
|
||
|
||
\begin{frame}{Code pour les simulations}
|
||
Le code pour les simulations, nos implémentations, le rapport et cette
|
||
présentation est disponible sur notre dépôt GitHub : \\
|
||
|
||
\begin{center}
|
||
\url{https://github.com/Polarolouis/anova-phylogenetique-projet-msv/}
|
||
\end{center}
|
||
\end{frame}
|
||
|
||
\begin{frame}{Concernant les fautes d'orthographes}
|
||
Après relecture du rapport, nous avons pu constater que celui-ci contenait de
|
||
nombreuses coquilles. Nous vous présentons nos excuses.
|
||
\end{frame}
|
||
\begin{frame}{Ornstein-Uhlenbeck}
|
||
Une autre possibilité de processus stochastique est le processus
|
||
d'Ornstein-Uhlenbeck (\emph{mean-reverting process}) décrit par l'EDS suivante:
|
||
|
||
$$\dd r_t = -\theta(r_t - \mu)+\sigma\dd W_t$$
|
||
|
||
\begin{figure}
|
||
\includegraphics[scale=0.20]{OrnsteinUhlenbeck3.png}
|
||
\caption{Exemple de trajectoires de processus d'Ornstein-Uhlenbeck (tiré de Wikipédia)}
|
||
\end{figure}
|
||
\note{
|
||
\begin{itemize}
|
||
\item Ce processus tend avec le temps vers la valeur $\mu$ et peut donc
|
||
modéliser l'existence d'un optimum pour le trait considéré.
|
||
\item C'est le processus qui sous-tend le modèle EVE
|
||
\item Et l'idée derrière que ce n'est pas le processus qui saute, mais
|
||
l'optimum et il modélise donc deux niches différentes.
|
||
\end{itemize}
|
||
}
|
||
\end{frame}
|
||
|
||
\begin{frame}{Lien modèle mixte}
|
||
\begin{align}
|
||
Y &= X\beta + Zu + \epsilon \text{, } \quad u \sim \mathcal{N}_q(0, G) \text{,} \quad \epsilon \sim \mathcal{N}_n(0, \sigma^2I_n) \notag\\
|
||
&\text{où Z de taille n $\times$ q est la matrice de design pour les effets aléatoires } \notag\\
|
||
&\text{et }u \perp \epsilon \notag
|
||
\end{align}
|
||
\end{frame}
|
||
|
||
\begin{frame}{Obtention des estimateurs}
|
||
En connaissant l'arbre et le processus stochastique, le modèle s'écrit:
|
||
$$Y = X\beta + \sigma E, E\sim\normal(0_n, K)$$
|
||
\begin{columns}
|
||
\begin{column}{0.49\textwidth}
|
||
\begin{align*}
|
||
K &= LL^T\\
|
||
\Var [L^{-1}E] &= L^{-1} V [L^{-1}]^T = I
|
||
\end{align*}
|
||
Et on peut alors écrire la régression décorrélée :
|
||
$$L^{-1}Y = (L^{-1}X)\beta + \sigma E', E'\sim\normal(0_n,I)$$
|
||
\end{column}
|
||
\begin{column}{0.49\textwidth}
|
||
De là par les formules classiques on a les estimateurs :
|
||
\begin{align*}
|
||
\hat{\beta} &= (X^T V^{-1} X)^{-1} X^T V^{-1} Y\\
|
||
\hat{\sigma^2} &= \frac{1}{n-p} (Y - X\hat{\beta})^T V^{-1} (Y - X\hat{\beta})\\
|
||
&= \frac{1}{n-p} \|Y - X\hat{\beta}\|_{V^{-1}}^2
|
||
\end{align*}
|
||
\end{column}
|
||
\end{columns}
|
||
\note{
|
||
La F stat
|
||
\begin{align*}
|
||
F &=\frac{\text{variance entre groupes}}{\text{variance intra-groupes}}\\
|
||
\text{variance extra} &= \frac{\sum_{i=1}^p n_i (\bar{Y_{i,.}} - \bar{Y})^2}{p-1}\\
|
||
\text{variance intra} &= \frac{\sum_{i=1}^p \sum_{j=1}^{n_i} (Y_{i,j} - \bar{Y_{i,.}})^2}{n-p}\\
|
||
\end{align*}
|
||
}
|
||
\end{frame}
|
||
|
||
% \begin{frame}[allowframebreaks]{questions posables}
|
||
% % - comment obtenir la stat de test pour anova phylo (Cholesky)
|
||
% % - en quoi c'est un modèle mixte pour Satterthwaite ?
|
||
% % - calcul de la Hessienne optim vs formule analytique, mettre formule analytique
|
||
% % - Le LRT un modèle emboité blabla ?
|
||
% % - sur quoi est basé EVEmodel ?
|
||
% % - Mettre la démo du calcul de la Hessienne
|
||
% % - Ornstein Uhleinbeck : qu'est ce que ça change par rapport au MB ? EVE dit optimum qui saute pas le processus qui saute
|
||
% % Modélise deux niches différentes.
|
||
% % Effet sur la moyenne masi ok , et sur la variance $K_\alpha$, ok pour satterthwaite masi prendre $\alpha$ en compte aussi
|
||
% % Modifie la structure de variance et ajoute un paramètre $\alpha$, $K(\alpha)$, un saut sur l'optima.
|
||
% % - données de comptage transformées donc ok de modéliser par MB
|
||
|
||
% % En écologie ne travaille pas sur autant de traits, spécificité de la RNA-seq
|
||
% % des milliers de données.
|
||
|
||
% % LIMMA pour le cas non phylogénétique.
|
||
% % Pour le cas phylogénétique \texttt{phylolimma}.
|
||
|
||
% % Méthodes d'amélioration essayer de faire quelque chose qui prennent en compte
|
||
% % plusieurs gènes à la fois
|
||
% % - Est ce q'on pourrait faire une méthode comme LIMMA et faire Satterthwaite ?
|
||
% % - c'est bizarre d'utiliser des mesures
|
||
|
||
% % Questions Mélina :
|
||
% % - Qu’est qu’une ANOVA phylogénétique ? En quoi différent l’ ANOVA classique et l’ ANOVA phylogénétique ?
|
||
% % - Comment modéliser l’évolution d’un trait continu sur un arbre (choix du processus dans l’ANOVA phylogénétique : savoir qu’il existe différentes manières de faire, soit on prend un brownien, soit on prend un OU … )
|
||
% % - Comment prendre en compte les erreurs de mesures dans l’anova phylogénétique ? (Car ici, dans le cadre de l’expression des gènes chez plusieurs espèces, on mesure plusieurs individus par espèce, on a donc une variabilité intra-espèce et une variabilité inter-espèces… il faut donc prendre en compte cela dans le modèle, et c’est d’ailleurs ce que fait EVE)
|
||
% % - Quel test effectuer pour tester si on a une différence d’expression significative entre différent groupes d’espèces ? (LRT ou test basé sur la stat de Fisher).
|
||
% % - Qu’est ce qu’un modèle mixte ? Comment estimer les paramètres dans un modèle mixte ? Quels tests stats ? Quel est le lien entre une anova phylo et un modèle mixte ?
|
||
% % - Pourquoi faire du REML au du ML classique ? Dans quel context?
|
||
% % - Pour l'analyse de données réelles, vous avez également été confrontés à un problème de tests multiples : puisque vous faites un test par gènes, et que vous avez des milliers de gènes, alors vous devez "corriger les p-values" pour extraire votre sous liste de "gènes différentiellement exprimés" (deux approches classiques : Bonferroni / BH )
|
||
% % Voici la reformulation de votre texte en utilisant l'environnement "itemize" pour qu'il soit copiable :
|
||
|
||
% \begin{itemize}
|
||
% \item Comment obtenir la stat de test pour ANOVA phylo (Cholesky)
|
||
% \item En quoi c'est un modèle mixte pour Satterthwaite ?
|
||
% \item Calcul de la Hessienne optim vs formule analytique, mettre formule analytique
|
||
% \item Le LRT un modèle emboîté blabla ?
|
||
% \item Sur quoi est basé EVEmodel ?
|
||
% \item Mettre la démo du calcul de la Hessienne
|
||
% \item Ornstein Uhleinbeck : qu'est-ce que ça change par rapport au MB ? EVE dit optimum qui saute pas le processus qui saute Modélise deux niches différentes. Effet sur la moyenne mais ok, et sur la variance $K_\alpha$, ok pour Satterthwaite mais prendre $\alpha$ en compte aussi Modifie la structure de variance et ajoute un paramètre $\alpha$, $K(\alpha)$, un saut sur l'optima.
|
||
% \item Données de comptage mais transformées donc ok de modéliser par MB
|
||
% \item En écologie, ne travaille pas sur autant de traits, spécificité de la RNA-seq des milliers de données.
|
||
% \item LIMMA pour le cas non phylogénétique. Pour le cas phylogénétique \texttt{phylolimma}.
|
||
% \item Méthodes d'amélioration essayer de faire quelque chose qui prennent en compte plusieurs gènes à la fois
|
||
% \begin{itemize}
|
||
% \item Est-ce qu'on pourrait faire une méthode comme LIMMA et faire Satterthwaite ?
|
||
% \item C'est bizarre d'utiliser des mesures
|
||
% \end{itemize}
|
||
% \item Questions Mélina :
|
||
% \begin{itemize}
|
||
% \item Qu’est-ce qu’une ANOVA phylogénétique ? En quoi diffère l’ANOVA classique et l’ANOVA phylogénétique ?
|
||
% \item Comment modéliser l’évolution d’un trait continu sur un arbre (choix du processus dans l’ANOVA phylogénétique : savoir qu’il existe différentes manières de faire, soit on prend un brownien, soit on prend un OU … )
|
||
% \item Comment prendre en compte les erreurs de mesures dans l’ANOVA phylogénétique ? (Car ici, dans le cadre de l’expression des gènes chez plusieurs espèces, on mesure plusieurs individus par espèce, on a donc une variabilité intra-espèce et une variabilité inter-espèces… il faut donc prendre en compte cela dans le modèle, et c’est d’ailleurs ce que fait EVE)
|
||
% \item Quel test effectuer pour tester si on a une différence d’expression significative entre différents groupes d’espèces ? (LRT ou test basé sur la stat de Fisher).
|
||
% \item Qu’est-ce qu’un modèle mixte ? Comment estimer les paramètres dans un modèle mixte ? Quels tests stats ? Quel est le lien entre une ANOVA phylo et un modèle mixte ?
|
||
% \item Pourquoi faire du REML au du ML classique ? Dans quel contexte?
|
||
% \item Pour l'analyse de données réelles, vous avez également été confrontés à un problème de tests multiples : puisque vous faites un test par gènes, et que vous avez des milliers de gènes, alors vous devez "corriger les p-values" pour extraire votre sous-liste de "gènes différentiellement exprimés" (deux approches classiques : Bonferroni / BH )
|
||
% \end{itemize}
|
||
% \end{itemize}
|
||
% \end{frame}
|
||
|
||
\end{document}
|