---
title: "Reconstruction d'états ancestraux"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ape)
library(phytools)
library(knitr)
library(kableExtra)


# If you want to turn off messages and warnings 
# in the whole document you can use this code chunk right after the yaml portion.
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
```

```{r, include=FALSE}
treefile <-"/home/quentin/Atelier_phylo/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 
strains_gc_file <-"/home/quentin/Atelier_phylo/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt"
tr <- read.tree(treefile)
ori.tr.label <- tr$tip.label
ori.tr <- tr
```

## Introduction

Si certaines combinaisons de caractères sont systématiquement associées à plusieurs espèces, cela pourrait suggérer que des **forces évolutives**, comme la sélection, ont façonné ces associations. 

Toutefois, les associations non aléatoires de certains traits chez certaines espèces peuvent être dues à l'**héritage commun* de leurs ancêtres et, par conséquent, les changements concomitants dans le temps ne peuvent être déduits.

Si les caractères ont évolué de façon aléatoire sans association, les espèces plus étroitement apparentées sont plus susceptibles d'être semblables que les autres, créant ainsi des **relations apparentes** entre les caractères.

Il est donc nécessaire de considérer les **relations phylogénétiques** entre les espèces lors de l'analyse de leurs caractères.

Deux questions peuvent être abordées lors de l'intégration de la phylogénie dans les données comparatives :

* prise en compte de la non-indépendance inter-espèces dans l'étude des caractères et de leurs relations,
* estimation des paramètres d'évolution des caractères. 

Ces deux questions sont étroitement liées. En effet l'impact de la phylogénie sur la distribution des caractères dépend non seulement de la phylogénie mais aussi de la façon dont ces caractères évoluent.

## Mise en oeuvre

Vérifiez que l'arbre tr est enraciné

```{r, include=FALSE}
is.rooted(tr)
```

Tracé de l'arbre

```{r, include=FALSE}
plot(tr)
```
Lecture des informations sur les souches 

```{r info souches, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
strains_info_file <-"/home/quentin/Atelier_phylo/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt"
strains_info <- read.table(file=strains_info_file, header=T,sep="\t", row.names=1)
strains_info <- strains_info[tr$tip.label,]
```

Changement des labels et annotation de l'arbre

```{r plot annotated tree, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Arbre espèce annoté."}
tr$tip.label <- as.character(strains_info$Strain)
tipcol <- c()
tipcol[1:4] <- 'red'
tipcol[5:16] <- 'blue'

plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol)
nodelabels(tr$node, bg='white', col='purple', frame='n', adj=c(1,-1), cex=0.5)
tiplabels(strains_info$Light, bg='white', col='black', frame='n', adj=c(-0.2, 0.5), cex=1)
add.scale.bar(length=0.1)
```

### GC et taille des génomes

```{r plot tailel_GC, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Taux de GC en fonction de la taille des génomes."}
strains_gc_file <-"/home/quentin/Atelier_phylo/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt"
gc_data <- read.table(file=strains_gc_file, header=F,sep="\t", row.names=1)
colnames(gc_data) <- c('Length', 'GC')
gc_data <- gc_data[ori.tr.label,]
plot(gc_data$Length, gc_data$GC, pch=20, ylim=c(0,1), xlab="taille des génomes", ylab="taux de G+C")

```

### Cartographie de l'évolution d'un caractère continue sur l'arbre

La fonction contMap (Map continuous trait evolution on the tree) trace l'arbre sur lequel est reporté l'évolution d'un caractère continu. La cartographie est réalisée en estimant les états aux nœuds internes à l'aide du ML avec la fonction fastAnc, et l'interpolation des états le long de chaque branche (Felsenstein, 1985).

```{r GC evolution, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Evolution du taux de GC sur les branches de l'arbre."}
gc <- gc_data$GC
names(gc) <- rownames(gc_data)
XX <- contMap(ori.tr, gc, sig=0, outline=F, lwd=5, lims=c(0.3,0.6))
```
```{r taille evolution, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Evolution de la taille des génomes sur les branches de l'arbre."}
ln <- gc_data$Length
names(ln) <- rownames(gc_data)

XX <- contMap(ori.tr, ln, sig=0, outline=F, lwd=5)
```

## Reconstruction par Maximum de vraisemblance

Il existe plusieurs fonctions R utilisant le ML pour estimer les états ancestraux aux nœuds internes pour les caractères continus:

* ace(method="ML"), 
* fastAnc, 
* anc.ML. 

Les trois méthodes effectuent l'estimation du ML en supposant un modèle brownien pour modéliser l'évolution du caractère. La seule différence est la méthode d'optimisation. 

* fastAnc utilise une méthode de ré-enracinement dans laquelle l'arbre est ré-enraciné à chaque nœud interne et l'algorithme de contraste est utilisé pour obtenir l'état ML pour ce noeud. 
* anc.ML utilise une optimisation numérique, mais initie la recherche avec les valeurs obtenues avec fastAnc.

```{r comparaisonML, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Comparaison des valeurs estimées aux noeuds par le trois méthodes qui effectuent l'estimation du ML."}
aceML <- ace(gc, ori.tr, type="continuous", method="ML")
fAnc <- fastAnc(ori.tr, gc, vars=TRUE,CI=TRUE)
ancML <- anc.ML(ori.tr, gc)
obj<-cbind(aceML$ace,fAnc$ace,ancML$ace)
colnames(obj)<-c("ace(method=ML)","fastAnc","anc.ML")
pairs(obj,pch=21,bg="grey",cex=1.5)

```

Tracé de l'arbre 

```{r arbreeaGC, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Etats ancestraux de la fréquence des bases aux noeuds de l'arbre obtenus avec aceML (rouge G+C et bleu A+T)."}
plotTree(tr)
nodelabels(pie=aceML$ace,cex=0.5)
```

### Lien entre caractères

Felsenstein a proposé une méthode qui prend en compte la phylogénie dans l'analyse des données comparatives (GC contre taille du génome). L'idée qui sous-tend sa méthode des **contrastes** est que, si nous supposons qu'un trait continu évolue de façon aléatoire dans n'importe quelle direction (modèle de mouvement brownien), le contraste entre deux espèces devrait avoir une **distribution normale** avec une moyenne nulle et une variance proportionnelle au temps écoulé depuis la divergence. Si les contrastes sont mis à l'échelle avec le temps écoulé depuis la divergence, alors ils ont une variance égale à un.

La fonction pic (phylogenetically independent contrasts)

```{r contrastGC, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Méthode des contrastes avec le taux de GC."}
pic.gc <- pic(gc, ori.tr)
plot(tr)
nodelabels(round(gc, 2), adj = c(0, -0.5), frame = "n")
tiplabels(round(gc, 2), adj = c(-1, -0.5), frame = "n")
tiplabels(round(gc, 2), adj = c(-1, -0.5), frame = "n")

#pic.gc <- pic(gc, ori.tr, rescaled.tree = TRUE)
#plot(pic.gc$rescaled.tree)
#nodelabels(round(gc, 2), adj = c(0, -0.5), frame = "n")
#tiplabels(round(gc, 2), adj = c(-1, -0.5), frame = "n")
#tiplabels(round(gc, 2), adj = c(-1, -0.5), frame = "n")
```

```{r contrastLen, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Méthode des contrastes avec les tailles de génomes."}
lns <- ln/10000000
pic.lns <- pic(lns, ori.tr)
plot(tr)
nodelabels(round(lns,2), adj = c(0, 1), frame = "n")
tiplabels(round(lns,2), adj = c(-1, 1), frame = "n")
```

### Lien entre les pairs de contrastes 

```{r contrastGCLen, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Lien entre les pairs de contrastes."}
plot(pic.lns, pic.gc, xlim=c(-0.15,0.15), ylim=c(-0.15,0.15), xlab="contraste sur le taux de GC", ylab="contraste sur la tailles des génomes")
abline(a = 0, b = 1, lty = 2)

cor(pic.gc, pic.lns)
lm(pic.lns ~ pic.gc)
```

Remarque, faire la régression entre les PICs en passant par l'origine est justifiée seulement si les caractères évoluent sous un modèle de mouvement brownien et s'il existe une relation linéaire entre eux.

### Autocorrélation phylogénétique

Comme les espèces ne sont pas indépendantes par leurs relations phylogénétiques, nous pouvons utiliser cette dépendance pour quantifier l'association entre les variables observées au niveau des espèce. Gittleman et Kot[109] propose une méthode basée sur une approche d'**autocorrélation**. 

Elle utilise l'indice d'autocorrélation I de Moran qui fait intervenir ''wij'' un poids qui quantifie la proximité phylogénétique entre espèces ''i'' et ''j''. 

Avec un arbre phylogénétique, nous avons 

''wij'' = 1/''dij'', 

où ''dij'' sont les distances mesurées sur un arbre, et ''wii'' = 0.

```{r autoGC, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
w <- 1/cophenetic(tr)
diag(w) <- 0
iMGC <- Moran.I(gc, w)
as.data.frame(iMGC)%>% kable(caption = "Table: indice d'autocorrélation de Moran pour le taux de GC.") %>% 
  kable_classic(full_width = F, html_font = "Cambria", font_size = 12)
```

```{r autoLen, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
w <- 1/cophenetic(tr)
diag(w) <- 0
iMlns <- Moran.I(lns, w)
as.data.frame(iMlns)%>% kable(caption = "Table: indice d'autocorrélation de Moran pour la taille des génomes.") %>% 
  kable_classic(full_width = F, html_font = "Cambria", font_size = 12)
```

Le résultat est une liste avec quatre éléments : 
* la valeur observée de I (observée),
* sa valeur attendue sous l'hypothèse nulle d'aucune corrélation,
* l'écart-type de l'I (sd) observé, 
* la P-valeur de l'hypothèse nulle.

La fonction **gearymoran** d'ade4 calcule le coefficient de Moran et test sa signification à l'aide d'une procédure de randomisation. L'option *nrepet* spécifie le nombre de réplications de la randomisation (999 par défaut).

```{r gearymoran, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
library(ade4)
gearymoran(w, data.frame(gc, ln), nrepet = 99999, alter="two-sided")
```
### Habitats

Nous allons maintenant nous intéresser à la profondeur à laquelle les organismes vivent. Cette information est absente pour la souche *Aaap*. Nous allons donc supprimer cette feuille à l'arbre.

```{r habitats, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Profondeur à laquelle les organismes vivent."}
strains_info %>% kable(caption = "Table: Propriétés des souches. Triées selon l'odre des feuilles dans l'arbre espèce.") %>% 
  kable_classic(full_width = F, html_font = "Cambria", font_size = 12)

to_remove <- strains_info[rownames(strains_info) == 'Aaap',1]
tr_d <- drop.tip(tr, as.character(to_remove))
strains_info_d <- strains_info[rownames(strains_info) != 'Aaap', ]
dp <- strains_info_d$Depth
names(dp) <- strains_info_d$Strain
XX <- contMap(tr_d, dp)
```

### Ecotypes

```{r ecotypes, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
ecotypes <- c()
ecotypes[grep('Syn', strains_info$Light)] <- 'Syn'
ecotypes[grep('LL', strains_info$Light)] <- 'LL'
ecotypes[grep('HL', strains_info$Light)] <- 'HL'
names(ecotypes)<-tr$tip.label
ecotypes %>% kable(caption = "Table: Ecotype des souches.") %>% 
  kable_classic(full_width = F, html_font = "Cambria", font_size = 12)
```

Reconstruction des états du caractère aux nœuds de l'arbre avec la fonction ace d'ape (type="discrete", model="ER" for equal-rates model). 

Les vitesses de transition entre les trois états du carctère sont identiques. 

```{r plot ecotypes, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Origine des ecotypes."}
eco_er <- ace(ecotypes, tr, type="discrete", CI=T, model="ER")
plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol)
tiplabels(strains_info$Light, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=eco_er$lik.anc,cex=0.5)
add.scale.bar(length=0.1)
```
```{r plot allecotypes, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Origine de tous les ecotypes."}
allecotypes <- strains_info$Light
eco_er <- ace(allecotypes, tr, type="discrete", CI=T, model="ER")
plot(tr, label.offset=0.2, show.node.label=F, cex=1, tip.color=tipcol)
tiplabels(strains_info$Light, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=eco_er$lik.anc,cex=0.5)
add.scale.bar(length=0.1)
```

## Profils des gènes

```{r lire matchtable, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
matchtable_file <- "/home/quentin/Atelier_phylo/work/ProchlorococcusSynechococcus/panoct/results//matchtable_0_1.txt"
genomes.list_file <- "/home/quentin/Atelier_phylo/work/ProchlorococcusSynechococcus/panoct/results//genomes.list"

matchtable <- read.table(matchtable_file, sep="\t", row.names=1)
genomes.list <- read.table(genomes.list_file)
colnames(matchtable) <- t(genomes.list)
#head(matchtable)
```

Transposer la matrice et ordonner les lignes comme les feuilles de l'arbre

```{r trnasposer, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
t_matchtable <- t(matchtable)
t_matchtable <- t_matchtable[ori.tr$tip.label,]
nb_strains <- nrow(t_matchtable)
```

Compiler les différents motifs de présence absence des gènes observés dans chaque groupe de gènes orthologues dans une liste

```{r compiler, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
motifs <- list()
for ( cluster in 1:ncol(t_matchtable) ) {
   pattern <- paste(t_matchtable[,cluster], sep='', collapse='')
   nb <- length(motifs[[pattern]]$occurences)
   if ( nb == 0 ) {
      motifs[[pattern]]$sring <- t_matchtable[,cluster]
   }
   motifs[[pattern]]$occurences[[nb+1]] <- cluster
   #cat(cluster, pattern, nb, sep="\t", "\n")
}

```
Définir un ensemble de motifs

Pour la mise en main des scripts, il peut être utile de travailler avec un sous ensemble de motifs! 

```{r sous ensemble, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
nb_start <- 1
nb_end <- 10
if ( nb_end > length(motifs)) nb_end <- length(motifs) 

nb_start <- 1
nb_end <- length(motifs) 

```

Taille des patterns trouvés

```{r plot taille pattern, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Occurence des motifs trouvés."}
nb <- nb_start
pat_size <- c()
pat_name <- c()
i <- 1

while (nb <= nb_end ) {
   pattern <- names(motifs)[nb]
   #cat(pattern, "\n")
   size <- length(motifs[[pattern]]$occurences)
   motifs[[pattern]]$size <- size
   if ( size > 10 ) {
      pat_size[i] <- size
      pat_name[i] <- pattern
      #cat(pattern, size, "\n")
      i <- i+1
   }
   nb <- nb+1
}
names(pat_size) <- pat_name
pat_sort <- sort(pat_size, decreasing=T)

par(mar=c(5,10,1,1))
barplot(pat_sort, horiz=T, cex.names = 0.6, las=1, cex=1)
```

### Inférence des états ancestraux

Les différentes fonctions utilisent des paramètres similaires.

Le nombre d'états est déduit des données.

Le modèle est spécifié avec une matrice d'entiers représentant les indices des paramètres : 1 représente le premier paramètre, 2 le second, et ainsi de suite.

Vous devrez choisir la matrice qui spécifie les probabilités de transition entre les états de votre caractère discret. Ace propose des raccourcis pour

* un modèle à un paramètre (ER) où tous les taux de transitions sont égaux,
* un modèle symétrique (SYM) dans lequel les taux de transitions dans les deux sens sont égaux,
* un modèle où tous les taux sont différents (ARD), toutes les transitions possibles entre les états reçoivent des paramètres distincts. 

Vous pouvez également spécifier votre propre matrice. Pour indiquer qu'une transition est impossible, un zéro doit être donné dans la cellule appropriée de la matrice.

La reconstruction la plus parcimonieuse

Cette fonction *MPR* (Most Parsimonious Reconstruction) réalise la reconstruction des caractères ancestraux par **parcimonie**, comme décrit dans Hanazawa et al. (1995) et modifié par Narushima et Hanazawa (1997). Ils ont utilisé l'approche proposée par Farris (1970) et Swofford and Maddison’s (1987) pour reconstruire les états ancestraux. Le caractère est supposé prendre des valeurs entières. 

L'algorithme recherche les ensembles de valeurs pour chaque nœud sous forme d'intervalles avec des valeurs inférieures et supérieures (ex. [1,1]).

L'arbre doit être non enraciné et entièrement dichotomique.

Exemple 
```{r exemple1, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Reconstruction des caractères ancestraux par parcimonie."}
nb <- 3
pattern <- names(motifs)[nb]
mrp <- MPR(motifs[[pattern]]$sring, unroot(ori.tr), "Aaao")

plot(unroot(tr), label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(paste("[", mrp[, 1], ",", mrp[, 2], "]", sep = ""))

```

Tracé MPR pour chaque motif

```{r exemple2, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Reconstruction des caractères ancestraux par parcimonie."}
nb <- nb_start

while (nb < 5 ) {

   pattern <- names(motifs)[nb]
   if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) {
      mrp <- MPR(motifs[[pattern]]$sring, unroot(ori.tr), "Aaao")

      plot(unroot(ori.tr), label.offset=0.2, show.node.label=F, cex=1, main=paste("MRP", nb))
      tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame="n", adj=c(-0.5, 0.5), cex=1)
      nodelabels(paste("[", mrp[, 1], ",", mrp[, 2], "]", sep = ""), bg="white", col="purple")

      #cat ("Press [enter] to continue")
      #line <- readline()

  }
   nb <- nb+1
}
nb_end <- length(motifs) 

```

Marquage stochastique des caractères 

La fonction make.simmap de phytools permet de réaliser un marquage stochastique des caractères aux nœuds de l'arbre.

* nsim: number of simulations.
* pi: gives the prior distribution on the root node of the tree, options are equal, estimated, or a vector with the frequencies. 

Exemple

```{r stochastique, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Marquage stochastique des caractères."}
nb <- 3
pattern <- names(motifs)[nb]
ms_er_trees <- make.simmap(ori.tr, motifs[[pattern]]$sring, model="ER", nsim=100, pi="estimated")

cols<-setNames(c("blue","red"),c(0,1))
plotSimmap(ms_er_trees[[1]], cols)
```

Calculons maintenant les fréquences d'état à partir des marquages stochastiques pour chaque noeud interne. Ce sont nos probabilités postérieures pour les noeuds.

Fonction pour calculer les états: 

```{r map_states, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
map_states <- function(x){
   y <- sapply(x$maps,function(x) names(x)[1])
   names(y) <- x$edge[,1]
   y <- y[as.character(length(x$tip)+1:x$Nnode)]
   return(y)
}

```

Apliquer la fonction à tous les arbres et tracer l'arbre et les probabilités postérieures aux nœuds. 

```{r probapost, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="probabilités postérieures aux nœuds"}
AA <- sapply(ms_er_trees, map_states)
piesA <- t(apply(AA,1,function(x,levels,Nsim) summary(factor(x,levels))/Nsim,levels=c(0,1), Nsim=100))

plot(ori.tr, label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=piesA, cex=0.5, piecol=cols)

```

### Maximum de vraisemblance

La commande *ace* de la librairie ape peut reconstruire des états ancestraux pour des caractères discrets en utilisant le **maximum de vraisemblance** (Pagel 1994 Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters. Proceedings of the Royal Society of London. Series B. Biological Sciences, 255, 37–45.).

Par défaut, la probabilité que les différents états ancestraux des caractères discrets sont calculés à l'aide d'une procédure d'estimation conjointe (**joint**) en utilisant une procédure semblable à celle décrite dans Pupko et al. (2000). Si "marginal = VRAI", une procédure d'estimation marginale est utilisée. Avec cette méthode, les valeurs de vraisemblance à un nœud donné sont calculées en utilisant uniquement les informations provenant des feuilles (et des branches) descendantes de ce nœud. Avec l'estimation conjointe, toutes les informations sont utilisées pour chaque nœud. 

La mise en œuvre actuelle de l'estimation conjointe utilise un algorithme "en deux passes" qui est beaucoup plus rapide que l'approche stochastique alors que les estimations des deux méthodes sont très proches.

Si l'option CI = TRUE est utilisée, alors la probabilité de chaque état ancestral est retournée pour chaque noeud dans une matrice appelée lik.anc.

obj<-anc.Bayes(tree,x,ngen=10000) # sample ancestral states print(obj,printlen=20) ## estimates

Remarque: la fonction fitDiscrete de geiger utilise la même syntaxe que ace et donne des résultats similaires, mais les valeurs de log-vraisemblance sont mises à l'échelle différemment cependant le test du rapport de vraisemblance reste le même.

Exemple motif 3 ("0000001100101111") avec:
* un modèle à un paramètre (ER),
* un modèle symétrique (SYM) ,
* un modèle où tous les taux sont différents (ARD).

```{r exemple3, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Maximum de vraisemblance"}
nb <- 5
pattern <- names(motifs)[nb]
cat(pattern, "\n")
er <- ace(motifs[[pattern]]$sring, ori.tr, type="discrete", model="ER")
sym <- ace(motifs[[pattern]]$sring, ori.tr, type="discrete", model="SYM")
ard <- ace(motifs[[pattern]]$sring, ori.tr, type="discrete", model="ARD")

plot(ori.tr, label.offset=0.2, show.node.label=F, cex=1)
tiplabels(motifs[[pattern]]$sring, bg='white', col='black', frame='n', adj=c(-0.5, 0.5), cex=1)
nodelabels(pie=er$lik.anc, cex=0.3, adj=c(0.5, 0.5))
nodelabels(pie=sym$lik.anc, cex=0.3, adj=c(0.53, 0.5))
nodelabels(pie=ard$lik.anc, cex=0.3, adj=c(0.56, 0.5))

```

La sortie du modèle ER correspond à celle du modèle SYM. Pour un caractère binaire, ces modèles sont identiques, bien qu'ils soient distincts pour les caractères à trois états ou plus. Pour un caractère à trois états, ER est un modèle à un paramètre, SYM un modèle à trois paramètres et ARD un modèle à six paramètres.

Le modèle ARD donne généralement la probabilité la plus élevée. Cependant, il inclut aussi plus de paramètres que le modèle ER. Il est bien connu que l'ajout de paramètres à un modèle augmente généralement sa probabilité. Pour déterminer si l'utilisation du modèle ARD est appropriée, vous devez exécuter un test de vraisemblance.

Test de vraisemblance

La différence de logarithme de probabilité entre deux modèles est distribuée sous forme de chi2, avec des degrés de liberté égaux au nombre de paramètres qui diffèrent entre les modèles. Ainsi, un test de vraisemblance de ces modèles demande si la différence de vraisemblance est suffisamment grande pour se situer dans la queue la plus à droite de la distribution du chi2.

La fonction pchisq(value, df) donne le pourcentage de la fonction de distribution cumulative pour chi2 qui se trouve à gauche de la valeur donnée pour un degré de liberté souhaité(df). Ainsi, 

```{r pchisq, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
 1-pchisq(2*abs(er$loglik - ard$loglik), 1)
```

### Boucle sur tous les motifs

Par curiosité nous allons réaliser le test pour tous les motifs et imprimer les valeurs les plus significatives. 

```{r pchisq significatifs, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
nbnodes <- length(ori.tr$node)
nb <- nb_start
nbpattern <- 1

pat <- c(0)
er <- c(0)
ard <- c(0)
pval <- c(0)

df <- data.frame(pat, er, ard, pval)
names(df) <- c("pattern", "er loglik", "ard loglik", "p-value")

while (nb <= nb_end ) {
   pattern <- names(motifs)[nb]
   if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) {
      #cat(nb, "pattern", pattern, "\n")
      er <- ace(motifs[[pattern]]$sring, ori.tr, type="discrete", model="ER")
      motifs[[pattern]]$er <- er
      #cat("er loglik", er$loglik, "\n")
      
      ard <- ace(motifs[[pattern]]$sring, ori.tr, type="discrete", model="ARD")
      motifs[[pattern]]$ard <- ard
      #cat("ard loglik", ard$loglik, "\n")
      
      #cat("pchisq", "\n")
      ml_test <- 1-pchisq(2*abs(er$loglik - ard$loglik), 1)
      if ( ml_test < 0.01 ) {
         #cat(pattern, er$loglik, ard$loglik, ml_test,"\n")
         df <- rbind(df, c(pattern, round(er$loglik,digit=3), round(ard$loglik,digit=3), round(ml_test,digit=5)))
      }
   }
   nb <- nb+1
}
df %>% kable(caption = "Table: Test de vraisemblance (pvalue < 0.01).") %>% 
  kable_classic(full_width = F, html_font = "Cambria", font_size = 12)
```

Fonction pour annoter les états des feuilles de l'arbre

Les inférences des états du caractère au nœuds sont calculées sur les noeuds interne de l'arbre, nous ajoutons ici les noeuds feuilles. 

```{r tip_states, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
tip_states  <- function(pattern) {
   leaves <- matrix(c(0), nrow=nb_strains, ncol=2)
   for ( i in 1:nb_strains ) {
      leaves[i,1] <- 0
      leaves[i,2] <- 0
      if (pattern[i] == 0 ) {
         leaves[i,1]  <- 1
      } else {
         leaves[i,2] <- 1
      }
      #cat(pattern[i], leaves[i,1], leaves[i,2], "\n")
   }
   return(leaves)
}

```

Exemple d'appel de la fonction

```{r include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
leaves <- tip_states(motifs[[pattern]]$sring)
```

Concaténation des deux matrices (feuilles et noeuds)

```{r include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
edge_states <- rbind(leaves, er$lik.anc)
```

Fonction pour annoter les transitions d'états sur les branches de l'arbre

Une transition d'état est identifiée sur une branche (edge) si l'état du caractère présent au nœud parent est différent de l'état du caractère présent au nœud fils.

* En entrée l'arbre tr et les états du caractère inférés aux noeuds (edge_states).
* En sortie une liste avec les listes de transitions (type de tansition, direction et une couleur).

```{r states_transitions, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
states_transitions <- function(tr, edge_states) {
   
   nb_edges <- nrow(tr$edge)
   direction <- c()
   transition <- c()
   transcolor <- c()

   for ( i in 1:nb_edges ) {
      child <- tr$edge[i,2]
      parent <- tr$edge[i,1]
      if ( edge_states[child, 1] > 0.5 ) {
         child_states <- 0 
      } else {
         child_states <- 1
      } 
      if ( edge_states[parent, 1] > 0.5 ) {
         parent_states <- 0 
      } else {
         parent_states <- 1
      } 
      if ( parent_states != child_states ) {
         transition <- c(transition, i)
         if ( parent_states == 0 ) {
            direction <- c(direction, '+')
            transcolor <- c(transcolor, 'blue')
         } else {
            direction <- c(direction, '-')
            transcolor <- c(transcolor, 'red')
        }
        cat(i, 'parent', parent , parent_states, '->', child_states, 'child', child, direction, "\n")     }
   }
   translist <- list(transition=transition, direction=direction, transcolor=transcolor)
   return(translist)
}
```

 Tracé pour chaque motif
    
```{r gain perte, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Gain perte de gènes (ER)"}
nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2)
mycolors <- c('red', 'blue')

nb <- nb_start
nbpattern <- 1

# while (nb <= 5 ) {
while (nb < nb_end ) {

   pattern <- names(motifs)[nb]
   if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) {
      leaves <- tip_states(motifs[[pattern]]$sring)
      edge_states <- rbind(leaves, motifs[[pattern]]$er$lik.anc)
      st <- states_transitions(tr, edge_states)
      if ( nb <= 5 ) {
        plot(tr, label.offset=0.2, show.node.label=F, cex=1)
        vec <- as.vector(motifs[[pattern]]$sring)
        bgcol <- unlist(lapply(vec, FUN= function(x) mycolors[x+1]))
        tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5)
        nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors)
        edgelabels(text=st$direction, edge=st$transition, frame="r", bg=st$transcolor, adj=c(0,-0.5),cex=0.2)
        add.scale.bar(length=0.1)
      #cat ("Press [enter] to continue")
      #line <- readline()
    }
  }
   nb <- nb+1
}

```

Flux de gènes

```{r gain perte total, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE}
nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2)
mycolors <- c('red', 'blue')

nb <- nb_start
nbpattern <- 1
gain <- vector (mode='numeric',length=nrow(tr$edge))
loss <- vector (mode='numeric',length=nrow(tr$edge))

while (nb < nb_end ) {
#while (nb < 100 ) {

   pattern <- names(motifs)[nb]
   if ( length(grep(1,motifs[[pattern]]$sring)) < nb_strains && length(grep(0,motifs[[pattern]]$sring)) < nb_strains ) {
      leaves <- tip_states(motifs[[pattern]]$sring)
      edge_states <- rbind(leaves, motifs[[pattern]]$er$lik.anc)
      st <- states_transitions(tr, edge_states)
      nb_trans <- length(st$transition)
      if ( nb_trans > 0 ) {
         for ( j in 1:nb_trans) {
           edge <- st$transition[j]
           #cat(nb_trans, j, edge, st$direction[j], "\n")
           if ( st$direction[j] == "+" ) {
            #cat('gain', nb_trans, j, st$direction[j], pattern, "motif size", motifs[[pattern]]$size, "\n")
            gain[edge] <- gain[edge] + motifs[[pattern]]$size
           } else if ( st$direction[j] == "-" ) {
            loss[edge] <- loss[edge] - motifs[[pattern]]$size
           }
         }
      } else {
          cat("no transition, unexpected!\n")
      }
  }
   nb <- nb+1
}

```

Gains et pertes de gènes

```{r plot gain perte, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Flux de gènes avec gains et pertes"}
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
flux <- gain + loss
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
#tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5)
#nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors)
edgelabels(text=gain, frame="n", col='blue', bg='white', adj=c(0,-0.5),cex=0.8)
edgelabels(text=loss, frame="n", col='red', bg='white', adj=c(0,1.5),cex=0.8)
add.scale.bar(length=0.1)

```

Flux de gènes

```{r bilan flux gènes, include=TRUE, eval=TRUE, cache=FALSE, echo=FALSE, fig.align='center', fig.width=8, fig.height=8, fig.cap="Bilan des flux de gènes avec gains et pertes"}
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
flux <- gain + loss
plot(tr, label.offset=0.2, show.node.label=F, cex=1)
#tiplabels(motifs[[pattern]]$sring, bg=bgcol, col='black', frame='c',cex=0.5)
#nodelabels(pie=motifs[[pattern]]$er$lik.anc,cex=0.5, piecol=mycolors)
edgelabels(text=flux, frame="n", col='black', bg='white', adj=c(0,-0.5),cex=0.8)
add.scale.bar(length=0.1)

```