Matrices PAM et BLOSUM

Prérequis

Cette section prolonge l'alignement global de Needleman-Wunsch et l'alignement local de Smith-Waterman : vous devez savoir ce qu'est une colonne d'alignement, un appariement et un gap. Les séquences se lisent comme en 1.1, et les matrices seront chargées avec l'API de Biopython présentée en 1.3. Aucune connaissance préalable des matrices de substitution n'est nécessaire ; nous partons du problème biologique pour arriver aux formules.

Introduction

Deux protéines homologues ne diffèrent pas par des substitutions choisies au hasard. Les remplacements entre résidus de taille et de propriétés chimiques voisines, comme l'isoleucine et la valine, ou l'acide aspartique et l'acide glutamique, sont plus fréquents que les remplacements entre une cystéine et un tryptophane. Un alignement qui traite toutes les substitutions de la même manière ignore cette information. Il peut préférer une série de remplacements improbables à une série de remplacements conservatifs, simplement parce que le nombre de colonnes identiques est légèrement plus grand. Les matrices de substitution corrigent ce biais en attribuant à chaque paire de résidus un score fondé sur les substitutions observées dans des protéines apparentées.

L'enjeu computationnel est de transformer ces observations en un modèle additif utilisable par la programmation dynamique. Un score d'alignement est la somme de scores de paires et de pénalités de gap ; si la matrice est mal choisie, l'optimum change, et avec lui l'alignement, la position des domaines et parfois la conclusion biologique. Les deux familles historiques, PAM et BLOSUM, répondent à la même question avec des constructions différentes. PAM extrapole un modèle de mutation à partir de protéines très proches ; BLOSUM compte directement les paires observées dans des blocs conservés, à un seuil d'identité donné. Comprendre cette différence évite de traiter une matrice comme une table de constantes universelles.

Théorie

Pourquoi un score log-odds

Notons $q_a$ la fréquence de fond du résidu $a$ et $p_{ab}$ la probabilité d'observer la paire $(a,b)$ dans une colonne d'alignement homologue. Sous indépendance, cette probabilité serait $q_a q_b$. Le rapport $p_{ab}/(q_a q_b)$ mesure donc la surreprésentation de la paire : supérieur à 1, il signale une paire favorisée ; inférieur à 1, une paire évitée. Le score log-odds est le logarithme en base 2 de ce rapport :

$S(a,b) = \log_2 \frac{p_{ab}}{q_a q_b}$.

Un score positif signifie que la paire est plus fréquente que sous indépendance, un score nul qu'elle est aussi fréquente, un score négatif qu'elle est moins fréquente. La base 2 exprime le score en bits : chaque unité correspond à un doublement du rapport de vraisemblance. Le logarithme transforme un produit de probabilités en une somme, si bien que le score d'un alignement est la somme des scores de ses colonnes. C'est ce modèle additif que la programmation dynamique optimise.

La matrice n'est ni une distance ni une probabilité : c'est une table de rapports de vraisemblance arrondis. Les valeurs dépendent de l'échantillon, des fréquences de fond et du seuil de divergence ; deux matrices ne décrivent pas le même modèle, et leurs scores bruts ne sont pas comparables.

PAM : extrapoler depuis les protéines proches

La première famille doit son nom à l'unité de Dayhoff : la mutation ponctuelle acceptée (PAM, accepted point mutation). Dayhoff et ses collaborateurs ont comparé des protéines très proches, souvent à plus de 85 % d'identité, et reconstruit les substitutions acceptées le long d'arbres phylogénétiques. Une PAM correspond à une substitution acceptée pour 100 résidus. PAM1 décrit donc un régime très proche, proche de 1 % de divergence, et non un alignement où 1 % des positions seraient observablement différentes.

Pour atteindre une divergence plus grande, Dayhoff extrapole PAM1 par un modèle de Markov : on élève la matrice de transition à une puissance. PAM250 vaut PAM1 puissance 250 et décrit un régime de divergence élevé, de l'ordre de 20 % d'identité résiduelle. Ce n'est pas « 250 % de divergence » : 250 est un nombre de mutations acceptées par 100 résidus, pas un pourcentage, et les mutations peuvent se superposer aux mêmes sites. L'extrapolation suppose des taux stationnaires, indépendants des sites et identiques entre lignées, hypothèses simplificatrices pour des protéines biaisées ou sélectionnées.

La numérotation PAM est inversée par rapport à l'intuition : plus le nombre est grand, plus la matrice décrit des séquences éloignées. PAM1 convient aux protéines très proches, PAM250 aux protéines lointaines. Les erreurs d'estimation de PAM1 sont amplifiées par les puissances successives. Les matrices PAM restent utiles pour raisonner en distance évolutive et pour certaines similarités lointaines, mais elles ne sont plus le défaut des aligneurs modernes.

BLOSUM : compter dans les blocs conservés

Les matrices BLOSUM (BLOcks SUbstitution Matrix) ont été proposées par Henikoff et Henikoff en 1992. Au lieu d'extrapoler un modèle à partir de protéines proches, elles comptent les paires de résidus observées dans des blocs conservés de familles protéiques : des régions alignées sans gap, jugées homologues. Pour éviter qu'une famille très représentée domine, les séquences d'un bloc sont regroupées lorsque leur identité dépasse un seuil. BLOSUM62 signifie que les séquences ont été regroupées à 62 % d'identité ou plus : les paires entre groupes distincts sont comptées, celles à l'intérieur d'un même groupe ne le sont pas. Le nombre 62 est un paramètre de construction, pas un seuil à atteindre pour qu'un alignement soit valide.

Le seuil agit sur la sensibilité. BLOSUM80 regroupe seulement les séquences très similaires et conserve beaucoup de paires proches ; la matrice est adaptée aux homologues proches. BLOSUM45 regroupe davantage de séquences et compte des paires plus divergentes ; la matrice est plus permissive et adaptée aux homologies lointaines. La numérotation BLOSUM est donc directe : plus le nombre est grand, plus la matrice vise des séquences proches. C'est l'inverse de PAM, une source de confusion fréquente.

La construction BLOSUM ne suppose pas d'arbre phylogénétique et n'extrapole pas par une puissance. Elle estime directement les fréquences de paires au seuil choisi, puis calcule les rapports log-odds et les arrondit. En contrepartie, les blocs conservés ne représentent pas toute la diversité des protéines : les régions variables, les protéines de faible complexité et les familles peu documentées y sont sous-représentées. BLOSUM62 est devenue le défaut de BLASTP parce qu'elle offre un bon compromis pour une divergence moyenne, pas parce qu'elle serait optimale partout.

Deux constructions, deux usages

Le tableau suivant résume les différences qui ont des conséquences pratiques.

La construction explique l'usage. PAM1 est estimée sur des alignements fiables puis extrapolée : elle convient à l'étude de la divergence et aux comparaisons évolutives. BLOSUM est estimée directement à un régime donné : elle convient à la recherche de similarité, où l'on veut une matrice calibrée pour un seuil d'identité observable. PAM250 pour des protéines proches ou BLOSUM80 pour des protéines très éloignées n'est pas interdit, mais ce n'est pas l'usage prévu.

Choisir selon la divergence et la longueur

Il n'existe pas de règle universelle, mais des règles explicites aident à raisonner. Pour des protéines très proches, BLOSUM80 préserve le signal des identités et des substitutions conservatrices. Pour une divergence moyenne, BLOSUM62 est un bon compromis et reste le défaut raisonnable. Pour des protéines très éloignées, BLOSUM45 ou PAM250 peut récupérer des similarités que BLOSUM62 pénalise trop, au prix d'un risque accru de faux positifs. Les seuils de ce cours sont : au moins 85 % d'identité, BLOSUM80 ; de 30 % à moins de 85 %, BLOSUM62 ; moins de 30 %, PAM250. Ce ne sont pas des lois biologiques, mais des conventions explicites et testables.

La longueur modifie l'interprétation. Un alignement court est dominé par quelques paires et sa variance est grande : un segment de 20 résidus peut scorer haut par hasard, surtout s'il est riche en résidus fréquents ou de faible complexité. Un alignement long accumule davantage de preuves, mais peut contenir des régions mal alignées. La E-value dépend de la longueur de la requête, de la taille de la base et du nombre d'alignements testés, pas seulement du score brut ; pour deux séquences, une permutation ou un modèle nul reste nécessaire.

Les pénalités de gap sont associées à la matrice. Le modèle affine s'écrit $g(k)=g_o+k g_e$, avec $g_o$ l'ouverture et $g_e$ l'extension. Ouvrir coûte cher, allonger coûte moins cher : cela favorise quelques gaps longs plutôt que de nombreux gaps d'un résidu. Les valeurs classiques de BLASTP sont 11/1 avec BLOSUM62 et plus élevées avec PAM250, mais les défauts varient selon l'outil et la version. Une pénalité trop faible fragmente l'alignement et peut relier des domaines non homologues ; une pénalité trop forte rate des insertions réelles. Matrice, gaps et filtrage forment un ensemble.

Faux sentiment de précision, normalisation et biais de composition

Les scores matriciels sont des entiers, ce qui donne une impression de précision. Ils sont pourtant issus d'estimations de fréquences, d'arrondis et d'un échantillon de protéines. Une différence d'un point n'est souvent pas significative. Les matrices dépendent des fréquences de fond, de la base utilisée et des hypothèses de construction : les présenter comme une vérité biologique revient à confondre un modèle avec la réalité.

Les scores bruts de matrices différentes ne sont pas comparables. BLOSUM62 et PAM250 n'ont pas la même échelle, les mêmes fréquences de fond ni le même arrondi ; un score de 95 avec l'une ne vaut pas 95 avec l'autre. Les aligneurs convertissent le score brut en bit score, puis en E-value. Le bit score normalise le score par les paramètres de la matrice, et la E-value est le nombre attendu d'alignements d'un score au moins égal dans une base de taille donnée. C'est cette normalisation qui permet de comparer recherches et matrices.

Le biais de composition est plus subtil. Les fréquences $q_a$ sont estimées sur un ensemble de protéines ; si la requête est riche en glycine, en proline ou en un résidu rare, le modèle d'indépendance est mal calibré. Les régions de faible complexité produisent des scores élevés pour des raisons de composition, pas d'homologie. Les matrices construites sur des protéines globulaires s'appliquent mal aux protéines membranaires, intrinsèquement désordonnées ou répétées. Des filtres et des statistiques de composition corrigent partiellement ce biais, sans remplacer l'inspection des alignements.

Exemple travaillé à la main

Prenons un modèle jouet à quatre résidus, A, C, D et E, avec les fréquences de fond $q_A = 0{,}30$, $q_C = 0{,}20$, $q_D = 0{,}30$ et $q_E = 0{,}20$. Les probabilités de paires sont choisies de sorte que les marges redonnent ces fréquences : $p_{AA}=0{,}12$, $p_{AC}=0{,}04$, $p_{AD}=0{,}10$, $p_{AE}=0{,}04$, $p_{CC}=0{,}12$, $p_{CD}=0{,}02$, $p_{CE}=0{,}02$, $p_{DD}=0{,}12$, $p_{DE}=0{,}06$ et $p_{EE}=0{,}08$. Les paires hors diagonale sont symétriques ; le tableau complet, qui compte les deux orientations, somme à 1.

Calculons quelques scores. Pour A-A, le rapport vaut $0{,}12/(0{,}30 \times 0{,}30)=1{,}333$, donc $S(A,A)=\log_2(1{,}333)=+0{,}42$ bit, arrondi à 0. Pour A-C, le rapport vaut $0{,}04/(0{,}30 \times 0{,}20)=0{,}667$, donc $S(A,C)=\log_2(0{,}667)=-0{,}58$ bit, arrondi à $-1$. Pour C-C, $0{,}12/(0{,}20 \times 0{,}20)=3$, donc $S(C,C)=+1{,}58$ bit, arrondi à $+2$. Pour C-D, $0{,}02/(0{,}20 \times 0{,}30)=0{,}333$, donc $S(C,D)=-1{,}58$ bit, arrondi à $-2$. Pour D-E, $0{,}06/(0{,}30 \times 0{,}20)=1$, donc $S(D,E)=0$. Pour E-E, $0{,}08/(0{,}20 \times 0{,}20)=2$, donc $S(E,E)=+1$ bit, arrondi à $+1$.

La petite matrice obtenue est donc : A-A 0, A-C $-1$, A-D 0, A-E $-1$, C-C $+2$, C-D $-2$, C-E $-1$, D-D 0, D-E 0, E-E $+1$. Avec une pénalité de gap de $-2$, l'alignement AACC face à ACCC vaut $0 + (-1) + 2 + 2 = 3$. L'alignement ACDE face à ADDE vaut $0 + (-2) + 0 + 1 = -1$ : trois résidus identiques sur quatre ne suffisent pas à compenser le remplacement C-D, ce qui illustre l'effet d'une substitution défavorable. Un gap dans A-CD face à AACD coûte $-2$, mais les deux C et le D apportent $+2$ et 0, soit un total de 0. Ces valeurs seront vérifiées par le code ci-dessous.

Implémentation

Les blocs suivants ont été exécutés avec Python 3.13.5 et Biopython 1.88 le 2026-11-01. log_odds_score refuse les probabilités hors de $]0,1]$, car le logarithme n'y est pas défini. score_alignment normalise la casse et accepte une matrice symétrique : si $(a,b)$ est absente, il essaie $(b,a)$.

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

La sortie vérifiée est :

AA: +0.42 bits -> +0
AC: -0.58 bits -> -1
AD: +0.15 bits -> +0
AE: -0.58 bits -> -1
CC: +1.58 bits -> +2
CD: -1.58 bits -> -2
CE: -1.00 bits -> -1
DD: +0.42 bits -> +0
DE: +0.00 bits -> +0
EE: +1.00 bits -> +1

Le passage aux entiers utilise l'arrondi au plus proche. En Python, round arrondit les demi-entiers vers le pair, ce qui peut surprendre. Les scores d'une matrice réelle sont eux-mêmes arrondis. La petite matrice ci-dessous reprend ces valeurs et ajoute une pénalité de gap de $-2$.

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

La sortie vérifiée est :

AACC / ACCC : 3.0
ACDE / ADDE : -1.0
A-CD / AACD : 0.0

La règle de choix des exercices est volontairement simple. Elle ne remplace pas une évaluation sur données réelles ; elle rend explicite le lien entre identité observée et matrice. La fonction lève une erreur hors de $[0,100]$.

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

La sortie vérifiée est :

   95 % -> BLOSUM80
   85 % -> BLOSUM80
 84.9 % -> BLOSUM62
   75 % -> BLOSUM62
   30 % -> BLOSUM62
 29.9 % -> PAM250
   20 % -> PAM250

Enfin, Biopython fournit les matrices réelles. L'objet Array s'indexe par une paire de lettres et renvoie un flottant. La fonction ci-dessous ne gère pas les gaps : elle additionne les scores de colonnes appariées, ce qui suffit pour comparer deux matrices sur un alignement déjà construit.

from Bio.Align import substitution_matrices
 
blosum62 = substitution_matrices.load("BLOSUM62")
pam250 = substitution_matrices.load("PAM250")
 
def score_ungapped(seq1, seq2, matrix):
    return sum(float(matrix[a, b]) for a, b in zip(seq1, seq2))
 
print("BLOSUM62 A-A :", blosum62["A", "A"])
print("PAM250 A-A :", pam250["A", "A"])
print("BLOSUM62 D-W :", blosum62["D", "W"])
print("PAM250 D-W :", pam250["D", "W"])

La sortie vérifiée est :

BLOSUM62 A-A : 4.0
PAM250 A-A : 2.0
BLOSUM62 D-W : -4.0
PAM250 D-W : -7.0

Ces quatre valeurs montrent une différence de calibrage : PAM250 donne un score plus faible à l'identité A-A et plus négatif au remplacement D-W. Cela ne signifie pas que PAM250 est plus sévère en général : les échelles diffèrent, et il faut comparer des alignements complets puis normaliser.

Étude de cas

Nous utilisons un exemple synthétique. La référence est un peptide de 20 résidus contenant les vingt acides aminés standards : ACDEFGHIKLMNPQRSTVWY. Le peptide proche remplace quatre résidus par des résidus voisins (F→Y, I→V, K→R, T→S), soit 80 % d'identité. Le peptide éloigné conserve cinq positions, soit 25 %. Un témoin négatif inverse la référence et ne partage aucune position.

Quelle matrice utiliser pour chaque cas, et que dit le score ? Nous appliquons la règle de cours, puis calculons les scores avec BLOSUM62 et PAM250. Le code est autonome.

from Bio.Align import substitution_matrices
 
blosum62 = substitution_matrices.load("BLOSUM62")
pam250 = substitution_matrices.load("PAM250")
 
def score_ungapped(seq1, seq2, matrix):
    return sum(float(matrix[a, b]) for a, b in zip(seq1, seq2))
 
def identity_percent(seq1, seq2):
    return 100.0 * sum(a == b for a, b in zip(seq1, seq2)) / len(seq1)
 
seq1 = "ACDEFGHIKLMNPQRSTVWY"
seq_close = "ACDEYGHVRLMNPQRSSVWY"
seq_far = "AVDDYGKVRMLQAEKTTIWF"
seq_control = seq1[::-1]
 
for label, seq2 in (("proche", seq_close), ("eloignee", seq_far), ("temoin", seq_control)):
    identity = identity_percent(seq1, seq2)
    choice = "BLOSUM80" if identity >= 85 else "BLOSUM62" if identity >= 30 else "PAM250"
    sb = score_ungapped(seq1, seq2, blosum62)
    sp = score_ungapped(seq1, seq2, pam250)
    print(f"{label:9s} identite={identity:5.1f} %  choix={choice:8s}  BLOSUM62={sb:6.1f}  PAM250={sp:6.1f}")

La sortie vérifiée est :

proche    identite= 80.0 %  choix=BLOSUM62  BLOSUM62= 105.0  PAM250= 111.0
eloignee  identite= 25.0 %  choix=PAM250    BLOSUM62=  54.0  PAM250=  73.0
temoin    identite=  0.0 %  choix=PAM250    BLOSUM62= -26.0  PAM250= -26.0

Pour le peptide proche, la règle choisit BLOSUM62. Les scores totaux sont voisins, 105 pour BLOSUM62 et 111 pour PAM250, soit environ 5,25 et 5,55 par position : le score brut ne départage pas les matrices. BLOSUM62 est calibrée pour cette divergence moyenne et reste le défaut de BLASTP ; PAM250 n'apporte pas d'avantage identifiable ici. Pour le peptide éloigné, la règle choisit PAM250. BLOSUM62 donne 54 (2,70 par position) et PAM250 73 (3,65 par position) : PAM250 est plus permissive pour les substitutions observées, ce qui est cohérent avec son usage pour les homologies lointaines. Le témoin reçoit $-26$ avec les deux matrices, soit $-1,30$ par position ; les scores positifs ne sont donc pas un artefact d'échelle.

La matrice choisie par une règle simple n'est pas celle qui maximise le score brut : PAM250 gagne ici, mais BLOSUM62 reste préférable pour le peptide proche par sa calibration et sa spécificité. Les scores bruts de deux matrices ne sont pas comparables : 105 et 111 ne mesurent pas la même échelle, et 54 contre 73 ne signifie pas que PAM250 est deux fois meilleure. Un score positif ne suffit pas à conclure à une homologie : il faut une E-value, un modèle nul ou une permutation.

La décision dépend du coût d'erreur : rechercher un domaine lointain et un orthologue proche n'ont pas les mêmes exigences.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-2/2.5-matrices-pam-blosum -q

python -m pytest exercises/module-2/2.5-matrices-pam-blosum -q

python -m pytest exercises/module-2/2.5-matrices-pam-blosum/projet -q

Résumé

• Une matrice de substitution est une table de scores log-odds, pas une mesure directe de la biologie.
• PAM1 décrit une mutation acceptée pour 100 résidus ; les PAM élevées sont extrapolées, et leur numérotation croît avec la divergence.
• BLOSUM compte des paires dans des blocs conservés, avec un seuil de regroupement ; sa numérotation croît avec la proximité.
• BLOSUM62 est un bon défaut pour la divergence moyenne, pas une vérité universelle : longueur, gaps et composition modifient le score.
• Les scores bruts de matrices différentes ne sont pas comparables ; utilisez le bit score et la E-value.
• Un exemple court se vérifie à la main, puis se contrôle avec un test reproductible.

Références

• Dayhoff, M. O., Schwartz, R. M. et Orcutt, B. C. (1978). A model of evolutionary change in proteins. Atlas of Protein Sequence and Structure, 5(3), 345–352.
• Henikoff, S. et Henikoff, J. G. (1992). Amino acid substitution matrices from protein blocks. PNAS, 89(22), 10915–10919.
• Altschul, S. F. (1991). Amino acid substitution matrices from an information theoretic perspective. Journal of Molecular Biology, 219(3), 555–565.
• Altschul, S. F., Gish, W., Miller, W., Myers, E. W. et Lipman, D. J. (1990). Basic local alignment search tool. Journal of Molecular Biology, 215(3), 403–410.
• Biopython. Bio.Align.substitution_matrices, documentation 1.88. https://biopython.org/docs/1.88/api/Bio.Align.substitution_matrices.html