Needleman-Wunsch

Prérequis

Cette section suppose connues les notions de séquence, d'alphabet et de complément inverse de la section 1.1, ainsi que la lecture d'un fichier FASTA vue en 1.2. Elle prolonge la section 2.2 : un alignement n'est pas une superposition de deux chaînes, c'est une correspondance entre positions, avec des insertions et des délétions possibles. Le module 2 situe Needleman-Wunsch parmi les autres algorithmes de comparaison. Aucune connaissance préalable de Biopython n'est nécessaire pour la première partie.

Introduction

Comparer deux séquences est l'une des opérations les plus fréquentes en bioinformatique. Deux gènes orthologues de deux bactéries proches ne diffèrent souvent que par quelques substitutions ; une protéine et sa paralogue peuvent avoir conservé des domaines entiers tout en ayant perdu ou gagné des segments. Pour décrire ces relations, il faut décider quelles positions se correspondent et quelles positions n'ont pas d'équivalent. Un alignement est précisément cette correspondance, enrichie d'un score qui agrège les appariements, les substitutions et les gaps.

Le problème computationnel est double. D'abord, le nombre d'alignements possibles entre deux séquences croît de façon exponentielle avec leur longueur : une recherche exhaustive est hors de portée dès quelques dizaines de résidus. Ensuite, la notion même de « bon » alignement dépend d'un modèle de score : identités, substitutions, pénalités d'insertion et de délétion. Needleman et Wunsch ont proposé en 1970 une méthode de programmation dynamique qui résout le premier problème de manière exacte pour un modèle additif donné. Leur algorithme reste la référence de l'alignement global, c'est-à-dire d'un alignement qui couvre les deux séquences de bout en bout.

L'enjeu biologique est de distinguer une ressemblance due à l'évolution d'une ressemblance fortuite. L'enjeu computationnel est de produire un score reproductible et un alignement inspectable, puis de savoir ce que ce score signifie. Nous construisons ici l'algorithme pas à pas, nous le vérifions sur un exemple court, nous le comparons à l'implémentation de Biopython et nous délimitons ce qu'il ne peut pas dire. La section 2.4 montrera l'alignement local, mieux adapté à la recherche de domaines isolés.

Théorie

Alignement global : définitions et intuition

Un alignement global (global alignment) de deux séquences $s = s_1 s_2 \dots s_n$ et $t = t_1 t_2 \dots t_m$ est une suite de colonnes, chaque colonne contenant soit un caractère de $s$ et un caractère de $t$, soit un caractère et un gap, soit un gap et un caractère. Aucune colonne ne contient deux gaps, et les caractères de chaque séquence apparaissent dans l'ordre. Les deux séquences alignées ont donc la même longueur, et l'on note $A$ l'alignement complet.

Un score additif attribue une valeur à chaque colonne. Le score de l'alignement est la somme des scores de ses colonnes :

où $\ell$ est le nombre de colonnes, $S(a,b)$ le score d'un appariement ou d'une substitution, et où chaque gap est compté $g$. Dans le modèle le plus simple, $S(a,b) = \mathrm{match}$ si $a=b$, et $S(a,b) = \mathrm{mismatch}$ sinon. Les matrices de substitution plus réalistes, comme BLOSUM ou PAM, remplacent ce modèle en 2.5 ; la récurrence reste la même.

L'idée de la programmation dynamique est de ne pas énumérer les alignements. On note $F(i,j)$ le score optimal de l'alignement global des préfixes $s_1 \dots s_i$ et $t_1 \dots t_j$. La propriété de sous-structure optimale affirme qu'un alignement optimal de ces préfixes se termine nécessairement par l'une de trois colonnes : $s_i$ apparié à $t_j$, $s_i$ apparié à un gap, ou un gap apparié à $t_j$. Sinon, on pourrait retirer la dernière colonne et améliorer le reste, ce qui contredirait l'optimalité. C'est cette observation qui donne la récurrence.

Initialisation

La première ligne et la première colonne représentent l'alignement d'un préfixe avec une séquence vide. Aligner $i$ caractères avec rien coûte $i$ gaps, et aligner $j$ caractères avec rien coûte $j$ gaps :

Avec $g=-1$, la première colonne vaut $0,-1,-2,\dots$ et la première ligne $0,-1,-2,\dots$. Cette initialisation est la marque de l'alignement global. Dans un alignement local, on initialiserait au contraire toute la première ligne et toute la première colonne à zéro, car on autorise l'alignement à commencer n'importe où. Confondre les deux revient à changer d'algorithme.

Récurrence

Pour $i \ge 1$ et $j \ge 1$, le score optimal vérifie :

Les trois termes correspondent aux trois dernières colonnes possibles. Le premier prolonge un alignement des préfixes $i-1$ et $j-1$ par l'appariement $s_i$ avec $t_j$. Le deuxième prolonge un alignement des préfixes $i-1$ et $j$ par $s_i$ face à un gap : on avance dans $s$ sans avancer dans $t$. Le troisième prolonge un alignement des préfixes $i$ et $j-1$ par un gap face à $t_j$ : on avance dans $t$ sans avancer dans $s$. Le maximum choisit la meilleure des trois possibilités, et la matrice se remplit ligne par ligne, de gauche à droite et de haut en bas.

La valeur $F(n,m)$ est le score optimal de l'alignement global. Elle est unique, car le maximum d'un ensemble fini est unique. En revanche, plusieurs chemins peuvent atteindre ce maximum : le score ne détermine pas l'alignement. C'est une distinction centrale, que nous retrouverons au moment du traceback.

Traceback : reconstruire l'alignement

La matrice seule donne le score. Pour obtenir un alignement, on part de la cellule $(n,m)$ et on remonte vers $(0,0)$ en suivant un prédécesseur qui a produit le maximum. À chaque pas, on ajoute une colonne à l'alignement, puis on inverse l'ordre à la fin. Trois mouvements sont possibles :

• diagonale : on ajoute $s_i$ et $t_j$, puis $i \leftarrow i-1$ et $j \leftarrow j-1$ ;
• vertical : on ajoute $s_i$ et un gap, puis $i \leftarrow i-1$ ;
• horizontal : on ajoute un gap et $t_j$, puis $j \leftarrow j-1$.

Une erreur classique consiste à tester l'égalité avec le mauvais terme ou à oublier le score de substitution dans la comparaison. Le traceback doit utiliser exactement les mêmes paramètres que le remplissage. Il faut aussi gérer les bords : lorsque $i=0$, seuls des mouvements horizontaux sont possibles ; lorsque $j=0$, seuls des mouvements verticaux le sont.

Départager les alignements optimaux

Lorsque plusieurs prédécesseurs atteignent le maximum, il faut choisir. Ce choix s'appelle le départage (tie-breaking). Il n'est pas anodin : il détermine l'alignement affiché, alors que le score reste identique. Une règle explicite et stable est indispensable pour que les tests et les comparaisons soient reproductibles.

Nous adoptons la règle suivante, également utilisée dans la section 2.4 : à score égal, on préfère la diagonale, puis le mouvement vertical (gap dans la seconde séquence), puis le mouvement horizontal (gap dans la première séquence). Cette règle est simple à implémenter et à tester. D'autres règles sont légitimes, par exemple préférer un appariement identique avant tout autre mouvement, ou préférer les gaps en fin de séquence. Ce qui compte est de la documenter et de ne pas la changer au milieu d'une analyse.

Sur l'exemple ci-dessous, trois alignements atteignent le score 0. Notre règle en sélectionne un ; Biopython, avec sa propre règle, en affiche un autre. Les deux sont optimaux. Un désaccord d'affichage n'est donc pas un bug tant que les scores coïncident et que chaque alignement se recalcule au bon score.

Score optimal, alignement optimal

Le score optimal $F(n,m)$ est une propriété de la matrice et des paramètres. Un alignement optimal est un chemin qui réalise ce score. Il peut y en avoir plusieurs, et leur nombre peut croître rapidement avec la longueur des séquences. Un programme qui ne renvoie que le score jette l'information de chemin ; un programme qui renvoie un alignement sans son score oblige l'utilisateur à recalculer. Les deux doivent être fournis et cohérents.

Il faut résister à deux tentations. La première est de croire qu'un score élevé identifie à lui seul une homologie : le score dépend de la longueur, de la composition et des paramètres. La seconde est de comparer des scores obtenus avec des pénalités de gap différentes. Un score de 120 avec $g=-1$ et un score de 90 avec $g=-4$ ne sont pas comparables sans conversion ni modèle statistique.

Complexité et mémoire

Le remplissage visite $(n+1)(m+1)$ cellules et effectue un nombre constant d'opérations par cellule. La complexité en temps est donc $\Theta(nm)$. La matrice complète occupe $\Theta(nm)$ cellules, soit $\Theta(nm)$ mémoire. Pour deux séquences de 10 000 caractères, cela représente environ $10^8$ cellules : hors de portée d'une implémentation Python naïve, même si la formule reste correcte.

Si seul le score importe, on peut ne conserver que deux lignes, ce qui ramène la mémoire à $O(m)$ et le temps reste $\Theta(nm)$. Pour reconstruire un alignement avec une mémoire linéaire, on utilise l'algorithme de Hirschberg, de même complexité temporelle. Les implémentations optimisées, comme celle de Biopython, sont écrites en C et restent quadratiques : l'optimisation change la constante, pas l'ordre de grandeur. Pour des séquences très proches, un alignement bandé restreint la zone explorée autour de la diagonale et peut réduire fortement le coût.

Exemple travaillé à la main

Prenons $s =$ GATTACA et $t =$ GCATGCU, avec $\mathrm{match}=+1$, $\mathrm{mismatch}=-1$ et $g=-1$. Les longueurs valent $n=7$ et $m=7$. La matrice complète est donnée dans le tableau 2.3.1 ; les lignes correspondent aux préfixes de $s$, les colonnes aux préfixes de $t$.

Vérifions quelques cellules. $F(1,1)=1$ car G et G s'apparient. $F(1,2)=0$ : le meilleur choix est de laisser C face à un gap, ce qui donne $1-1=0$. $F(3,4)=2$ : les préfixes GAT et GCAT s'alignent par G/G, gap/C, A/A et T/T, soit $1-1+1+1=2$. La cellule finale vaut $F(7,7)=0$.

Le traceback avec notre règle donne :

G-ATTACA
GCA-TGCU
| | | |

Recomptons le score colonne par colonne : G/G $+1$, gap/C $-1$, A/A $+1$, T/gap $-1$, T/T $+1$, A/G $-1$, C/C $+1$, A/U $-1$. La somme vaut $0$, ce qui confirme la cellule finale. L'alignement contient quatre appariements identiques, deux substitutions et deux gaps.

Deux autres alignements ont aussi le score 0 : G-ATTACA avec GCATG-CU, et G-ATTACA avec GCAT-GCU. Ils déplacent simplement un gap. Le nombre total de chemins optimaux vaut 3 ; la fonction count_optimal_alignments des exercices le vérifie sans énumérer les alignements. Cette multiplicité est normale et doit être documentée, pas masquée.

Implémentation

Le code ci-dessous est une version compacte et exécutable de l'algorithme. Il normalise les séquences en majuscules, remplit la matrice, puis remonte le chemin. Les indices Python sont 0-based, alors que les formules utilisent des positions 1-based : seq1[i-1] désigne le $i$-ième caractère.

Aucune sortie. Exécutez le code pour voir le résultat.

La sortie vérifiée le 2026-11-01 avec Python 3.13.5 est :

         -   G   C   A   T   G   C   U
   -     0  -1  -2  -3  -4  -5  -6  -7
   G    -1   1   0  -1  -2  -3  -4  -5
   A    -2   0   0   1   0  -1  -2  -3
   T    -3  -1  -1   0   2   1   0  -1
   T    -4  -2  -2  -1   1   1   0  -1
   A    -5  -3  -3  -1   0   0   0  -1
   C    -6  -4  -2  -2  -1  -1   1   0
   A    -7  -5  -3  -1  -2  -2   0   0
 
Score : 0
G-ATTACA
GCA-TGCU
| | | |

Trois points méritent attention. D'abord, l'initialisation doit être faite avant la double boucle ; si l'on oublie les bords, les premières cellules utilisent des zéros implicites et le score est faux. Ensuite, le traceback doit tester la diagonale avant les gaps pour respecter la règle de départage annoncée. Enfin, la comparaison F[i][j] == F[i-1][j-1] + s doit employer le même match et le même mismatch que le remplissage ; un décalage d'indice ou un score différent produit un chemin incohérent, voire une boucle infinie.

La version des exercices ajoute la normalisation des espaces, la vérification des types, le refus d'une pénalité de gap positive et une fonction alignment_score qui recalcule le score d'un alignement. Cette dernière fonction est un garde-fou précieux : elle permet de vérifier qu'un alignement renvoyé par un programme, le nôtre ou celui d'un autre, réalise bien le score annoncé.

Comparaison avec Bio.Align.PairwiseAligner

Biopython expose un aligneur global et local dans Bio.Align. Les paramètres doivent être posés explicitement : le mode, le score d'appariement, le score de substitution et les pénalités de gap. Le code suivant compare les deux implémentations sur l'exemple.

from Bio import Align
 
seq1, seq2 = "GATTACA", "GCATGCU"
aligner = Align.PairwiseAligner()
aligner.mode = "global"
aligner.match_score = 1
aligner.mismatch_score = -1
aligner.open_gap_score = -1
aligner.extend_gap_score = -1
 
alignments = aligner.align(seq1, seq2)
print("Score Biopython :", alignments.score)
print("Nombre d'alignements optimaux :", len(alignments))
print(alignments[0])

La sortie vérifiée est :

Score Biopython : 0.0
Nombre d'alignements optimaux : 3
target            0 G-ATTACA 7
                  0 |-||.-|. 8
query             0 GCATG-CU 7

Le score concorde avec le nôtre. L'alignement affiché diffère : Biopython place le gap après GCATG, nous le plaçons après GCA. Les deux ont le score 0. Le nombre d'alignements optimaux, 3, est une information que notre fonction principale ne renvoie pas, mais que count_optimal_alignments calcule.

La sémantique des gaps de Biopython mérite une attention particulière. Avec open_gap_score = -1 et extend_gap_score = -1, un gap de longueur $k$ coûte $-1 + (k-1)(-1) = -k$, ce qui correspond exactement à notre pénalité linéaire $g=-1$ par base. Si l'on règle open_gap_score = -2 et extend_gap_score = -0.5, on obtient un modèle affine, différent de celui de cette section. Comparer deux scores sans vérifier ces réglages n'a aucun sens.

Étude de cas

Exemple synthétique. Les deux fragments ci-dessous ne proviennent pas d'un génome public : ils ont été construits pour isoler deux événements simples, une substitution et une délétion. Les séquences font 21 et 20 nucléotides et ne prétendent pas représenter un gène réel. L'étiquette « exemple synthétique » fait partie du résultat : elle interdit de tirer une conclusion biologique de ce qui suit.

Nous comparons deux paires. La première oppose ATGGCTGAACGTAAACCGTAA à ATGGCCGAACGTAAACCGTAA : une seule position diffère, T contre C. La seconde oppose la même référence à ATGGCTGAAGTAAACCGTAA, qui a perdu un C interne. La question est double : notre implémentation retrouve-t-elle le score de Biopython, et le score suffit-il à distinguer une substitution d'une délétion ?

Le script suivant importe le corrigé des exercices, construit un aligneur Biopython avec les mêmes paramètres, puis affiche les deux scores, les deux alignements et leur écart. Il s'exécute depuis la racine du dépôt.

from sys import path
 
path.insert(0, "exercises/module-2/2.3-needleman-wunsch")
 
from Bio import Align
 
from solution_2_3 import alignment_score, needleman_wunsch
 
 
def biopython(seq1, seq2):
    aligner = Align.PairwiseAligner()
    aligner.mode = "global"
    aligner.match_score = 1
    aligner.mismatch_score = -1
    aligner.open_gap_score = -1
    aligner.extend_gap_score = -1
    alignments = aligner.align(seq1, seq2)
    best = alignments[0]
    return float(alignments.score), (str(best[0]), str(best[1])), len(alignments)
 
 
pairs = {
    "substitution": ("ATGGCTGAACGTAAACCGTAA", "ATGGCCGAACGTAAACCGTAA"),
    "deletion": ("ATGGCTGAACGTAAACCGTAA", "ATGGCTGAAGTAAACCGTAA"),
}
 
for label, (a, b) in pairs.items():
    ours = needleman_wunsch(a, b)
    bio_score, bio_alignment, n = biopython(a, b)
    print(label)
    print("  NW     :", ours["score"], ours["aligned1"], "/", ours["aligned2"])
    print("  Bio    :", bio_score, bio_alignment[0], "/", bio_alignment[1])
    print("  difference :", ours["score"] - bio_score)
    print("  chemins:", n)
    print("  score de l'alignement NW :", alignment_score(ours["aligned1"], ours["aligned2"]))

La sortie vérifiée le 2026-11-01 est :

substitution
  NW     : 19 ATGGCTGAACGTAAACCGTAA / ATGGCCGAACGTAAACCGTAA
  Bio    : 19.0 ATGGCTGAACGTAAACCGTAA / ATGGCCGAACGTAAACCGTAA
  difference : 0.0
  chemins: 1
  score de l'alignement NW : 19
deletion
  NW     : 19 ATGGCTGAACGTAAACCGTAA / ATGGCTGAA-GTAAACCGTAA
  Bio    : 19.0 ATGGCTGAACGTAAACCGTAA / ATGGCTGAA-GTAAACCGTAA
  difference : 0.0
  chemins: 1
  score de l'alignement NW : 19

Les deux paires donnent un écart nul : notre implémentation et Biopython sont d'accord sur le score et sur l'alignement, et le score recalculé par alignment_score confirme la valeur 19. Pour la substitution, l'alignement contient 20 appariements identiques et une substitution, soit $20-1=19$. Pour la délétion, il contient 19 appariements identiques et un gap, soit $19-1=19$. Le même score agrège donc deux histoires différentes.

C'est la leçon principale de l'étude de cas : un score d'alignement global n'est pas un diagnostic. Il ne dit pas si un nucléotide a été remplacé ou perdu ; il dit seulement que le modèle de score préfère cette explication à toutes les autres. L'alignement, lui, montre l'événement. Pour décider si cette ressemblance est biologiquement significative, il faudrait un modèle nul, des séquences de contrôle, une matrice adaptée et, pour une recherche dans une base, une E-value. Aucun de ces éléments n'est fourni par Needleman-Wunsch.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-2/2.3-needleman-wunsch -q

python -m pytest exercises/module-2/2.3-needleman-wunsch -q

python -m pytest exercises/module-2/2.3-needleman-wunsch/projet -q

Résumé

• Needleman-Wunsch calcule un alignement global exact pour un score additif donné.
• L'initialisation globale vaut $F(i,0)=ig$ et $F(0,j)=jg$ ; des bords à zéro changent d'algorithme.
• La récurrence compare trois prédécesseurs : diagonale, gap dans la seconde séquence, gap dans la première.
• Le score optimal est unique, mais plusieurs alignements peuvent l'atteindre ; le départage doit être documenté.
• Le traceback reconstruit un chemin, puis il faut inverser l'ordre des colonnes.
• Le temps et la mémoire sont en $\Theta(nm)$ pour la matrice complète ; des variantes réduisent la mémoire.
• Biopython reproduit le score si le mode et les pénalités de gap sont réglés explicitement.
• Un score brut n'est pas une preuve d'homologie : il manque un modèle statistique et un contrôle.

Références

• Needleman, S. B. et Wunsch, C. D. (1970). A general method applicable to the search for similarities in the amino acid sequences of two proteins. Journal of Molecular Biology, 48(3), 443–453.
• Durbin, R., Eddy, S. R., Krogh, A. et Mitchison, G. (1998). Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press.
• Eddy, S. R. (2004). What is dynamic programming? Nature Biotechnology, 22(7), 909–910.
• Hirschberg, D. S. (1975). A linear space algorithm for computing maximal common subsequences. Communications of the ACM, 18(6), 341–343.
• Smith, T. F. et Waterman, M. S. (1981). Identification of common molecular subsequences. Journal of Molecular Biology, 147(1), 195–197.
• Cock, P. J. A. et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11), 1422–1423.
• Biopython. Bio.Align.PairwiseAligner documentation, version 1.88. https://biopython.org/docs/1.88/api/Bio.Align.html