Utiliser latex pour expliquer la programmation dynamique

24 avril 2010

Vous voulez faire découvrir la programmation dynamique à des lycéens ? Rien de mieux qu’une petite présentation… En bon geek, vous allez la faire en latex. Mais pour ça, il va falloir connaître quelques ficelles, surtout quand on veut pouvoir représenter la façon dont on remplit la matrice, pas à pas:

Pour que ce soit lisible, j’ai mis le code ici, sur gist.github. Une fois que vous l’avez récupérer, lancez la commande « pdflatex », vous obtiendrez un document comme ça: dynamic_programming.
\documentclass{beamer}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{verbatim}
\usepackage{url}</code>

\setbeamertemplate{navigation symbols}{}

\setbeamertemplate{footline}[page number]

\title{About optimal sequence alignment}
\subtitle{A short glimpse into bioinformatics}

\begin{document}

\begin{frame}
\titlepage
\end{frame}

\begin{frame}
\frametitle{Pairwise sequence alignment}
Assumptions:
\begin{itemize}
\item sequences $S_1$ and $S_2$ are homologous, they share a common ancestor;
\item differences between them are due to only two kinds of events, substitutions and insertion-deletions.
\end{itemize}

Strategy:
\begin{itemize}
\item choose a scoring matrix (reward for match, penalty for mismatch and gap);
\item compute the editing distance (number of matches, mismatches and gaps) to go from one sequence to the other;
\item keep the alignment with the highest score.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Needleman-Wunsch algorithm}
\begin{itemize}
\item Aim: find the optimal global alignment of sequences $S_1$ and $S_2$
\item Recursion rule:
\begin{align}
D(i,j) = max
\begin{cases}
D(i-1,j-1) + score( S_1[i], S_2[j] )\\
D(i-1,j) + gap\\
D(i,j-1) + gap
\end{cases}
\end{align}
\item Scoring scheme: identity=0 transition=-2 transversion=-5 gap=-10
\item Sequences: $S_1$=TTGT $S_2$=CTAGG
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Fill the matrix}
\begin{tabular}{l|cccccccccccc}
&amp; &amp; &amp; C &amp; &amp; T &amp; &amp; A &amp; &amp; G &amp; &amp; G \\ \hline
&amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; \\
&amp; \visible{0} &amp; \visible{$\rightarrow$} &amp; \visible{-10} &amp; \visible{$\rightarrow$} &amp; \visible{-20} &amp; \visible{$\rightarrow$} &amp; \visible{-30} &amp; \visible{$\rightarrow$} &amp; \visible{-40} &amp; \visible{$\rightarrow$} &amp; \visible{-50} \\
&amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \\
T &amp; \visible{-10} &amp; \visible{$\rightarrow$} &amp; \visible{-2} &amp; \visible{$\rightarrow$} &amp; \visible{-10} &amp; \visible{$\rightarrow$} &amp; \visible{-20} &amp; \visible{$\rightarrow$} &amp; \visible{-30} &amp; \visible{$\rightarrow$} &amp; \visible{-40} \\
&amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \\
T &amp; \visible{-20} &amp; \visible{$\rightarrow$} &amp; \visible{-12} &amp; \visible{$\rightarrow$} &amp; \visible{-2} &amp; \visible{$\rightarrow$} &amp; \visible{-12} &amp; \visible{$\rightarrow$} &amp; \visible{-22} &amp; \visible{$\rightarrow$} &amp; \visible{-32} &amp; \\
&amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \\
G &amp; \visible{-30} &amp; \visible{$\rightarrow$} &amp; \visible{-22} &amp; \visible{$\rightarrow$} &amp; \visible{-12} &amp; \visible{$\rightarrow$} &amp; \visible{-4} &amp; \visible{$\rightarrow$} &amp; \visible{-12} &amp; \visible{$\rightarrow$} &amp; \visible{-22} &amp; \\
&amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \visible{$\searrow$} &amp; \visible{$\downarrow$} &amp; \\
T &amp; \visible{-40} &amp; \visible{$\rightarrow$} &amp; \visible{-32} &amp; \visible{$\rightarrow$} &amp; \visible{-22} &amp; \visible{$\rightarrow$} &amp; \visible{-14} &amp; \visible{$\rightarrow$} &amp; \visible{-9} &amp; \visible{$\rightarrow$} &amp; \visible{-17} &amp;
\end{tabular}
\end{frame}

\begin{frame}
\frametitle{Traceback}
\begin{tabular}{l|cccccccccccc}
&amp; &amp; &amp; C &amp; &amp; T &amp; &amp; A &amp; &amp; G &amp; &amp; G \\ \hline
&amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; &amp; \\
&amp; \color{red}{0} &amp; $\rightarrow$ &amp; -10 &amp; $\rightarrow$ &amp; -20 &amp; $\rightarrow$ &amp; -30 &amp; $\rightarrow$ &amp; -40 &amp; $\rightarrow$ &amp; -50 \\
&amp; $\downarrow$ &amp; \color{red}{$\searrow$} &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; \\
T &amp; -10 &amp; $\rightarrow$ &amp; \color{red}{-2} &amp; $\rightarrow$ &amp; -10 &amp; $\rightarrow$ &amp; -20 &amp; $\rightarrow$ &amp; -30 &amp; $\rightarrow$ &amp; -40 \\
&amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; \color{red}{$\searrow$} &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; \\
T &amp; -20 &amp; $\rightarrow$ &amp; -12 &amp; $\rightarrow$ &amp; \color{red}{-2} &amp; \color{red}{$\rightarrow$} &amp; \color{red}{-12} &amp; $\rightarrow$ &amp; -22 &amp; $\rightarrow$ &amp; -32 &amp; \\
&amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; \color{red}{$\searrow$} &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; \\
G &amp; -30 &amp; $\rightarrow$ &amp; -22 &amp; $\rightarrow$ &amp; -12 &amp; $\rightarrow$ &amp; -4 &amp; $\rightarrow$ &amp; \color{red}{-12} &amp; $\rightarrow$ &amp; -22 &amp; \\
&amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; $\searrow$ &amp; $\downarrow$ &amp; \color{red}{$\searrow$} &amp; $\downarrow$ &amp; \\
T &amp; -40 &amp; $\rightarrow$ &amp; -32 &amp; $\rightarrow$ &amp; -22 &amp; $\rightarrow$ &amp; -14 &amp; $\rightarrow$ &amp; -9 &amp; $\rightarrow$ &amp; \color{red}{-17} &amp;
\end{tabular}
\end{frame}

\begin{frame}[fragile]
\frametitle{Output}
Plot the optimal alignment:
\begin{verbatim}
CTAGG
*| |*
TT-GT
\end{verbatim}

Score: -17

Complexity in time: $O(nm)$

Complexity in memory: $O(nm)$
\end{frame}

\begin{frame}
\frametitle{Acknowledgments}
Bellman; Levenstein; Needleman and Wunsch; Sankoff and Sellers; Hirschberg; Smith and Waterman; Gotoh; Ukkonen, Myers and Fickett; and many others...

\vspace{1cm}
Want to know more? start reading!
\url{http://lectures.molgen.mpg.de/online_lectures.html}
\end{frame}

\end{document}


Significativité d’un alignement pour BLAST

12 juillet 2009

On a vu il y a quelque temps comment estimer la significativité d’un alignement en utilisant la formule de Bayes. Mais pratiquement, ce n’est pas la méthode la plus utilisée. Il est temps d’introduire BLAST…

Derrière cet acronyme percutant se cache l’algorithme le plus utilisé en alignement de séquences biologiques. Il correspond à une heuristique d’alignement local deux-à-deux présenté dans l’article « Basic local alignment search tool » dont il tire son nom. Cet article d’Altschul, Gish, Miller, Myers et Lipman, paru en 1990, a fait date: à ce jour il a été cité dans au moins 30000 publications, un record! Il est très utilisé parce qu’il est rapide et permet de trouver les alignements significatifs entre une sequence inconnue et une (très) grande banque de séquences connues.

Cet algorithme fonctionne en plusieurs étapes. Étant donné deux séquences à comparer, il commence par rechercher des segments de longueur égale, sans gap,  identiques entre eux. On appelle ça les seeds, ce sont des k-mers avec k=11 pour des séquences nucléotidiques. Puis l’algorithme essaie d’étendre ces seeds en autorisant les gaps: il obtient alors des HSPs (pour high-scoring segment pairs). Finalement, pour chaque HSP (donc pour chaque alignement local entre les deux séquences),  l’algorithme choisit ceux dont le score est considéré comme étant significatif.

Optimisation of the Smith-Waterman algorithm

Supposons que l’on aligne une séquence de fonction inconnue, la query, de taille m, avec une séquence de fonction connue, la subject, de taille n, via BLAST. On obtient N HSPs, chacun ayant un score H. Pour estimer la significativité d’un HSP, on a besoin de savoir à quel point c’est probable d’obtenir un HSP de score H par hasard, sans que les deux séquences aient quoi que ce soit à voir entre elles. Cela revient à calculer la probabilité que le score maximal des N HSPs obtenus dépasse un seuil S donné.

Maintenant considérons que les H_i sont des variables aléatoires indépendamment distribuées selon une densité de probabilité f(H).  Si on prend leur maximum X = \max_i ( H_i ), on peut écrire:

P_N ( X \le S ) = P(H_1 \le S) ... P(H_N \le S)

P_N ( X \le S ) = \left( \int_{-\infty}^{S} p(H) \, dH \right)^N

P_N ( X \le S ) = \left( 1 - \int_{S}^{+\infty} p(H) \, dH \right)^N

Si les deux séquences ont beaucoup d’HSPs entre elles (N grand) et comme S est vraisemblablement élevé, il est dans la queue de la distribution et donc la surface sous l’intégrale est petite. On peut donc appliquer l’approximation bien connue (1-x)^k \approx \exp^{-kx} quand x est petit:

P_N ( X \le S ) \approx \exp \left( - N \int_{S}^{+\infty} p(H) \, dH \right)

Maintenant, supposons pour l’instant que p(H) tombe exponentiellement sur ces extrêmes (les queues de distribution, quand S est suffisamment grand):

\int_{S}^{+\infty} p(H) dH = \int_{S}^{+\infty} a \exp^{-\lambda H} dH = \frac{a}{\lambda} \exp^{-\lambda S}

Ainsi: P_N ( X \le S ) \approx \exp \left( - \frac{Na}{\lambda} \exp^{-\lambda S} \right)

Finalement, la probabilité que le score maximal X des N HSPs obtenus entre nos deux séquences dépasse un seuil S donné vaut:

P_N ( X > S ) \approx 1 - \exp \left( - \frac{Na}{\lambda} \exp^{-\lambda S} \right)

Cette formule nous permet donc de dire quand le score X tombe dans les 1% ou 5% de la distribution. C’est bien ce qu’on cherchait, non ?

Pour la petite histoire, X suit une distribution dite « des valeurs extrêmes ». Vous pouvez trouver ici l’article de Gumbel, en français, décrivant les premiers résultats. On se sert de cette théorie dans de nombreux domaines, notamment en finance et assurance, pour estimer la variabilité d’un système et donc les risques d’atteindre des valeurs extrêmes. Les lois des valeurs extrêmes ont l’une de leur queues de distribution plus épaisse que la loi normale, donc la probabilité d’occurrence d’une valeur extrême, c’est-à-dire tombant dans cette queue de distribution plus épaisse, est plus forte que pour une loi normale. Mais certains aiment quand même bien la loi normale, et un prix Nobel ça inspire confiance, quoique après ça, certains vont réviser leur modèles…

Loi Normale et loi de GumbelLa 1e image est tirée de « The estimation of statistical parameters for local alignment score distribution », Althschul et al., Nucleic Acids Research, 2001. Quant à la 2e image, voici le code R pour la générer:


require( distrEx )
require( Cairo )
Cairo( "normal_vs_gumbel.png", width=700, height=500, bg="white" )
x <- seq(-4,6,0.01)
par( mar=c(5,5,3,2), font=2, font.axis=2, font.lab=2, cex=1.5, lwd=2 )
plot( x, dnorm(x), type="l", col="black", main="Loi Normale et loi de Gumbel", ylab="densité de probabilité", ylim=c(0.0,0.5), xlim=c(-4,6) )
points( x, dgumbel(x), type="l", col="red" )
abline( v=0, lty=2, lwd=1 )
legend( 0.2, 0.5, c( "densité de la loi Normale", "densité de la loi de Gumbel" ), col=c("black","red"), lwd=2, bty="n" )
dev.off()

Ajout: un tutoriel sur BLAST vient de sortir dans PLoS Biology, avec une figure qui résume très bien les différentes étapes.


Significativité d’un alignement en bayésien

14 mai 2009

On a vu précédemment comment aligner (globalement) deux séquences l’une avec l’autre et comment calculer le score de l’alignement résultant étant donné un scoring scheme (match, mismatch, gap). Mais bien que l’on sache obtenir un alignement optimal, cela ne nous dit pas si notre alignement est significatif, c’est-à-dire s’il est biologically meaningful

Replaçons nous dans le contexte; nous venons d’obtenir (par des moyens expérimentaux) une nouvelle séquence. Pour en savoir plus sur l’information qu’elle code, on se dit qu’on peut essayer de l’aligner avec une séquence déjà connue. C’est ce que nous faisons et on trouve l’alignement optimal A_1 de score S_1. C’est bien, mais pas suffisant. En effet, nous n’avons (pour l’instant) aucun moyen de savoir si cet alignement est dû au hasard ou bien si ces deux séquences sont homologues, c’est-à-dire si elles ont un ancêtre commun.

Et oui, on peut très bien avoir un alignement A_2 entre notre séquence inconnue et une 2e séquence que l’on connaît déjà, de score S_2 celui-là. Or on aimerait bien avoir un moyen pour juger que deux séquences sont homologues puisque dans ce cas-là on se permettra (grossièrement) d’inférer que la fonction de notre séquence inconnue est la même que la fonction de la séquence qu’on connaît déjà (à peu de choses près).  Dans notre cas, quel alignement doit-on utiliser pour inférer la fonction ? A_1 ou A_2 ?

Et comme toujours, il existe plusieurs approches pour (tenter de) répondre à cette question, l’approche bayésienne (voir ci-dessous) et la distribution du maximum de scores (prochainement dans un autre billet).

On a vu que le score S d’un alignement était un log de ratio de vraisemblances:

S = log( \frac{P(x,y/M)}{P(x,y/R)} )

Mais en fait, au lieu de P(x,y/M), la probabilité d’obtenir les séquences x et y sachant qu’elles ont un ancêtre commun (le match modèle M), ce que l’on veut vraiment, c’est calculer la probabilité que les deux séquences aient un ancêtre commun connaissant ces deux séquences: P(M/x,y). Grâce à la formule de Bayes, on peut calculer cette probabilité:

P(M/x,y) = \frac{P(x,y/M)P(M)}{P(x,y)}

Cependant, au préalable, il faut spécifier les prior: P(M) la probabilité a priori que les deux séquences sont reliées (par un ancêtre commun) et P(R) la probabilité a priori du modèle aléatoire:  P(R) = 1 - P(M).

Ainsi:

P(M/x,y) = \frac{P(x,y/M)P(M)}{P(x,y/M)P(M)+P(x,y/R)P(R)}

P(M/x,y) = \frac{P(x,y/M)P(M) / P(x,y/R)P(R)}{1+P(x,y/M)P(M) / P(x,y/R)P(R)}

Finalement, on pose:

S' = S + log( \frac{P(M)}{P(R)} )

On obtient alors:

P(M/x,y) = \sigma(S') avec \sigma(x)=\frac{\exp^x}{1+\exp^x}

\sigma est la fonction logistique. C’est une sigmoide tendant vers 1 quand x tend vers \infty, et vers 0 quand x tend vers -\infty. Cette fonction est bien utile pour convertir des scores faits de sommes en probabilités.

Une fois qu’on a fait tout ça, on peut comparer la valeur finale de P(M/x,y) à 0 pour savoir si nos séquences sont reliées entre elles ou non. Sans rentrer dans les détails, il faut être bien sûr de travailler avec des probabilités (leur somme sur toutes les paires possibles de séquences doit faire 1), ce qui n’est pas toujours le cas dans un scoring scheme construit de façon ad hoc. De plus, les priors prennent beaucoup d’importance quand on cherche un alignement significatif entre une séquence inconnue et une banque de séquences connues. Dans ce cas-là, il faut prendre en compte la taille de la banque puisque la probabilité d’obtenir, par chance, un bon alignement entre deux séquences non reliées entre elles augmente avec la taille de la banque.


L’alignement global deux-à-deux

29 mars 2009

Le titre du billet ne vous évoque probablement rien, mais peut-être en sera-t-il autrement après avoir lu la suite, enfin j’espère parce que l’article est long et qu’il introduit des concepts majeurs, indispensables à la compréhension des méthodes utilisées dans les recherches portant sur l’évolution des génomes !

Comme vu précédemment (ici), l’alignement de séquences biologiques est une tâche fondamentale à l’heure de la génomique. En effet, à notre époque, on peut séquencer le génome d’organismes vivants « à haut débit », de manière quasi-industrielle (le « quasi » est même superflu dans certains cas). Mais obtenir toutes ces séquences n’a pas pour seul but de les stocker dans un ordinateur bien sûr: on séquence de l’ADN dans le but de répondre à des questions scientifiques. Supposons que l’on ait deux séquences d’ADN à notre disposition, ces fragments d’ADN codant chacun pour une protéine permettant de fixer l’oxygène. Or ces deux fragments sont légèrement différents: l’un code pour l’hémoglobine (dans le sang) et l’autre pour la myoglobine (dans les muscles). Comment fait-on pour « voir » ces différences ? Où sont-elles le long de la séquence ? Ces différences au niveau de l’ADN ont-elles des conséquences au niveau de la séquence d’acides aminés ? Ces différences sont-elles importantes pour la fonction de ces protéines ? …

On se rend bien compte qu’il nous faut trouver un moyen pour comparer deux séquences d’ADN. On peut commencer par dire qu’il y a x% de A dans la première séquence alors qu’il y en a seulement y% dans la deuxième, mais ça ne nous avancera pas à grand chose puisque c’est l’enchaînement des nucléotides qui est important. De plus, on sait qu’au cours de l’évolution des mutations apparaissent le long du génome, par exemple un A remplacé par erreur par un T. Il faut donc essayer de repérer le long de la séquence où on eu lieu ces mutations. Dans le cas A->T, on parle de substitutions, c’est une mutation ponctuelle, par opposition aux insertions-délétions (appelés « indels« ) de plusieurs nucléotides comme AATGTCA->AACA où TGT a disparu.

Et bien figurez-vous que ce type de question, les informaticiens se les étaient déjà posées, avant même que l’on sache que l’ADN était le support de l’information génétique. En 1950, Hamming proposait une distance pour quantifier les substitutions entre deux chaines de symboles de même longueur: d_H=2 entre AATGC et ATTGA (A_2->T_2 et C_5->A_5). En 1965, Levenshtein étendait la mesure aux cas de séquences de longueur différente pour prendre en compte les insertions-délétions. Le problème revient alors à trouver le nombre de changements qu’il faut faire subir à une séquence pour la transformer en l’autre, par exemple AATGC -> ATTAGCA. Ainsi, on cherche à calculer la distance d’édition, le nombre d’opérations parmi « substitution », « insertion » et « délétion », qui fait passer d’une chaine de caractères (string en anglais) à une autre. En informatique, ce problème revient à trouver la plus longue séquence en commun entre deux chaines de caractères (LCS pour longest common substring).

Le problème LCS a deux caractéristiques majeures. Il a une sous-structure optimale: le problème peut être décomposé en sous-problèmes plus simples, qui eux-mêmes peuvent être décomposés en sous-sous-problèmes également plus simples, et ainsi de suite jusqu’à ce que la solution soit triviale à trouver. De plus, il a ce qu’on appelle des sous-problèmes chevauchants: la solution d’un sous-problème dépend des solutions aux sous-problèmes de niveau inférieur. De tels problèmes sont résolus grâce à une technique appelée la « programmation dynamique« : une solution optimale du problème peut être construite efficacement à partir de solutions optimales à des sous-problèmes. Pour cela on a besoin d’une procédure appelée « mémoisation » (et non pas « mémorisation »!): sauveguarder dans une table les solutions à un niveau de sous-problèmes de telle sorte qu’elles soient disponibles pour le niveau suivant de sous-problèmes. La programmation dynamique a été inventée par Richard Bellman dans les années 1940-1950. Cette méthode est plus rapide que les méthodes naïves de comparaison de deux séquences, et elle est exacte: contrairement aux méthodes heuristiques qui font des approximations pour aller plus vite, la programmation dynamique trouve la solution optimale du problème et rapidement. Et pour ces raisons, ça déboîte !

Levenshtein, quand il a généralisé la distance de Hamming, a également proposé un algorithme de programmation dynamique pour la calculer. Enfin, on y arrive, c’est en 1970 que Needleman et Wunsch vont comparer deux séquences de protéines à l’aide des méthodes décrites ci-dessus (ça marche pareil pour des séquences d’ADN). Et le reste du billet va justement montrer comment on fait, en pratique.

On suppose que les deux séquences sont homologues, c’est-à-dire qu’elles ont évoluées à partir d’un ancêtre commun, et que les différences entre elles ne sont dues qu’à deux types d’événements, des indels et des substitutions. Pour les aligner l’une avec l’autre, on applique le principe de parcimonie: l’alignement nécessitant le moins d’événements est dit « optimal », c’est la solution. Comme les deux types d’événements arrivent à des fréquences différentes, on utilise une parcimonie pondérée. Par exemple, à une certaine position, un match (deux lettres identiques) vaut +1, un mismatch (deux lettres différentes) vaut 0, et un gap vaut -1: c’est le scoring scheme. Needleman et Wunsch ont donc cherché l’alignement qui maximisait un score (la similarité entre les deux séquences), celui-ci étant la somme des matches, mismatches et pénalités de gap le long de l’alignement.

Quelques notations: on cherche à aligner la séquence s_1 de longueur l_1 avec la séquence s_2  de longueur l_2. Le i-ème élément de s_1 est noté s_1[i] et s_1[3:7] correspond à la sous-séquence de s_1 allant de l’élément 3 à l’élément 7. Le score est ici la similarité entre les deux séquences, notée S_{ij} pour l’alignement de s_1[1:i] avec s_2[1:j].

Comment marche la programmation dynamique? On commence par créer une matrice F dont l’entrée F_{i,j} correspond au score  de l’alignement entre s_1[1:i] et s_2[1:j]. Ensuite, on remplit la matrice récursivement en commençant en haut à gauche: remplir F(i,j) ne dépend que de 3 cases remplies précédemment F(i-1,j-1), F(i,j-1) et F(i-1,j). Bien sûr, il y a des conditions au bornes: dans la première ligne et la première colonne, on ne met que des gaps.

matrixf

Pour calculer F(i,j), on commence par  déterminer le score de s_1[i] aligné avec s_2[j] (match ou mismatch), puis on applique la règle de décision suivante:

equationfillf

Une fois que la matrice est remplie, on fait ce qu’on appelle le traceback. On part d’en bas à droite de la matrice, case avec le score le plus élevé (similarité de l’alignement optimal), et on remonte en suivant le meilleur chemin: ce chemin dans la matrice correspond à l’alignement optimal. En fin de compte, on dispose bien d’une méthode pour comparer quantitativement deux séquences d’ADN.

Prenons un exemple simple: AAT et ACAT.

1. on définit le scoring scheme:

match=2 ; mismatch=0 ; gap=-1

2. on remplit la première ligne et la première colonne avec des gaps:

fillf_row1col1

3. on remplit le reste de la matrice avec la formule donnée précédemment:

fillf1

4. on fait le traceback:

traceback

5. on montre l’alignement obtenu, avec un score de 5:

alignJe précise que cette méthode marche aussi pour les protéines, simplement, au lieu que l’alphabet soit de 4 lettres {A;T;G;C}, il est de 20 lettres, une par acides aminés. De plus, la comparaison de deux séquences peut avoir plusieurs alignements optimaux. Et pour finir d’expliquer le titre, on parle ici d’alignement global parce qu’on essaie d’aligner les deux séquences sur toute leur longueur. Vous vous en doutez, il existe aussi l’alignement local…

Maintenant qu’on dispose d’un algorithme pour comparer deux séquences, on peut s’amuser un peu à changer les paramètres (oui, on s’amuse de peu en ce bas monde). Avec les paramètres match=2, mismatch=0 et gap=-1, on obtient l’alignement suivant, de score 67:

aligngap-1

En passant la pénalité de gap à -2, on obtient un alignement tout différent, de score 58:

aligngap-2

Et en mettant le gap à -3, on obtient encore un autre alignement, de score 52 celui-là:

aligngap-3

Alors, lequel est le meilleur ? Chacun est optimal, mais cette optimalité est fonction du scoring scheme choisi, donc garder celui de score le plus élevé n’est pas forcément le plus judicieux. On peut aussi garder celui qui est le plus plaisant à l’oeil (vous rigolez mais au début les chercheurs faisaient comme ça !). Une autre possibilité est se fixer un scoring scheme qui semble relativement correct et de garder le même tout le temps. Enfin bon, tout cela n’est pas très satisfaisant, vous en conviendrez. Le vrai problème derrière tout ça consiste à écrire un modèle d’évolution des séquences et de l’utiliser lors de l’alignement: on garderait dans ce cas l’alignement ayant la vraisemblance la plus élevée (au sens statistique du terme: l’alignement qui maximise la vraisemblance d’observer ces séquences sous un modèle d’évolution donné). Mais ce sera pour d’autres billets…


%d blogueurs aiment cette page :