Bio.PDB : hiérarchie et géométrie

Prérequis

Cette section prolonge 4.1, consacrée au format PDB et à ses colonnes fixes. Elle suppose que vous savez déjà ouvrir un fichier texte, lire une ligne ATOM et repérer les champs de coordonnées. Les objets de 1.3 vous sont familiers, notamment la distinction entre un conteneur Biopython et une chaîne Python. Le module 4 fournit le fil rouge structural : la crambine 1CRN et l'hémoglobine 4HHB. Biopython 1.88 et NumPy 2.x sont installés.

Introduction

La crambine est une petite protéine de réserve de Crambe hispanica, longue de 46 résidus et riche en cystéines. Sa structure cristallographique, déposée sous le code 1CRN, a été déterminée à 1,5 Å et contient 327 atomes de protéine. Trois ponts disulfure relient les six cystéines et stabilisent un repliement compact, hydrophobe et remarquablement résistant. La position des atomes et la géométrie des ponts expliquent une partie de cette stabilité.

Pour le calcul, la difficulté n'est pas de lire trois nombres par atome, mais de conserver la structure qui les organise. Un fichier PDB empile des enregistrements ATOM et HETATM, plusieurs modèles pour une structure RMN, plusieurs chaînes, des résidus numérotés avec des trous et parfois des codes d'insertion. Une liste plate de coordonnées perd les appariements atome-résidu-chaîne et rend les sélections fausses.

Biopython expose cette organisation sous forme d'arbre : Structure, Model, Chain, Residue, Atom. Bio.PDB fournit les parseurs, les sélections et un index spatial, NeighborSearch, pour éviter de comparer tous les atomes deux à deux. Cette section apprend à parcourir cet arbre sans supposer que les identifiants sont des indices, à calculer des distances et des contacts, puis à détecter des ponts disulfure. L'étude de cas mesure ces grandeurs sur 1CRN et discute ce que l'absence d'un pont signifierait.

Théorie

Un arbre, pas un tableau

Un fichier PDB se lit ligne par ligne, mais une structure biologique n'est pas une liste. Biopython reprend le modèle SMCRA (Structure-Model-Chain-Residue-Atom) : une Structure contient des Model, un modèle contient des Chain, une chaîne contient des Residue, un résidu contient des Atom. Chaque objet possède un identifiant, des enfants et un lien vers son parent. En remontant la chaîne des parents avec get_parent(), un atome de carbone alpha retrouve son résidu, sa chaîne et son modèle.

Cette hiérarchie n'est pas une commodité d'affichage. Elle porte les questions biologiques : un contact entre deux atomes n'a de sens que si l'on sait à quels résidus et à quelles chaînes ils appartiennent ; un pont disulfure ne relie pas deux atomes isolés mais deux cystéines. Elle porte aussi les pièges : deux chaînes peuvent numéroter leurs résidus de 1 à 100, et deux modèles peuvent décrire la même chaîne avec des coordonnées différentes.

Un Structure peut contenir plusieurs Model. Une structure cristallographique n'en a généralement qu'un ; une structure RMN en compte souvent vingt, qui représentent un ensemble de conformations compatibles avec les données. Choisir structure[0] est un acte d'analyse, pas une formalité. Pour un calcul de contacts, on précise toujours le modèle utilisé.

Parcourir et sélectionner

Deux styles coexistent. L'indexation directe descend l'arbre : structure[0]["A"][3]["CA"] désigne l'atome CA du résidu 3 de la chaîne A du premier modèle. Elle est lisible lorsque les identifiants sont simples, mais elle échoue si la chaîne, le résidu ou l'atome n'existe pas, et elle ne dit rien des autres atomes.

L'itération est plus robuste : for chain in model, puis chain.get_residues() et residue.get_atoms(). Structure.get_chains() traverse tous les modèles ; pour un décompte par modèle, partez de model. L'ordre d'itération suit l'ordre du fichier, pas un tri garanti : triez explicitement si votre rapport doit être reproductible.

L'indexation d'une chaîne par un entier, chain[3], vise le résidu polymérique sans code d'insertion. Dès qu'un code d'insertion existe, utilisez l'identifiant complet, par exemple chain[(" ", 100, "A")]. Cette différence est une source classique de résidus manqués dans une sélection. Pour une sélection réutilisable, Bio.PDB.Selection.unfold_entities(model, "A") aplatit l'arbre en liste d'atomes, que NeighborSearch accepte directement.

PDBParser, MMCIFParser et identifiants

PDBParser(QUIET=True).get_structure("1crn", path) renvoie une Structure. QUIET=True masque les avertissements de construction, pas les erreurs de fichier absent. Le paramètre PERMISSIVE vaut True par défaut dans Biopython 1.88 : un enregistrement mal formé peut être ignoré et produire une structure incomplète sans exception. Pour un contrôle strict, passez PERMISSIVE=False. L'en-tête est disponible dans structure.header : nom, méthode expérimentale, résolution, date de dépôt.

MMCIFParser(QUIET=True).get_structure("1crn", path) lit le format PDBx/mmCIF, devenu le format d'archivage principal. Ses options auth_chains=True et auth_residues=True sont actives par défaut : elles exposent les identifiants d'auteur (auth_asym_id, auth_seq_id) plutôt que les identifiants canoniques (label_asym_id, label_seq_id). Les deux jeux peuvent différer, et label_seq_id est strictement croissant, sans code d'insertion. Un script qui compare des résidus entre PDB et mmCIF doit donc documenter le jeu choisi.

Le format PDB historique code la chaîne sur un caractère, le numéro de résidu sur quatre chiffres et le numéro d'atome sur cinq chiffres. Il ne représente pas commodément les très grands complexes, les longues chaînes et certains identifiants. mmCIF lève ces limites et sépare les identifiants de l'auteur de ceux de l'archive.

L'identifiant d'un résidu est un triplet (hetflag, resseq, icode). Pour un acide aminé polymérique, hetflag vaut " ". Pour un hétéro-résidu, il vaut "H_" suivi du nom, par exemple "H_SO4" ; pour une molécule d'eau, "W". resseq est le numéro de l'auteur, qui peut commencer à un autre nombre que 1, comporter des trous et se répéter d'une chaîne à l'autre. icode est le code d'insertion, souvent vide, parfois A, B ou C dans les anticorps. Un résidu est donc identifié par la chaîne et par ce triplet, jamais par son rang dans l'itération.

Les atomes possèdent aussi un nom (CA, N, O, SG) et parfois un indicateur de conformation alternative (altloc). Deux atomes de même nom peuvent coexister dans un résidu désordonné ; Biopython en sélectionne un, mais get_altloc() et l'occupation permettent de vérifier lequel. Un atome nommé CA peut être un carbone alpha ou un calcium : le nom du résidu et l'élément lèvent l'ambiguïté, pas la position dans la ligne.

Distance, contacts et index spatial

La distance euclidienne entre deux atomes s'écrit

Elle se calcule directement avec NumPy, ou avec l'opérateur de soustraction d'un Atom Biopython, qui renvoie déjà une distance. Le calcul d'une paire est en $O(1)$ ; le calcul de toutes les paires est en $O(N^2)$, ce qui devient coûteux pour un complexe de plusieurs dizaines de milliers d'atomes.

NeighborSearch construit une grille tridimensionnelle et ne compare que les atomes de cellules voisines. La construction est linéaire en nombre d'atomes, et search_all(rayon) renvoie les paires réellement proches. Pour un rayon de 4 Å, on obtient les contacts atomiques ; il reste à décider ce qu'est un contact entre résidus. Une définition reproductible précise le rayon, les atomes retenus (lourds ou tous), l'exclusion des paires d'un même résidu et l'exclusion des voisins immédiats de la chaîne.

Sans exclusion, les contacts $i, i+1$ et $i, i+2$ dominent le décompte et n'apportent presque rien : ils sont imposés par la géométrie du squelette. Pour 1CRN, search_all(4.0) renvoie 1 960 paires d'atomes, qui correspondent à 172 paires de résidus ; en retirant les résidus séparés de moins de trois positions dans une même chaîne, il reste 90 paires. Ce filtrage doit être écrit dans la légende de tout tableau de contacts.

Ponts disulfure

Un pont disulfure est une liaison covalente entre les atomes de soufre de deux cystéines. La distance $d(\mathrm{SG}_1,\mathrm{SG}_2)$ vaut environ 2,05 Å. En pratique, on détecte les paires de cystéines dont les atomes SG sont à moins de 2,5 Å, puis on vérifie la cohérence avec les enregistrements SSBOND et CONECT du fichier.

Le seuil est un filtre, pas une définition. Une distance de 2,4 Å peut être un pont contraint ou une simple proximité ; une distance de 2,6 Å peut être un pont mal affiné ou une conformation alternative. La présence d'un atome SG manquant, d'un altloc ou d'un modèle différent peut faire disparaître un pont réel. Inversement, deux cystéines proches sans liaison covalente ne doivent pas être annotées comme pont disulfure sur la seule foi d'un seuil.

Liaisons hydrogène

Une liaison hydrogène relie un donneur $D$ portant un hydrogène à un accepteur $A$. On retient souvent une distance $d(D,A) \le 3{,}5$ Å et un angle $\theta = \angle(D-H\cdots A) \ge 120^\circ$, parfois $150^\circ$ selon les conventions. La distance seule est un critère nécessaire mais insuffisant : dans une structure à 1,5 Å, les hydrogènes ne sont généralement pas observés, et un couple N···O proche peut être mal orienté.

Sur 1CRN, une recherche de couples N···O de la chaîne principale, à moins de 3,5 Å et séparés d'au moins trois résidus, trouve 41 candidats. Le couple N11-O7 mesure 2,894 Å ; l'angle Cα-N···O, pris comme proxy de l'angle au donneur, vaut 121,2°. C'est un candidat plausible, pas une preuve : l'angle réel utilise l'hydrogène, et un couple à 2,884 Å peut avoir un proxy de 114,0° et échouer au critère d'angle. Pour une annotation de production, on utilise un programme spécialisé comme DSSP ou un serveur qui reconstruit les hydrogènes et traite les bifurcations.

Modèles, désordre et résidus manquants

len(structure) donne le nombre de modèles. Pour une structure RMN, analysez chaque modèle séparément ou choisissez-en un en le justifiant ; ne mélangez jamais les coordonnées de deux modèles. Pour une structure cristallographique, il peut exister des conformations alternatives dans un même modèle : atom.is_disordered() et atom.get_altloc() décrivent la situation.

Les résidus manquants n'existent pas dans l'arbre. Biopython ne fabrique pas de coordonnées pour un résidu absent du fichier. Une chaîne de 46 résidus peut donc être décrite par 44 résidus observés, avec un trou dans la numérotation. Le fichier PDB liste souvent ces absences dans REMARK 465 ; mmCIF possède des catégories dédiées. Comparez toujours le nombre de résidus observés au nombre annoncé par SEQRES ou par l'entité mmCIF avant de conclure à une protéine tronquée. Pour 1CRN, les 46 résidus sont présents et la numérotation est continue.

Exemple travaillé

Le fichier 1CRN contient un modèle, une chaîne A, 46 résidus et 327 atomes. Les 46 carbones alpha sont présents. La distance entre les CA des résidus 1 et 2 vaut 3,794 Å ; entre les CA des résidus 3 et 40, elle vaut 5,687 Å. Les distances CA-CA consécutives varient de 3,762 à 3,863 Å, avec une moyenne de 3,816 Å, ce qui est cohérent avec une chaîne polypeptidique repliée.

Les six cystéines sont aux résidus 3, 4, 16, 26, 32 et 40. Les trois paires de SG les plus proches mesurent 2,004 Å, 2,035 Å et 2,051 Å. Elles correspondent aux ponts 3-40, 4-32 et 16-26. Toutes les cystéines sont engagées : il n'y a pas de cystéine libre dans ce modèle. Nous retrouverons ces valeurs dans l'étude de cas.

Implémentation

Les sorties ci-dessous ont été exécutées le 2026-11-01 avec Biopython 1.88, NumPy 2.3.1 et Python 3.13. Le fichier est public/pdb/1CRN.pdb.

Descendre l'arbre

from pathlib import Path
 
from Bio.PDB import PDBParser
 
path = Path("public/pdb/1CRN.pdb")
structure = PDBParser(QUIET=True).get_structure("1crn", str(path))
model = structure[0]
 
print(structure.id, len(structure), len(list(structure.get_chains())))
print(len(list(model.get_residues())), len(list(model.get_atoms())))
residue = model["A"][3]
print(residue.id, residue.resname, residue["SG"].coord)

1crn 1 1
46 327
(' ', 3, ' ') CYS [15.885 12.426 10.016]

Le triplet (' ', 3, ' ') se lit : résidu polymérique, numéro 3, pas de code d'insertion. residue["SG"].coord est un tableau NumPy de trois nombres ; il ne faut pas le confondre avec l'identifiant du résidu.

Extraire les CA et mesurer

import numpy as np
 
 
def atom_distance(atom1, atom2):
    return float(np.linalg.norm(atom1.coord - atom2.coord))
 
 
ca = [residue["CA"] for residue in model["A"] if "CA" in residue]
print(len(ca))
print(round(atom_distance(ca[0], ca[1]), 3))
print(round(atom_distance(ca[2], ca[39]), 3))
consecutive = [atom_distance(ca[i], ca[i + 1]) for i in range(len(ca) - 1)]
print(
    round(min(consecutive), 3),
    round(max(consecutive), 3),
    round(sum(consecutive) / len(consecutive), 3),
)

46
3.794
5.687
3.762 3.863 3.816

La liste ca est construite par compréhension : l'ordre est celui de la chaîne, et l'indice 39 correspond au 40e résidu. Cette correspondance n'est valable que parce que la numérotation de 1CRN est continue. Dans un fichier avec des trous, indexez par identifiant, pas par position.

Détecter les ponts disulfure

from itertools import combinations
 
sg_atoms = []
for residue in model["A"]:
    if residue.resname == "CYS" and "SG" in residue:
        sg_atoms.append((residue.id[1], residue["SG"]))
 
for (seq1, atom1), (seq2, atom2) in combinations(sg_atoms, 2):
    distance = atom_distance(atom1, atom2)
    if distance <= 2.5:
        print(f"CYS {seq1:2d} - CYS {seq2:2d} : {distance:.2f} A")

CYS  3 - CYS 40 : 2.00 A
CYS  4 - CYS 32 : 2.03 A
CYS 16 - CYS 26 : 2.05 A

Le test <= 2.5 inclut les trois paires et aucune autre. Le format :2d aligne les numéros ; :.2f arrondit à deux décimales. Pour un rapport, conservez la valeur non arrondie et indiquez le seuil.

Chercher les contacts avec NeighborSearch

from Bio.PDB.NeighborSearch import NeighborSearch
 
atoms = list(model.get_atoms())
search = NeighborSearch(atoms)
pairs = search.search_all(4.0, level="A")
print(len(pairs))
 
best = {}
for atom1, atom2 in pairs:
    residue1 = atom1.get_parent()
    residue2 = atom2.get_parent()
    if residue1 is residue2:
        continue
    if residue1.get_parent().id == residue2.get_parent().id:
        if abs(residue1.id[1] - residue2.id[1]) < 3:
            continue
    key = tuple(sorted([residue1.id[1], residue2.id[1]]))
    distance = atom_distance(atom1, atom2)
    best[key] = min(best.get(key, distance), distance)
 
print(len(best))
for (seq1, seq2), distance in sorted(best.items(), key=lambda item: item[1])[:3]:
    print(seq1, seq2, round(distance, 3))

1960
90
3 40 2.004
4 32 2.035
16 26 2.051

search_all travaille au niveau des atomes (level="A"), puis le dictionnaire agrège par paire de résidus. Le filtre abs(...) < 3 retire les voisins immédiats d'une même chaîne. Sans lui, le décompte serait de 172 paires, gonflé par les contacts de squelette. Ce bloc suppose une chaîne unique ; pour plusieurs chaînes, incluez l'identifiant de chaîne dans la clé.

Un candidat de liaison hydrogène

from Bio.PDB.vectors import calc_angle
 
donor = model["A"][11]
acceptor = model["A"][7]
distance = atom_distance(donor["N"], acceptor["O"])
angle = calc_angle(
    donor["CA"].get_vector(),
    donor["N"].get_vector(),
    acceptor["O"].get_vector(),
)
print(round(distance, 3), round(float(np.degrees(angle)), 1))

2.894 121.2

L'angle est calculé au niveau de l'azote, entre le CA et l'oxygène accepteur. Ce n'est pas l'angle $D-H\cdots A$ : l'hydrogène manque. Le résultat sert à filtrer des candidats, pas à certifier une liaison hydrogène.

Étude de cas

La question est simple : les six cystéines de la crambine forment-elles trois ponts disulfure, et que signifierait l'absence de l'un d'eux ? Le fichier local ne contient que la protéine, sans eau ; cette absence de coordonnées interdit toute conclusion sur l'hydratation.

Première étape : inventorier la structure. Un modèle, une chaîne A, 46 résidus, 327 atomes et six atomes SG. La numérotation va de 1 à 46 sans trou. Les 46 CA sont présents, ce qui permet de calculer les distances de squelette sans traitement des manquants.

Deuxième étape : calculer toutes les distances SG-SG. Les trois paires engagées mesurent 2,004 Å, 2,035 Å et 2,051 Å. Le fichier confirme ces ponts par trois enregistrements SSBOND (Cys 3-Cys 40, Cys 4-Cys 32, Cys 16-Cys 26) et par les enregistrements CONECT correspondants. Détection géométrique et annotation sont cohérentes.

Troisième étape : tester la sensibilité au seuil. Avec 2,00 Å, aucune paire n'est retenue, car la plus courte mesure 2,004 Å. Avec 2,02 Å, une seule ; avec 2,04 Å, deux ; avec 2,06 Å, trois. Un seuil « raisonnable » de 2,5 Å fonctionne ici, mais cette expérience montre qu'un seuil trop strict produit des faux négatifs. À l'inverse, un seuil de 3,0 Å pourrait inclure des contacts soufrés non covalents.

Quatrième étape : discuter la présence et l'absence. Les six cystéines sont appariées, sans cystéine libre. Si une cystéine n'avait aucun partenaire à moins de 2,5 Å, trois explications seraient possibles : un thiol réellement libre, un atome SG absent ou désordonné, ou un pont réel mais éloigné du seuil. L'annotation SSBOND et la connectivité CONECT permettent de trancher, pas la distance seule. Dans 1CRN, aucune ambiguïté : les trois ponts sont annotés et géométriquement nets.

Cinquième étape : replacer les ponts dans les contacts. À 4 Å, les trois paires de cystéines figurent parmi les contacts les plus courts, avec 6, 6 et 8 paires d'atomes. Le pont 3-40 est même le contact non local le plus court du modèle. Une liaison covalente S-S est plus courte qu'un contact de van der Waals.

L'interprétation doit rester prudente : la structure montre trois ponts, pas leur nécessité. Une mutation, une réduction chimique ou une expression en conditions réductrices sont des expériences différentes. Le calcul fournit une hypothèse géométrique, pas une mesure de stabilité.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-4/4.2-biopython-pdb -q

python -m pytest exercises/module-4/4.2-biopython-pdb -q

python -m pytest exercises/module-4/4.2-biopython-pdb/projet -q

Résumé

• Bio.PDB organise une entrée en Structure > Model > Chain > Residue > Atom ; chaque objet connaît ses enfants et son parent.
• PDBParser et MMCIFParser produisent la même hiérarchie, mais mmCIF sépare identifiants d'auteur (auth_*) et identifiants canoniques (label_*).
• Un résidu est identifié par (hetflag, resseq, icode) : le numéro de l'auteur n'est ni un rang ni un index de tableau.
• La distance euclidienne se calcule en $O(1)$ par paire ; NeighborSearch évite le coût $O(N^2)$ d'une comparaison exhaustive.
• Un contact doit être défini avant d'être compté : rayon, atomes lourds, exclusion des voisins immédiats et des mêmes résidus.
• Un pont disulfure se détecte par une distance SG-SG d'environ 2,05 Å, puis se confirme avec SSBOND et CONECT.
• Une liaison hydrogène exige une distance donneur-accepteur et un angle ; sans hydrogène, le critère d'angle reste un proxy.
• Les modèles multiples, les conformations alternatives et les résidus manquants doivent être diagnostiqués avant toute conclusion biologique.

Références

• Berman, H. M. et al. (2000). The Protein Data Bank. Nucleic Acids Research, 28(1), 235–242.
• Cock, P. J. A. et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11), 1422–1423.
• RCSB PDB. Entry 1CRN : Water structure of a hydrophobic protein at atomic resolution. https://www.rcsb.org/structure/1CRN
• Biopython. Bio.PDB : Structure, Model, Chain, Residue, Atom et documentation de Bio.PDB. https://biopython.org/docs/latest/api/Bio.PDB.html
• wwPDB. PDBx/mmCIF format et spécifications. https://mmcif.wwpdb.org/