silico.biotoul.fr
 

Atelier Phylogénomique Etats ancestraux

From silico.biotoul.fr

Revision as of 16:05, 15 October 2021 by Quentin (Talk | contribs)
Jump to: navigation, search

Contents

Liens

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.

Répertoire

Nous allons travailler dans un nouveau répertoire:

mkdir ~/work/ProchlorococcusSynechococcus/Ancestral_Characters 
cd ~/work/ProchlorococcusSynechococcus/Ancestral_Characters
cp /home/formation/work/ProchlorococcusSynechococcus/Ancestral_Characters/Strains_info.txt ~/work/ProchlorococcusSynechococcus/Ancestral_Characters

Copier l'arbre espèce de référence

Choisissez l'arbre que vous allez utiliser comme référence et copiez le dans le nouveau répertoire.

Comme alternative, vous pouvez aussi copier celui-ci:

cp /home/formation//work/ProchlorococcusSynechococcus/OG/alignment/alignment_0.4_70_31_nuc_GTR+F+R4_rooted.treefile ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph

Composition en GC des génomes

Comme caractère continue, nous allons étudier l'évolution de la composition en G+C des génomes. Le taux de GC est calculé avec gc_freq.pl.

for i in ~/work/Prochlorococcus/prokka/Aaa*/Aaa*.fna
do      
~/work/scripts/gc_freq.pl --file $i
done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt

for i in ~/work/Synechococcus/prokka/Aaa*/Aaa*.fna
do      
~/work/scripts/gc_freq.pl --file $i
done > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Synechococcus_gc.txt

cat ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Prochlorococcus_gc.txt ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/Synechococcus_gc.txt > ~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt

Arbre espèces

Nous allons utiliser le logiciel R pour réaliser nos analyses et les librairies ape et phytools.

R

library(ape)
treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 
strains_gc_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt"
tr <- read.tree(treefile)

Vérifiez que l'arbre tr est enraciné:

is.rooted(tr)

plot(tr)
plot(tr <- root(tr, interactive = TRUE)) # pour changer la racine de façon interactive

Informations sur les souches

strains_info_file <-"~/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 de l'arbre

tr$tip.label <- as.character(strains_info$Strain)
tipcol <- c()
tipcol[1:4] <- 'red'
tipcol[5:16] <- 'blue'

pdf(file="~/work/ProchlorococcusSynechococcus/images/species_tree.pdf", paper="a4r")
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)
dev.off()
evince ~/work/ProchlorococcusSynechococcus/images/species_tree.pdf

GC et taille des génomes

Lire les fichiers arbre et données:

library(phytools)

treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 
tr <- read.tree(treefile)

strains_gc_file <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/all_gc.txt"
data <- read.table(file=strains_gc_file, header=F,sep="\t", row.names=1)

Ajout des noms des colonnes et réordonnement des lignes de strains_gc en fonction des feuilles de l'arbre tr.

pdf(file="~/work/ProchlorococcusSynechococcus/images/length_GC.pdf", paper="a4r")
colnames(data) <- c('Length', 'GC')
data<- data[tr$tip.label,]
plot(data$Length, data$GC, pch=20, ylim=c(0,1), main="Length/GC")
dev.off()

Number of genes

Question 6.1:
Existe-t-il un lien entre la taille du génome et son contenu en GC? Commentez.
Existe-t-il un lien entre la taille du génome et son contenu en gènes?

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

La fonction 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).

gc <- data$GC
names(gc) <- rownames(data)
ln <- data$Length
names(ln) <- rownames(data)

pdf(file="~/work/ProchlorococcusSynechococcus/images/length_GC_length.pdf", paper="a4r")
XX <- contMap(tr, gc, sig=2, outline=F, lwd=5, lims=c(0.3,0.6))
XX <- contMap(tr, ln, sig=2, outline=F, lwd=5)
dev.off()
Question 6.2:
Commentez les figures obtenues.

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, et 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, par exemple, 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.

Comparaison des trois méthodes sur notre exemple (pour chaque méthode '<method>$ace' est un vecteur contenant les états des nœuds internes)

aceML <- ace(gc,tr,type="continuous",method="ML")
fAnc <- fastAnc(tr,gc,vars=TRUE,CI=TRUE)
ancML <- anc.ML(tr,gc)
obj<-cbind(aceML$ace,fAnc$ace,ancML$ace)
colnames(obj)<-c("ace(method=\"ML\")","fastAnc","anc.ML")
pdf(file="~/work/ProchlorococcusSynechococcus/images/ML_fastAnc_anc.pdf", paper="a4r")
pairs(obj,pch=21,bg="grey",cex=1.5)
dev.off()
Question 6.3:
Que peut on conclure de cette analyse?

Tracé de l'arbre

pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_aceML.pdf", paper="a4r")
plotTree(tr)
nodelabels(pie=aceML$ace,cex=0.5)
dev.off()

Lien entre caractères

Felsenstein a proposé une méthode qui prend en compte la phylogénie dans l'analyse des données comparatives. 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 ces derniers, alors ils ont une variance égale à un.

La fonction pic (phylogenetically independent contrasts)

lns <- ln/10000000
pic.gc <- pic(gc, tr)
pic.lns <- pic(lns, tr)

pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_contrasts.pdf", paper="a4r")
plot(tr)
nodelabels(round(gc, 2), adj = c(0, -0.5), frame = "n")
nodelabels(round(lns,2), adj = c(0, 1), frame = "n")
tiplabels(round(gc, 2), adj = c(-1, -0.5), frame = "n")
tiplabels(round(lns,2), adj = c(-1, 1), frame = "n")
dev.off()

Lien entre les pairs de contrastes

plot(pic.lns, pic.gc)
pdf(file="~/work/ProchlorococcusSynechococcus/images/tree_pairs_contrasts.pdf", paper="a4r")
plot(pic.lns, pic.gc, xlim=c(-0.15,0.15), ylim=c(-0.15,0.15))
abline(a = 0, b = 1, lty = 2) # x = y line

cor(pic.gc, pic.lns)
lm(pic.lns ~ pic.gc)
dev.off()
Question 6.4:
Existe-t-il un lien entre la taille du génome et son contenu en GC (coefficients de la régression)? Commentez.

Remarque, faire la régression entre les PICs en passant par l'origine est justifiée 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

masquée

Habitats

Question 6.6:
Rappelez d'ou proviennent les différents isolats analysés et expliquez comment ont été défini les écotypes.

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.

tr_d <- drop.tip(tr, "Aaap")
strains_info_d <- strains_info[rownames(strains_info) != 'Aaap', ]
dp <- strains_info_d$Depth
names(dp) <- rownames(strains_info_d)

pdf(file="~/work/ProchlorococcusSynechococcus/images/Habitats_profondeur.pdf", paper="a4r")
XX <- contMap(tr_d, dp)
dev.off()
Question 6.7:
Existe-t-il un lien entre profondeur et l'évolution des souches (voir Biller et al., 2015)? Commentez.

Ecotypes

Codage des ecotypes HL et LL

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

Reconstruction des états du caractère aux nœuds de l'arbre:

pdf(file="~/work/ProchlorococcusSynechococcus/images/Ecotypes.pdf", paper="a4r")
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)
dev.off()

Vous pouvez utiliser toutes l'information des ecotype et refaire l'inférence des états du caractère aux nœuds de l'arbre.

ecotypes <- strains_info$Light
Question 6.8:
Sur quelles branches de l'arbre placez vous les transitions d'écotypes? Commentez.

Profils des gènes

R
library(ape)

Lecture de l'arbre espèces

treefile <-"~/work/ProchlorococcusSynechococcus/Ancestral_Characters/species_tree.ph" 
tr <- read.tree(treefile)

Matrice de présence/absence

Nous allons utiliser les groupes de gènes orthologues calculés par panOCT.

matchtable_file <- "~/work/ProchlorococcusSynechococcus/panoct/results//matchtable_0_1.txt"
genomes.list_file <- "~/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

t_matchtable <- t(matchtable)
t_matchtable <- t_matchtable[tr$tip.label,]

nb_strains <- nrow(t_matchtable)

Compiler les différents motifs observés dans une liste

Plusieurs groupes de gènes orthologues peuvent avoir le même profil phylogénétique. Il n'est donc pas nécessaire de faire les calculs sur la totalité.

motifs[[pattern]]$sring : motif associé aux souches
motifs[[pattern]]$occurences : occurence du motif 
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!

nb_start <- 1
nb_end <- 10
if ( nb_end > length(motifs)) nb_end <- length(motifs) 

ou

nb_start <- 1
nb_end <- length(motifs) 

Taille des patterns 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)

pdf(file="~/work/ProchlorococcusSynechococcus/images/taille_des_motifs.pdf", paper="a4r")
par(mar=c(5,10,1,1))
barplot(pat_sort, horiz=T, cex.names = 0.6, las=1, cex=1)
dev.off()
Question 6.9:
Quel sera l'usage de ce calcul?

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 (Most Parsimonious Reconstruction our MPR) 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.

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

Exemple

pdf(file="~/work/ProchlorococcusSynechococcus/images/plus_parcimonieuse.pdf", paper="a4r")
nb <- 3
pattern <- names(motifs)[nb]
mrp <- MPR(motifs[[pattern]]$sring, unroot(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 = ""))
dev.off()

Tracé MPR pour chaque motif

pdf(file="~/work/ProchlorococcusSynechococcus/images/plus_parcimonieuse_all.pdf", paper="a4r")
nb <- nb_start
nb_end < 10
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 ) {
      mrp <- MPR(motifs[[pattern]]$sring, unroot(tr), "Aaao")

      plot(unroot(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
}
dev.off()
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

motif 3 ("0000001100101111")
pdf(file="~/work/ProchlorococcusSynechococcus/images/marquage_stochastique.pdf", paper="a4r")
nb <- 3
pattern <- names(motifs)[nb]
ms_er_trees <- make.simmap(tr, motifs[[pattern]]$sring, model="ER", nsim=100, pi="estimated")

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

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:

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.

pdf(file="~/work/ProchlorococcusSynechococcus/images/probabilites_posterieures.pdf", paper="a4r")
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(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)
dev.off()

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. For discrete characters only maximum likelihood estimation is available

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")

pdf(file="~/work/ProchlorococcusSynechococcus/images/ML3.pdf", paper="a4r")
nb <- 3
pattern <- names(motifs)[nb]
er <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ER")
sym <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="SYM")
ard <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ARD")

plot(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))
dev.off()

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,

 1-pchisq(2*abs(er$loglik - ard$loglik), 1)

donne le pourcentage de la distribution du chi2 qui se trouve à droite de la différence de vraisemblance observée pour les modèles ER et ARD (qui diffèrent d'un paramètre).

Boucle sur tous les motifs

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

nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2, )
nb <- nb_start
nbpattern <- 1
#while (nb <= length(motifs) ) {
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 ) {
      er <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ER")
      motifs[[pattern]]$er <- er
      ard <- ace(motifs[[pattern]]$sring, tr, type="discrete", model="ARD")
      motifs[[pattern]]$er <- ard
      ml_test <- 1-pchisq(2*abs(er$loglik - ard$loglik), 1)
      if ( ml_test < 0.01 ) {
         cat(nb, pattern, er$loglik, ard$loglik, ml_test, "\n")
      }
   }
   nb <- nb+1
}

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.

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

leaves <- tip_states(motifs[[pattern]]$sring)

Concaténation des deux matrices (feuilles et noeuds)

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).
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)
}

Exemple d'appel de la fonction

states_transitions(tr, edge_states)

Tracé pour chaque motif

nbnodes <- length(tr$node)
nodesflux <- matrix(c(0), nrow=nbnodes, ncol=2)
mycolors <- c('red', 'blue')

nb <- nb_start
nbpattern <- 1
#while (nb <= length(motifs) ) {
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)

      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

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 ) {

   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], "\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
pdf(file="~/work/ProchlorococcusSynechococcus/images/gains_pertes_genes.pdf", paper="a4r")
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)
dev.off()
evince ~/work/ProchlorococcusSynechococcus/images/gains_pertes_genes.pdf
Flux de gènes
pdf(file="~/work/ProchlorococcusSynechococcus/images/flux_genes.pdf", paper="a4r")
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)
dev.off()
Question 6.10:
Commentez les figures obtenues.
Comparez ces résultats avec ceux de Sun et Blanchard, 2014)
evince ~/work/ProchlorococcusSynechococcus/images/flux_genes.pdf

Inférence Bayésienne

BayesTraits

BayesTraits is a computer package for performing analyses of trait evolution among groups of species for which a phylogeny or sample of phylogenies is available. This new package incoporates our earlier and separate programes Multistate, Discrete and Continuous. BayesTraits can be applied to the analysis of traits that adopt a finite number of discrete states, or to the analysis of continuously varying traits. Hypotheses can be tested about models of evolution, about ancestral states and about correlations among pairs of traits.

BayesTraitsV3.Manual.pdf

La mise en œuvre de ce logiciel est complexe et ne sera pas traitée dans le cadre de cd TP.