Site WWW de Laurent Bloch
Slogan du site

ISSN 2271-3905
Cliquez ici si vous voulez visiter mon autre site, orienté vers des sujets moins techniques.

Pour recevoir (au plus une fois par semaine) les nouveautés de ce site, indiquez ici votre adresse électronique :

Alignement de séquences
Algorithme de Needleman et Wunsch
Article mis en ligne le 27 mai 2008
dernière modification le 24 avril 2022

Cet article a une suite qui, par le retour sur trace (backtracking), mène à l’alignement.

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 (ou similitude, les deux termes sont ici équivalents) 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 [1], qui mettent en œuvre 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(kn) ; 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 nœuds 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 nucléotides, la séquence n° 2 n=7 nucléotides.

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 nucléotide de rang i de la première séquence avec le nucléotide de rang j de le seconde séquence sera Si,j = 1 si les deux nucléotides sont identiques, sinon :

 Si,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 M1,j) et la première colonne (les Mi,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
G 0
A 0
T 0
C 0
G 0
A 0

Remplissage de la matrice

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

Mi,j = maximum de :

 Mi-1,j-1 + Si,j (concordance ou discordance dans la diagonale) ;
 Mi,j-1 + w (gap dans la séquence n° 1) ;
 Mi-1,j + w (gap dans la séquence n° 2).

Ce que l’on peut représenter par un schéma ainsi :

Nous voyons que pour calculer Mi,j il faut (et il suffit de) connaître Mi-1,j, Mi,j-1 et Mi-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 nucléotide G (concordance), S1,1 = 1. Nous avons posé par hypothèse w = 0. Donc :

M1,1 = Max[M0,0 +1, M1,0 +0, M0,1 +0]
= Max[1,0,0]
= 1

Nous pouvons donc inscrire un 1 en M1,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 :

M2,1 = Max[M1,0 +0, M2,0 +0, M1,1 +0]
= Max[0,0,1]
= 1

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(kn)). 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 Mm,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 :

A
¦
A

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


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

  1. (module nw:lb
  2.    (main main)
  3.    (import nw:matrices)
  4.    (import nw:chains)
  5.    (import nw:alignment))
  6.  
  7. (define (nw-2 s1 s2 match-bonus gap-penalty)
  8.   (let ((ncol (+ (chain-length s1) 2))
  9.         (nlin (+ (chain-length s2) 2)))
  10.     (let ((T (make-matrix nlin ncol 0)))
  11.        (matrix-margins T s1 s2) ;;
  12.        (do ((j 2 (+ j 1)))
  13.            ((= j ncol))
  14.            (matrix-set! T 1 j (* (- j 1) gap-penalty)))
  15.        (do ((i 2 (+ i 1)))
  16.            ((= i nlin) T)
  17.            (matrix-set! T i 1 (* (- i 1) gap-penalty))
  18.            (do ((j 2 (+ j 1)))
  19.                ((= j ncol))
  20.                (let ((val
  21.                       (max
  22.                        (+ (matrix-ref T (- i 1) (- j 1))
  23.                           (if (char=? (matrix-ref T i 0)
  24.                                      (matrix-ref T 0 j))
  25.                               match-bonus 0))
  26.                        (+ (matrix-ref T (- i 1) j)
  27.                           gap-penalty)
  28.                        (+ (matrix-ref T i (- j 1))
  29.                           gap-penalty))))
  30.                   (matrix-set! T i j val)))))))
  31.  
  32. ;; on suppose qu'une séquence est dans un fichier fasta
  33. ;; read-fasta lit ce fichier et en rend la première séquence
  34. (define (read-fasta port)
  35.   (let ((titre (read-line port)))
  36.     (if (or (eof-object? titre)
  37.             (zero? (string-length titre))
  38.             (not (char=? (string-ref titre 0)
  39.                          #\>)))
  40.         (error 'read-fasta "not a fasta file" port)
  41.         (let loop ((str ""))
  42.           (let ((lu (read-line port)))
  43.             (if (or (eof-object? lu)
  44.                     (char=? #\> (string-ref lu 0)))
  45.                 (cons titre str)
  46.                 (begin (print lu)
  47.                        (loop (string-append str lu)))))))))
  48.  
  49. (define (usage) ; pour corriger une erreur d'invocation
  50.   (print "nw fichier-1 fichier-2 match-bonus gap-penalty")
  51.   (exit 1))
  52.  
  53. (define (main argv)
  54.   (if (not (= (length argv) 5))
  55.       (usage)
  56.       (let ((f1 (cadr argv))
  57.             (f2 (caddr argv))
  58.             (match-bonus
  59.              (string->number (cadddr argv)))
  60.             (gap-penalty
  61.              (string->number (cadddr (cdr argv)))))
  62.         (let* ((s1 (make-chain
  63.                     (cdr (let ((port (open-input-file f1)))
  64.                             (read-fasta port)))))
  65.                (s2 (make-chain
  66.                     (cdr (let ((port (open-input-file f2)))
  67.                             (read-fasta port)))))
  68.                (the-score-matrix
  69.                 (nw-2 s1 s2 match-bonus gap-penalty)))
  70.            (matrix-print the-score-matrix)
  71.            (matrix-print
  72.             (alignment  the-score-matrix
  73.                         match-bonus gap-penalty))))))

Télécharger

Pour compiler un programme constitué de plusieurs modules, donc répartis dans plusieurs fichiers, Bigloo a besoin d’un fichier access file qui donne la liste des modules et des fichiers dans lesquels ils se trouvent, ainsi :

  1. ((nw:lb "NW-lb.scm")
  2.  (nw:matrices "NW-matrices.scm")
  3.  (nw:chains "NW-chains.scm")
  4.  (nw:alignment "NW-alignment.scm"))

Télécharger

Ce fichier, par défaut, se nommera .afile, si on veut lui donner un autre nom il faut le préciser dans la commande d’invocation du compilateur (drapeau -afile).

Le module de manipulation de matrices

Voici les deux séquences au format FASTA :

Chaque fichier ne peut comporter qu’une séquence, mais le texte de la séquence peut s’étendre sur plusieurs lignes. Vous pouvez visiter au NCBI le site de référence du format FASTA.

Le module de matrices :

  1. (module nw:matrices
  2.    (export
  3.     (make-matrix n m . fill)
  4.     (matrix? obj)
  5.     (matrix-ref T i j)
  6.     (matrix-set! T i j val)
  7.     (matrix-nlines T)
  8.     (matrix-ncols T)
  9.     (matrix-margins M s1 s2)
  10.     (matrix-print T))
  11.    (import nw:chains))
  12.    
  13. (define matrix-tag "*MATRIX*")
  14.  
  15. (define (make-matrix lin col . fill)
  16.    (let ((the-table
  17.           (vector matrix-tag
  18.                   (make-vector lin #f))))
  19.       (do ((i 0 (+ i 1)))
  20.           ((= i lin))
  21.           (vector-set! (vector-ref the-table 1)
  22.                        i
  23.                        (if (null? fill)
  24.                            (make-vector col)
  25.                            (make-vector col (car fill)))))
  26.       the-table))
  27.  
  28. ;; un prédicat d'appartenance, pour vérifier qu'un
  29. ;; objet appartient bien au type :
  30.  
  31. (define (matrix? obj)
  32.    (and (vector? obj)
  33.         (string=? (vector-ref obj 0) matrix-tag)
  34.         (vector?  (vector-ref obj 1))))
  35.  
  36. ;; un mutateur, pour modifier un objet du type en
  37. ;; affectant une valeur à un élément du tableau :
  38.  
  39. (define (matrix-set! T i j val)
  40.    (if (matrix? T)
  41.        (vector-set!
  42.         (vector-ref (vector-ref T 1) i) j val)))
  43.  
  44. ;; un sélecteur, pour accéder à un élément du tableau :
  45.  
  46. (define (matrix-ref T i j)
  47.    (if (matrix? T)
  48.        (vector-ref
  49.           (vector-ref (vector-ref T 1) i) j)))
  50.  
  51. (define (matrix-margins M s1 s2) ;; pour remplir les marges
  52.    (let ((nlin (matrix-nlines M));; du tableau avec les
  53.          (ncol (matrix-ncols M)));; textes des séquences
  54.       (do ((j 2 (+ j 1))        
  55.            (c (chain-ref s1 1)  
  56.               (chain-ref s1
  57.                           (min j (- ncol 2)))))
  58.           ((= j ncol) 'fait)
  59.           (matrix-set! M 0 j c))
  60.       (do ((i 2 (+ i 1))
  61.            (c (chain-ref s2 1)
  62.               (chain-ref s2
  63.                           (min i (- nlin 2)))))
  64.           ((= i nlin) 'fait)
  65.           (matrix-set! M i 0 c))
  66.       (matrix-set! M 0 0 #\space)
  67.       (matrix-set! M 0 1 #\space)
  68.       (matrix-set! M 1 0 #\space)))
  69.  
  70. ;; diverses procédures utilitaires dont la fonction se
  71. ;; comprend d'elle-même :
  72.  
  73. (define (matrix-nlines T)
  74.    (if (matrix? T)
  75.        (vector-length
  76.         (vector-ref T 1))))
  77.  
  78. (define (matrix-ncols T)
  79.    (if (matrix? T)
  80.        (vector-length
  81.         (vector-ref
  82.          (vector-ref T 1) 0))))
  83.  
  84. (define (matrix-print T)
  85.    (if (matrix? T)
  86.        (let ((n (matrix-nlines T))
  87.              (m (matrix-ncols T)))
  88.           (do ((i 0 (+ 1 i)))
  89.               ((= i n))
  90.               (let ((this-line
  91.                      (vector-ref
  92.                       (vector-ref T 1) i)))
  93.                  (do ((j 0 (+ 1 j)))
  94.                      ((= j m))
  95.                      (display
  96.                       (vector-ref this-line j))
  97.                      (display " "))
  98.                  (newline))))))

Télécharger

Le module de chaînes :

  1. (module nw:chains
  2.    (export
  3.     (make-chain s)
  4.     (chain-ref s i)
  5.     (chain-set! s i c)
  6.     (chain-length s)))
  7. ;; il nous faut des chaînes numérotées de 1 à longueur
  8. ;; de chaîne pour éviter trop de gymnastique mentale :
  9. (define (make-chain s) ; prend une string et rend
  10.   s)                   ; un chaîne à partir de 1
  11.  
  12. (define (chain-ref s i)
  13.   (string-ref s (- i 1)))
  14.  
  15. (define (chain-set! s i c)
  16.   (string-set! s (- i 1) c))
  17.  
  18. (define (chain-length s)
  19.   (string-length s))

Télécharger

Le Makefile :