Programmation dynamique

Prérequis

Cette section prolonge la comparaison de séquences abordée en 2.1 et la lecture de séquences vue en 1.1. Vous devez savoir extraire une séquence d'un fichier FASTA concaténé, comme dans 1.2, et distinguer un motif exact d'un alignement avec erreurs, comme dans 1.4. Les formules utilisent des indices 1-based, alors que le code Python travaille en 0-based ; chaque décalage est explicité. Aucune bibliothèque d'alignement n'est nécessaire pour comprendre la matrice.

Introduction

Deux séquences apparentées ne sont pas identiques : au fil des générations, des substitutions remplacent des bases, des insertions ajoutent des caractères et des délétions en retirent. Retrouver l'histoire de ces événements derrière deux séquences actuelles est l'un des problèmes centraux de la bioinformatique. Un alignement propose une correspondance position par position : les résidus alignés sont supposés homologues, les gaps représentent des événements d'insertion ou de délétion. Cette correspondance sert ensuite à inférer une fonction, à construire un arbre phylogénétique, à dessiner des amorces ou à interpréter un variant.

Le problème computationnel est redoutable si on l'attaque frontalement. Pour deux séquences de longueur 4 et 5, il existe déjà 681 alignements possibles ; pour deux séquences de longueur 10, plus de 8 millions ; pour deux séquences de longueur 20, environ $2,6 \times 10^{14}$. Énumérer ces candidats et les scorer un par un devient impossible dès que les séquences dépassent quelques dizaines de résidus, et absurde pour un gène de plusieurs kilobases.

La programmation dynamique (dynamic programming) transforme cette explosion en un calcul polynomial. Au lieu de considérer chaque alignement comme un objet global, elle construit la solution à partir de solutions de préfixes plus courts. Le nombre de sous-problèmes distincts est petit, même si le nombre de chemins qui les traversent est immense. Cette idée, formulée par Bellman dans les années 1950, s'applique à l'alignement global, à l'alignement local, aux gaps affines et à bien d'autres problèmes de séquence.

Théorie

Un nombre d'alignements qui explose

Un alignement global de $s$ (longueur $n$) et $t$ (longueur $m$) est un chemin dans une grille de $(n+1)(m+1)$ nœuds. Chaque pas diagonal aligne deux résidus, chaque pas vertical aligne un résidu de $s$ avec un gap, chaque pas horizontal fait l'inverse. Le nombre de chemins est le nombre de Delannoy $D(n,m)$ : $D(4,5)=681$, $D(10,10)=8\,097\,453$ et $D(20,20)=260\,543\,813\,797\,441$. L'énumération exhaustive est donc hors de portée dès quelques dizaines de résidus. La programmation dynamique ne parcourt pas les chemins : elle calcule une valeur par nœud.

Sous-structure optimale

La sous-structure optimale (optimal substructure) est la propriété qui rend ce calcul possible. Si un alignement optimal de $s$ et $t$ aligne $s_n$ avec $t_m$, alors l'alignement des préfixes $s_1 \dots s_{n-1}$ et $t_1 \dots t_{m-1}$ doit être optimal : sinon, remplacer ce préfixe par un meilleur alignement améliorerait le score global, contradiction. Le même argument vaut si le dernier résidu est aligné à un gap. Un alignement optimal se décompose donc en un sous-alignement optimal et une dernière colonne.

Cette propriété dépend du barème. Avec des pénalités affines, la dernière colonne ne suffit plus : il faut savoir si un gap est déjà ouvert, d'où des états supplémentaires.

Sous-problèmes chevauchants

Le calcul naïf de $F(i,j)$ appelle $F(i-1,j-1)$, $F(i-1,j)$ et $F(i,j-1)$, qui rappellent des préfixes déjà visités. L'arbre d'appels est exponentiel, mais il n'existe que $(n+1)(m+1)$ sous-problèmes distincts. C'est le chevauchement des sous-problèmes (overlapping subproblems). Mémoriser chaque résultat fait passer le temps à $O(nm)$.

Mémorisation ou tabulation

La mémorisation (memoization) est descendante : une fonction récursive met ses résultats en cache, par exemple avec functools.lru_cache. Elle est lisible et ne calcule que les états atteints, mais paie les appels récursifs et peut dépasser la limite de récursion.

La tabulation est ascendante : on remplit les bords, puis les cellules par lignes et colonnes croissantes. Chaque cellule trouve ses trois dépendances déjà calculées. C'est la version habituelle pour l'alignement, car elle donne accès à toute la matrice, donc au traceback. Les deux approches doivent produire le même score ; sinon, un bord ou un indice est faux.

La récurrence générale

Soit $F(i,j)$ le score du meilleur alignement global des préfixes $s_1 \dots s_i$ et $t_1 \dots t_j$, et $S(a,b)$ le score de substitution. Alors

avec les bords

La réponse est $F(n,m)$. Les trois candidats sont la diagonale (deux résidus alignés), le gap dans $t$ et le gap dans $s$. La diagonale ne signifie pas « identité » : c'est $S(s_i,t_j)$ qui distingue match et mismatch. Un gap n'est pas un caractère biologique, mais une colonne vide introduite pour aligner la suite.

La matrice de score

La matrice de score donne $S(a,b)$. Ici : match $+1$, mismatch $-1$, gap $-1$. Elle est symétrique. Les matrices BLOSUM et PAM, vues plus tard, encodent des préférences observées. Un score n'est ni une probabilité ni une grandeur universelle : il dépend du barème et croît avec la longueur. Un score brut n'est donc comparable qu'avec un barème et une longueur comparables.

Pénalités de gap : linéaire et affine

Une pénalité linéaire facture $k \cdot g$ pour un gap de longueur $k$ : un gap de trois bases coûte trois fois un gap d'une base, ce qui traite un indel unique comme trois événements.

Une pénalité affine sépare ouverture et extension :

Avec $g_o=-5$ et $g_e=-1$, un gap de 1, 2 et 3 bases coûte $-6$, $-7$ et $-8$. Ouvrir est cher, allonger l'est moins : c'est le comportement attendu d'un indel unique. Certaines bibliothèques facturent $g_o+(k-1)g_e$ ; un gap d'une base coûte alors $-5$, et la différence vaut $g_e$ par gap. Cette convention doit être écrite dans le code et la documentation. Avec l'affine, la récurrence linéaire ne suffit plus : il faut suivre l'état « gap ouvert », donc trois matrices couplées. Le temps reste $O(nm)$, mais la constante augmente.

Conditions aux bords

En global, les préfixes vides paient tous leurs gaps : la première colonne vaut $0,g,2g,\dots,ng$ et la première ligne $0,g,2g,\dots,mg$. Seule $(0,0)$ vaut zéro. En local, toute la première ligne et toute la première colonne valent zéro, car on peut commencer n'importe où. Initialiser les bords à zéro en global surestime le score et ignore les extrémités.

La matrice est de dimensions $(n+1) \times (m+1)$ : une matrice $5 \times 6$ correspond à des séquences de longueur 4 et 5. Dans le code, $F(i,j)$ est matrix[i][j] et $s_i$ est seq[i-1]. Ce décalage est une source classique d'erreurs.

Traceback

Le score est en $(n,m)$ ; pour retrouver l'alignement, on remonte vers $(0,0)$ en suivant la décision qui a produit le maximum. Un pas diagonal aligne $s_i$ avec $t_j$, un pas vertical aligne $s_i$ avec un gap, un pas horizontal aligne un gap avec $t_j$. On teste la diagonale, puis le vertical, puis l'horizontal : cet ordre départage les égalités. Le score est unique, l'alignement ne l'est pas forcément. Les colonnes sont accumulées à l'envers et doivent être inversées.

Complexité et réduction mémoire

Le remplissage visite $(n+1)(m+1)$ cellules en temps constant : $O(nm)$ temps et $O(nm)$ mémoire pour la matrice complète. Pour le score seul, deux lignes suffisent, donc $O(m)$ mémoire : chaque cellule a besoin de la diagonale et de la verticale de la ligne précédente, et de l'horizontale de la ligne courante. Mais le traceback devient impossible sans la matrice.

L'algorithme de Hirschberg divise le problème, calcule les scores direct et inverse avec deux lignes, puis récurse : $O(nm)$ temps et $O(n+m)$ mémoire. L'alignement en bande limite la matrice à une diagonale de largeur $k$ : $O(kn)$ temps et mémoire pour des séquences proches. Enfin, le type entier compte : un tableau int16 plafonne à 32 767, insuffisant pour un chromosome de plusieurs millions de bases ; utilisez int32 ou int64.

Exemple travaillé à la main

Prenons $s=$ GATT, $t=$ GCATT, avec match $+1$, mismatch $-1$, gap $-1$. La matrice est $5 \times 6$. Les bords valent $0,-1,-2,-3,-4,-5$ et $0,-1,-2,-3,-4$.

Pour $F(1,1)$ (G/G), le candidat diagonal vaut $0+1=1$, le vertical $-1-1=-2$, l'horizontal $-1-1=-2$ : maximum 1 par la diagonale. Pour $F(1,2)$ (G/C), le candidat diagonal vaut $-1$, le vertical $-3$, l'horizontal $1-1=0$ : maximum 0 par la gauche. Pour $F(2,3)$ (A/A), le candidat diagonal vaut $0+1=1$, le vertical $-2$, l'horizontal $-1$ : maximum 1 par la diagonale.

Le tableau complet, avec la direction retenue, est :

La dernière cellule compare T et T : $F(3,4)+1=2+1=3$, $F(3,5)-1=0$, $F(4,4)-1=0$. Le score optimal est 3. Le traceback suit $(4,5)\to(3,4)\to(2,3)\to(1,2)\to(1,1)\to(0,0)$. Les deux premiers pas sont diagonaux ; le pas $(2,3)\to(1,2)$ est horizontal et aligne un gap avec C. On obtient :

s : G - A T T
t : G C A T T

Le score se relit : $1-1+1+1+1=3$, avec quatre identités et un gap. Le départage place le gap après le premier G ; d'autres paires de séquences peuvent avoir plusieurs alignements optimaux.

Implémentation

Les blocs ont été exécutés avec Python 3.13.5 et NumPy 2.3.1 le 2026-11-01 ; les sorties sont réelles. Le code utilise des indices 0-based et laisse la validation de l'alphabet au corrigé. Le premier bloc remplit la matrice complète ; le score se lit dans le coin inférieur droit.

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

La sortie vérifiée est :

[0, -1, -2, -3, -4, -5]
[-1, 1, 0, -1, -2, -3]
[-2, 0, 0, 1, 0, -1]
[-3, -1, -1, 0, 2, 1]
[-4, -2, -2, -1, 1, 3]
score = 3

Chaque ligne correspond à un préfixe de GATT, chaque colonne à un préfixe de GCATT. Le score 3 est en matrix[4][5] : les dimensions incluent les préfixes vides.

Le bloc suivant calcule le même score avec deux lignes, donc en $O(m)$ mémoire.

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

La sortie vérifiée est :

3
4
-4
6

Une séquence vide contre ACGT vaut quatre gaps, soit $-4$. Avec match $+2$ et gap $-2$, G-ATT vaut $2-2+2+2+2=6$. La matrice n'existe plus : cette fonction ne peut pas produire de traceback.

Le troisième bloc reconstruit l'alignement et compte identités et gaps.

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

La sortie vérifiée est :

G-ATT
GCATT
{'identities': 4, 'mismatches': 0, 'gaps': 1, 'aligned_length': 5, 'identity_rate': 0.8}

Le traceback recalcule la substitution au lieu de stocker des pointeurs ; si la matrice ne correspond pas au barème, aucun mouvement ne convient et la fonction lève une erreur. Le taux 0,8 vaut 4 identités sur 5 colonnes, gaps inclus.

Le quatrième bloc chiffre la pénalité affine de cette section, $g_o+k g_e$.

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

La sortie vérifiée est :

[0, -6, -7, -8, -9]
[0, -2, -3, -4, -5]
[0, -2, -2, -2, -2]

Avec $g_o=-5$, un gap d'une base coûte $-6$ ; avec $g_o=-1$, $-2$. La troisième ligne montre une ouverture payante sans extension. Un gap de longueur nulle coûte zéro ; le corrigé refuse les coûts positifs.

Le dernier bloc donne l'ordre de grandeur de la mémoire d'une matrice $10\,000 \times 10\,000$.

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

La sortie vérifiée est :

int16 200 Mo
int32 400 Mo
float64 800 Mo
max int16 : 32767

Une matrice de 100 millions de cellules pèse 200 Mo en int16, 400 Mo en int32 et 800 Mo en float64, hors objets Python. Le maximum de int16 est 32 767 : un alignement plus long déborde. La réduction à deux lignes reste la meilleure option sans traceback.

Étude de cas

Les deux fragments ci-dessous forment un exemple synthétique, pas une région réelle : le second a été obtenu à partir du premier en supprimant le triplet GTA et en remplaçant le dernier T par C. Ils servent à observer l'effet du barème sur un cas contrôlé, sans prétendre décrire un locus d'Escherichia coli.

• frag_a : GTTGACATGCGTAAAGGCTTACGT, 24 nt.
• frag_b : GTTGACATGCAAGGCTTACGC, 21 nt.

La question est simple : quel alignement global obtenons-nous, et combien le score change-t-il lorsque le barème change ? Le script suivant remplit la matrice, reconstruit l'alignement, compte les identités et les gaps, puis compare trois barèmes. Il ajoute un contrôle par permutation : la même composition en bases, dans un ordre aléatoire, est alignée avec frag_a.

from random import Random
 
def align(seq1, seq2, match=1, mismatch=-1, gap=-1):
    n, m = len(seq1), len(seq2)
    matrix = [[0] * (m + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        matrix[i][0] = i * gap
    for j in range(1, m + 1):
        matrix[0][j] = j * gap
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            substitution = match if seq1[i - 1] == seq2[j - 1] else mismatch
            matrix[i][j] = max(
                matrix[i - 1][j - 1] + substitution,
                matrix[i - 1][j] + gap,
                matrix[i][j - 1] + gap,
            )
    i, j = n, m
    top, bottom = [], []
    while i > 0 or j > 0:
        if i > 0 and j > 0:
            substitution = match if seq1[i - 1] == seq2[j - 1] else mismatch
            if matrix[i][j] == matrix[i - 1][j - 1] + substitution:
                top.append(seq1[i - 1])
                bottom.append(seq2[j - 1])
                i -= 1
                j -= 1
                continue
        if i > 0 and matrix[i][j] == matrix[i - 1][j] + gap:
            top.append(seq1[i - 1])
            bottom.append("-")
            i -= 1
            continue
        top.append("-")
        bottom.append(seq2[j - 1])
        j -= 1
    return matrix[-1][-1], "".join(reversed(top)), "".join(reversed(bottom))
 
frag_a = "GTTGACATGCGTAAAGGCTTACGT"
frag_b = "GTTGACATGCAAGGCTTACGC"
score, top, bottom = align(frag_a, frag_b)
print(score)
print(top)
print(bottom)
 
for params in [(1, -1, -1), (2, -1, -2), (1, 0, -1)]:
    score, top, bottom = align(frag_a, frag_b, *params)
    identities = sum(a == b for a, b in zip(top, bottom) if a != "-" and b != "-")
    gaps = sum(a == "-" or b == "-" for a, b in zip(top, bottom))
    print(params, "score =", score, "identities =", identities, "gaps =", gaps)
 
chars = list(frag_b)
Random(20261101).shuffle(chars)
shuffled = "".join(chars)
print(shuffled, align(frag_a, shuffled)[0])

La sortie vérifiée est :

16
GTTGACATGCGTAAAGGCTTACGT
GTTGACATGC---AAGGCTTACGC
(1, -1, -1) score = 16 identities = 20 gaps = 3
(2, -1, -2) score = 33 identities = 20 gaps = 3
(1, 0, -1) score = 17 identities = 20 gaps = 3
CGAACTGGTCCCAGATATGGT -1

Le meilleur alignement place un gap de trois bases dans frag_b, juste après GTTGACATGC. Les 21 colonnes appariées comprennent 20 identités et une substitution finale (T contre C) ; le taux d'identité est $20/24 \approx 0{,}833$. Le gap correspond exactement au triplet supprimé, ce qui est rassurant : le score et la structure de l'alignement racontent la même histoire.

Les trois barèmes sélectionnent le même alignement, mais produisent des scores très différents : 16, 33 et 17. Le barème $+2/-1/-2$ double les matches et les pénalités de gap, d'où un score 33. Le barème avec mismatch nul gagne un point, car la substitution finale ne coûte plus rien. Aucun de ces scores n'est « le vrai » : ils mesurent une hypothèse de correspondance, pas une distance physique.

Le contrôle par permutation est encore plus instructif. La séquence mélangée CGAACTGGTCCCAGATATGGT a la même composition que frag_b, mais un score de $-1$. L'ordre des bases porte donc l'essentiel du signal. Attention : un seul mélange ne donne pas une p-value. Pour tester la significativité, il faudrait répéter de nombreuses permutations, estimer la distribution des scores sous un modèle de fond et corriger pour les tests multiples. Le score brut reste une statistique descriptive.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-2/2.2-programmation-dynamique -q

python -m pytest exercises/module-2/2.2-programmation-dynamique -q

python -m pytest exercises/module-2/2.2-programmation-dynamique/projet -q

Résumé

• La programmation dynamique remplace l'énumération exponentielle des alignements par le remplissage de $(n+1)(m+1)$ cellules.
• La sous-structure optimale permet de décomposer un alignement en un sous-alignement optimal et une dernière colonne.
• Les sous-problèmes chevauchants justifient la mémorisation ou la tabulation : chaque état n'est calculé qu'une fois.
• La récurrence globale compare diagonale, gap dans la seconde séquence et gap dans la première ; les bords paient les gaps.
• Le score se lit en $F(n,m)$ et le traceback remonte depuis ce coin pour reconstruire les colonnes.
• Le temps est $O(nm)$, la mémoire $O(nm)$ ; deux lignes suffisent pour le score, mais pas pour le traceback.
• Une pénalité affine $g_o + k g_e$ distingue l'ouverture de l'extension ; la convention doit être documentée.
• Un score brut n'est ni une identité ni une preuve d'homologie : il dépend du barème, de la longueur et de la composition.

Références

• Needleman, S. B. et 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.
• Smith, T. F. et Waterman, M. S. (1981). Identification of common molecular subsequences. Journal of Molecular Biology, 147(1), 195–197.
• Bellman, R. (1957). Dynamic Programming. Princeton University Press.
• Durbin, R., Eddy, S., Krogh, A. et Mitchison, G. (1998). Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press.
• Hirschberg, D. S. (1975). A linear space algorithm for computing maximal common subsequences. Communications of the ACM, 18(6), 341–343.
• Biopython. Pairwise sequence alignment — Bio.Align.PairwiseAligner. Documentation officielle : https://biopython.org/docs/latest/api/Bio.Align.html
• Python Software Foundation. functools — lru_cache. Documentation officielle : https://docs.python.org/3/library/functools.html#functools.lru_cache