Ici nous appliquons les méthodes implémentées sur l'arbre de \cite{chenQuantitativeFrameworkCharacterizing2019}.

% TODO Décrire les données en détails

Les données compilées par \cite{chenQuantitativeFrameworkCharacterizing2019} 
sont des données de RNA-seq, c'est-à-dire des données quantifiant l'expression 
des gènes, par le biais du transcriptome, parmi les différentes espèces du bout 
de l'arbre.

Le but est alors d'identifier les gènes différentiellement exprimés, au sens de
nombre d'ARN par gène différent entre les espèces. 

<< '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))
@

\subsection{Modalités des tests}
<< 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))
}
@

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 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]
    \centering
    << graphique_all_pvalues, echo = FALSE>>=
    all_plots <- plot_spacer()
    for (test in unique(pvalues_dataframe$test_method)) {
        all_plots <- all_plots + ggplot(pvalues_dataframe %>% filter(test_method == test)) +
            aes(x = reorder(gene, pvalue), 
                y = pvalue, color = as.factor(selected)) +
            labs(color = "Selected", x = "Genes", y = "P-values")+
            # geom_bar(stat = "identity", position = "dodge") +
            geom_point()+
            geom_hline(yintercept = 0.05) +
            facet_grid(cols = vars(test_method)) +
            theme(axis.text.x=element_blank(), axis.ticks.x = element_blank())
    }
    all_plots + patchwork::plot_layout(guides = "collect", 
        axis_titles = "collect", tag_level = "new") + 
        plot_annotation(title = "Selected genes by tested methods")
    @
    \caption{\emph{p-values} ordonnées pour les différents tests}
    \label{fig:pval-methods}
\end{figure}

Pour la suite de cette analyse, nous allons appliquer un ajustement des p-values
pour les test multiples, nommément la correction de \cite{benjaminiControllingFalseDiscovery1995}.

<< 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) %>%
        data.frame()
@

Une fois ces corrections appliquées, nous allons comparer les gènes sélectionnés,
c'est-à-dire différentiellement exprimés. 

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 == "ANOVA Phylo Satterthwaite Ajustée",]$selected)}$ 
comme étant différentiellement exprimés.

Ce résultat n'étant pas biologiquement crédible, nous préférons ne pas 
l'afficher dans le \emph{UpSet diagram}, 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}

<<'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)
@

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 pour cette problématique.

Elle détecte ici \Sexpr{sum(evemodel_dataframe[evemodel_dataframe$test_method == "EVEAdj",]$selected)}
gènes différentiellement exprimés.

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, des optimisations n'ont pas convergées. Ces gènes sont
présentés dans le tableau~\ref{tab:na-evemodel}.

\subsection*{Toutes les méthodes}

Nous allons ici comparer toutes les méthodes dans un \emph{UpSet diagram} (figure~\ref{fig:venn-all-methods-eve}) afin de 
voir les gènes sélectionnés en commun et les éventuelles différences entre les 
méthodes.

<< 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]
    << full_plot, 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}

\paragraph{Analyse des résultats} Nous pouvons voir que la méthode la plus 
parcimonieuse est celle utilisant le LRT, qui sélectionne
$\Sexpr{sum(pvalues_adj_dataframe[pvalues_adj_dataframe$test_method == "LRT Ajusté",]$selected)}$
gènes qui sont eux-mêmes \textbf{sélectionnés par toutes les méthodes}.
Cette unanimité sur ces gènes nous invite à penser qu'ils sont en effet 
différentiellement exprimés.

La seconde méthode sélectionnant le moins de gènes est l'ANOVA Phylogénétique 
avec REML. Elle sélectionne $\Sexpr{sum(pvalues_adj_dataframe[pvalues_adj_dataframe$test_method == "ANOVA Phylo REML Ajustée",]$selected)}$
gènes. Ces sélections se décompose en plusieurs sous ensembles

TODO Ici nous avons supposé un mouvement brownien comme processus sous-jacent de
l'arbre mais ce n'est peut-être pas le meilleur modèle et un OU pourrait être 
intéressant. Intéressant pour l'ouverture.

% \begin{table}
%     <<'table-anova-phylo-reml'>>=
%     kable()
%     @
% \end{table}