Matrice de points

Prérequis

Cette section suppose que vous savez lire un enregistrement FASTA et distinguer une séquence d'ADN d'une séquence protéique, comme dans 1.2. Les conventions de notation du cours restent valables : une séquence est notée $s = s_1 s_2 \dots s_n$, les positions sont 1-based dans le texte et 0-based dans le code. Une familiarité avec les tableaux NumPy n'est pas nécessaire : les quelques opérations utilisées sont expliquées au fur et à mesure. La section 1.4 fournit un utile contrepoint : une expression régulière cherche un motif exact, alors qu'une matrice de points compare toutes les positions deux à deux.

Introduction

Comparer deux séquences est l'une des opérations les plus fréquentes en bioinformatique. Deux gènes peuvent partager un ancêtre commun, deux protéines peuvent remplir la même fonction, ou une région du génome peut contenir une copie d'un segment distant. L'évolution accumule des substitutions, des insertions, des délétions, des duplications et des inversions. Avant de calculer un alignement, il est souvent utile de regarder la structure globale de la comparaison : où se trouvent les zones qui se ressemblent, comment elles sont décalées, et si plusieurs copies d'un même segment apparaissent. La matrice de points, ou dot plot, répond à cette question avec une image.

L'idée est simple. On place une séquence sur l'axe horizontal, l'autre sur l'axe vertical, et on remplit une case lorsque les deux positions comparées portent des caractères identiques ou suffisamment similaires. Une région conservée dessine une diagonale ; une insertion décale la diagonale ; une duplication fait apparaître une diagonale parallèle ; une inversion produit une anti-diagonale. Cette lecture visuelle précède l'alignement et aide à choisir un modèle, une pénalité de gap ou une fenêtre d'analyse. Elle ne remplace pas un alignement : elle montre toutes les correspondances possibles, y compris celles qui ne participent pas à un chemin cohérent.

Le défi computationnel est double. D'une part, la matrice complète contient $n \times m$ cases, ce qui devient prohibitif pour deux génomes entiers. D'autre part, le bruit de fond est important : dans deux séquences d'ADN aléatoires, environ un quart des cases sont remplies par simple coïncidence. Sans fenêtre, sans filtre de faible complexité et sans seuil, le signal biologique se noie dans les points isolés. Cette section construit une matrice avec NumPy, apprend à la lire, à la lisser et à en extraire les diagonales dominantes. Le fil rouge reste Escherichia coli K-12 MG1655, accession NC_000913.3, mais les expériences quantitatives de la section sont des exemples synthétiques explicitement étiquetés.

Théorie

Comparer sans aligner

Un alignement répond à la question : quel est le meilleur appariement global ou local entre deux séquences, compte tenu d'un schéma de score et de pénalités de gap ? Needleman et Wunsch ont formulé ce problème en programmation dynamique en 1970, et Smith et Waterman en ont proposé la version locale en 1981. Une matrice de points répond à une autre question : pour chaque paire de positions, les deux caractères se ressemblent-ils ? Elle ne choisit pas un chemin, elle affiche le paysage complet des correspondances possibles. C'est une étape d'exploration, pas une preuve d'homologie.

Formellement, pour deux séquences $a = a_1 \dots a_n$ et $b = b_1 \dots b_m$, la matrice de points d'identité est définie par

La matrice a $n$ lignes et $m$ colonnes : la ligne $i$ correspond à $a_i$, la colonne $j$ à $b_j$. Cette convention doit être fixée avant d'écrire du code, car transposer la matrice inverse la lecture des insertions et des délétions. Pour des protéines, on remplace souvent l'égalité stricte par un seuil de similarité :

où $\mathrm{sim}$ peut être un score de substitution comme BLOSUM62 et $\tau$ un seuil, par exemple zéro pour ne garder que les appariements non défavorables. On parle alors de matrice de similarité plutôt que d'identité : une case remplie signale une paire intéressante, pas un alignement validé.

La diagonale principale, d'offset $d = j - i = 0$, correspond aux positions alignées sans décalage. Un offset positif signifie que la colonne est en avance sur la ligne : la séquence $b$ contient une insertion, ou $a$ une délétion, avant la région considérée. Un offset négatif correspond au cas inverse ; la matrice n'est symétrique que si $a = b$.

Fenêtres, bruit et faible complexité

Sur de l'ADN uniforme, la probabilité qu'une case soit remplie par identité stricte vaut $1/4$. Pour deux séquences de 1 000 pb, on attend donc environ $1\,000\,000 / 4 = 250\,000$ points, soit un quart de la matrice. Avec une fenêtre de $w$ caractères, une case n'est remplie que si les $w$ caractères concordent :

La probabilité d'une correspondance fortuite devient $4^{-w}$ pour l'ADN, et $20^{-w}$ pour des protéines si les vingt acides aminés sont équiprobables. Pour 1 000 pb et $w = 8$, l'espérance tombe à $1\,000\,000 / 65\,536 \approx 15$ points. La fenêtre exige une correspondance plus longue, donc plus rare au hasard. Elle a un coût : un exon court, un domaine de quelques résidus ou une région très divergente peuvent disparaître si la fenêtre est trop grande. Le choix de $w$ est un compromis à documenter, et la matrice perd $w-1$ lignes et colonnes : ses coordonnées désignent des débuts de fenêtres, pas des bases isolées.

Les répétitions et la faible complexité trompent encore plus qu'un alphabet uniforme. Une région riche en A et T, ou un microsatellite ATATAT, produit de nombreuses correspondances locales et des motifs en peigne. Une séquence de poly-A contre elle-même remplit une diagonale entière, même sans aucune homologie au sens évolutif. Les familles de répétitions créent des diagonales parallèles ; les répétitions inversées créent des anti-diagonales. Un filtre de faible complexité, fondé par exemple sur l'entropie

ou un masquage des régions peu complexes, réduit ces faux positifs. La composition biaisée suffit à produire des points, sans conservation évolutive.

Diagonales, inversions et duplications

La forme du nuage de points est une signature de l'événement. Le tableau suivant résume les motifs les plus utiles, avec la convention $d = j - i$.

Une duplication interne crée souvent une « marche d'escalier » : la diagonale principale se coupe, reprend avec un offset égal à la longueur insérée, et une diagonale parallèle apparaît pour la copie. Une inversion, elle, ne se lit pas sur une diagonale nord-ouest → sud-est : le segment inversé produit une anti-diagonale, c'est-à-dire une ligne où $i + j$ est constant. Un algorithme qui ne parcourt que les offsets $j - i$ ne la verra pas ; il faut soit inverser le complément de la séquence, soit parcourir explicitement les anti-diagonales.

Complexité et mémoire

Remplir la matrice demande $n \times m$ comparaisons, soit un temps $O(nm)$. La mémoire est également $O(nm)$ si l'on conserve la matrice complète. Un tableau booléen NumPy occupe un octet par case ; pour deux séquences de 10 000 pb, la matrice fait $10^8$ cases, soit environ 100 Mo. Pour deux génomes de 5 millions de pb, elle ferait $25 \times 10^{12}$ cases, soit environ 25 To : le calcul complet est hors de portée. La vectorisation NumPy accélère les constantes, pas l'asymptote.

Plusieurs stratégies évitent la matrice complète. On peut stocker uniquement les positions des correspondances, mais elles restent nombreuses sur de l'ADN. On peut indexer les $k$-mers de l'une des séquences et ne comparer que les positions partageant un $k$-mer, avec un coût attendu proche de $O(n+m)$ pour la construction et la recherche. On peut aussi se limiter à une bande autour d'une diagonale attendue, ce qui donne $O((n+m)w)$ pour une largeur $w$. Les aligneurs modernes combinent index de $k$-mers, minimizers et programmation dynamique bandée. Le dot plot complet reste un outil d'exploration pour de petites régions, pas un passage obligé pour deux génomes entiers.

Exemple travaillé à la main

Prenons $A =$ ACGTA et $B =$ TACGT, deux séquences synthétiques de longueur 5. On compare chaque caractère de $A$ (lignes) à chaque caractère de $B$ (colonnes). La matrice obtenue est :

     T  A  C  G  T
A  [ 0  1  0  0  0 ]
C  [ 0  0  1  0  0 ]
G  [ 0  0  0  1  0 ]
T  [ 1  0  0  0  1 ]
A  [ 0  1  0  0  0 ]

Il y a 25 comparaisons et 6 cases remplies. Quatre points forment une diagonale d'offset $+1$ : les positions 1 à 4 de $A$ (ACGT) correspondent aux positions 2 à 5 de $B$ (ACGT). L'identité sur cette diagonale est de 100 % sur 4 bases. Les deux autres points, en (4,1) et (5,2) en coordonnées 1-based, sont des coïncidences isolées : un T contre un T, puis un A contre un A. Avec min_length = 3, ils disparaissent et seule la diagonale de longueur 4 est retenue : le filtre sépare un signal de quatre bases de deux coïncidences.

Implémentation

La construction tient en une ligne grâce au broadcasting de NumPy. On encode chaque séquence en un tableau d'octets, on ajoute un axe aux deux tableaux, puis on compare toutes les paires. Le résultat est un tableau booléen de forme (len(seq_a), len(seq_b)). Le code suivant est volontairement autonome : il ne lit aucun fichier et n'utilise que NumPy.

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

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

[[0 1 0 0 0]
 [0 0 1 0 0]
 [0 0 0 1 0]
 [1 0 0 0 1]
 [0 1 0 0 0]]
matches : 6
diagonales >= 3 : [(1, 0, 1, 4)]
diagonales >= 1 : [(1, 0, 1, 4), (-3, 3, 0, 2)]

Trois points méritent attention. D'abord, np.frombuffer évite de construire une liste Python de caractères ; l'encodage ASCII est explicite et une séquence non ASCII doit être rejetée ou normalisée. Ensuite, le broadcasting compare un tableau de forme $(n,1)$ à un tableau de forme $(1,m)$ et produit une matrice $(n,m)$ sans boucle Python. Enfin, diagonal_runs parcourt chaque offset $d$, découpe les suites de True et filtre par longueur. Les indices renvoyés sont 0-based ; pour afficher des coordonnées biologiques, il faut ajouter 1. Une anti-diagonale n'est jamais visitée par cette boucle, ce qui explique pourquoi une inversion n'apparaît pas dans le résultat.

La visualisation utilise Matplotlib avec le backend Agg, qui écrit un fichier sans ouvrir de fenêtre. Le script compare trois variantes d'une séquence synthétique de 60 pb : une insertion, une duplication interne et une inversion. La ligne du haut montre la matrice brute, dominée par le bruit ; la ligne du bas applique une fenêtre de 5.

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
 
def dotplot(seq_a, seq_b):
    a = np.frombuffer(seq_a.upper().encode("ascii"), dtype=np.uint8)
    b = np.frombuffer(seq_b.upper().encode("ascii"), dtype=np.uint8)
    return a[:, None] == b[None, :]
 
def window_dotplot(seq_a, seq_b, window):
    a = np.frombuffer(seq_a.upper().encode("ascii"), dtype=np.uint8)
    b = np.frombuffer(seq_b.upper().encode("ascii"), dtype=np.uint8)
    n = len(a) - window + 1
    m = len(b) - window + 1
    matrix = np.zeros((max(0, n), max(0, m)), dtype=bool)
    if n <= 0 or m <= 0:
        return matrix
    matrix[:] = True
    for k in range(window):
        matrix &= a[k:k + n, None] == b[None, k:k + m]
    return matrix
 
reference = "TAAGCCACGCTTGTGATGACATGTCGGTAGGTGCCCACAATGACCGTACCTCAACGCCGC"
complement = str.maketrans("ACGT", "TGCA")
reverse_complement = lambda s: s.translate(complement)[::-1]
 
pairs = [
    ("Insertion de 6 pb", reference[:20] + "AAGCTT" + reference[20:]),
    ("Duplication interne", reference[:30] + reference[10:20] + reference[30:]),
    ("Inversion de 12 pb", reference[:34] + reverse_complement(reference[34:46]) + reference[46:]),
]
 
fig, axes = plt.subplots(2, 3, figsize=(12, 8), constrained_layout=True)
for column, (label, query) in enumerate(pairs):
    raw = dotplot(reference, query)
    windowed = window_dotplot(reference, query, 5)
    axes[0, column].imshow(raw, cmap="Greys", origin="upper", aspect="equal", interpolation="nearest")
    axes[0, column].set_title(f"{label} : brut, {int(raw.sum())} points", fontsize=10)
    axes[1, column].imshow(windowed, cmap="Greys", origin="upper", aspect="equal", interpolation="nearest")
    axes[1, column].set_title(f"fenêtre 5 : {int(windowed.sum())} points", fontsize=10)
    for row in (0, 1):
        axes[row, column].set_xlabel(f"référence, {len(reference)} pb")
        axes[row, column].set_ylabel(f"requête, {len(query)} pb")
 
fig.suptitle("Matrice de points : bruit brut et lissage par fenêtre", fontsize=13)
fig.savefig("dotplot_synthetique.png", dpi=150)
plt.close(fig)
print("Figure enregistrée : dotplot_synthetique.png (1800 x 1200 px)")

La sortie vérifiée est Figure enregistrée : dotplot_synthetique.png (1800 x 1200 px). La figure produite le 2026-11-01 avec Matplotlib 3.10.7 montre trois panneaux bruts à environ 25 % de cases noires, puis trois panneaux fenêtrés à 1–2 % de points. Dans le panneau « insertion », la diagonale se coupe et reprend décalée vers la droite ; dans « duplication », trois segments parallèles apparaissent ; dans « inversion », une anti-diagonale traverse le carré. L'axe origin="upper" place la première base en haut à gauche, et interpolation="nearest" évite d'inventer des valeurs entre les cases.

Étude de cas

Nous utilisons trois séquences synthétiques dérivées d'une même référence de 60 pb. Elles ne proviennent pas de NC_000913.3 : le génome réel sert de fil rouge conceptuel, mais aucune coordonnée réelle n'est nécessaire ici, et aucune n'est inventée. La question est la suivante : quelle signature visuelle et quelles diagonales dominantes produit chacun des trois événements, et que peut-on conclure sans alignement ?

La procédure est la suivante. On calcule d'abord la matrice brute, puis on compte les points. On applique ensuite une fenêtre de 5 et on extrait les trois diagonales de plus grand support avec la fonction dotplot_report du mini-projet. Enfin, on compare les motifs attendus à ceux du tableau 2.1.1.

Tableau 2.1.1 — Résultats mesurés sur les séquences synthétiques, avec NumPy 2.3.1 et une fenêtre de 5. Les pourcentages entre parenthèses sont la densité de points de la matrice.

Pour l'insertion, la variante B contient AAGCTT après la position 20 de la référence. La diagonale dominante a un offset $+6$, un support de 36 fenêtres et une identité de 100 % sur 40 pb : c'est exactement la partie située après l'insertion. La diagonale d'offset 0 conserve un support de 17 fenêtres, limité au début de la séquence. L'interprétation est directe : six bases ont été insérées, et le reste de la séquence est inchangé. Le dot plot suggère l'événement, mais il ne dit rien sur son origine ni sur sa fonction.

Pour la duplication, la variante B reprend les positions 11 à 20 et les insère après la position 30. La diagonale principale est coupée en deux : un premier segment d'offset 0 (support 26) et un second d'offset $+10$ (support 26). La copie insérée produit une troisième diagonale d'offset $+20$ (support 6). Le rapport choisit l'offset 0 comme meilleure diagonale, car les deux premières ont le même support et le tri départage par offset croissant. Cette égalité n'est pas un défaut : elle traduit le fait que la duplication coupe la correspondance principale en deux parties de longueur comparable. La diagonale parallèle est la signature de la copie.

Pour l'inversion, la variante B remplacer les positions 35 à 46 par leur complément inverse. La diagonale principale garde un support de 42 fenêtres et une identité de 86,67 % sur 60 pb, mais elle est interrompue : les suites nord-ouest → sud-est mesurent 35 et 15 bases. Le segment inversé forme une anti-diagonale, visible dans la figure mais absente de la liste des diagonales dominantes, car dotplot_report ne parcourt que les offsets $j - i$. L'identité de 86,67 % ne signifie pas que 13 % des bases ont été substituées : une inversion n'est pas une série de substitutions indépendantes, et les coïncidences de bases peuvent modifier le chiffre. Pour confirmer l'inversion, il faut aligner la séquence avec son complément inverse ou tester explicitement les anti-diagonales.

Le dot plot sert à formuler une hypothèse sur la structure de la comparaison ; l'alignement, l'annotation et un test statistique servent ensuite à la vérifier. Une diagonale courte, une anti-diagonale ou un point isolé ne suffisent jamais à conclure. Sur un génome réel comme NC_000913.3, les répétitions, les séquences d'insertion et les régions de faible complexité produiraient de nombreux motifs secondaires qu'il faudrait confronter à l'annotation.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-2/2.1-matrice-de-points -q

python -m pytest exercises/module-2/2.1-matrice-de-points -q

python -m pytest exercises/module-2/2.1-matrice-de-points/projet -q

Résumé

• Une matrice de points compare toutes les paires de positions et révèle la structure globale d'une comparaison.
• Les lignes correspondent à la première séquence, les colonnes à la seconde ; la convention doit être documentée.
• Une fenêtre de $w$ caractères réduit le bruit de fond, mais efface les correspondances plus courtes et modifie les coordonnées.
• Une insertion décale la diagonale, une duplication crée une diagonale parallèle, une inversion une anti-diagonale.
• Le coût d'une matrice complète est $O(nm)$ en temps et en mémoire ; les grands génomes exigent un index ou une bande.
• Un point ou une diagonale courte ne prouve pas l'homologie : il faut un score, un contrôle de faible complexité et souvent un alignement.
• Les exemples de cette section sont synthétiques ; NC_000913.3 sert de fil rouge, mais aucune coordonnée réelle n'a été inventée.

Références

• 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.
• Smith, T. F. & Waterman, M. S. (1981). Identification of common molecular subsequences. Journal of Molecular Biology, 147(1), 195–197.
• 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. et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11), 1422–1423.
• NumPy Developers. NumPy documentation : broadcasting, ndarray, frombuffer. https://numpy.org/doc/stable/
• Matplotlib Development Team. Matplotlib documentation : imshow, backend Agg, sauvegarde des figures. https://matplotlib.org/stable/
• NCBI RefSeq. Escherichia coli str. K-12 substr. MG1655, complete genome, accession NC_000913.3. https://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3