Il existe toute une série d’outils écrits en Python et destinés à lire des fichiers PDB, en mode filtre ou parseur dont le plus notable est celui du module Bio.PDB du projet Biopython. Ce dernier est généraliste, d’autres se focalisent sur des aspects particuliers comme extraire des métadonnées de l’entête PDB, calculer des propriétés structurales, statistiques, physicochimiques … en fait il y a de la place pour tous types de codes. En marge du développement du paquetage buildez.pdb je vais donner quelques exemples obtenus avec le parseur de Biopython, qui est stable, disponible, assez documenté.

D’après un article sur buildblog.buildez.net – Mis à jour en Avril 2025.

Nous prendrons comme exemple la structure 1ADO [Blom_1997] qui correspond à une Aldolase et dont la structure du fichier PDB a été examinée dans l’article [ Comprendre un fichier PDB ] sur ce blog. Il est recommandé de lire cet article avant de passer à la suite, rappelons simplement ici que la structure de cette enzyme se présente sous la forme d’une unité asymétrique comportant 4 chaines polypeptidiques qui sont étiquetées A, B, C, D et qui ligandent deux molécules de 13P (1,3-dihydroxyacetonephosphate, chaînes A et B) et deux ions sulfate (SO4, chaînes C et D), ainsi que des molécules d’eau (toutes les chaînes).

1. Parser l’entête PDB

Dans un premier temps nous allons nous intéresser à l’entête (header) du fichier PDB, avec le code suivant pour extraire différentes informations :

Nous chargeons le fichier PDB (1ado.pdb) dans une chaîne de caractères s que l’on décompose (en utilisant le caractère de fin de ligne en tant que séparateur) en lignes de texte, qui sont collectées dans une liste Python l. Cette forme va être acceptée par le parseur p qui va renvoyer un objet structure contenant les données.
La troisième partie consiste à exploiter les données de l’objet structure. Nous pouvons ainsi afficher différentes parties de l’entête (structure.header) pour peu que l’on connaisse les mots clés utilisés (voir le module parse_pdb_header.py du paquetage Bio.PDB) par ce qui est un dictionnaire Python.

Ce qui donnera le résultat suivant:

[data\1ado.pdb] found 17564 lines ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain A is discontinuous at line 14182. warnings.warn( ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain B is discontinuous at line 14206. warnings.warn( ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain C is discontinuous at line 14218. warnings.warn( ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain D is discontinuous at line 14223. warnings.warn( ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain A is discontinuous at line 14228. warnings.warn( ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain B is discontinuous at line 14990. warnings.warn( ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain C is discontinuous at line 15672. warnings.warn( ...StructureBuilder.py:100: PDBConstructionWarning: WARNING: Chain D is discontinuous at line 16625. warnings.warn( resolution = 1.9 name = fructose 1,6-bisphosphate aldolase from rabbit muscle head = lyase deposition_date = 1996-12-02 release_date = 1997-12-24 structure_reference = ['n.s.blom,j.sygusch enhanced electron density envelopes by extended solvedefinition to be published '] journal_reference = n.blom,j.sygusch product binding and role of the c-terminal region in d-fructose 1,6-bisphosphate aldolase. nat.struct.biol. v. 4 36 1997 issn 1072-8368 8989320 10.1038/nsb0197-36 author = N.S.Blom,J.Sygusch compound = {'1': {'misc': '', 'molecule': 'aldolase', 'chain': 'a, b, c, d', 'ec_number': '4.1.2.13', 'ec': '4.1.2.13'}}

Nous obtenons, après quelques warnings, le contenu recherché de l’entête du fichier PDB: résolution , nom, dates, référence de structure, référence de la publication (qui risque d’être interprétée et signalée si vous avez une extension type Click and Read dans votre navigateur) et des informations sur la structure avec la ligne structure.header['compound'].

Une évolution: des warnings

Dans les versions anciennes de Bio.PDB il n’y avait pas ces warnings, le parseur fonctionnait et perdait de l’information. Première remarque: il indique des lignes de la liste l au lieu d’interpréter les objets, ce qui aurait pu être plus parlant. Si nous examinons le premier warning  (Chain A is discontinuous at line 14182) en nous référant à la ligne signalée dans le fichier, nous avons :

Nous constatons qu’il s’agit du passage de la fin de la protéine (chaine D, TYR363, lignes ATOM) au début du premier ligand (composé 13P, lignes HETATM) de la chaine A. Pourtant il y avait une ligne TER, donc ce warning ne nous apprends pas grand chose. Les autres non plus car ils marquent les ‘frontières’ entre les blocs HETATM qui changent de chaine, pour les ligands : 13C, SO4 et H0H.

2. Détection de vraies exceptions (délétions)

Ceci dit, puisqu’il y a eu une évolution, il peut être intéressant de tester le parseur avec des chaines présentant des défauts objectifs: pseudo délétions, manque de chaines latérales … Nous allons prendre quelques protéines tests et résumer les résultats dans le tableau :

structure	délétions	résidus	warnings
1HXD [Weaver_2001]	2 délétions A: MET211 à GLN221, B: MET211 à GLN221	–	–
4JIK [Meng_2013]	3 délétions A: ALA19 à GLY21, A: VAL40 à ASN51, A: GLU76 à  ILE79	5 résidus incorrects [PRO4, ALA19, ASN51, GLU76, GLY272]	–
1OIM:A [Zechel_2003]	1 délétion A: LYS232 à GLU234	1 résidu incorrect [LYS232]	–

Donc visiblement il n’y a pas de warnings à l’exécution du code.  Prenons le cas de 4JIK qui comporte 3 délétions, si nous déroulons la molécule (voir les sections suivantes) nous pouvons avoir le nom de chaque résidu, par la propriété resname de l’objet Bio.PDB.Residue.Residue. Mais nous n’avons pas de propriété qui serait équivalente au numéro de résidu (par exemple RESSEQ). Donc nous n’avons pas de moyen simple pour détecter une continuité dans les résidus qui marquerait cette délétion au cours de la construction hiérarchique (itérateur) de la protéine.

Contournement

Pour faire une détection, nous devons utiliser le module Bio.PDB.Polypeptide avec un code tout simple qui va afficher la séquence de chaque polypeptide trouvé dans la chaine :

Ce qui va nous donner (j’ai ajouté une étoile en fin de polypeptide) :

Nous avons effectivement 3 étoiles donc 3 délétions pour une seule chaine. Les résidus correspondent bien, ce qui donne les moyens de détecter les problèmes. A condition de savoir combien de chaines nous avons, ce qui se compte en ajoutant des instructions dans l’itérateur, car nous allons avoir les polypeptides de chaque chaine sans délimiteurs de chaine.
Et avec cette méthode nous n’aurons que des numéros d’index (1 à n), et non de vrais numéros RESSEQ (la protéine ne commence pas forcément avec un résidu noté 01). C’est le cas avec 4JIK, qui commence avec une PRO4 :

Nous remarquons que cette PRO4 est incomplète, il manque plusieurs atomes à ce résidu, sans warning du parseur.

3. Et coté ligands ?

L’exemple montre ce qu’il est à peu près possible de réaliser avec ce type de parseur, qui est adapté plus à des biologistes structuraux qu’a des chimistes médicinaux. Dans le second cas, les ligands vont nous intéresser au moins autant que la structure, cette information est disponible dans l’entête sous différentes déclinaisons :

Mais elle n’est pas collectée par le parseur et stockée dans le dictionnaire final (par exemple compound['2'] renvoie une erreur).

4. Parser les coordonnées

Dans le cas des coordonnées moléculaires l’approche va être différente, il va falloir implémenter une boucle qui va utiliser un itérateur Python et faire remonter tous les éléments de la structure: atomes -> résidus -> chaînes -> modèles.
Mais d’abord nous allons définir deux fonctions get_atom_info et get_res_info destinées à alléger le code de l’exemple. Ces fonctions permettent de sélectionner les attributs des atomes et des résidus que l’on veut renvoyer vers une autre structure de données pendant l’exécution de la boucle.

Nous constatons que l’on retrouve des champs (colonnes) du format PDB qui sont classiques des blocs ATOM ou HETATM et qui sont distribués dans les attributs des résidus ou des atomes. La structure étant dans l’objet structure il suffit de dérouler la boucle suivante:

Pour les biochimistes ces boucles sont auto explicatives.
Nous notons que i le numéro de résidu est 363 et le nom de résidu TYR on va dérouler l’ensemble des atomes du résidu. Ces actions sont réalisées pour toutes les chaînes de la structure et éventuellement tous les modèles. Pour compter chaînes et modèles on utilise deux compteurs model_n et chain_n qui vont s’incrémenter au fur et à mesure, nous oublions RESSEQ. Le code précédent va donner le résultat (des lignes ont été masquées par ‘…’ pour garder un peu de visibilité) pour la première chaîne (il ne s’agit pas d’une structure multi modèles) :

Nous notons la valeur de residue.disordered qui passe de zéro à un dans le cas du résidu 13P : ['13P', 1, ' '] 364. En effet, une des deux molécules de 13P présente deux conformations alternatives dans la chaîne A avec un taux d’occupation équivalent a 50% pour chaque configuration. Cela ne concerne que la chaîne A, et si on fait dérouler la boucle pour chain_n=2 nous aurons le résultat ['13P', 0, ' '] compatible avec l’examen attentif des données du fichier PDB.
L’affichage des atomes ne présente pas de surprises, j’ai simplement masqué l’affichage de atom.vector qui, par exemple, aurait donné <Vector 25.33, 54.01, 39.47> pour le premier atome du résidu A:13P.

Par contre on remarquera que le numéro d’atome n’est pas implémenté dans get_atom_info, pour cela il suffit d’ajouter la ligne info.append(self.get_serial_number()) dans la fonction. Il faut se reporter au module Atom.py du paquetage Bio.PDB pour détecter d’autres attributs utilisables.

Et coté ligand ?

Cette fois nous pouvons afficher les coordonnées atomiques du ligand (residue.resname='13P' ou res_n=364) nous aurons  :

Nous identifions les coordonnées moléculaires dans une structure de type array (Numpy) et nous voyons que le ligand est bien doté d’atomes désordonnées, avec des facteurs d’occupation à 0.5 pour chaque atome du 13P.

5. Conclusion

Il s’agit ici d’un panorama limité mais qui donne déjà une idée précise des possibilité du parseur PDB inclus dans Biopython. Mais même à la date de révision de l’article il ne correspond pas encore à nos besoins orientés drug design, il manque tous ces outils qui nous permettent d’évaluer la qualité d’une structure et de la préparer pour des calculs, d’une manière automatisée.

C’est ce qui a motivé la conception d’un parseur sur des principes simples. Une fois que la structure de la protéine (hiérarchique par nature) est ‘sérialisée’ dans le fichier PDB, il suffit de s’en servir et d’alimenter des structures de données spécialisées: liste chainée, matrice, dictionnaire …. Cette approche n’est pas choisie par principe mais parce qu’elle est la plus optimale en matière de simplification algorithmique ou de calculs. Si nous avons besoin de travailler sur les résidus, nous générons une resmap, s’il s’agit d’interactions nous générons une proxmap … Différentes matrices adaptées à des utilisations type, qui nous rendent la vie plus facile que de reconstruire à chaque fois les infos issues d’un itérateur. Le parseur étant un objet de niveau supérieur qui regroupe ces structures de données dans une construction cohérente.

Par contre, PDB.Bio possède des fonctions qui sont très intéressantes et faciles à programmer : rotations et translations de coordonnées moléculaires, alignement structuraux, superpositions, calcul de surfaces MSMS.