Atelier Phylogénomique Etats ancestraux
From silico.biotoul.fr
m (→Liens) |
m (→Liens) |
||
Line 1: | Line 1: | ||
==Liens== | ==Liens== | ||
*[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Reconstruction_d.27.C3.A9tats_ancestraux Etats ancestraux] | *[http://silico.biotoul.fr/p/Atelier_Phylog%C3%A9nomique#Reconstruction_d.27.C3.A9tats_ancestraux Etats ancestraux] | ||
+ | |||
+ | * R files: /home/formation/work/ProchlorococcusSynechococcus/Reconstruction_etats_ancestraux | ||
+ | |||
==Introduction== | ==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 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. |
Current revision as of 13:24, 8 December 2022
Liens
- R files: /home/formation/work/ProchlorococcusSynechococcus/Reconstruction_etats_ancestraux
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.
La mise en œuvre de ce logiciel est complexe et ne sera pas traitée dans le cadre de cd TP.