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. Nous réanalysons les données, en utilisant les méthodes développées
et testées ci-dessus.

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.

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} La barre indiquant 3681 gènes comptabilise les
gènes qui n'ont été sélectionnés par aucune méthode. Le nombre de gènes 
sélectionnés par méthode est présenté dans le tableau~\ref{tab:data-genes-selectionnes}. 

La méthode EVE détecte ici
\Sexpr{sum(evemodel_dataframe[evemodel_dataframe$test_method == "EVEAdj",]$selected)}
gènes différentiellement exprimés. Étant la méthode état de l'art nous pouvons
nous en servir comme référence. 

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 bel et 
bien 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. Des méthodes 
que nous avons utilisées c'est celle-ci qui semble s'en sortir le mieux, elle
donne des résultats semblables à EVE.


\begin{table}[H]
    \centering
    <<'table-anova-phylo-reml'>>=
    kable(colSums(pvalueseve_dataframe_wide[,-1]), col.names = c("Nombre de gènes sélectionnés"), booktabs = TRUE)
    @
    \caption{Nombre de gènes sélectionnés par méthode}
    \label{tab:data-genes-selectionnes}
\end{table}