Savoir intégrer un alignement de séquence dans un processus structural n’est pas une option. Par exemple UCSF Chimera/Matchmaker utilise ce type d’approche pour faciliter les alignements de structures. A notre niveau, et même si nous utilisons des structures déjà alignées (superposables en 3D) c’est aussi parfois le cas. Je vais évoquer les étapes préalables à un calcul de RMSD (résidu à résidu) dont la finalité sera la caractérisation quantitative de fluctuations structurales. Nous pouvons en avoir besoin pour de la classification (in fine pour un choix éclairé des structures à docker).

D’après un article publié dans buildblog.buidez.net en 2014 – Mise à jour mars 2025.

1. Séquences ‘structurales’ vs. séquences

Nous devrons d’abord prendre en compte des particularités dans les structures PDB et c’est la que l’alignement de séquences va nous aider. Pour une même protéine, nous sommes d’accord que les séquences peuvent être différentes si elles sont issues d’organismes différents (ex: mutations). Mais dans le cas d’une même protéine, issue du même organisme, tissu, sous-ensemble cellulaire, les séquences doivent être identiques. Par contre, si elles sont issues non pas d’une séquence primaire mais déduites d’une structure 3D (en particulier une cristallographie RX) ce que j’appelle des séquences structurales, elles peuvent être différentes.
En effet, il est est tout à fait possible que certaines parties des deux protéines n’apparaissent pas dans les structures PDB. Par exemple, des portions de boucles (plus ou moins mobiles ou qui ont besoin d’interactions pour se stabiliser). Ce n’est jamais prévisible et cela dépends de facteurs à la protéine ou aux processus expérimentaux. Nous trouverons donc des régions apparemment sans chaine et résidus sur l’une et/ou l’autre des structures, que nous allons appeler des ‘pseudos’ délétions. Il ne s’agit pas de vraies délétions, les résidus existent dans le cristal, mais il n’apparaissent pas dans la structure 3D. Sans statistiques précises, cela dépends des séries protéiques, j’estime que les deux tiers de la PDB pourraient correspondre à des structures avec une ou plusieurs pseudo délétions.
Dans l’entête d’un fichier PDB ces informations peuvent être présentes ou déductibles car la séquence primaire est incluse. Premier problème, lorsqu’on travaille sur des fichiers qui ne contiennent que des coordonnées structurales, ce qui est l’essentiel de nos données, nous n’avons pas d’entête et faire suivre des annotations serait trop lourd. Second problème, il ne suffit pas de prendre chaque structure par un bout et de la dérouler pour identifier les pseudo délétions. En effet, les structures ne débutent pas forcément à partir du même résidu dans la séquence, pour des raisons expérimentales (bioingéniérie) et parce que les extrémités des chaines polypeptidiques sont souvent flexibles ou (misfolded) mal repliées. Donc certains acides aminés d’une séquence vont se trouver sans contrepartie puis décalés dans l’autre. Finalement c’est l’alignement de ces séquences ‘pseudo’ tronquées qui va nous permettre de remettre en vis à vis les résidus (identiques, alignés et conservés) que l’on va comparer dans le calcul de RMSD.

2. Un exemple d’alignement (pairwise2)

Tout d’abord prenons cas général d’alignement de séquences (non identiques) issues d’organismes différents. Il s’agit d’un cas d’école ancien décrit dans blog de bio-informatique, qui utilise une partie de p53 (facteur de transcription) en comparant une séquence humaine et une séquence de souris. Le code suivant permet de lancer un alignement global, il utilise un appel direct à pairwise2 que l’on trouve dans le paquetage Biopython :

La première partie consiste à définir les chaines de caractères r_seq (p53 humain) et q_seq (p53 souris) qui définissent la séquence de référence et la séquence de test (q pour query). Puis on sélectionne la matrice blosum62, les valeurs de pénalités pour les gap (gap_open, gap_extend), puis l’alignement est réalisé au moyen de la fonction pairwise2.align.globalds, ce qui donne:

*** Sequence global [GLOBAL] alignment [DS] -> Reference sequence (60) residues MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP -> Query sequence (60) residues MEESQSDISLELPLSQETFSGLWKLLPPEDILPSPHCMDDLLLPQDVEEFFEGPSEALRV -> Set matrix to {blosum62} -> matrix loaded -> Set gap penalties: gap_open=-10.000000, gap_extend=-0.500000 -> Doing algnment Found [3] alignments (r_aln, q_aln, score, begin, end)

On constate que 3 alignements optimaux ont été trouvés, en effet la fonction pairwise2.align.globalds renvoie une liste d’alignements possibles et correspondant à un score donné. Ceci explique la boucle d’affichage à la suite du calcul, ce qui nous donnera le résultat :

====================================================================== MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP---- MEESQSDISLELPLSQETFSGLWKLLPPEDIL-PSP-HCMDDLLL-PQDVEEFF-EGPSEALRV 141.5 0 64 MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP------ MEESQSDISLELPLSQETFSGLWKLLPPEDIL-PSP-HCMDDLLL-PQDVEEFFE---GPSEALRV 141.5 0 66 MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP------ MEESQSDISLELPLSQETFSGLWKLLPPEDIL-PSP-HCMDDLLL-PQDVEEFF---EGPSEALRV 141.5 0 66 ====================================================================== Score [141.500000] begin [0] end [64] : for first alignment

On constate que les alignements ont le même score et diffèrent par le nombre de gaps et la position des gaps, en particulier à la fin de la séquence. On remarquera aussi le nom plus ou moins bizarre de la fonction qui est appelée, en particulier les deux dernières lettres (pairwise2.align.globalds) qui sont utilisées pour spécifier un mode de fonctionnement (d: A dictionary returns the score of any pair of characters, s: Same open and extend gap penalties for both sequences).

Limitations
En ce qui concerne le type (algorithme) d’alignement optimal utilisé par pairwise2, sous la forme Needleman-Wunsch (global) ou Smith-Waterman (local) il existe des limitations en termes de lenteur et de bugs. Dans un contexte de calculs intensifs sur les séquences, il faudra peut être utiliser un autre module ou un autre outil. Mais dans le contexte ‘aide à des comparaisons de structures‘ (et non de séquences) c’est suffisant, avec un avantage, c’est qu’il n’y a ni logiciel externe, ni passage par un fichier intermédiaire.

3. Un cas réel (structures InhA)

Une fois que l’on a vu comment les choses fonctionnent, je vais utiliser les fonctions d’alignement de pairwise2 au moyen de deux fonctions d’interface du paquetage buildez.bio: bio_interface_2str2aln_globalds pour l’alignement global et bio_interface_2str2aln_localds pour un alignement local. Cette fois le jeu de séquences (référence et query) correspond à deux structures PDB d’une même protéine InhA :

## taken from inha protein r_seq = "TGLLDGKRILSGIITDSSIAFHIARVAQEQGAQLVLTGFDRLQVLRRLTDRLPAKAPLLELDLQNEEHLASLAGRVTEAIGA" q_seq = "TGLLDGKRILVSGIITDSSIAFHIAQGAQLVLTGFDRLRLIQRITDRLPAKAPLLELDVQNEEHLASLAGRVTEAIGA"

Ces séquences identifient un motif structural, d’un numéro de résidu à un autre et elles sont issues des structures cristallographiques. Nous voyons qu’elles n’ont pas la même longueur. Cette partie de la séquence d’InhA n’est pas due au hasard, elle encadre la SBL (Substrate Binding Loop) dont on veut observer les fluctuations et qui comporte parfois ces pseudos délétions. Donc nous allons faire un alignement de ces séquences.

Calculs (global vs. local)
Nous allons utiliser un code légèrement différent, pour les deux types d’alignement:

Ce qui nous donnera comme résultat les deux alignements suivants:

*** Sequence alignments (same parameters) using wrappers -> bio_interface_2str2aln_globalds Score [141.500000] begin [0] end [64] : for first alignment MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP---- MEESQSDISLELPLSQETFSGLWKLLPPEDIL-PSP-HCMDDLLL-PQDVEEFF-EGPSEALRV -> bio_interface_2str2aln_localds Score [154.000000] begin [0] end [60] : for first alignment MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP------ MEESQSDISLELPLSQETFSGLWKLLPPEDIL-PSP-HCMDDLLL-PQDVEEFFE---GPSEALRV

Nous constatons qu’il y a quelques différences minimes entre alignement local et alignement global. Les fonctions d’interface analysent les résultats de l’alignement et renvoient un seul alignement, le premier dans la liste retournée par pairwise2.align.globalds ou pairwise2.align.localds. On peut programmer les fonction d’interface pour différents choix, mais celui ci est le plus simple. De toutes manières, l’objectif est de fournir un cadre aux fonctions qui vont calculer les RMSD.

Cet exemple montre comment intégrer un processus d’alignement de séquence au milieu d’un code de biologie structurale, pour peu que nous manipulions des séquences sous la forme de chaines de caractères simples. Maintenant nous avons deux chaines de caractères alignées qui vont nous servir de guide pour les calculs, il suffira de s’en servir pour annoter une resmap, pour d’utiliser la fonction pdb_resmap2rmsd de buildez.pdb.localrms pour faire les calculs. Les résidus qui n’ont pas de vis à vis (gap) afficheront un RMSD égal à zéro (par convention).