# ★ Parser un fichier PDB en python facilement et efficacement

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.

### Articles peut-être en rapport

- ★ De l'OpenData au LinkedData : exemple de nosdonnees.fr
- ★ Pourquoi Python et Django
- ★ Django-ROA, pour une architecture orientée ressources

### Commentaires

Martin le 18/04/2006 :

Bonsoir
Merci beaucoup pour partager ce "bout" de code... Maintenant, moi la seule chose que je veux est 'spliter' un PDB en ces chaines (A, B, C, etc). Si je prends vos classes et je met le print suivant:
print PDB['ATOM']['C'].get_coordinates()[:10]
j'obtient:
#####
$ ./PDBparser.py 1RJD.pdb
Traceback (most recent call last):
File "./PDBparser.py", line 64, in ?
print PDB['ATOM']['C'].get_coordinates()[:10]
NameError: name 'PDB' is not defined
#####
des idees ? merci encore!
martin
[mgrana_at_gmail.com]

David, biologeek le 19/04/2006 :

Il faut d'abord créer une instance de PDBParser de la manière suivante dans le code :
PDB = PDBParser(open("1SDK.pdb"))
Bon je suis d'accord c'est pas très pratique (surtout qu'à la dernière mise à jour les archives ont sautées :/) et il faut que je documente beaucoup mieux tout ça. À l'époque je ne savais pas le faire mais maintenant oui :)
Je mettrais à jour prochaînement le billet avec de la documentation agile.

Jerome Pansanel le 30/04/2007 :

Est-ce que vous avez pensé à utiliser cette classe pour écrire un outil permettant de préparer une protéine pour un docking ?

David, biologeek le 30/04/2007 :

Non, cette bibliothèque est un peu à l'abandon en fait... mais bon le format pdb évolue peu donc elle est toujours utilisable.
Si vous avez des suggestions d'améliorations, n'hésitez pas, je peux les publier ici.

Yann le 21/06/2009 :

Salut,

je suis un développeur Python qui dispose de pas mal de temps libre devant lui et j'aimerai bien jouer un peu avec les protéines. Léger problème : je manque cruellement de culture générale (pour l'instant) à ce sujet.

Par où (quelle documentation) pensez-vous que je devrais commencer pour, à terme, comprendre quelle protéïne fait quoi et pourquoi.

Merci :)

David, biologeek le 23/06/2009 :

Hello Yann,

Question difficile :)

En fait ça dépend un peu de ta façon d'apprendre, il y a 3 options à mon avis :

* une formation, coûteux mais si les formateurs sont intéressants ça peut faire gagner un temps fou.
* un livre, je sais que "le Campbell" (du nom de son auteur) était notre bible en prépa. Mais il doit sûrement y avoir un peu plus accessible pour commencer.
* internet, les articles de wikipedia sont à mon avis un bon point de départ.

Il y aurait peut-être une 4ème option qui serait de rencontrer des biologistes qui sont souvent intarissables sur le sujet et qui combineraient les 3 ;-).

Bon courage !

Maamri sarra le 26/10/2009 :

bonjour ,
est ce que je peu avoir votre aide pour utilisé la classification CATCH dans le site PDB :et comment je lit les resultats

David, biologeek le 26/10/2009 :

Si le problème est suffisamment bien décrit pourquoi pas, puis ça te permettra peut-être de trouver la solution par toi-même d'essayer de le formuler :-)