6. Eclatement de structures
C’est une fonction dont on a très souvent besoin en drug design. Donc nous éclatons une structure en chaines (A, B …) puis nous alignons les chaines entre elles pour les comparer, puis nous classifions. Une fois le processus répété et croisé sur plusieurs structures PDB, nous pouvons définir des chaines canoniques, représentatives d’un état structural, puis nous allons en choisir une, qui va être utilisée pour des processus in silico. Mais à la base, il faut éclater (split) les structures en chaînes, avec plusieurs modalités opératoires :
- Avec un éditeur de texte, en supprimant les lignes correspondant aux autres chaînes et en répétant le processus sur toutes les chaines.
- Avec du logiciel, j’utilise DSV de Discovery Studio à cause de sa présentation hiérarchisée de la structure et de ses possibilités étendues de couper-copier-coller.
- En mode automatique, avec un protocole naïf, par exemple un programme type GREP qui va filtrer les lignes avec la colonne ‘ A’ (pour la chaine A). Mais nous allons arriver rapidement aux limite de la méthode, car il y aura des exceptions non gérées, liées à l’écriture du fichier PDB.
Quel est le problème ?
En général les ligands sont affectés de la même lettre s’ils sont dans la même chaine, une lettre qui apparaît à la colonne 5 (index) des lignes ATOM ou HETATM. Sauf que … les ligands (ni même les molécules d’eau) n’ont pas toujours cette affectation (colonne CHAINID) de la chaine. C’est un problème de qualité au niveau des fichiers PDB, qui a tendance à s’atténuer mais qui existe. Pour vous en convaincre, il suffit de regarder dans les exemples précédents et de cherche l’identifiant A ou B, pour les lignes (par exemple HETATM) correspondant à la biotine (ligand BTN).
Contournement
D’autres approches (parser based) plus robustes et permettant plus de possibilités sont intégrées dans le paquetage buildez.pdb mais en mode chaine nous pouvons quand même faire quelque chose … car nous disposons de fonctions d’extraction de blocs.
Le principe est simple, imaginons que nous avons une structure avec ceux chaines A et B, un ligand dans chaque chaine, non identifié. Il suffira d’extraire le premier bloc protéine et le premier bloc ligand, pour avoir la première chaine complète, protéine + ligand. Même si ces blocs ne sont pas consécutifs dans un fichier PDB.
D’où l’intérêt d’une fonction pour compter les blocs comme pdb_coordsblock_count qui n’apparaissait pas évidente à première vue, mais qui va nous permettre de diriger le calcul.
La séquence suivante montre comment on peut extraire de manière automatisée la chaine A de la structure protéique incluant ou non le ligand correspondant, sans que le ligand soit affecté d’un identifiant de chaine. Avec une petite gymnastique pour extraire les données, sinon il aurait suffi de sortir directement les lignes HETATM et A qui ne correspondent pas à de l’eau (HOH ici) ou qui correspondent à BTN.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
print("=> check for chain 'A' block in coordinates") bn = pdb_coordsblock_count(s,"ATOM",pdb_field_coords_id("CHAINID"),"A","=") print("number of blocks for this request %d" % (bn)) print("=> check for HETATMs blocks in coordinates") bh = pdb_coordsblock_count(s,"HETATM",pdb_field_coords_id("RESSEQ"),"0",">") print("number of blocks for this request %d" % (bh)) print("=> get the first block (ATOM and HETATM) corresponding to 'A'") s1 = pdb_str_coordsblock_get(s, "ATOM HETATM", pdb_field_coords_id("CHAINID"), "A", "=", 0) bl = len(pdb_str_pgrep_get(s1, "HETATM")) print("=> number of HETATMs in this block %d" % (bl)) print("=> get the residue number 500") s2 = pdb_str_coordsblock_get(s,"HETATM",pdb_field_coords_id("RESSEQ"),"500","=", 0) print("=> combine the two strings to get the first subunit with ligand\n") s = s1 + "\\n" + s2 print(s) |
Ce qui va nous donner un joli bloc protéine + ligand que nous n’aurions pas pu atteindre avec une approche basée sur un GREP externe + identifiant de chaine.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
ATOM 1 N ASN A 4 66.039 30.785 50.117 1.00 88.31 N ATOM 2 CA ASN A 4 66.049 29.799 51.200 1.00 85.31 C ATOM 3 C ASN A 4 67.440 29.525 51.712 1.00 73.80 C ... ATOM 2354 CZ ARG A 317 43.198 -17.198 63.688 1.00100.00 C ATOM 2355 NH1 ARG A 317 43.450 -18.165 62.799 1.00100.00 N ATOM 2356 NH2 ARG A 317 42.253 -17.475 64.605 1.00100.00 N HETATM 4715 C11 BTN 500 55.016 4.107 53.277 1.00 24.55 C HETATM 4716 O11 BTN 500 53.849 4.252 52.886 1.00 21.58 O HETATM 4717 O12 BTN 500 55.291 3.685 54.427 1.00 19.36 O HETATM 4718 C10 BTN 500 56.168 4.363 52.354 1.00 48.06 C HETATM 4719 C9 BTN 500 56.855 5.658 52.754 1.00 10.86 C HETATM 4720 C8 BTN 500 57.925 6.100 51.752 1.00 2.02 C HETATM 4721 C7 BTN 500 58.527 7.477 52.110 1.00 5.93 C HETATM 4722 C2 BTN 500 59.794 7.809 51.391 1.00 13.36 C HETATM 4723 S1 BTN 500 60.680 9.282 52.000 1.00 22.04 S HETATM 4724 C6 BTN 500 61.626 9.423 50.461 1.00 13.06 C HETATM 4725 C5 BTN 500 60.674 9.025 49.349 1.00 42.69 C HETATM 4726 N1 BTN 500 59.967 10.188 48.871 1.00 11.08 N HETATM 4727 C3 BTN 500 58.738 10.251 49.397 1.00 83.40 C HETATM 4728 O3 BTN 500 58.041 11.266 49.478 1.00 18.36 O HETATM 4729 N2 BTN 500 58.385 9.033 49.882 1.00 12.65 N HETATM 4730 C4 BTN 500 59.522 8.140 49.938 1.00 60.80 C |
A noter qu’il nous manque quelque chose : les blocs CONECT à la fin du fichier, qui sont indispensables pour identifier quels atomes de ligand sont liés entre eux. En effet, sur ces lignes, il n’y a pas d’information de chaîne, ni de marqueur de terminaison pour chaque ligand. C’est un bloc unique à numérotation souvent contiguë et qui nécessite une analyse sur la base des numéros d’atomes (de ligand) que l’on va récupérer dans les lignes HETATM, donc au delà d’un filtre texte.