davidbgk
/
larlet-fr-david

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&#39;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

<p>L'une des <del>galère</del> 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 <a href="http://www.rcsb.org/pdb/">Protein Data Bank</a>, 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 «&nbsp;squelette » de parser de PDB en python, par «&nbsp;squelette » j'entend que vous allez devoir coder les parties spécifiques à vos besoins.</p>

<h2>Structure d'un fichier PDB</h2>


<p>Voila en gros à quoi ressemble un fichier PDB&nbsp;:</p>


<pre>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</pre>

<p>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 <a href="#">celui de la protéine 1SDK</a> (sinon elle est aussi disponible dans l'archive proposée en téléchargement). Là vous vous dites comme n'importe quel geek&nbsp;: <em>«&nbsp;Les bio-informaticiens ne connaissent pas le XML ?! »</em> 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...</p>


<h2>Le template PDBParser.py</h2>


<p>Le code est vraiment court mais il est très important de le comprendre pour continuer&nbsp;:</p>


<pre>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)</pre>


<p>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.</p>


<p>Par exemple PDB['COMPND'] contient&nbsp;:</p>


<pre>['COMPND    MOL_ID: 1;                                                            ',
 'COMPND   2 MOLECULE: HEMOGLOBIN A;                                              ',
 'COMPND   3 CHAIN: A, B, C, D;                                                   ',
 'COMPND   4 OTHER_DETAILS: CROSS-LINKED, CARBONMONOXY                            ']</pre>


<p>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.</p>


<p>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&nbsp;! Elle est assez simple à comprendre&nbsp;: si une classe intitulée <strong>SixpremierscaractèresdelaligneParser</strong> 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'.</p>


<p>Et c'est tout&nbsp;? 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 :)</p>


<h2>Les exemples d'application</h2>


<p>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).</p>


<h3>Récupérer les coordonnées atomiques de la chaîne C des acides aminés leucine (LEU) et des carbones alpha (CA)</h3>


<p>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&nbsp;:</p>


<pre>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]</pre>


<p>Cette classe permet d'obtenir les résultats suivants&nbsp;:</p>


<pre>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']]</pre>


<p>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.</p>


<p>Bon maintenant comment ça marche&nbsp;?</p>


<p>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&nbsp;: 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 <strong>get_coordinates()</strong> permet de récupérer les coordonnées atomiques des atomes considérés sous forme de liste sous la forme [x, y, z].</p>


<p>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 ;-)</p>


<h3>Quelle est la position des hélices alpha de la chaîne B&nbsp;?</h3>


<p>On passe maintenant aux structures secondaires et plus particulièrement aux lignes préfixées par 'HELIX'.</p>


<pre>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)</pre>


<p>Cette classe permet d'obtenir les résultats suivants&nbsp;:</p>


<pre>PDB['HELIX']['A']['GLY'].__repr__()
['15', '18']</pre>


<p>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.</p>


<p>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é.</p>


<h3>Je voudrais les premier et dernier auteurs s'il vous plaît</h3>


<p>Bon là vous commencez vraiment à comprendre le mécanisme, un dernier exemple sans grande utilité pour être sûr.</p>


<pre>class AuthorParser(DefaultClass):
    def __init__(self, authors):
        self.authors = authors[0][6:-9].strip().split(',')
        DefaultClass.__init__(self, self.authors)</pre>


<p>Cette classe permet d'obtenir les résultats suivants&nbsp;:</p>


<pre>print PDB['AUTHOR'][0], PDB['AUTHOR'][-1]
M.A.SCHUMACHER R.G.BRENNAN</pre>


<p>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.</p>


<h2>Performances</h2>


<p>Tout à un prix quand même et rien ne sera jamais plus rapide qu'un parsing «&nbsp;à 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&nbsp;:</p>


<pre>PDB = PDBParser(open("%s.pdb" % pdbid))
coordinates.append(PDB['ATOM'][chain]['LEU']['CA'].get_coordinates())</pre>


<p>contre 52 secondes pour&nbsp;:</p>


<pre>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']</pre>


<p>Bon premier constat&nbsp;: vive les list-comprehension, second constat&nbsp;: on n'obtient pas les mêmes résultats en fait puisque le premier <strong>coordinates</strong> 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&nbsp;:</p>


<pre>_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)</pre>


<p>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.</p>


<p><strong>[edit]</strong>&nbsp;: je n'ai pas fait de tests plus poussés mais il est certain que l'appel à <strong>globals()</strong> dans la fonction __getitem__ de la classe <strong>PDBParser</strong> 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.</p>


<h2>Téléchargement</h2>


<p>Vous pouvez télécharger le code commenté, les exemples et un fichier PDB contenus <del><a href="#">dans cette archive</a></del> (mise à jour cf. ci-dessous).</p>


<h2>En conclusion</h2>


<p>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&nbsp;!</p>


<h2>Mise à jour</h2>


<p>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 <strong>get_resolution()</strong> renvoie un float)&nbsp;:</p>


<pre>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</pre>


<p>Et voici le <del><a href="#">fichier complet mis à jour</a></del> (mise à jour cf. ci-dessous).</p>


<h2>Mise à jour du 25/04/06</h2>


<p>Bon je voulais pas que ça devienne une usine à gaz mais c'est pas une raison pour ne pas faire une documentation potable :-).</p>


<p>Voici donc la <del><a href="#">version 3 du parser</a></del>. 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.</p>


<h2>Mise à jour du 20/01/08</h2>


<p>Un grand merci à François Vallée qui a remonté ses améliorations&nbsp;: 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 <a href="/static/david/biologeek/realisations/PDBParser/PDBParser4.zip">télécharger la version 4</a>.</p>