Estimation du likelihood pour un arbre phylogénétique fixé

Estimation du likelihood pour un arbre phylogénétique fixé Étudiants: Nicolo.DeCoi@unil.ch Larisa.Espinoza@unil.ch Josip.Mikulic@unil.ch Professeur:Salamin Nicolas Assistante:Maryam Zaheri

Buts principaux • Likelihood (L) estimation for a simple tree of nucleotides sequences • Understand the mathematicals processes under this estimation • Programming a general fonction in R to calculate L for a specific tree

Intérêt et pertinence du projet • Utilisation et interprétation des données phylogénétiques • Etude et modélisation de l’évolution des espèces à partir des séquences trouvées • La méthode du ML permet d’optimiser les paramètres pour trouver le meilleur arbre phylogénétique • Détecter la pression sélective pour étudier le rôle d’un gène, les gènes importants sont très conservés

Méthodologie • la distance entre 2 séquences s’estime par le nombre de nucléotides différents à chaque site • Problème: sous estimation des substitutions (substitutions multiples, hidden substitutions) • Markov chain (fig. 1) • Continuous-time Markov process: est un modèle probabilistique qui décrit tous les possibles substitutions qui peuvent se passer à un site on se basant sur la situation actuelle (markov property) et on assumant que chaque site évolue indépendamment Figure 1: représentation des substitutions multiples à un site (Yang, Computational molecular evolution, Oxford universtiy press, 2006)

Markov chain model JC69: • Modèle simple pour nucléotides (T,C,A,G), bon modèle pour des petites distances • Q={qij} est la matrice des taux de substitution: • qij est le taux instantané de substitution du nucléotide i au nucléotide j • la somme de chaque colonne vaut 0 •  a la même valeur pour chaque nucléotide (ex =1) • qij∆t donne la probabilité de substitution dans un petit intervalle de temps • P(t)={pij(t)} est la matrice des probabilités de transition: • pij(t) est la probabilité qu’un nucléotide i devient j après un temps t • cette matrice considère et calcule tous les possible voies par lesquelles est passé le processus évolutif • même dimensions de la matrice Q, la somme de chaque colonne vaut 1 • se calcule ainsi:

Likelihood method: • méthode pour calculer la distance entre 2 séquences • la topologie de l’arbre est fixé • on assume que chaque site évolue indépendamment des autres • the nesting rule : on somme les états du nœud ancestrale après avoir fait la somme pour chaque nœud descendant • utilise une matrice data de dimension s x n Pour un nœud: Pour un site de l’arbre: Pour tous les sites de l’arbre:

Résultats Topologie de l’arbre phylogénétique • The Newick tree format est une écriture utilisé pour représenter des arbres tr <- read.tree() (((S1:0.2,S2:0.2)node1:0.1,S3:0.2)node3:0.1,(S4:0.2,S5:0.2)node2:0.1)node4:0.1; plot.phylo(tr,root.edge=TRUE,underscore=TRUE,no.margin=TRUE,font=2)nodelabels(c("node4","node3","node1","node2"),frame="c",bg="white")

site <- function(seq,Q){ dna<-c(-3,1,1,1,1,-3,1,1,1,1,-3,1,1,1,1,-3) Q<-matrix(dna,nrow=4,ncol=4,dimnames=list("bases"=c("T","C","A","G"), "bases"=c("T","C","A","G"))) seq <- read.dna("/Users/Deketek/Scuola/Mathematique/mc.paml", format="sequential",as.character=TRUE) L <- numeric(length=length(seq[1,])) t1 <- 0.2 t2 <- 0.2 t3 <- 0.2 t4 <- 0.2 t5 <- 0.2 t6 <- 0.1 t7 <- 0.1 t8 <- 0.1 p <- c(0.25,0.25,0.25,0.25) x <- matrix(data=NA,nrow=5,ncol=4) for(hin1:length(seq[1,])){ for(kin1:length(seq[,1])){ if(seq[k,h]=="t"){ x[k,] <- c(1,0,0,0)} else{if(seq[k,h]=="c"){ x[k,] <- c(0,1,0,0)} else{if(seq[k,h]=="a"){ x[k,] <- c(0,0,1,0)} else{if(seq[k,h]=="g"){ x[k,] <- c(0,0,0,1) }}}}} P1<-expm(Q*t1) l <- length(P1[1,]) c <- length(P1[,1]) sum1 <- c(0,0,0,0) C1 <- matrix(data=0,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C1[i,j] <- P1[i,j]*x[1,j] sum1[i] <- C1[i,j]+sum1[i]}} P2<-expm(Q*t2) l <- length(P2[1,]) c <- length(P2[,1]) sum2 <- c(0,0,0,0) C2 <- matrix(data=0,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C2[i,j] <- P2[i,j]*x[2,j] sum2[i] <- C2[i,j]+sum2[i]}} node1 <- sum1*sum2 P3<-expm(Q*t4) l <- length(P1[1,]) c <- length(P1[,1]) sum3 <- c(0,0,0,0) C3 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C3[i,j] <- P3[i,j]*x[4,j] sum3[i] <- C3[i,j]+sum3[i]}} P4<-expm(Q*t5) l <- length(P2[1,]) c <- length(P2[,1]) um4 <- c(0,0,0,0) C4 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C4[i,j] <- P4[i,j]*x[5,j] sum4[i] <- C4[i,j]+sum4[i]}} node2 <- sum3*sum4 P5<-expm(Q*t7) l <- length(P1[1,]) c <- length(P1[,1]) sum5 <- c(0,0,0,0) C5 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C5[i,j] <- P5[i,j]*node1[j] sum5[i] <- C5[i,j]+sum5[i]}} P6<-expm(Q*t2) l <- length(P2[1,]) c <- length(P2[,1]) sum6 <- c(0,0,0,0) C6 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C6[i,j] <- P6[i,j]*x[3,j] sum6[i] <- C6[i,j]+sum6[i]}} node3 <- sum5*sum6 P7<-expm(Q*t6) l <- length(P1[1,]) c <- length(P1[,1]) sum7 <- c(0,0,0,0) C7 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C7[i,j] <- P7[i,j]*node3[j] sum7[i] <- C7[i,j]+sum7[i]}} P8<-expm(Q*t8) l <- length(P2[1,]) c <- length(P2[,1]) sum8 <- c(0,0,0,0) C8 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C8[i,j] <- P8[i,j]*node2[j] sum8[i] <- C8[i,j]+sum8[i]}} node4 <- sum7*sum8 L[h]<- log(sum(node4*p)) } LT <- sum(L) return(LT) }

La matrice Q Les séquences à analyser > Q bases bases T C A G T -3 1 1 1 C 1 -3 1 1 A 1 1 -3 1 G 1 1 1 -3 > seq[,1:5] [,1] [,2] [,3] [,4] [,5] S1 "t" "g" "t" "a" "g" S2 "t" "t" "g" "g" "c" S3 "g" "t" "g" "t" "t" S4 "t" "g" "t" "a" "g" S5 "a" "g" "g" "c" "a" Résultats Les inputs de notre fonction: site <- function(seq,Q){

site <- function(seq,Q){ dna<-c(-3,1,1,1,1,-3,1,1,1,1,-3,1,1,1,1,-3) Q<matrix(dna,nrow=4,ncol=4,dimnames=list("bases"=c("T","C","A","G"), "bases"=c("T","C","A","G"))) seq <- read.dna("/Users/Deketek/Scuola/Mathematique/mc.paml", format="sequential",as.character=TRUE) L <- numeric(length=length(seq[1,])) t1 <- 0.2 t2 <- 0.2 t3 <- 0.2 t4 <- 0.2 t5 <- 0.2 t6 <- 0.1 t7 <- 0.1 t8 <- 0.1 p <- c(0.25,0.25,0.25,0.25) x <- matrix(data=NA,nrow=5,ncol=4) for(hin1:length(seq[1,])){ for(kin1:length(seq[,1])){ if(seq[k,h]=="t"){ x[k,] <- c(1,0,0,0)} else{if(seq[k,h]=="c"){ x[k,] <- c(0,1,0,0)} else{if(seq[k,h]=="a"){ x[k,] <- c(0,0,1,0)} else{if(seq[k,h]=="g"){ x[k,] <- c(0,0,0,1) }}}}} P1<-expm(Q*t1) l <- length(P1[1,]) c <- length(P1[,1]) sum1 <- c(0,0,0,0) C1 <- matrix(data=0,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C1[i,j] <- P1[i,j]*x[1,j] sum1[i] <- C1[i,j]+sum1[i]}} P2<-expm(Q*t2) l <- length(P2[1,]) c <- length(P2[,1]) sum2 <- c(0,0,0,0) C2 <- matrix(data=0,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C2[i,j] <- P2[i,j]*x[2,j] sum2[i] <- C2[i,j]+sum2[i]}} node1 <- sum1*sum2 P3<-expm(Q*t4) l <- length(P1[1,]) c <- length(P1[,1]) sum3 <- c(0,0,0,0) C3 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C3[i,j] <- P3[i,j]*x[4,j] sum3[i] <- C3[i,j]+sum3[i]}} P4<-expm(Q*t5) l <- length(P2[1,]) c <- length(P2[,1]) um4 <- c(0,0,0,0) C4 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C4[i,j] <- P4[i,j]*x[5,j] sum4[i] <- C4[i,j]+sum4[i]}} node2 <- sum3*sum4 P5<-expm(Q*t7) l <- length(P1[1,]) c <- length(P1[,1]) sum5 <- c(0,0,0,0) C5 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C5[i,j] <- P5[i,j]*node1[j] sum5[i] <- C5[i,j]+sum5[i]}} P6<-expm(Q*t2) l <- length(P2[1,]) c <- length(P2[,1]) sum6 <- c(0,0,0,0) C6 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C6[i,j] <- P6[i,j]*x[3,j] sum6[i] <- C6[i,j]+sum6[i]}} node3 <- sum5*sum6 P7<-expm(Q*t6) l <- length(P1[1,]) c <- length(P1[,1]) sum7 <- c(0,0,0,0) C7 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C7[i,j] <- P7[i,j]*node3[j] sum7[i] <- C7[i,j]+sum7[i]}} P8<-expm(Q*t8) l <- length(P2[1,]) c <- length(P2[,1]) sum8 <- c(0,0,0,0) C8 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C8[i,j] <- P8[i,j]*node2[j] sum8[i] <- C8[i,j]+sum8[i]}} node4 <- sum7*sum8 L[h]<- log(sum(node4*p)) } LT <- sum(L) return(LT) }

Résultats • x <- matrix(data=NA,nrow=5,ncol=4) • for(h in 1:length(seq[1,])){ • for(k in 1:length(seq[,1])){ • if(seq[k,h]=="t"){ • x[k,] <- c(1,0,0,0) • } • else{ • if(seq[k,h]=="c"){ • x[k,] <- c(0,1,0,0) • } • else{ • if(seq[k,h]=="a"){ • x[k,] <- c(0,0,1,0) • } • else{ • if(seq[k,h]=="g"){ • x[k,] <- c(0,0,0,1) • }}}}} > xsite1 bases seq T C A G S1 1 0 0 0 S2 0 0 1 0 S3 0 0 1 0 S4 0 0 0 1 S5 1 0 0 0 

Node 1 sum2 sum1 X t2 t1 P1<-expm(Q*t1) sum1 <- c(0,0,0,0) l <- length(P1[1,]) c <- length(P1[,1]) C1 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C1[i,j] <- P1[i,j]*x[1,j] sum1[i] <- C1[i,j]+sum1[i] }} P2<-expm(Q*t2) sum2 <- c(0,0,0,0) l <- length(P2[1,]) c <- length(P2[,1]) C2 <- matrix(data=NA,nrow=l,ncol=c) for(iin1:l){ for(jin1:c){ C2[i,j] <- P2[i,j]*x[2,j] sum2[i] <- C2[i,j]+sum2[i] }} node1 <- sum1*sum2 X[1,j] X[2,j] Exemple avec t = 0.2 > P 4 x 4 Matrix of class "dsyMatrix" bases bases T C A G T 0.5869967 0.1376678 0.1376678 0.1376678 C 0.1376678 0.5869967 0.1376678 0.1376678 A 0.1376678 0.1376678 0.5869967 0.1376678 G 0.1376678 0.1376678 0.1376678 0.5869967 Résultats Exemple pour le premier noeud

Résultats Vecteur vide à remplir avec le likelihood pour chaque site L <- numeric(length=length(seq[1,])) Loop pour calculer le likelihood de chaque site for(hin1:length(seq[1,])){ Définir la valeur qui entre dans le vecteur L et l’output de la fonction L[h]<- log(sum(node4*p)) } LT <- sum(L) return(LT) }

Perspectives • Optimisation des paramètres de l’arbre • Passage du modèle JC69 à des modèles pour nucléotides plus complexes • Passage au Codon-based model: • La chaîne de Markov décrit les substitutions d’un codon i en j • Ne considère pas les délétions/insertions • Les mutations se passent indépendamment au 3 sites des codons • Les paramètres de Q: {qij} = • 0 si i et j diffèrent à 2 ou 3 positions du codon • jsi i et j diffèrent d’une transversion synonyme • jsi i et j diffèrent d’une transition synonyme • jsi i et j diffèrent d’une transversion non-synonyme • jsi i et j diffèrent d’une transition non-synonyme • Estimation du taux de substitutions nonsynonyme/ synonyme (ω= dn/ds) est un important indicateur de la pression sélective: • ω = 1 -> neutral mutations • ω < 1 -> purifyng selection • ω > 1 -> diversifyng positive selection

Bibliographie Livres: Yang Z., Computational Molecular Evolution, Oxford university press, 2006 Paradis E., Analysis of Phylogenetics and Evolution with R, Springer, 2006 Masatoshi N./Sudhir K., Molecular Evolution and Phylogenetics, Oxford university press, 2000 Articles: • Goldman N./Yang Z., "A Codon-based Model of Nucleotide Substitution for Protein-coding DNA Sequences", The University of Chicago, 1994 • Yang & all, "Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites", Genetics nb 155, p. 431-449, 2000

Merci pour l’attention!! Questions?

Estimation du likelihood pour un arbre phylogénétique fixé