anova-phylogenetique-projet.../prez.Rnw
Louis Lacoste e669d7e99c Premier gros travail slides
Co-authored-by: AzGeffroy <AzGeffroy@users.noreply.github.com>
2024-03-22 17:29:12 +01:00

422 lines
16 KiB
Text
Raw Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

\documentclass[10pt]{beamer}
\usetheme[progressbar=frametitle]{metropolis}
\usecolortheme{aggie}
\usepackage{tikz}
\usepackage[compatibility=false]{caption}
\usepackage{subcaption}
\usepackage{graphicx}
\graphicspath{{img/}}
\usepackage{appendixnumberbeamer}
\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}
\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
intro/contexte: biologique avec l'exemple de Chen (mettre l'arbre) + figure de l'article ? -> trouver les gènes différentiellement exprimés
Il existe déjà des méthodes statistiques pour cette problématique (EVEmodel ?
State of the Art)
Transition avec le pourquoi du projet, trouver d'autres méthodes statistiques,
adaptées de méthodes classiques qui pourraient bien marcher
Méthode pas par nous : 1 slide par tiret
- Reprendre la forme matricielle de l'ANOVA phylo (mettre en rouge les diffs)
- Présenter le MB qui évolue sur l'arbre + lien matrice K
- Mettre la statistique de test (mettre en rouge la projection (donc diffs))
Transition vers notre travail
- 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$
Méthode par nous :
- 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
Simulations :
- les 2 arbres avec les groupes
- Modalités de simulations, bien préciser que l'idée de simuler c'est pour voir erreur de type I et puissance
- Les résultats de simulations: pour les résultats
Mettre ANOVA , ANOVA phylo Satterthwaite LRT
Applications aux données réelles :
- Rappel du type de données, RNA-seq sur pleins de gènes (éventuellement un extrait du tableau ?)
- 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
Conclusions/Ouvertures:
Conclusions:
- Récap du projet sur son contenu scientifique
Ouvertures :
- Utiliser un autre processus stochastique Ornstein-Uhlenbeck
- Comprendre ppourquoi Sattertwhaite a sur-sélectionné dans l'application: mausvaise implémentation ? évaluer l'impact de l'approx
- Prendre un autre arbre ou ré-échantillonner les groupes dans les simus
- Agrandir le cadre de simulations
- Appliquer les méthodes à d'autre données
- modèle qui fait gène par gène: imaginer en prenant tous les gènes : Limma
\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}
intro/contexte: biologique avec l'exemple de Chen + figure de l'article ?
trouver les gènes différentiellement exprimés
\end{frame}
\begin{frame}{Mouvement brownien}
Présenter le MB qui évolue sur l'arbre + lien matrice K
\end{frame}
% \begin{frame}{partie maths}
% % 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}
\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
\]
% Au tableau pas oublier de dire \beta1 = \mu1, \beta2 = \mu1 - \mu2
\end{frame}
\begin{frame}{Statistique de test}
Pour $l=\begin{bmatrix}0 \\1 \end{bmatrix}$ :
\[ H_0 : \beta_2 =0 \Leftrightarrow l^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}
\end{frame}
\begin{frame}{ANOVA phylogénétique et erreur de mesure}
\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}
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{frame}
\section{Calculs}
\begin{frame}{Calcul avec approximation de Satterthwaite}
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
\end{frame}
\section{Simulations}
<<'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, allowframebreaks]{Modalités de simulations}
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.
\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}
<<'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_{measure} = \Sexpr{total_variance}$$
Et alors les paramètres du modèles sont :
\begin{align*}
h \in (\Sexpr{heri}),&\\
\sigma^2_{phylo} &= h\times v_{tot},\\
\sigma^2_{measure} &= (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.
\end{figure}
% 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}
\section{Application aux données réelles}
\begin{frame}{Rappel des données}
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{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}[allowframebreaks]{Code pour les simulations}
Le code pour les simulations 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 coquille. Nous vous présentons nos excuses.
\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 :
- Quest quune ANOVA phylogénétique ? En quoi différent l ANOVA classique et l ANOVA phylogénétique ?
- Comment modéliser lévolution dun trait continu sur un arbre (choix du processus dans lANOVA phylogénétique : savoir quil 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 lanova phylogénétique ? (Car ici, dans le cadre de lexpression 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 cest dailleurs ce que fait EVE)
- Quel test effectuer pour tester si on a une différence dexpression significative entre différent groupes despèces ? (LRT ou test basé sur la stat de Fisher).
- Quest ce quun 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 )
\end{frame}
\end{document}