article.md 19KB
title: ★ Parser un fichier PDB en python facilement et efficacement slug: parser-un-fichier-pdb-en-python-facilement-et-efficacement date: 2005-10-17 16:21:03 type: post vignette: images/logos/python_nouveau.png contextual_title1: ★ De l'OpenData au LinkedData : exemple de nosdonnees.fr contextual_url1: 20101130-de-lopendata-au-linkeddata-exemple-de-nosdonneesfr contextual_title2: ★ Pourquoi Python et Django contextual_url2: 20091211-pourquoi-python-et-django contextual_title3: ★ Django-ROA, pour une architecture orientée ressources contextual_url3: 20090526-django-roa-pour-une-architecture-orientee-ressources
L'une des galère tâche quotidienne d'un bio-informaticien s'intéressant un tant soit peu à la biologie structurale et à la position d'une protéine dans l'espace est de parser des fichiers PDB. C'est fichiers, issus de la Protein Data Bank, contiennent une foultitude d'informations plus ou moins pertinentes en fonction de son sujet de recherche. Or, les parsers actuels que je connais sont pour la plupart de véritables usines à gaz souvent inadaptés et lourds pour le traitement que j'ai à faire. Voici donc un « squelette » de parser de PDB en python, par « squelette » j'entend que vous allez devoir coder les parties spécifiques à vos besoins.
Structure d'un fichier PDB
Voila en gros à quoi ressemble un fichier PDB :
HEADER OXYGEN TRANSPORT 26-FEB-96 1SDK TITLE CROSS-LINKED, CARBONMONOXY HEMOGLOBIN A COMPND MOL_ID: 1; COMPND 2 MOLECULE: HEMOGLOBIN A; COMPND 3 CHAIN: A, B, C, D; COMPND 4 OTHER_DETAILS: CROSS-LINKED, CARBONMONOXY SOURCE MOL_ID: 1; SOURCE 2 ORGANISM_SCIENTIFIC: HOMO SAPIENS; SOURCE 3 ORGANISM_COMMON: HUMAN KEYWDS HEME, OXYGEN TRANSPORT, RESPIRATORY PROTEIN, ERYTHROCYTE, KEYWDS 2 DISEASE MUTATION, POLYMORPHISM EXPDTA X-RAY DIFFRACTION AUTHOR M.A.SCHUMACHER,M.M.DIXON,R.KLUGER,R.T.JONES,R.G.BRENNAN REVDAT 1 01-AUG-96 1SDK 0 [SNIP] HELIX 1 1 PRO A 4 GLY A 15 1 12 HELIX 2 2 GLY A 18 SER A 35 5 18 HELIX 3 3 PRO A 37 PHE A 43 5 7 HELIX 4 4 ALA A 53 ALA A 71 1 19 HELIX 5 5 VAL A 73 HIS A 89 5 17 HELIX 6 6 PRO A 95 HIS A 112 1 18 HELIX 7 7 PRO A 119 TYR A 140 1 22 HELIX 8 8 PRO B 5 LYS B 17 1 13 HELIX 9 9 VAL B 23 VAL B 34 1 12 HELIX 10 10 PRO B 36 PHE B 42 5 7 [SNIP] ATOM 1 N VAL A 1 36.846 49.802 -14.919 1.00 89.28 N ATOM 2 CA VAL A 1 37.907 49.146 -14.140 1.00 89.20 C ATOM 3 C VAL A 1 37.655 49.392 -12.632 1.00 84.23 C ATOM 4 O VAL A 1 38.531 49.598 -11.742 1.00 82.54 O ATOM 5 CB VAL A 1 39.361 49.471 -14.649 1.00 93.91 C ATOM 6 CG1 VAL A 1 39.441 49.501 -16.172 1.00 93.10 C ATOM 7 CG2 VAL A 1 39.890 50.778 -14.094 1.00 94.26 C ATOM 8 N LEU A 2 36.364 49.205 -12.405 1.00 68.64 N ATOM 9 CA LEU A 2 35.719 49.461 -11.152 1.00 62.16 C ATOM 10 C LEU A 2 35.826 50.891 -10.976 1.00 69.86 C ATOM 11 O LEU A 2 36.748 51.514 -10.445 1.00 70.82 O ATOM 12 CB LEU A 2 35.859 48.504 -10.004 1.00 56.68 C ATOM 13 CG LEU A 2 34.686 47.515 -10.109 1.00 45.20 C ATOM 14 CD1 LEU A 2 34.868 46.618 -11.293 1.00 33.18 C ATOM 15 CD2 LEU A 2 34.615 46.658 -8.886 1.00 57.98 C ATOM 16 N SER A 3 35.120 51.235 -12.022 1.00 61.69 N ATOM 17 CA SER A 3 34.707 52.459 -12.504 1.00 57.56 C ATOM 18 C SER A 3 33.527 52.737 -11.624 1.00 50.93 C
Les morceaux présentés n'ont pas été choisis par hasard et vont servir dans la suite de la présentation. Si vous souhaitez obtenir un fichier PDB complet, voici celui de la protéine 1SDK (sinon elle est aussi disponible dans l'archive proposée en téléchargement). Là vous vous dites comme n'importe quel geek : « Les bio-informaticiens ne connaissent pas le XML ?! » en fait si, ce format est aussi disponible mais malheureusement pour une même protéine, on atteint 5,1 Mo contre 440 Ko pour la même en texte brut. Sachant qu'il faut généralement commencer par les télécharger pour travailler dessus et qu'on compte en milliers de protéines, vous voyez aisément le problème sans parler de l'espace disque qu'il faut alors considérer...
Le template PDBParser.py
Le code est vraiment court mais il est très important de le comprendre pour continuer :
class PDBParser:
def __init__(self, file_in):
self.PDB = {}
self.properties = []
for line in file_in.readlines():
try:
self.PDB[line[0:6].strip()].append(line[0:-1])
except KeyError:
self.PDB[line[0:6].strip()] = [line[0:-1]]
self.properties.append(line[0:6].strip())
def __getitem__(self, key):
try:
parserName = "%sParser" % key.capitalize()
parserClass = globals()[parserName]
return parserClass(self.PDB[key])
except KeyError:
return DefaultClass(self.PDB[key])
def __getslice__(self, begin, end):
return [line for property in self.properties \
for line in self.PDB[property]][begin:end]
def __repr__(self):
return [line for property in self.properties \
for line in self.PDB[property]]
def __str__(self):
return "
".join([line for property in self.properties \
for line in self.PDB[property]])
class DefaultClass:
def __init__(self, defaults): self.defaults = defaults
def __getslice__(self, begin, end): return self.defaults[begin:end]
def __getitem__(self, key): return self.defaults[key]
def __repr__(self): return self.defaults
def __str__(self): return "
".join(self.defaults)
Dans la classe PDBParser, __init__ et __getitem__ sont les méthodes spéciales de classe importantes, la première crée le dictionnaire PDB qui va contenir, pour chaque clé correspondant aux 6 premiers caractères d'une ligne, une liste de valeurs correspondant aux lignes en question.
Par exemple PDB['COMPND'] contient :
['COMPND MOL_ID: 1; ', 'COMPND 2 MOLECULE: HEMOGLOBIN A; ', 'COMPND 3 CHAIN: A, B, C, D; ', 'COMPND 4 OTHER_DETAILS: CROSS-LINKED, CARBONMONOXY ']
La liste properties permet de conserver l'ordre des clés de façon à pouvoir restituer le PDB complet en cas de besoin. Par exemple pour un affichage au moyen des méthodes __getslice__, __repr__ ou __str__. Par exemple PDB[:10] renverra les dix premières lignes du fichier PDB sous forme de liste.
La méthode __getitem__ est vraiment LA méthode clé de la classe, c'est elle qui rend le template modulable et donc si puissant ! Elle est assez simple à comprendre : si une classe intitulée SixpremierscaractèresdelaligneParser existe, elle est instanciée, sinon c'est la classe DefaultClass qui va l'être. DefaultClass contenant les méthodes spéciales de classes usuelles, elle permet d'afficher tout ou partie des données, par exemple PDB['COMPND'].__repr__() renvoie la liste des lignes préfixées par 'COMPND'.
Et c'est tout ? Pour l'instant oui, après c'est à vous de coder ce dont vous avez besoin, allez je suis sympa voici quelques exemples concrêts :)
Les exemples d'application
Pas évident de voir à première vue comment fonctionne le template alors voila quelques classes d'exemples qui pourront peut-être même vous servir (probablement à adapter).
Récupérer les coordonnées atomiques de la chaîne C des acides aminés leucine (LEU) et des carbones alpha (CA)
Il est assez courant de ne vouloir que quelques atomes ou une seule chaîne ou qu'un type d'acides aminés, cette classe permet de récupérer simplement ces informations dans votre PDB, puisqu'on va traiter les lignes préfixées par 'ATOM', notre classe doit s'appeler AtomParser :
class AtomParser(DefaultClass):
def __init__(self, atoms, *par):
if len(par) == 0:
self.atoms = atoms
elif len(par[0]) == 1:
self.atoms = [atom for atom in atoms if atom[21] == par[0]]
elif len(par[0]) == 2:
self.atoms = [atom for atom in atoms if atom[13:15] == par[0]]
elif len(par[0]) == 3:
self.atoms = [atom for atom in atoms if atom[17:20] == par[0]]
DefaultClass.__init__(self, self.atoms)
def __getitem__(self, key): return AtomParser(self.atoms, key)
def get_coordinates(self):
return [[atom[29:38].strip(),atom[38:46].strip(),atom[46:54].strip()]\
for atom in self.atoms]
Cette classe permet d'obtenir les résultats suivants :
print PDB['ATOM']['C']['LEU']['CA'] ATOM 2293 CA LEU C 2 56.540 41.891 -11.769 1.00 31.34 C ATOM 2483 CA LEU C 29 53.070 39.062 -36.319 1.00 5.06 C ATOM 2530 CA LEU C 34 49.611 32.028 -40.266 1.00 9.46 C ... print PDB['ATOM']['LEU']['C '].get_coordinates()[:10] [['35.826', '50.891', '-10.976'], ['28.628', '24.722', '-5.569'], ['27.926', '19.794', '-11.999'], ['28.129', '17.183', '-1.981'], [SNIP], ['32.035', '40.188', '-0.663']]
C'est donc assez pratique, l'ordre des paramètres n'a pas d'importance, PDB['ATOM']['C']['LEU']['CA'] est équivalent à PDB['ATOM']['C']['CA']['LEU'] et on peut omettre une information au besoin, par exemple pour récupérer les mêmes informations pour toutes les chaînes, PDB['ATOM']['LEU']['CA'] suffit.
Bon maintenant comment ça marche ?
En fait l'astuce est de considérer ici la taille du paramètre passé en argument, par exemple pour le choix de la chaîne, sa taille est de 1 et l'on filtre les lignes composant atoms en conséquence. On initialise ensuite la classe DefaultClass parente avec les données appropriées ce qui permet d'hériter de ses méthodes spéciales. Deuxième astuce : la classe est récursive, ce qui permet d'enchaîner ou non les arguments et de les traiter un par un. Enfin la fonction get_coordinates() permet de récupérer les coordonnées atomiques des atomes considérés sous forme de liste sous la forme [x, y, z].
Vous l'aurez compris, la première astuce est limitante pour désigner par exemple certains atomes qui sont nommés par plus de 2 lettres (par exemple CD1), si vous avez besoin spécifiquement de ce type d'atomes (très rare à mon avis), une solution rapide est de considérer alors que les atomes ont une taille de 4 lettres. Encore une fois, ce n'est qu'un exemple à adapter ;-)
Quelle est la position des hélices alpha de la chaîne B ?
On passe maintenant aux structures secondaires et plus particulièrement aux lignes préfixées par 'HELIX'.
class HelixParser(DefaultClass):
def __init__(self, helixs, *par):
if len(par) == 0:
self.helixs = helixs
elif len(par[0]) == 1:
self.helixs = [helix for helix in helixs if helix[19] == par[0]]
elif len(par[0]) == 3:
self.helixs = []
for helix in helixs:
if helix[15:18] == par[0]:
self.helixs.append(helix[20:26].strip())
elif helix[27:30] == par[0]:
self.helixs.append(helix[32:38].strip())
DefaultClass.__init__(self, self.helixs)
def __getitem__(self, key): return HelixParser(self.helixs, key)
Cette classe permet d'obtenir les résultats suivants :
PDB['HELIX']['A']['GLY'].__repr__() ['15', '18']
Ici aussi on considère la taille de l'argument (1 pour la chaîne et 3 pour l'acide aminé) ce qui permet d'obtenir une liste des positions des acides aminés particuliers impliqués dans une hélice alpha. Lorsque l'on considère un acide aminé, il faut aller vérifier à deux endroits différents étant donné qu'il y a deux réponses par ligne.
Vous remarquerez qu'on appel de PDB['HELIX']['A'] ne renvoie pas une liste de positions mais la liste contenant les lignes 'HELIX' de la chaîne A, si vous voulez obtenir une liste dans ce cas là, il faut un peu coder car je n'ai pas encore eu cette utilité.
Je voudrais les premier et dernier auteurs s'il vous plaît
Bon là vous commencez vraiment à comprendre le mécanisme, un dernier exemple sans grande utilité pour être sûr.
class AuthorParser(DefaultClass):
def __init__(self, authors):
self.authors = authors[0][6:-9].strip().split(',')
DefaultClass.__init__(self, self.authors)
Cette classe permet d'obtenir les résultats suivants :
print PDB['AUTHOR'][0], PDB['AUTHOR'][-1] M.A.SCHUMACHER R.G.BRENNAN
Tout simple, on commence par créer une liste d'auteurs à l'initialisation et le fait d'initialiser la classe parente DefaultClass permet de disposer de toutes ses méthodes spéciales. On demande ensuite pour l'exemple d'afficher les premier et dernier auteurs.
Performances
Tout à un prix quand même et rien ne sera jamais plus rapide qu'un parsing « à la main » par exemple en testant sur un panel de 2300 protéines, on obtient 1 minute 20 secondes (avec un Celeron M 1,4GHz, disque dur 5400tr/min et 768Mo de ram) pour :
PDB = PDBParser(open("%s.pdb" % pdbid))
coordinates.append(PDB['ATOM'][chain]['LEU']['CA'].get_coordinates())
contre 52 secondes pour :
coordinates = [[atom[29:38].strip(),atom[38:46].strip(),atom[46:54].strip()] \
for atom in open("%s.pdb" % pdbid).readlines() \
if atom.startswith('ATOM') and atom[21] == chain and \
atom[13:15] == 'CA' and atom[17:20] == 'LEU']
Bon premier constat : vive les list-comprehension, second constat : on n'obtient pas les mêmes résultats en fait puisque le premier coordinates est une liste de listes alors que le second n'est qu'une liste, si l'on réitère les tests de façon à obtenir les mêmes résultats :
_coord = []
for atom in open("%s.pdb" % pdbid).readlines():
if atom.startswith('ATOM') and atom[21] == chain and \
atom[13:15] == 'CA' and atom[17:20] == 'LEU':
_coord.append([atom[29:38].strip(),atom[38:46].strip(),atom[46:54].strip()])
coordinates.append(_coord)
On obtient le même résultat en 59 secondes... ce qui fait environ 25% de temps en plus lorsque l'on se sert du parser. Ça peut paraître énorme mais si l'on n'a à parser que quelques milliers de protéines ça peut devenir rentable en terme de temps d'implémentation/lisibilité du code. Sans compter la modularité vous permettant d'obtenir plusieurs informations pour un coût pratiquement équivalent. À vous encore de faire le bon compromis.
[edit] : je n'ai pas fait de tests plus poussés mais il est certain que l'appel à globals() dans la fonction __getitem__ de la classe PDBParser nuit aux performances du script donc si vous n'avez qu'une seule classe de parsing appelée, il vaut mieux la spécifier explicitement dans le script.
Téléchargement
Vous pouvez télécharger le code commenté, les exemples et un fichier PDB contenus dans cette archive (mise à jour cf. ci-dessous).
En conclusion
Ce n'est pas révolutionnaire mais c'est bien pratique et plutôt élégant, je laisse votre imagination et surtout vos besoins compléter ces classes :-). Toute remarque/suggestion/correction/contribution est la bienvenue !
Mise à jour
On a souvent besoin des résolutions dans lesquelles ont été cristallisées les protéines, le problème est qu'il n'y a pas vraiment de standard pour placer cette information dans les PDB donc ça finit toujours en vrac dans les remarques... j'ai essayé de faire au mieux, ça passe pour 3000 protéines je vous laisse vérifier sur votre set (la fonction get_resolution() renvoie un float) :
class RemarkParser(DefaultClass):
def __init__(self, remarks):
self.remarks = remarks
DefaultClass.__init__(self, self.remarks)
def get_resolution(self):
__angstrom_lines = [line for line in self.remarks if 'ANGSTROM' in line]
# new standard (?)
if 'REMARK 2 RESOLUTION.' in __angstrom_lines[0]:
return float(__angstrom_lines[0][23:27])
# old random remarks...
else:
try:
return float(__angstrom_lines[0].split('ANGSTROM')[0][-5:].split()[-1])
except:
# print wrong line if needed (uncomment)
#print angstrom_lines[0]
return None
Et voici le fichier complet mis à jour (mise à jour cf. ci-dessous).
Mise à jour du 25/04/06
Bon je voulais pas que ça devienne une usine à gaz mais c'est pas une raison pour ne pas faire une documentation potable :-).
Voici donc la version 3 du parser. Vous pouvez maintenant lancer le programme en ligne de commande suivi d'une liste de fichiers PDB, je vous laisse voir les options dans la documentation. Bien sûr la classe est toujours utilisable par simple import.
Mise à jour du 20/01/08
Un grand merci à François Vallée qui a remonté ses améliorations : il y a dorénavant une classe pour les lignes HETATM et une méthode (get_biomt) qui récupère les BIOMT des REMARK. Vous pouvez télécharger la version 4.