Superposition et RMSD

Prérequis

Cette section suppose que vous savez lire un fichier PDB et parcourir la hiérarchie Structure > Model > Chain > Residue > Atom, comme présenté en 4.2. Le format lui-même est décrit en 4.1. Elle utilise les tableaux NumPy, la norme euclidienne et le produit scalaire, ainsi que la lecture d'un fichier PDB avec Biopython vue en 1.3. La comparaison de séquences de 2.3 sera utile pour construire une correspondance entre résidus homologues. Le module 4 replace cette opération parmi les autres analyses structurales.

Introduction

Comparer deux structures de protéines est une question biologique courante. Deux copies d'une même enzyme cristallisées dans des conditions différentes, deux conformations d'un canal ionique, ou deux homologues issus d'espèces éloignées partagent souvent un même repliement tout en différant par des boucles, des domaines mobiles ou des interactions cristallines. Pour décrire ces différences, il faut d'abord décider quels atomes se correspondent, puis trouver le déplacement rigide qui les rapproche le mieux. La superposition et le RMSD (root mean square deviation, écart quadratique moyen) sont les réponses classiques à ce problème.

L'enjeu biologique est de distinguer une différence de conformation d'un simple changement de repère. Une protéine n'a pas d'orientation intrinsèque : la comparer à une autre sans alignement revient à comparer deux objets après une rotation arbitraire. Une fois le repère fixé, le RMSD résume l'écart entre les positions des atomes homologues. Il permet de suivre une trajectoire de dynamique moléculaire, de valider un modèle prédit contre une structure expérimentale, ou de quantifier l'effet d'une mutation.

L'enjeu computationnel est double. Il faut résoudre un problème d'optimisation : trouver la rotation et la translation qui minimisent une somme de distances au carré. Il faut aussi construire une correspondance crédible entre atomes, car le RMSD n'a de sens que pour des paires définies. Nous construirons la solution de Kabsch pas à pas, en montrant le rôle du centrage, de la matrice de covariance et du déterminant. Nous l'appliquerons à la crambine 1CRN, puis nous verrons pourquoi un RMSD global peut cacher une flexion entre deux domaines parfaitement rigides.

Théorie

Correspondance d'atomes homologues

Une superposition met en relation deux ensembles ordonnés de points. On note $x_i$ les coordonnées d'un atome de la structure cible et $y_i$ celles de l'atome correspondant dans la structure mobile, pour $i = 1, \dots, N$. La correspondance est une bijection entre les deux sélections : chaque atome de la cible est apparié à exactement un atome mobile, et réciproquement. Sans cette bijection, la formule du RMSD n'a pas de sens.

En pratique, la correspondance se construit à partir de l'annotation. Pour une même protéine, on peut apparier les atomes par nom de chaîne, numéro de résidu, nom de résidu et nom d'atome. Pour deux homologues, il faut d'abord aligner les séquences ou les structures, puis reporter les positions alignées sur les atomes. On choisit souvent les carbones alpha (Cα), plus robustes aux mouvements des chaînes latérales, ou le squelette N, Cα, C, O. Comparer tous les atomes lourds est plus sensible aux conformations latérales et exige une correspondance atomique complète.

Le choix de la sélection fait partie du résultat. Un RMSD calculé sur 46 Cα, sur 300 atomes de squelette ou sur tous les atomes lourds ne mesure pas la même chose. Il faut donc toujours indiquer le nombre de paires $N$, le type d'atomes, les chaînes et les résidus inclus, ainsi que le traitement des atomes manquants et des conformations alternatives. Un chiffre de RMSD sans ces précisions est difficile à interpréter, voire trompeur.

Transformation rigide

Une transformation rigide déplace un objet sans le déformer. Elle s'écrit, pour des vecteurs colonnes dans $\mathbb{R}^D$, $x_i \approx R y_i + t$, où $R$ est une matrice de rotation $D \times D$ et $t$ un vecteur de translation. En trois dimensions, elle possède six degrés de liberté : trois pour la rotation et trois pour la translation. La matrice $R$ doit être orthogonale, c'est-à-dire vérifier $R^{\top} R = I$, et son déterminant doit valoir $+1$. Ces deux conditions garantissent que les distances et les angles sont conservés, et que l'orientation de l'espace n'est pas inversée.

Une matrice orthogonale de déterminant $-1$ est une réflexion, parfois appelée rotation impropre. Elle conserve les distances mais transforme un objet en son image miroir. Une protéine est chirale : ses acides aminés sont presque tous de configuration L, et une réflexion produirait une molécule qui n'existe pas. La distinction entre rotation propre et réflexion est donc une contrainte physique, pas une subtilité algébrique.

Définition du RMSD

Après avoir appliqué une transformation, on mesure l'écart résiduel par le RMSD : $\mathrm{RMSD} = \sqrt{\frac{1}{N}\sum_{i=1}^{N} \lVert x_i - R y_i - t \rVert^2}$. La somme parcourt les $N$ paires homologues, et la norme est la norme euclidienne. Le résultat s'exprime dans l'unité des coordonnées, généralement l'ångström. Le RMSD est nul si et seulement si chaque atome mobile atteint exactement son homologue cible. Il croît avec l'amplitude des écarts, mais il s'agit d'une moyenne quadratique : un petit nombre d'atomes très déplacés peut dominer le chiffre.

Le RMSD n'est pas une distance maximale. Deux structures peuvent avoir un RMSD faible avec un ou deux atomes très mal placés, ou un RMSD plus élevé avec un décalage réparti sur tous les atomes. Il ne dit rien non plus sur la direction des mouvements ni sur leur signification biologique. C'est un résumé, utile mais incomplet, qui doit être accompagné d'une inspection des écarts par résidu ou par domaine.

Centrage et élimination de la translation

Pour une rotation $R$ fixée, la translation optimale se calcule analytiquement. La fonction à minimiser est une somme de carrés, donc convexe en $t$. En annulant le gradient par rapport à $t$, on obtient $t = \bar{x} - R \bar{y}$, où $\bar{x}$ et $\bar{y}$ sont les barycentres des deux nuages. Cette relation a une lecture simple : après rotation, on translate le barycentre mobile sur le barycentre cible. Le problème se réduit alors à trouver la rotation qui aligne les coordonnées centrées $x_i' = x_i - \bar{x}$ et $y_i' = y_i - \bar{y}$.

Le centrage préalable n'est pas une étape cosmétique. Si l'on oublie de centrer, la matrice de covariance contient un terme de translation et la rotation calculée est fausse. Si l'on centre seulement l'un des deux nuages, le biais persiste. La bonne pratique est de centrer les deux, de calculer la rotation, puis de reconstruire la translation avec les barycentres d'origine.

Algorithme de Kabsch par SVD

Une fois les nuages centrés, on cherche $R$ qui minimise $\sum_i \lVert x_i' - R y_i' \rVert^2$. En développant le carré, les termes $\sum_i \lVert x_i' \rVert^2$ et $\sum_i \lVert y_i' \rVert^2$ ne dépendent pas de $R$. Le problème équivaut donc à maximiser $\sum_{i=1}^{N} x_i'^{\top} R y_i' = \operatorname{tr}(R H)$, où $H$ est la matrice de covariance croisée $H = \sum_{i=1}^{N} y_i' x_i'^{\top}$. Avec des tableaux NumPy de forme $(N, D)$ où chaque ligne est un atome, on calcule directement $H = M^{\top} T$, avec $M$ les coordonnées mobiles centrées et $T$ les coordonnées cibles centrées. La décomposition en valeurs singulières s'écrit $H = U S V^{\top}$, avec $U$ et $V$ orthogonales et $S$ diagonale positive. La rotation qui maximise la trace est alors $R = V U^{\top}$. Cette formule est le cœur de la méthode de Kabsch. Les valeurs singulières mesurent la qualité de l'accord selon les axes principaux du nuage ; les matrices $U$ et $V$ alignent ces axes. En dimension trois, la décomposition porte sur une matrice $3 \times 3$, donc son coût est négligeable devant le parcours des atomes.

Déterminant, réflexion et chiralité

La matrice $V U^{\top}$ est orthogonale, mais son déterminant peut valoir $+1$ ou $-1$. Si le déterminant vaut $+1$, c'est une rotation propre et le problème est résolu. S'il vaut $-1$, la décomposition a trouvé une réflexion, qui n'est pas acceptable pour une molécule chirale. La correction consiste à inverser le signe de la dernière colonne de $V$ associée à la plus petite valeur singulière, puis à recalculer $R = V U^{\top}$. On obtient ainsi la rotation propre la plus proche de la solution orthogonale.

Cette correction ne rend pas le RMSD nul si les données sont réellement incompatibles avec une rotation. Elle choisit simplement la meilleure rotation possible. Un déterminant négatif doit donc alerter : il peut signaler des coordonnées miroir, une erreur de signe dans un axe, ou un problème de construction du modèle. Accepter la réflexion donnerait une superposition visuellement plausible mais biologiquement fausse.

Complexité

Pour $N$ atomes en dimension $D$, la construction de $H$ coûte $O(ND^2)$, la décomposition $O(D^3)$ et la transformation $O(ND^2)$. En dimension fixe $D = 3$, le coût total est linéaire en $N$, et la mémoire est $O(ND)$. La méthode est donc très rapide, même pour des milliers d'atomes. Le coût dominant vient souvent de la construction de la correspondance et du parsing du fichier, pas de la décomposition elle-même.

Exemple travaillé à la main

Prenons trois points mobiles $y_1 = (-1, 0)$, $y_2 = (1, 0)$ et $y_3 = (0, 1)$, et trois points cibles $x_1 = (3, 3)$, $x_2 = (3, 5)$ et $x_3 = (2, 4)$. Les barycentres valent $\bar{y} = (0, 1/3)$ et $\bar{x} = (8/3, 4)$. Les coordonnées centrées sont $y_1' = (-1, -1/3)$, $y_2' = (1, -1/3)$, $y_3' = (0, 2/3)$ d'une part, et $x_1' = (1/3, -1)$, $x_2' = (1/3, 1)$, $x_3' = (-2/3, 0)$ d'autre part. La matrice de covariance croisée vaut $H = \sum_i y_i' x_i'^{\top} = \begin{pmatrix} 0 & 2 \\ -2/3 & 0 \end{pmatrix}$.

On calcule $H^{\top} H = \operatorname{diag}(4/9, 4)$, donc les valeurs singulières sont $2$ et $2/3$. Une décomposition possible est $U = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$ et $V = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$. La rotation vaut alors $R = V U^{\top} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$, soit une rotation de $90^\circ$ dans le sens direct. Son déterminant vaut $+1$. La translation se déduit des barycentres : $t = \bar{x} - R \bar{y} = (8/3, 4) - (-1/3, 0) = (3, 4)$.

Vérifions un point : $R y_1 + t = (0, -1) + (3, 4) = (3, 3) = x_1$. Les deux autres points donnent de même $x_2$ et $x_3$, donc le RMSD vaut exactement $0$. Cet exemple montre que la rotation et la translation se lisent directement dans les données, sans recherche numérique.

RMSD global et RMSD local

Le RMSD global est calculé après une seule transformation appliquée à tous les atomes sélectionnés. Il résume l'accord d'ensemble, mais il peut compenser des écarts opposés : un domaine peut être déplacé vers la gauche et un autre vers la droite, avec un RMSD global modéré. Le RMSD local peut désigner deux choses qu'il faut distinguer. Après la superposition globale, on peut calculer le RMSD résiduel de chaque région : il mesure alors l'écart de cette région dans le repère commun. On peut aussi superposer chaque région indépendamment, ce qui donne le RMSD local au sens d'un ajustement propre à la région.

Ces deux quantités ne répondent pas à la même question. L'écart résiduel local après ajustement global indique où le modèle global se trompe. Le RMSD après ajustement local indique si la région est rigide en elle-même. Une région peut avoir un écart résiduel global élevé et un RMSD local nul : elle est rigide, mais mal placée par la transformation commune. C'est exactement le signal d'une articulation entre domaines. À l'inverse, un écart résiduel global faible avec un RMSD local élevé suggère une déformation interne compensée par l'ajustement global.

Limites du RMSD

Trois limites doivent rester visibles. La première est la correspondance : un RMSD calculé sur des paires arbitraires ou sur des sélections différentes n'a pas de sens biologique, même si le calcul est exact. La deuxième est la flexibilité : une seule transformation rigide ne décrit pas un mouvement concerté, et un RMSD global peut mélanger des domaines qui bougent différemment. La troisième est la sensibilité aux valeurs aberrantes : comme le RMSD élève les écarts au carré, quelques atomes très déplacés, une boucle mal résolue ou un ligand mal apparié peuvent dominer la moyenne. Des variantes robustes, des RMSD par résidu et une inspection visuelle sont alors nécessaires.

Implémentation

Le code ci-dessous reprend les quatre fonctions de la section. Il valide les dimensions, centre les nuages, calcule la rotation par SVD, corrige le déterminant, puis reconstruit la translation. Les indices NumPy sont 0-based, alors que les résidus d'un fichier PDB sont numérotés selon l'annotation.

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 :

Rotation :
[[ 0. -1.]
 [ 1.  0.]]
Translation : [3. 4.]
RMSD : 0.0
Points transformés :
[[3. 3.]
 [3. 5.]
 [2. 4.]]

Quatre points méritent attention. D'abord, center renvoie un nouveau tableau et ne modifie pas l'entrée ; cela évite les effets de bord. Ensuite, kabsch_rotation centre elle-même les deux nuages, ce qui la rend sûre à appeler directement, mais superimpose recalcule les barycentres d'origine pour la translation. Troisièmement, la correction de déterminant agit sur la dernière valeur singulière, qui est la plus petite : c'est la direction la moins contrainte par les données. Enfin, le RMSD est calculé sur les points transformés, pas sur une formule théorique, ce qui vérifie la cohérence de toute la chaîne.

Comparaison avec Bio.PDB.Superimposer

Biopython fournit une implémentation éprouvée dans Bio.PDB.Superimposer. Elle s'utilise avec des objets Atom, et non avec de simples tableaux. La méthode set_atoms(fixed, moving) calcule la transformation qui amène moving sur fixed, puis expose rms et rotran. Le code suivant superpose une structure sur elle-même, ce qui doit donner un RMSD nul.

from Bio.PDB import PDBParser, Superimposer
 
parser = PDBParser(QUIET=True)
structure = parser.get_structure("1CRN", "public/pdb/1CRN.pdb")
atoms = [res["CA"] for res in structure[0]["A"] if "CA" in res]
 
superimposer = Superimposer()
superimposer.set_atoms(atoms, atoms)
print("RMSD :", superimposer.rms)
print("Rotation :")
print(superimposer.rotran[0])

La sortie donne un RMSD de l'ordre de $10^{-15}$ Å, soit zéro à la précision machine. Un détail de convention mérite d'être signalé : rotran[0] est organisée pour une multiplication à droite, mobile @ rotation, alors que notre fonction renvoie une rotation destinée à mobile @ rotation.T. Les deux matrices sont donc transposées l'une de l'autre, même si elles représentent la même transformation. Le RMSD, lui, est identique. Cette différence de convention est une source classique de confusion lorsqu'on compare les sorties de deux bibliothèques.

Étude de cas

Nous utilisons la crambine, structure PDB 1CRN, dont le fichier local contient 46 résidus. Nous sélectionnons les 46 atomes Cα, du résidu THR 1 au résidu ASN 46. Les coordonnées sont réelles. Pour tester la méthode, nous fabriquons deux copies synthétiques de cette structure : une copie rigide, obtenue par une rotation de $25^\circ$ autour de l'axe $z$ et de $10^\circ$ autour de l'axe $x$, suivie d'une translation de $(12, -5, 8)$ Å ; puis une copie pliée, où la queue C-terminale (résidus 31 à 46) subit en plus une rotation de $20^\circ$ autour du Cα 31. La déformation est synthétique, mais elle s'applique à des coordonnées expérimentales réelles.

La question est double. Notre implémentation retrouve-t-elle exactement la transformation rigide de la première copie ? Et que signifie un RMSD global non nul lorsque chaque sous-région de la seconde copie est en fait rigide ? Nous procédons en cinq étapes : parser le fichier, construire les deux copies, superposer avec notre fonction, calculer les écarts résiduels par groupe, puis comparer le RMSD global à celui de Bio.PDB.Superimposer.

from math import cos, radians, sin
import sys
 
import numpy as np
from Bio.PDB import PDBParser, Superimposer
 
sys.path.insert(0, "exercises/module-4/4.3-superposition-rmsd")
from solution_4_3 import rmsd, superimpose
 
 
def rotation_z(angle):
    c, s = cos(angle), sin(angle)
    return np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]])
 
 
def rotation_x(angle):
    c, s = cos(angle), sin(angle)
    return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]])
 
 
parser = PDBParser(QUIET=True)
structure = parser.get_structure("1CRN", "public/pdb/1CRN.pdb")
chain = next(next(structure.get_models()).get_chains())
residues = [res for res in chain if res.id[0] == " " and "CA" in res]
atoms = [res["CA"] for res in residues]
target = np.array([atom.coord for atom in atoms], dtype=float)
 
q = rotation_z(radians(25.0)) @ rotation_x(radians(10.0))
offset = np.array([12.0, -5.0, 8.0])
mobile_rigid = target @ q.T + offset
print("atomes CA :", len(atoms))
print("RMSD copie rigide :", f"{superimpose(mobile_rigid, target)['rmsd']:.6e}")
 
hinge = target[30].copy()
tail = (target[30:] - hinge) @ rotation_z(radians(20.0)).T + hinge
mobile_flex = mobile_rigid.copy()
mobile_flex[30:] = tail @ q.T + offset
 
fit = superimpose(mobile_flex, target)
print("RMSD copie pliée :", f"{fit['rmsd']:.6f}")
for name, start, stop in (("coeur", 0, 30), ("queue", 30, len(target))):
    indices = np.arange(start, stop)
    delta = fit["transformed"][indices] - target[indices]
    global_residual = np.sqrt(np.mean(np.sum(delta * delta, axis=1)))
    local_rmsd = rmsd(mobile_flex[indices], target[indices])
    print(f"{name} ({len(indices)} atomes) : global={global_residual:.6f} local={local_rmsd:.6f}")
 
structure_b = parser.get_structure("1CRN-mobile", "public/pdb/1CRN.pdb")
chain_b = next(next(structure_b.get_models()).get_chains())
atoms_b = [res["CA"] for res in chain_b if res.id[0] == " " and "CA" in res]
for atom, coord in zip(atoms_b, mobile_flex):
    atom.set_coord(coord)
superimposer = Superimposer()
superimposer.set_atoms(atoms, atoms_b)
print("Bio.PDB.Superimposer :", f"{superimposer.rms:.6f}")

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

atomes CA : 46
RMSD copie rigide : 4.791942e-15
RMSD copie pliée : 1.374813
coeur (30 atomes) : global=1.345387 local=0.000000
queue (16 atomes) : global=1.428352 local=0.000000
Bio.PDB.Superimposer : 1.374813

La copie rigide donne un RMSD de l'ordre de $10^{-15}$ Å, soit du bruit de virgule flottante : la rotation et la translation sont retrouvées à la précision machine. Ce contrôle valide l'implémentation. La copie pliée donne un RMSD global de 1,3748 Å. Pourtant, le cœur et la queue ont chacun un RMSD local nul : pris séparément, chaque bloc est parfaitement rigide. Le RMSD global ne mesure donc pas une déformation interne ; il mesure le compromis entre deux blocs séparés par une articulation.

L'écart résiduel global du cœur (1,3454 Å) et celui de la queue (1,4284 Å) sont du même ordre. La transformation commune répartit l'erreur entre les deux domaines au lieu de choisir l'un ou l'autre. Bio.PDB.Superimposer retrouve exactement le même RMSD global, ce qui confirme la cohérence des deux implémentations. Biopython utilise simplement une convention de matrice transposée pour rotran, sans effet sur la valeur du RMSD.

L'interprétation biologique doit rester prudente. La déformation est synthétique : elle ne prouve pas que la crambine possède une charnière à cet endroit. Elle montre en revanche qu'un RMSD global unique peut masquer une flexibilité localisée. Pour une analyse réelle, on calculerait les RMSD par résidu, on identifierait les zones d'écart résiduel élevé, on testerait plusieurs découpages en domaines et on confronterait le résultat aux facteurs B, à la présence de ligands et à la résolution des structures. Un RMSD de 1,37 Å sur 46 résidus n'a pas la même signification que le même chiffre sur 500 résidus.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-4/4.3-superposition-rmsd -q

python -m pytest exercises/module-4/4.3-superposition-rmsd -q

python -m pytest exercises/module-4/4.3-superposition-rmsd/projet -q

Résumé

• Une superposition compare des atomes homologues, pas deux listes dans un ordre arbitraire.
• Une transformation rigide combine une rotation propre et une translation ; une réflexion est exclue.
• Le centrage élimine la translation et réduit le problème à l'orientation relative.
• Kabsch maximise la trace $R H$ grâce à la SVD de la matrice de covariance croisée.
• Le déterminant de la rotation doit valoir $+1$ ; sinon, on corrige la décomposition.
• Le RMSD global résume l'accord d'ensemble, mais peut compenser des écarts locaux.
• Les RMSD locaux après ajustement global et après ajustement indépendant répondent à des questions différentes.
• Un RMSD sans sélection d'atomes, sans $N$ et sans analyse locale est difficile à interpréter.

Références

• Kabsch, W. (1976). A solution for the best rotation to relate two sets of vectors. Acta Crystallographica Section A, 32(5), 922–923.
• Theobald, D. L. (2005). Rapid calculation of RMSDs using a quaternion-based characteristic polynomial. Acta Crystallographica Section A, 61(4), 478–480.
• 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.PDB.Superimposer documentation, version 1.88. https://biopython.org/docs/1.88/api/Bio.PDB.Superimposer.html
• Biopython. Bio.PDB.PDBParser documentation, version 1.88. https://biopython.org/docs/1.88/api/Bio.PDB.PDBParser.html