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}
& & & C & & T & & A & & G & & G \\ \hline
& & & & & & & & & & & & \\
& \visible{0} & \visible{$\rightarrow$} & \visible{-10} & \visible{$\rightarrow$} & \visible{-20} & \visible{$\rightarrow$} & \visible{-30} & \visible{$\rightarrow$} & \visible{-40} & \visible{$\rightarrow$} & \visible{-50} \\
& \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \\
T & \visible{-10} & \visible{$\rightarrow$} & \visible{-2} & \visible{$\rightarrow$} & \visible{-10} & \visible{$\rightarrow$} & \visible{-20} & \visible{$\rightarrow$} & \visible{-30} & \visible{$\rightarrow$} & \visible{-40} \\
& \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \\
T & \visible{-20} & \visible{$\rightarrow$} & \visible{-12} & \visible{$\rightarrow$} & \visible{-2} & \visible{$\rightarrow$} & \visible{-12} & \visible{$\rightarrow$} & \visible{-22} & \visible{$\rightarrow$} & \visible{-32} & \\
& \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \\
G & \visible{-30} & \visible{$\rightarrow$} & \visible{-22} & \visible{$\rightarrow$} & \visible{-12} & \visible{$\rightarrow$} & \visible{-4} & \visible{$\rightarrow$} & \visible{-12} & \visible{$\rightarrow$} & \visible{-22} & \\
& \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \visible{$\searrow$} & \visible{$\downarrow$} & \\
T & \visible{-40} & \visible{$\rightarrow$} & \visible{-32} & \visible{$\rightarrow$} & \visible{-22} & \visible{$\rightarrow$} & \visible{-14} & \visible{$\rightarrow$} & \visible{-9} & \visible{$\rightarrow$} & \visible{-17} &
\end{tabular}
\end{frame}
\begin{frame}
\frametitle{Traceback}
\begin{tabular}{l|cccccccccccc}
& & & C & & T & & A & & G & & G \\ \hline
& & & & & & & & & & & & \\
& \color{red}{0} & $\rightarrow$ & -10 & $\rightarrow$ & -20 & $\rightarrow$ & -30 & $\rightarrow$ & -40 & $\rightarrow$ & -50 \\
& $\downarrow$ & \color{red}{$\searrow$} & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & \\
T & -10 & $\rightarrow$ & \color{red}{-2} & $\rightarrow$ & -10 & $\rightarrow$ & -20 & $\rightarrow$ & -30 & $\rightarrow$ & -40 \\
& $\downarrow$ & $\searrow$ & $\downarrow$ & \color{red}{$\searrow$} & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & \\
T & -20 & $\rightarrow$ & -12 & $\rightarrow$ & \color{red}{-2} & \color{red}{$\rightarrow$} & \color{red}{-12} & $\rightarrow$ & -22 & $\rightarrow$ & -32 & \\
& $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & \color{red}{$\searrow$} & $\downarrow$ & $\searrow$ & $\downarrow$ & \\
G & -30 & $\rightarrow$ & -22 & $\rightarrow$ & -12 & $\rightarrow$ & -4 & $\rightarrow$ & \color{red}{-12} & $\rightarrow$ & -22 & \\
& $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & $\searrow$ & $\downarrow$ & \color{red}{$\searrow$} & $\downarrow$ & \\
T & -40 & $\rightarrow$ & -32 & $\rightarrow$ & -22 & $\rightarrow$ & -14 & $\rightarrow$ & -9 & $\rightarrow$ & \color{red}{-17} &
\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}
Publié par walrus 
-mers avec 
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.
, avec une séquence de fonction connue, la subject, de taille 
, via BLAST. On obtient 
HSPs, chacun ayant un score 
. 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 
donné.
sont des variables aléatoires indépendamment distribuées selon une densité de probabilité 
. Si on prend leur maximum 
, on peut écrire:
quand 
est petit:
tombe exponentiellement sur ces extrêmes (les queues de distribution, quand 
des 
de score 
. 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.
entre notre séquence inconnue et une 2e séquence que l’on connaît déjà, de score 
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 ? 
, la probabilité d’obtenir les séquences 
sachant qu’elles ont un ancêtre commun (le match modèle 
), 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: 
. Grâce à la formule de Bayes, on peut calculer cette probabilité:
la probabilité a priori que les deux séquences sont reliées (par un ancêtre commun) et 
la probabilité a priori du modèle aléatoire: 
.
avec 
est la fonction logistique. C’est une sigmoide tendant vers 1 quand 
, et vers 0 quand 
. Cette fonction est bien utile pour convertir des scores faits de sommes en probabilités.
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 
), 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.
entre AATGC et ATTGA (
et 
). 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).
de longueur 
avec la séquence 
de longueur 
. Le i-ème élément de ![s_1[i]](http://s0.wp.com/latex.php?latex=s_1%5Bi%5D&bg=ffffff&fg=333333&s=0 "s_1[i]")
et ![s_1[3:7]](http://s0.wp.com/latex.php?latex=s_1%5B3%3A7%5D&bg=ffffff&fg=333333&s=0 "s_1[3:7]")
correspond à la sous-séquence de 
pour l’alignement de ![s_1[1:i]](http://s0.wp.com/latex.php?latex=s_1%5B1%3Ai%5D&bg=ffffff&fg=333333&s=0 "s_1[1:i]")
avec ![s_2[1:j]](http://s0.wp.com/latex.php?latex=s_2%5B1%3Aj%5D&bg=ffffff&fg=333333&s=0 "s_2[1:j]")
.
correspond au score de l’alignement entre 
![s_2[j]](http://s0.wp.com/latex.php?latex=s_2%5Bj%5D&bg=ffffff&fg=333333&s=0 "s_2[j]")
(match ou mismatch), puis on applique la règle de décision suivante:
Je 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…
