Algorithme de Needleman et Wunsch

Sommaire

Principes de l'algorithme
Initialisation
Remplissage de la matrice
Déterminer l'alignement optimal
L'algorithme
Le programme
Le module de manipulation de matrices

Principes de l'algorithme

Dans un autre article de ce site sont présentés des algorithmes de recherche d’un mot dans un texte, notamment celui de Knuth-Morris-Pratt (KMP). Ces algorithmes sont dévolus au problème de la recherche exacte : il s’agit de trouver, si elle existe, la première occurrence exacte de ce mot dans ce texte.

Nous allons maintenant étudier, parce que c’est un problème central en bioinformatique, une recherche approximative : il s’agit de savoir si deux mots se ressemblent, quel est leur degré de ressemblance, ou de trouver, dans un ensemble de mots, celui qui ressemble le plus à un mot-cible. Et nous allons voir que ce problème relève de solutions très différentes de celles qui valent pour la recherche exacte.

Notons d’abord que la ressemblance est une notion imprécise : la plupart des algorithmes utilisés proposent différents paramètres pour ajuster les facteurs de ressemblance aux caractéristiques du problème traité.

Les algorithmes utilisés fournissent en général deux résultats :

pour chaque comparaison de deux chaînes, un score de ressemblance, qui permet ensuite de trouver la meilleure ressemblance parmi un ensemble de comparaisons ;
un alignement des deux chaînes (qui n’ont pas forcément la même longueur) selon la configuration qui procure le meilleur score ; on dit bien un alignement, et non pas l’alignement, parce qu’en effet, comme nous le verrons plus loin, le problème peut admettre plusieurs solutions conduisant au même score.

Le plus caractéristique de cette famille d’algorithmes est peut-être celui de Needleman et Wunsch, que nous étudierons ici ; il réalise un alignement global de deux séquences (chaînes de caractères).

Calculer un alignement global peut être coûteux si les séquences à aligner sont longues, ou s’il y en a beaucoup. D’autres algorithmes, qui ressemblent à celui-ci, ont été conçus pour limiter la taille du problème en ne réalisant l’alignement que pour des régions « intéressantes ». La détermination des régions intéressantes est bien sûr en elle-même un problème intéressant. Citons l’algorithme de Smith et Waterman, qui réalise des alignements locaux, et le logiciel BLAST¹, qui mettent en oeuvre des méthodes similaires à celles de Needleman et Wunsch, après des optimisations éventuellement complexes.

Le problème de la comparaison de séquences est exponentiel, la solution est en O(kⁿ) ; ces algorithmes sont susceptibles d’une multiplicité de solutions ; une des techniques les plus généralement utilisées pour en réduire la complexité est la programmation dynamique, qui fait l’objet d’un autre article sur ce site.

La programmation dynamique résout des problèmes en combinant des solutions de sous-problèmes. (Thomas Cormen, Charles Leiserson, Ronald Rivest et Clifford Stein, Introduction à l’algorithmique)

L’idée de la programmation dynamique est de mémoriser les résultats de calculs intermédiaires qui seront probablement répétés. La programmation dynamique est par exemple souvent un bon choix lorsque l’on aura besoin, après les avoir calculées, des valeurs stockées dans tous les noeuds d’un arbre ou dans toutes les cases d’un tableau. Parfois aussi cette conservation des résultats intermédiaires est imposée par un problème tel que le calcul d’une valeur se fait en fonction de toutes les précédentes. L’art algorithmique consiste à chercher des solutions qui évitent ce type de contrainte mais c’est parfois impossible. Et puis il y a des problèmes intrinsèquement récursifs pour lesquels n’existe pas d’algorithme itératif.
Nous allons donc chercher des procédés pour associer un algorithme qui calcule des valeurs successives avec une structure de données qui les archive.

Cette section doit beaucoup au travail préalable de William Saurin pour cet enseignement, ainsi qu’à une page créée par Eric C. Rouchka, de l’université Washington à Saint-Louis, et reprise par Per Kraulis au Stockholm Bioinformatics Center :

http://www.sbc.su.se/~pjk/molbioinf...

Supposons que nous souhaitions calculer un alignement global de deux séquences :
séquence n° 1 : G A A T T C A G T T A
séquence n° 2 : G G A T C G A
La séquence n° 1 a m=11 résidus, la séquence n° 2 n=7 résidus.

Nous allons ici étudier l’algorithme avec des paramètres particulièrement simples, peut-être même simplistes : pénalité nulle pour les trous (gaps) et les discordances (substitutions), une pénalité négative, ou prime, égale à 1 pour les concordances (matches). Le but est d’acquérir une vue d’ensemble de l’architecture de la solution, qui permettra au lecteur d’envisager ensuite des exemples plus compliqués, avec des formules de calcul plus élaborées pour les scores et pour les pénalités de gap.

Le principe de pondération que nous adopterons sera le suivant :

la « prime de score » pour la comparaison du résidu de rang i de la première séquence avec le résidu de rang j de le seconde séquence sera S_i,j = 1 si les deux résidus sont identiques, sinon :
S_i,j = 0 (score de discordance) ;
w = 0 (pénalité de gap).

L’algorithme opère en trois étapes :

initialisation ;
calcul des scores et remplissage de la matrice ;
calcul de l’alignement en « remontant » dans la matrice.

Initialisation

Création d’une matrice M de m+2=13 colonnes et n+2=9 lignes : la ligne et la colonne de rangs 0 contiendront les textes des séquences, la seconde ligne (de rang 1, les M_1,j) et la première colonne (les M_i,1) de M sont remplies de 0 parce que nous avons posé qu’il n’y avait pas de pénalités pour des gaps initiaux ou finals.

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0

A 0

T 0

C 0

G 0

A 0

Remplissage de la matrice

À chaque position M_i,j de la matrice M (i est le numéro de ligne, j le numéro de colonne) le score se calcule ainsi :

M_i,j = Maximum de :

M_i−1,j−1 + S_i,j (concordance ou discordance

dans la diagonale)

M_i,j−1 + w (gap dans la séquence n° 1)

M_i−1,j + w (gap dans la séquence n° 2)

Nous voyons que pour calculer M_i,j il faut (et il suffit de) connaître M_i−1,j, M_i,j−1 et M_i−1,j−1 ; de ce point de vue le problème est assez analogue à ceux posés par Fibonacci ou par le triangle de Pascal.

Ainsi, comme chaque séquence commence par le résidu G (concordance), S_1,1=1. Nous avons posé par hypothèse w=0. Donc :

M_1,1 = Max[M_0,0+1, M_1,0+0, M_0,1+0] (1)

= Max[1,0,0] (2)

= 1 (3)

Nous pouvons donc inscrire un 1 en M_1,1 :

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1

G 0

A 0

T 0

C 0

G 0

A 0

Ceci fait, toujours parce que w=0, nous pouvons facilement remplir la ligne 1 et la colonne 1 avec des 1 ; ainsi :

M_2,1 = Max[M_1,0+0, M_2,0+0, M_1,1+0] (4)

= Max[0,0,1] (5)

= 1 (6)

soit :

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1 1 1 1 1 1 1 1 1 1 1

G 0 1

A 0 1

T 0 1

C 0 1

G 0 1

A 0 1

Finalement :

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1 1 1 1 1 1 1 1 1 1 1

G 0 1 1 1 1 1 1 1 2 2 2 2

A 0 1 2 2 2 2 2 2 2 2 2 3

T 0 1 2 2 3 3 3 3 3 3 3 3

C 0 1 2 2 3 3 4 4 4 4 4 4

G 0 1 2 2 3 3 4 4 5 5 5 5

A 0 1 2 3 3 3 4 5 5 5 5 6

Nous avons signalé ci-dessus que le problème général de la comparaison de séquences était exponentiel (O(kⁿ)). L’utilisation de la programmation dynamique, avec le graphe représenté par ce tableau, permet de le réduire à un problème quadratique (O(m × n), m et n étant les longueurs respectives des séquences). En effet, il y a m × n valeurs dans la table, et le calcul de chacune s’effectue en temps constant.

Déterminer l'alignement optimal

L’étape précédente nous a déjà permis de savoir que le score d’alignement maximum pour nos deux séquences est 6. Souvent, cette information est suffisante, parce que l’on cherche en fait les meilleurs scores parmi une collection de séquences à comparer à la cible. Mais peut être aussi intéressant de connaître un alignement qui donne ce score.

Nous allons maintenant déterminer l’alignement effectif qui donne ce résultat.

Pour cela, on considère la case du tableau qui contient le score maximum, qui est M_m,n, et on la compare à ses voisines. Ici toutes les voisines contiennent la valeur 5. Comme la différence de scores est 1 dans tous les cas, et que le seul moyen d’avoir un accroissement de 1 est une concordance (match) (toutes les autres situations donnent un accroissement nul), c’est que la case précédente était la voisine en diagonale :

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1 1 1 1 1 1 1 1 1 1 1

G 0 1 1 1 1 1 1 1 2 2 2 2

A 0 1 2 2 2 2 2 2 2 2 2 3

T 0 1 2 2 3 3 3 3 3 3 3 3

C 0 1 2 2 3 3 4 4 4 4 4 4

G 0 1 2 2 3 3 4 4 5 5 5 5

A 0 1 2 3 3 3 4 5 5 5 5 6

Ce qui nous donne un alignement :

Maintenant nous considérons la case courante et cherchons celle qui la précède : c’est la voisine avec le score maximum, soit celle de la même ligne.

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1 1 1 1 1 1 1 1 1 1 1

G 0 1 1 1 1 1 1 1 2 2 2 2

A 0 1 2 2 2 2 2 2 2 2 2 3

T 0 1 2 2 3 3 3 3 3 3 3 3

C 0 1 2 2 3 3 4 4 4 4 4 4

G 0 1 2 2 3 3 4 4 5 5 5 5

A 0 1 2 3 3 3 4 5 5 5 5 6

Cet alignement correspond à un gap dans la séquence n° 2 :

T A

_ A

Encore une fois, le prédécesseur immédiat donne un gap dans la séquence n° 2 :

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1 1 1 1 1 1 1 1 1 1 1

G 0 1 1 1 1 1 1 1 2 2 2 2

A 0 1 2 2 2 2 2 2 2 2 2 3

T 0 1 2 2 3 3 3 3 3 3 3 3

C 0 1 2 2 3 3 4 4 4 4 4 4

G 0 1 2 2 3 3 4 4 5 5 5 5

A 0 1 2 3 3 3 4 5 5 5 5 6

T T A

_ _ A

Au bout du compte :

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1 1 1 1 1 1 1 1 1 1 1

G 0 1 1 1 1 1 1 1 2 2 2 2

A 0 1 2 2 2 2 2 2 2 2 2 3

T 0 1 2 2 3 3 3 3 3 3 3 3

C 0 1 2 2 3 3 4 4 4 4 4 4

G 0 1 2 2 3 3 4 4 5 5 5 5

A 0 1 2 3 3 3 4 5 5 5 5 6

G A A T T C A G T T A

| | | | | |

G G A _ T C _ G _ _ A

Il y a une autre solution possible :

G A A T T C A G T T A

0 0 0 0 0 0 0 0 0 0 0 0

G 0 1 1 1 1 1 1 1 1 1 1 1

G 0 1 1 1 1 1 1 1 2 2 2 2

A 0 1 2 2 2 2 2 2 2 2 2 3

T 0 1 2 2 3 3 3 3 3 3 3 3

C 0 1 2 2 3 3 4 4 4 4 4 4

G 0 1 2 2 3 3 4 4 5 5 5 5

A 0 1 2 3 3 3 4 5 5 5 5 6

qui donne l’alignement suivant :

G _ A A T T C A G T T A

| | | | | |

G G _ A _ T C _ G _ _ A

L'algorithme

Algo : NW Données : s1, s2, S et gap ; des chaînes numérotées de ; 1 à longueur(s1) et de 1 à longueur(s2), ; S le score de concordance entre caractères ; et gap un coût de gap. Résultat : la matrice de calcul des scores Créer C une matrice à longueur(s1) + 2 colonnes et à longueur(s2) + 2 lignes pour j allant de 2 à longueur(s1)+2 faire C[1, j] <- gap * j-2 fait pour i allant de 2 à longueur(s2)+2 faire C[i, 1] <- gap * i-2 pour j allant de 2 à longueur(s2)+2 faire C[i, j] = max(C[i-1,j-1] + si s1[i] = s2[j] alors 1 sinon 0, C[i-1,j] + gap, C[i,j-1] + gap) fait fait retourner C
Nous verrons dans un prochain article l’algorithme de remontée dans le graphe (backtracking) pour trouver un alignement optimal.
Au sujet de ces algorithmes on consultera avec profit le livre de Maxime Crochemore, Christophe Hancart et Thierry Lecroq, Algorithmique du texte, chez Vuibert.

Le programme

(module nw:lb (main main) (import nw:matrices) (import nw:chains) (import nw:alignment)) (define (nw-2 s1 s2 match-bonus gap-penalty) (let ((ncol (+ (chain-length s1) 2)) (nlin (+ (chain-length s2) 2))) (let ((C (make-matrix nlin ncol 0))) (matrix:margins C s1 s2) ;; (do ((j 2 (+ j 1))) ((= j nlin) 'fait) (matrix:set! C j 1 (* j gap-penalty))) (do ((i 2 (+ i 1))) ((= i ncol) 'fait) (matrix:set! C 1 i (* gap-penalty i))) (do ((i 2 (+ i 1))) ((= i ncol) C) (do ((j 2 (+ j 1))) ((= j nlin) 'fait) (let ((val (max (+ (matrix:ref C (- j 1) (- i 1)) (if (char=? (matrix:ref C 0 i) (matrix:ref C j 0)) match-bonus 0)) (+ (matrix:ref C j (- i 1)) gap-penalty) (+ (matrix:ref C (- j 1) i) gap-penalty)))) (matrix:set! C j i val))))))) ;; là on suppose qu'une séquence est dans un fichier fasta ;; read-fasta lit un fichier fasta et rend la premiere séquence trouvée (define (read-fasta port) (let ((titre (read-line port))) (if (or (eof-object? titre) (zero? (string-length titre)) (not (char=? (string-ref titre 0) #\>))) (error 'read-fasta "not a fasta file" port) (let loop ((str "")) (let ((lu (read-line port))) (if (or (eof-object? lu) (char=? #\> (string-ref lu 0))) (cons titre str) (begin (print lu) (loop (string-append str lu))))))))) (define (f-concorde base1 base2 match-bonus) (if (char=? base1 base2) match-bonus 0)) ;; (define (usage) (print "nw fichier-1 fichier-2 match-bonus gap-penalty") (exit 1)) (define (main argv) (if (not (= (length argv) 5)) (usage) (let ((f1 (cadr argv)) (f2 (caddr argv)) (match-bonus (string->number (cadddr argv))) (gap-penalty (string->number (cadddr (cdr argv))))) (let* ((s1 (make-chain (cdr (let ((port (open-input-file f1))) (read-fasta port))))) (s2 (make-chain (cdr (let ((port (open-input-file f2))) (read-fasta port))))) (the-score-matrix (nw-2 s1 s2 match-bonus gap-penalty))) (matrix:print the-score-matrix) (matrix:print (alignment the-score-matrix match-bonus gap-penalty))))))

Le module de manipulation de matrices

Voici les deux séquences au format FASTA :

> La séquence 1 GAATTCAGTTA > La séquence 2 GGATCGA
Le module de matrices :

(module nw:matrices (export (make-matrix n m . fill) (matrix? obj) (matrix:ref T i j) (matrix:set! T i j val) (matrix:nlines T) (matrix:ncols T) (matrix:margins M s1 s2) (matrix:print T)) (import nw:chains)) ;; il nous faut des matrices (define matrix:tag "*MATRIX*") (define (make-matrix lin col . fill) (let ((the-table (vector "*MATRIX*" (make-vector lin #f)))) (do ((i 0 (+ i 1))) ((= i lin)) (vector-set! (vector-ref the-table 1) i (if (null? fill) (make-vector col) (make-vector col (car fill))))) the-table)) ;; un prédicat d'appartenance, pour vérifier qu'un ;; objet appartient bien au type : (define (matrix? obj) (and (vector? obj) (string=? (vector-ref obj 0) "*MATRIX*") (vector? (vector-ref obj 1)))) ;; un mutateur, pour modifier un objet du type en ;; affectant une valeur à un élément du tableau : (define (matrix:set! T i j val) (if (matrix? T) (vector-set! (vector-ref (vector-ref T 1) i) j val))) ;; un sélecteur, pour accéder à un élément d'un ;; tableau du type : (define (matrix:ref T i j) (if (matrix? T) (vector-ref (vector-ref (vector-ref T 1) i) j))) (define (matrix:margins M s1 s2) (let ((nlin (matrix:nlines M)) (ncol (matrix:ncols M))) (do ((j 2 (+ j 1)) (c (chain-ref s1 1) (chain-ref s1 (min j (- ncol 2))))) ((= j ncol) 'fait) (matrix:set! M 0 j c)) (do ((i 2 (+ i 1)) (c (chain-ref s2 1) (chain-ref s2 (min i (- nlin 2))))) ((= i nlin) 'fait) (matrix:set! M i 0 c)) (matrix:set! M 0 0 #\space) (matrix:set! M 0 1 #\space) (matrix:set! M 1 0 #\space))) ;; diverses procédures utilitaires dont la fonction se ;; comprend d'elle-même : (define (matrix:nlines T) (if (matrix? T) (vector-length (vector-ref T 1)))) (define (matrix:ncols T) (if (matrix? T) (vector-length (vector-ref (vector-ref T 1) 0)))) (define (matrix:print T) (if (matrix? T) (let ((n (matrix:nlines T)) (m (matrix:ncols T))) (do ((i 0 (+ 1 i))) ((= i n)) (let ((this-line (vector-ref (vector-ref T 1) i))) (do ((j 0 (+ 1 j))) ((= j m)) (display (vector-ref this-line j)) (display " ")) (newline))))))

Le module de chaînes :

(module nw:chains (export (make-chain s) (chain-ref s i) (chain-set! s i c) (chain-length s))) ;; il nous faut des chaîne numérotées de 1 à longueur de chaîne: (define (make-chain s) ; prend une string et rend un chaîne à partir de 1 s) (define (chain-ref s i) (string-ref s (- i 1))) (define (chain-set! s i c) (string-set! s (- i 1) c)) (define (chain-length s) (string-length s))

Le Makefile :

BIGLOO = bigloo AFILE = NW-afile.scm BGL_FLAGS = -afile $(AFILE) -Obench -farithmetic -static-bigloo CIBLE = nw REPERT_CIBLE = . %.o: %.scm @ $(BIGLOO) $(BGL_FLAGS) -c $*.scm -o $*.o OBJECTS = NW-lb.o NW-matrices.o NW-chaines.o NW-alignement.o SOURCES = NW-lb.scm NW-matrices.scm NW-chaines.scm NW-alignement.scm all: $(REPERT_CIBLE)/$(CIBLE) $(REPERT_CIBLE)/$(CIBLE): $(OBJECTS) @ echo "Edition de liens..." @ $(BIGLOO) $(BGL_FLAGS) $(OBJECTS) \ -o $(REPERT_CIBLE)/$(CIBLE) @ echo "$(REPERT_CIBLE)/$(CIBLE) construit." @ echo "-------------------------------" $(REPERT_CIBLE): @ mkdir -p $(REPERT_CIBLE) install: cp $(REPERT_CIBLE)/$(CIBLE) $(REPERT_CGI) clean: -rm -f $(OBJECTS) $(SOURCES_C) -rm -f *~ Src/*~ Src/*.o Src/*.mco @ echo "nettoyage fait..." @ echo "-------------------------------" # destruction aussi des binaires : cleanall: clean -rm -f $(REPERT_CIBLE)/$(CIBLE)

1: BLAST (Basic Local Alignment Search Tool) est un logiciel écrit par S.F. Altschul, W. Gish, W. Miller, Gene Myers et D.J. Lipman. Il est l'outil de travail quotidien de tous ceux qui travaillent en biologie moléculaire (cf. http://en.wikipedia.org/wiki/BLAST).

This document was translated from L^AT_EX by H^EV^EA.