mirror of
https://github.com/Polarolouis/anova-phylogenetique-projet-msv.git
synced 2026-06-17 18:25:25 +02:00
✨ Add some frames and figures
This commit is contained in:
AzGeffroy
parent
59a8a45c96
commit
a5136dee9c
2 changed files with 333 additions and 40 deletions
373
prez.Rnw
373
prez.Rnw
|
|
@ -96,47 +96,53 @@ source(here("R","utils.R"))
|
|||
|
||||
\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
|
||||
\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}
|
||||
|
||||
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
|
||||
|
||||
|
||||
|
||||
|
|
@ -362,6 +368,7 @@ Et alors les paramètres du modèles sont :
|
|||
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
|
||||
|
||||
|
|
@ -384,6 +391,14 @@ Et alors les paramètres du modèles sont :
|
|||
% ANOVA phylo vs Sattertwhaite
|
||||
% Satterthwaite vs LRT
|
||||
\end{frame}
|
||||
\begin{frame}{Résultats: Erreur de type I}
|
||||
Add erreur de type 1 pour comparer que ça avec Satterthwaite REML, ANOVA et ANOVA phylo REML
|
||||
|
||||
\end{frame}
|
||||
\begin{frame}{Résultats: puissance}
|
||||
Comparer les puissances des méthodes principales
|
||||
|
||||
\end{frame}
|
||||
|
||||
|
||||
\section{Application aux données réelles}
|
||||
|
|
@ -393,6 +408,284 @@ Et alors les paramètres du modèles sont :
|
|||
RNA-seq chez 17 espèces et de l'arbre phylogénétique présenté
|
||||
figure~\ref{fig:arbre-chen2019}.
|
||||
\end{frame}
|
||||
\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 != "ANOVA Phylo Satterthwaite 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
|
||||
|
||||
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))
|
||||
))
|
||||
|
||||
# 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}
|
||||
\end{frame}
|
||||
\section{Conclusions et ouvertures}
|
||||
\begin{frame}{Conclusions scientifiques}
|
||||
Récap du projet sur son contenu scientifique
|
||||
\end{frame}
|
||||
\begin{frame}{Problèmes rencontrés}
|
||||
\begin{itemize}
|
||||
\item Utilisation de l'approx de la Hessienne -> expression analytique
|
||||
\end{itemize}
|
||||
\end{frame}
|
||||
\begin{frame}{Ouvertures}
|
||||
|
||||
\begin{itemize}
|
||||
\item Utilisation du processus d'Ornstein-Uhlenbeck
|
||||
\item Pourquoi Satterthwaite a surselectionné ? Creuser
|
||||
\item Prendre un autre arbre, autres données, ou ré-échantillonner les groupes dans les simus
|
||||
\item Modèle qui fait gène par gène: imaginer en prenant tous les gènes -> méthode LIMMA
|
||||
\end{itemize}
|
||||
\end{frame}
|
||||
|
||||
\begin{frame}{Remerciements}
|
||||
|
||||
|
|
|
|||
BIN
prez.pdf
BIN
prez.pdf
Binary file not shown.
Loading…
Add table
Reference in a new issue