L’alignement global deux-à-deux

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…

Laisser un commentaire

Entrez vos coordonnées ci-dessous ou cliquez sur une icône pour vous connecter:

Logo WordPress.com

Vous commentez à l'aide de votre compte WordPress.com. Déconnexion / Changer )

Image Twitter

Vous commentez à l'aide de votre compte Twitter. Déconnexion / Changer )

Photo Facebook

Vous commentez à l'aide de votre compte Facebook. Déconnexion / Changer )

Photo Google+

Vous commentez à l'aide de votre compte Google+. Déconnexion / Changer )

Connexion à %s

%d blogueurs aiment cette page :