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", "UpSetR", "evemodel", 
    "compcodeR", "mvSLOUCH")

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")
} else {
    require(phylotools)
    require(phytools)
    require(phylolm)
    require(limma)
    require(edgeR)
    require(here)
    require(ggplot2)
    require(dplyr)
    require(tidyr)
    require(UpSetR)
    require(evemodel)
    require(compcodeR)
    require(mvSLOUCH)
    require(patchwork)
}

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*{Vanilla (ML et REML), Satterthwaite (ML et REML), LRT}
<< 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, REML = TRUE)
    })

    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( "VanillaAdj", "VanillaAdjREML",
        "SatterthwaiteAdj", "LRTAdj",
        "SatterthwaiteAdjREML"), 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( "Vanilla", "VanillaREML",
        "Satterthwaite", "LRT",
        "SatterthwaiteREML"), 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és dans le 
code.

Ci-dessous la figure~\ref{fig:pval-methods} présente les p-values des 
différentes méthodes. Il est important de noter que ce graphique présente les
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} 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 %>% 
        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. 


Ces résultats sont présentés dans le diagramme de Venn (figure~\ref{fig:venn-all-methods}).


On peut voir que la méthode de Satterthwaite sans REML a sélectionné énormément
de gènes, 
$\Sexpr{sum(pvalues_adj_dataframe[pvalues_adj_dataframe$test_method == "SatterthwaiteAdj",]$selected)}$ 
comme étant différentiellement exprimés. 
% TODO Comprendre la sur-sélection de Satterthwaite et pas de SatterthwaiteREML

\begin{figure}[H]
    << upset_selection, echo = FALSE>>=
    upset(pvalues_adj_dataframe_wide, 
        nsets = 5,
        mainbar.y.label = "Nombre de gènes en commun",
        sets.x.label = "Nombre de gènes sélectionnés")
    @
    \caption{Diagramme de Venn comparant les gènes sélectionnés selon les méthodes}
    \label{fig:venn-all-methods}
\end{figure}

\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.

\emph{Remarque :} La méthode a produit des \texttt{NA} pour certains gènes,
d'après le message d'erreur, une optimisation n'a pas convergé. Ces gènes sont
présentés dans le tableau~\ref{tab:na-evemodel}.

\subsection*{Toutes les méthodes}

Nous allons ici comparer toutes les méthodes dans un diagramme de Venn (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  %>% 
    pivot_wider(id_cols = gene,
    names_from = test_method, 
    values_from = selected, values_fill = 0) %>%
    data.frame()
@

\begin{figure}[H]
    << full_plot, echo = FALSE>>=
    upset(pvalueseve_dataframe_wide, 
        nsets = 10,
        mainbar.y.label = "Nombre de gènes en commun",
        sets.x.label = "Nombre de gènes sélectionnés")
    @
    \caption{Diagramme de Venn de toutes les méthodes en incluant la méthode EVE}
    \label{fig:venn-all-methods-eve}
\end{figure}