Dans cette page, on s'intéresse au problème de l'alignement de
gènes, qui sont des séquences d'ADN codant une protéine ou
un caractère. De telles séquences peut être modélisée par
une suite de bases azotées notées A, C, T et G (que nous appelerons «lettres» dans la suite).
Pour établir les relations de parenté entre un ensemble
d'espèces, on prend un même gène chez ces différentes
espèces, et on étudie les mutations survenues entre ces
gènes chez ces différentes espèces. On considère que deux
espèces sont d'autant plus proches parentes qu'il y a peu
de mutations entre leurs gènes.
Lorsque l'on compare une séquence $(1)$ à une séquence $(2)$, le terme mutation peut
se référer (entre autres) à trois choses différentes :
Une substitution correspond au remplacement d'une lettre par une autre.
Une délétion correspond à la disparition d'une lettre.
Une insertion correspond à l'apparition d'une lettre.
Alignement de deux séquences
Considérons les deux séquences jouet suivantes :
ACACTTACAGCA (1)
ACACTACAGGA (2)
Il est naturel de les aligner de la façon suivante :
ACACTTACAGCA (1)
ACACTA_CAGGA (2)
Cet alignement met en évidence le fait que l'on a pu
passer de la séquence $(1)$ à la $(2)$ avec une délétion
(du 3ème A de (1) - matérialisée par un tiret bas (_) dans (2)) et d'une substitution (du
dernier C de (1), substitué par un G). Le coût de
cet alignement est $1c_s+1c_d$, où $c_s > 0$ désigne le coût d'une
substitution, et $c_d>0$ celui d'une délétion. La question qui se pose alors est comment
calculer l'alignement ayant le coût le plus bas ? Ce
coût minimal sera considéré comme la distance entre les
deux séquences de départ (il s'agit de ce qu'on appelle la distance de Levenshtein entre les deux séquences). À noter qu'il existe dans la littérature de nombreuses autres façons de définir le coût. On peut donner des coûts différents aux différentes substitutions, donner un coût non-linéaire à des délétions successives (en pratique, $2$ délétions successives sont beaucoup plus probables que $2$ délétions séparées), etc. Dans la suite de cette page, on reste sur le modèle le plus simple décrit ci-dessus.
La résolution de ce problème peut se ramener à une
recherche de plus court chemin dans un graphe tel que
celui apparaissant dans l'application ci-dessous.
Les deux
séquences $(1)$ et $(2)$ à aligner apparaissent
respectivement en abscisses et en ordonnées.
Chaque sommet
du graphe correspond à un couple constitué d'un début de (1) (on
parle plutôt de préfixe) et d'un préfixe de $(2)$.
On obtient le préfixe de $(1)$ en regardant les lettres de
$(1)$ à gauche du sommet, et le préfixe de $(2)$ en
regardant les lettres de $(2)$ au dessus du sommet. Ainsi,
le sommet tout en haut à gauche (sommet de départ)
correspond aux préfixes vides des deux séquences. Le
sommet tout en bas à droite (sommet d'arrivée) correspond
aux séquences entières.
Les arcs entre les sommets ont différentes significations selon leur orientation.
Les arcs horizontaux correspondent à une délétion (et ont une longueur égal au coût d'une insertion/délétion)
Les arcs verticaux correspondent à une insertion (et ont une longueur égal au coût d'une insertion/délétion)
Les arc diagonaux sont de deux types :
si les lettres correspondantes dans les séquences $(1)$ et $(2)$ sont égales, l'arc correspond au cas où la lettre n'a pas muté, et est de longueur 0.
si les lettres correspondantes dans les séquences $(1)$ et $(2)$ sont différentes, l'arc correspond à une substitution (et a donc un longueur égal au coût d'une substitution)
Les longueurs sont matérialisée par l'épaisseur du trait.
Un chemin du sommet de départ vers le sommet
d'arrivée correspond à un alignement de (1) et (2). Sa
longueur est par construction le coût de l'alignement.
Le plus court chemin correspond donc à
l'alignement de coût minimal entre les deux séquences.
Séquence (1) : Séquence (2) :
Algorithme de Dijkstra
Notons $n$ et $m$ la longueur des deux séquences à
aligner. Le graphe obtenu à $O(nm)$ sommets, et $O(nm)$
arcs. Calculer un plus court chemin en utilisant l'algorithme de Dijkstra est
alors de complexité $O(nm\log(nm))$. Si les deux séquences
sont de même longueur $m\approx n$, l'algorithme obtenu
pour leur alignement est donc de complexité
$O(n^2\log(n))$ en temps (et $O(n^2)$ en espace).
Algorithme A*
On peut améliorer la méthode précédente en utilisant
l'algorithme A* -
avec pour heuristique le coût d'une insertion/délétion
multiplié par le nombre d'insertions/délétions nécessaires
pour ramener les deux séquences à la même longueur
(indépendamment donc des substitutions). On peut montrer
qu'il s'agit d'une heuristique monotone, qui mène
empiriquement à un gain de la moitié des itérations
environ.
Algorithme de Wagner et Fischer
Il est en fait possible de calculer un plus court chemin
plus rapidement, sans faire appel à des algorithmes
généraux tels que Dijkstra.
On créé un tableau tel que celui ci-dessous - dans chaque
case $(i,j)$, on va mettre la distance $d_{i,j}$ entre les
suffixes correspondants à cette case (i.e. les parties des
deux séquences respectivement situées à droite et en
dessous de cette case).
Dans les explications suivantes, on suppose que coût de
substitution vaut 3, et que le coût d'insertion/délétion
vaut 4.
On initialise la case tout en bas à droite à 0 (la
distance entre la chaîne vide et elle-même).
On initialise la dernière colonne et la dernière ligne. À chaque fois, il s'agit de comparer la chaîne vide à une chaîne d'une certaine longueur $l$ leur distance est donc $lc_d$.
On remplit les cases de la droite vers la gauche
et du bas vers le haut, en utilisant la formule suivante :
$$c_{(i,j)}\longleftarrow \min(c_{(i+1,j)}+c_d,\ c_{(i,j+1)}+c_d,\
c_{(i+1,j+1)}+\epsilon_{i,j}c_s),$$ où $\epsilon_{i,j}$ vaut 0 si les
caractères $i+1$ de la première séquence et $j+1$ de la seconde sont
égaux, et 1 sinon.
Pour obtenir le chemin correspondant à l'alignement
optimal, on part alors de la case en haut à gauche (qui
contient donc la distance entre les deux séquences), et
l'on se déplacer de cellule en cellule en suivant à chaque
fois une cellule dont est originaire le minimum calculé
dans la formule ci-dessus.
La complexité de cet algorithme pour deux séquences de
longueur $n$ est $O(n^2)$. Il est en pratique beaucoup
plus efficaces que ceux présentés ci-dessus.
On peut également avec un peu plus de travail faire
tomber la complexité spatiale à $O(n)$ - le tableau de
taille $O(n^2)$ posant un problème de mémoire dès que les
séquences sont longues (pour $n=10\;000$ - ce qui est dans
la moyenne d'un gène humain, le tableau a déjà environ
$10^{8}$ cases. Et il existe des gènes de plus de 2 millions de lettres…)
Alignement multi-séquences
On souhaite maintenant aligner $p\geq 3$ séquences de
longueur $\approx n$. Considérons par exemple les séquences suivantes :
ACCTAAGCTT (1)
AGCAAGCTT (2)
ACCTAAGGCTA (3)
et un alignement correspondant :
ACCTAAG_CTT (1)
AGC_AAG_CTT (2)
ACCTAAGGCTA (3)
On définit le coût de cet alignement comme étant égal au
sommes des coûts des alignements des séquences prises deux
à deux. Sur l'exemple ci-dessus, avec les notations introduites précédemment,
Le coût de l'alignement de $(1)$ et $(2)$ est $c_s+c_d$
Le coût de l'alignement de $(1)$ et $(3)$ est $c_s+c_d$
Le coût de l'alignement de $(2)$ et $(3)$ est $2c_s+2c_d$
Le coût total de l'alignement des trois séquences est donc
$4c_s+4c_d$. La question est ici aussi de trouver un
alignement de coût minimal.
On peut reprendre l'approche présentée pour l'alignement
de deux séquences. On construit un graphe dont chaque
sommet correspond à un $p$-uplet de préfixes des $p$
séquences à aligner. Le graphe obtenu est alors de taille
$O(n^p)$, et la recherche de plus court chemin à l'aide de
l'algorithme de Dijkstra de complexité temporelle
$O(n^p\log(n^p))$. L'espace mémoire utilisé est rapidement
prohibitif, ainsi que la complexité temporelle. Pour un
alignement de 4 séquences de $10^4$ lettres (ce qui est
raisonnable au regard des besoins de biologistes), on a déjà
$n^p=10^{16}$…
L'idée est donc de ne pas précalculer l'intégralité du
graphe, et d'utiliser des algorithmes d'exploration guidés. L'application ci-dessous permet de tester notamment A* (et également Best First et Beam Search, qui sont nettement moins rapide, et donnent des résultats fortement sous-optimaux pour ce problème), avec les heuristiques suivantes :
«longueurs» il s'agit de l'heuristique décrite ci-dessus - adaptée le cas échéant à plus de deux séquences.
«paires». On commence par utiliser l'algorithme
de Wagner et Fischer pour précalculer les tableaux
de distances entre les chaque paire de séquences. On
prend alors comme heuristique la somme (sur chaque
paire de séquences) des distances des suffixes
(i.e. des séquences restant à aligner à partir de la
position courante dans le graphe).
Cette heuristique semble très efficace en
pratique, et permet de trouver dans des temps très
raisonnables des alignements des 4 séquences de
plusieurs milliers de lettres, difficilement
abordables avec les autres méthodes codées dans
l'application… (typiquement, sur 4 séquences de longueur 1000 avec 3% de mutations, elle réduit d'un facteur 1 million l'espace exploré).