Smith-Waterman

Prérequis

Cette section prolonge la programmation dynamique de 2.2 et l'alignement global de Needleman-Wunsch étudié en 2.3 : vous devez savoir remplir une matrice, distinguer un déplacement diagonal d'un gap et remonter un chemin. Les séquences se lisent comme dans 1.1 et 1.2, et les positions 0-based du code seront opposées aux positions 1-based du texte. Aucune connaissance préalable de Biopython n'est nécessaire : PairwiseAligner est introduit progressivement.

Introduction

Les protéines sont des mosaïques. Un domaine de liaison à l'ADN, un site actif de kinase ou un motif de type zinc finger peut apparaître dans deux protéines qui ne partagent presque rien ailleurs. Une région courte conservée porte souvent l'information fonctionnelle, même quand les extrémités divergent. La question biologique n'est donc pas toujours « ces deux séquences sont-elles globalement apparentées ? » mais « quelle est la meilleure région commune, et où se trouve-t-elle ? ». Un alignement global répond mal à cette seconde question : il force l'appariement de toutes les positions et pénalise les extrémités non homologues.

Le même problème se pose avec les lectures de séquençage. Une lecture de 150 pb peut recouvrir un exon, une jonction d'épissage ou une région intergénique, sans couvrir le gène entier. L'aligner de bout en bout contre un transcrit complet accumule des gaps et des mésappariements dans les parties qui ne sont pas réellement comparables. Ce que l'on veut, c'est le meilleur segment local, celui qui correspond à l'origine biologique de la lecture. Les aligneurs de lectures utilisent des index et des heuristiques, mais le cœur du problème reste le même : trouver une sous-séquence commune de score maximal.

L'algorithme de Smith-Waterman, publié en 1981, résout ce problème par programmation dynamique. Il ajoute au Needleman-Wunsch deux idées simples : une cellule ne descend jamais sous zéro, et le chemin optimal peut commencer n'importe où dans la matrice. Le score final est le maximum de toute la matrice, et le traceback part de ce maximum pour s'arrêter dès qu'il rencontre un zéro. Ces deux modifications transforment un alignement global en un alignement local, avec des conséquences profondes sur l'interprétation.

La difficulté n'est pas seulement algorithmique. Un alignement local est, par construction, le meilleur parmi des millions de candidats : le maximum est donc biaisé vers le haut, et deux séquences sans lien produisent presque toujours un petit alignement positif. Un score brut ne suffit pas à conclure à une homologie.

Théorie

Un alignement qui ne couvre pas tout

Un alignement global de deux séquences $s$ et $t$ de longueurs $n$ et $m$ apparie tous les résidus, quitte à insérer des gaps. Un alignement local apparie deux sous-séquences $s[i..i']$ et $t[j..j']$ et laisse les préfixes et les suffixes non alignés. Formellement, on cherche :

où le score d'un appariement vide vaut zéro. Cette définition autorise un alignement vide, ce qui est cohérent : si aucune région ne score positivement, la meilleure réponse est de ne rien aligner.

Le choix local se justifie lorsque les séquences ont des histoires différentes. Dans ce cas, l'alignement global aligne de force des régions non homologues et dilue le signal du domaine. L'alignement local, lui, peut ignorer les extrémités sans coût. C'est aussi ce qui le rend utile pour les lectures, les fragments de gènes et les séquences chimériques.

Il faut résister à une tentation : croire que l'alignement local est « gratuit ». Il ne l'est pas. Le score maximal est choisi après examen de tous les départs et de toutes les fins possibles, ce qui multiplie les occasions de trouver une ressemblance fortuite.

La récurrence locale

Notons $a_i$ le $i$-ième résidu de $s$ et $b_j$ le $j$-ième résidu de $t$. La matrice $H$ stocke le meilleur score d'un alignement local qui se termine exactement en $a_i$ et $b_j$. Pour $i \ge 1$ et $j \ge 1$ :

avec $H(i,0)=H(0,j)=0$ pour tout $i$ et tout $j$. La pénalité $g$ est négative ou nulle ; dans le code de cette section, $g=-2$. Le terme diagonal apparie $a_i$ et $b_j$ ; le terme $H(i-1,j)+g$ aligne $a_i$ avec un gap dans $t$ ; le terme $H(i,j-1)+g$ aligne un gap dans $s$ avec $b_j$. Le zéro est le point de redémarrage : si les trois possibilités sont négatives, aucun alignement positif ne se termine dans cette cellule.

L'initialisation à zéro est la différence structurelle avec Needleman-Wunsch. Un alignement global initialise la première ligne et la première colonne avec des pénalités cumulées : aligner les $j$ premiers résidus de $t$ contre un préfixe vide coûte $j \times g$. Un alignement local n'a pas à commencer au début : n'importe quel préfixe peut être laissé de côté gratuitement, donc toutes les cellules de bord valent zéro. Mettre des pénalités cumulées sur les bords revient à interdire les départs tardifs et à sous-estimer les domaines internes.

La valeur $H(i,j)$ représente le meilleur alignement local qui finit en $(i,j)$, pas l'alignement des préfixes complets. Un alignement local peut donc commencer après le début des séquences et se terminer avant la fin ; c'est cette liberté qui donne son sens au maximum global de la matrice.

Le coût est $O(nm)$ en temps : chaque cellule combine trois valeurs en temps constant. En mémoire, la matrice complète occupe $O(nm)$ cellules, ce qui est nécessaire pour le traceback. Si l'on ne veut que le score, deux lignes suffisent et la mémoire tombe à $O(\min(n,m))$. Les pénalités affines utilisent trois matrices et restent en $O(nm)$ ; c'est le modèle par défaut de Biopython.

Traceback, bornes et coordonnées

Une fois la matrice remplie, on repère la cellule de score maximal. S'il y a plusieurs maxima, on choisit le premier dans l'ordre de lecture (ligne croissante, puis colonne croissante) ; ce choix doit être documenté, car il détermine l'alignement retourné. Depuis cette cellule, on remonte les déplacements enregistrés pendant le remplissage. Un déplacement diagonal émet la paire $(a_i,b_j)$ ; un déplacement vertical émet $a_i$ contre un gap ; un déplacement horizontal émet un gap contre $b_j$. On s'arrête dès que la cellule courante vaut zéro : le zéro n'est pas aligné, il marque la frontière gauche du segment.

Le point d'arrêt donne les coordonnées de départ. Si le traceback s'arrête à la cellule $(i_0,j_0)$ de valeur zéro, le premier résidu aligné est $a_{i_0+1}$ et $b_{j_0+1}$ en notation 1-based. En Python, les positions renvoyées sont 0-based : start1 = i0 et start2 = j0. La longueur consommée dans chaque séquence est le nombre de résidus non-gap de la chaîne alignée. Les coordonnées 1-based d'un segment de longueur $k$ commençant à l'indice 0-based $i$ sont donc $i+1$ à $i+k$, fin incluse.

L'ordre de priorité en cas d'égalité a des conséquences visibles. Si le score diagonal et le score vertical sont égaux, préférer la diagonale produit un alignement avec moins de gaps ; préférer le gap produit un alignement différent mais de même score. Aucun des deux n'est « plus vrai » : la matrice ne contient pas assez d'information pour trancher. L'essentiel est que le choix soit stable et que le score de l'alignement reconstruit soit égal au score de la cellule de départ.

Exemple travaillé à la main

Prenons $s =$ ACGT et $t =$ TACGTA, avec match $+2$, mismatch $-1$ et gap $-2$. La matrice $H$ complète est :

Vérifions quelques cellules. En $(1,1)$, $A$ contre $T$ est un mismatch : le terme diagonal vaut $-1$, le vertical et l'horizontal valent $-2$, donc $H(1,1)=0$. En $(1,2)$, $A$ contre $A$ : le terme diagonal vaut $0+2=2$, et les deux gaps sont négatifs, donc $H(1,2)=2$. En $(2,3)$, $C$ contre $C$ : le diagonal vaut $H(1,2)+2=4$. En $(3,4)$, $G$ contre $G$ : $H(2,3)+2=6$. Enfin, en $(4,5)$, $T$ contre $T$ : $H(3,4)+2=8$, qui est le maximum de la matrice.

Le traceback part de $(4,5)$ et suit quatre diagonales : $(4,5) \to (3,4) \to (2,3) \to (1,2) \to (0,1)$. La cellule $(0,1)$ vaut zéro, donc le segment commence après elle : start1 = 0 et start2 = 1 en 0-based. En 1-based, le segment couvre les positions 1 à 4 de $s$ et 2 à 5 de $t$. Les deux chaînes alignées sont ACGT et ACGT, pour un score de 8.

Comparons avec un alignement global des mêmes séquences et des mêmes paramètres. Le meilleur alignement global est -ACGT- contre TACGTA : quatre matches valent 8, deux gaps coûtent 4, soit un score de 4. L'alignement local gagne 4 points parce qu'il ne paie pas les deux extrémités non homologues. Ce n'est pas une petite différence de convention : c'est le signal du domaine qui devient visible.

Ce que le score ne dit pas

Le score local est un maximum. Même pour deux séquences aléatoires, il existe presque toujours un petit segment qui score positivement, surtout avec une pénalité de gap faible et un barème favorable aux matches. Le score dépend de la longueur des séquences, de leur composition, de la matrice de substitution et des pénalités. Un score de 35 n'a pas le même sens pour deux protéines de 30 résidus que pour deux protéines de 500 résidus, ni pour deux séquences riches en résidus fréquents.

L'évaluation statistique replace le score dans une distribution de référence. La E-value (expectation value) est le nombre attendu d'alignements de score au moins égal à $S$ dans une recherche de cette taille, sous un modèle nul. Plus elle est petite, plus le score est surprenant. Le bit score normalise le score brut par les paramètres de la matrice et permet de comparer des recherches différentes. Pour deux séquences seulement, on peut approcher la distribution nulle par permutations : on mélange l'ordre des résidus de l'une des séquences, on recalcule le score, et on répète. La proportion de permutations qui atteignent le score observé donne une p-valeur empirique.

Cette p-valeur a ses limites. Les permutations préservent la composition mais détruisent la structure locale, et elles ne remplacent pas une recherche contre une banque, où l'espace de recherche est bien plus grand. Pour les domaines protéiques, les profils HMM de Pfam ou HMMER sont souvent plus sensibles, car ils exploitent la conservation position par position. Smith-Waterman reste la brique de base et un bon vérificateur.

Implémentation

Les blocs suivants ont été exécutés avec Python 3.13.5 et Biopython 1.88 le 2026-11-01. La première fonction n'utilise que la bibliothèque standard. Elle normalise les séquences en majuscules, supprime les espaces, refuse une pénalité de gap positive et remplit une matrice de taille $(n+1) \times (m+1)$. Elle stocke aussi la direction choisie dans chaque cellule, ce qui rend le traceback indépendant d'une seconde comparaison de scores.

Trois choix structurent le code. La diagonale est prioritaire en cas d'égalité, puis le gap dans la seconde séquence, puis dans la première. Le maximum est retenu avec une comparaison stricte : le premier dans l'ordre de lecture gagne. Les positions de départ restent 0-based ; la conversion 1-based est faite à l'affichage.

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

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

8 ACGT ACGT
0 1
[0, 0, 0, 0, 0, 0, 0]
[0, 0, 2, 0, 0, 0, 2]
[0, 0, 0, 4, 2, 0, 0]
[0, 0, 0, 2, 6, 4, 2]
[0, 2, 0, 0, 4, 8, 6]
12 ACGTACGT ACG-ACGT

Le premier appel retrouve le segment ACGT et les positions attendues. Le second montre un gap interne : la séquence ACGACGT est alignée sur ACGTACGT avec un - après le G, et le score vaut 12. Le gap coûte 2, mais les sept matches rapportent 14. Le score de l'alignement reconstruit est bien égal à la valeur de la cellule de départ ; c'est cette cohérence que les tests vérifient.

Biopython expose le même algorithme dans Bio.Align.PairwiseAligner. En mode local, il suffit de régler le mode et les scores. Le bloc suivant est autonome et compare les deux exemples précédents. Le paramétrage open_gap_score = extend_gap_score = -2 reproduit une pénalité linéaire de 2 par résidu de gap ; avec des valeurs différentes, Biopython applique une pénalité affine.

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

La sortie vérifiée est :

ACGT TACGTA 8.0
[[0, 4], [1, 5]]
ACGTACGT ACGACGT 12.0
[[0, 3, 4, 8], [0, 3, 3, 7]]
Smith-Waterman

Les scores coïncident avec ceux de notre fonction. Les coordonnées de Biopython sont 0-based et décrivent les blocs alignés : pour le second exemple, [0, 3, 4, 8] côté cible et [0, 3, 3, 7] côté requête signifient un premier bloc de trois résidus, puis un gap d'un résidu dans la requête, puis un bloc de quatre résidus. aligner.align peut renvoyer plusieurs alignements optimaux ; [0] choisit le premier, et ce choix peut différer du nôtre en cas d'égalité. C'est une raison de tester le score et les invariants plutôt que de comparer des chaînes au caractère près avec une autre bibliothèque.

Étude de cas

Nous travaillons sur deux protéines courtes synthétiques, construites pour la démonstration. La première, de 33 résidus, contient le segment AYIAKQRQISFVKSHFSRQ ; la seconde, de 30 résidus, contient le même segment avec une substitution S en T à la position 17 du domaine. Les régions flanquantes sont volontairement différentes. Il ne s'agit pas d'un jeu public et les séquences ne correspondent à aucune protéine réelle : l'exemple est étiqueté comme synthétique.

seq1  MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQ
seq2  GGSAYIAKQRQISFVKSHFTRQDILDLWIY

La question est de localiser le domaine conservé, de vérifier que l'alignement local ignore les flancs et de mesurer si le score observé est rare par rapport à un modèle nul. Nous utilisons match $+2$, mismatch $-1$ et gap $-2$, les paramètres de la section.

Le script suivant utilise PairwiseAligner pour que le bloc soit autonome ; notre implémentation, appelée avec les mêmes paramètres, renvoie le même score et les mêmes bornes. Il calcule l'alignement local, puis un alignement global avec les mêmes scores, puis une distribution nulle par 500 permutations de la seconde séquence avec une graine fixée.

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

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

target            3 AYIAKQRQISFVKSHFSRQ 22
                  0 ||||||||||||||||.|| 19
query             3 AYIAKQRQISFVKSHFTRQ 22
 
score local 35.0
coords [[3, 22], [3, 22]]
score global 24.0
moyenne nulle 5.134
maximum nul 10.0
p empirique 0.002

L'alignement local couvre exactement les positions 3 à 21 en 0-based, soit 4 à 22 en 1-based, dans les deux séquences. Les 19 colonnes contiennent 18 identités et une substitution ; le score vaut $18 \times 2 - 1 = 35$. Le domaine est retrouvé sans les flancs, alors que l'alignement global, qui doit absorber les extrémités et un décalage de longueur, tombe à 24. L'écart de 11 points mesure le coût des régions non homologues.

La distribution nulle est nettement séparée : moyenne 5,134, maximum 10 sur 500 permutations, aucun score n'atteint 35. La p-valeur empirique vaut $(0+1)/(500+1) \approx 0,002$. Cette valeur est la plus petite que ce test puisse produire avec 500 permutations ; elle signifie « aucun des mélanges n'a fait aussi bien », pas « probabilité exacte de 0,002 ». Elle dépend de la composition de seq2, de la longueur des séquences et du barème. Une protéine courte a moins de puissance statistique qu'une protéine longue, et un domaine de 19 résidus est un échantillon modeste.

L'interprétation doit rester prudente. Le domaine a été inséré volontairement : l'expérience vérifie que l'algorithme le retrouve, pas qu'il découvre une homologie. Dans une vraie recherche, on interrogerait une banque avec BLASTP ou un profil HMM, en regardant la E-value, la couverture et la conservation des résidus. Un score de 35 contre une moyenne de 5 est un signal fort dans ce cadre synthétique ; il ne prouve pas qu'un domaine réel de même score soit biologiquement pertinent.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-2/2.4-smith-waterman -q

python -m pytest exercises/module-2/2.4-smith-waterman -q

python -m pytest exercises/module-2/2.4-smith-waterman/projet -q

Résumé

• Smith-Waterman cherche le meilleur alignement de deux sous-séquences, pas l'alignement de deux séquences entières.
• La récurrence combine un plancher à zéro, le diagonal avec $S(a_i,b_j)$ et deux gaps de pénalité $g$ ; les bords valent zéro.
• Le traceback part de la cellule maximale et s'arrête au premier zéro ; les positions de départ sont 0-based dans le code.
• Needleman-Wunsch initialise les bords avec des pénalités cumulées et remonte depuis le coin inférieur droit ; Smith-Waterman fait l'inverse.
• Le coût est $O(nm)$ en temps et en mémoire pour la matrice complète ; le score seul peut tenir en $O(\min(n,m))$.
• PairwiseAligner en mode local reproduit le score et fournit des coordonnées 0-based, mais peut choisir un autre alignement optimal.
• Un score local est un maximum biaisé : E-value, bit score ou permutations sont nécessaires pour écarter les faux positifs.
• Les domaines protéiques et les lectures fragmentaires justifient l'alignement local ; pour l'homologie distante, les profils HMM sont souvent plus sensibles.

Références

• Smith, T. F., & Waterman, M. S. (1981). Identification of common molecular subsequences. Journal of Molecular Biology, 147(1), 195–197.
• Needleman, S. B., & Wunsch, C. D. (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of Molecular Biology, 48(3), 443–453.
• Altschul, S. F., Gish, W., Miller, W., Myers, E. W., & Lipman, D. J. (1990). Basic local alignment search tool. Journal of Molecular Biology, 215(3), 403–410.
• Henikoff, S., & Henikoff, J. G. (1992). Amino acid substitution matrices from protein blocks. Proceedings of the National Academy of Sciences, 89(22), 10915–10919.
• Cock, P. J. A., Antao, T., Chang, J. T., Chapman, B. A., Cox, C. J., Dalke, 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 officielle : https://biopython.org/docs/latest/api/Bio.Align.html
• Biopython. Aligning sequences with Biopython. Documentation officielle : https://biopython.org/docs/latest/Tutorial/chapter_pairwise.html