Il y a toujours un moment ou nous avons besoin d’aligner des séquences protéiques, cela va de la modélisation par homologie à un simple calcul de fluctuation/écart entre deux structures protéiques. Heureusement il existe toute une série d’outils, dont certains sont prévus pour être utilisés avec un code Python. C’est le cas de Biopython (paquetage Bio) qui permet de réaliser des alignements simples ou multiples en utilisant différents algorithmes, en propre (partie intégrante de la bibliothèque) ou à travers des interfaces en ligne de commande. Dans ce cas, un binaire du logiciel de comparaison de séquences doit être installé (ex: NCBI standalone BLAST). Le paquetage  peut travailler avec d’autres paquetages tels que ReportLab ou des SGBD relationnels.

D’après un article publié dans buildblog.buidez.net en 2014.

1. Deux séquences protéiques

Nous allons comparer deux séquences protéiques, donc à priori 2 chaines de caractères. Si les séquences proviennent d’une base de séquences, nous restons dans le cas standard. Si les séquences proviennent de deux structures 3D de protéines, il faudra soit les télécharger à partir d’un site (PDB, Uniprot), soit les extraire à la main à partir du fichier PDB, soit via un logiciel de visualisation moléculaire s’il permet un copier-coller ou un export. En mode full Python, Bio permettra d’obtenir la séquence à partir de fichiers, avec des raffinements liés au formats, dans le cas le plus simple (exemple au format FASTA) :

from Bio import SeqIO for record in SeqIO.parse("example.fasta", "fasta"):     print(record.id)

Dans le contexte buildez, le paquetage buildez.pdb offre plusieurs structures de données, la plus facile à utiliser étant la resmap, car elle peut être annotée (ce qui peut être utile pour la suite des calculs). Dans ce cas on peut déléguer à Bio la partie alignement puis reprendre la main avec notre propre code et l’annotation d’un structure de données telle que la resmap qui est positionnée entre séquence et structure, nous fournira le lien entre les deux paquetages. A cet effet, un paquetage d’interface buildez.bio est disponible avec un module buildez.bio.aln pour les alignements.

Mais une fois que la séquence est disponible, il faut la valider, car si on travaille à partir d’une structure, des résidus pouvant être non standard (par exemple des sélénométhionines). Puis enfin il faudra transformer l’alphabet a 3 lettres/résidu en séquence protéique (1 lettre/résidu).

2. Gestion des matrices de similarité

Qui dit alignement de séquences, dit matrice de similarité (substitution matrix). Le module bio.aln gère les matrices de Biopython, qui sont identifiées sous forme de chaines de caractères dont on peut connaître la liste, par exemple le code suivant :

Donnera le résultat suivant, grace aux fonctions bio_interface_getsubstmatrix_list (renvoie la liste des matrices), bio_interface_getsubstmatrix_info (renvoie True ou False si la matrice existe) et bio_interface_getsubstmatrix (qui renvoie la matrice) :

Biopython a suivi des évolutions, par exemple les noms de matrices pouvaient être en minuscules et il pouvait y avoir plus de matrices disponibles. Dans une ancienne version du module j’avais :

[benner6, benner22, benner74, blosum100, blosum30, blosum35, blosum40, blosum45, blosum50, blosum55, blosum60, blosum62, blosum65, blosum70, blosum75, blosum80, blosum85, blosum90, blosum95, feng, fitch, genetic, gonnet, grant, ident, johnson, levin, mclach, miyata, nwsgappep, pam120, pam180, pam250, pam30, pam300, pam60, pam90, rao, risler, structure]

C’est comme cela, les codes de référence évoluent et c’est justement pour cela que l’on fait des paquetages d’interface qui nous permettent d’assurer une stabilité dans le temps, le code appelant ne change pas, le code d’interface change en fonction de la librairie qui est appelée. Le format de la matrice pouvait changer aussi, par exemple la fonction bio_interface_getsubstmatrix('blosum62') renvoyait un dictionnaire :

{('B', 'N'): 3, ('W', 'L'): -2, ('G', 'G'): 6, ('X', 'S'): 0, ('X', 'D'): -1, ('K', 'G'): -2, ('S', 'E'): 0, ('X', 'M'): -1, ('Y', 'E'): -2, ... ('Y', 'Y'): 7, ('T', 'K'): -1, ('Z', 'I'): -3, ('T', 'P'): -1, ('V', 'L'): 1, ('F', 'I'): 0, ('G', 'Q'): -2, ('L', 'A'): -1, ('M', 'I'): 1}

Dans tous les cas ces structures de données sont simple et on arrive toujours à les convertir, soit pour continuer à les utiliser, soit pour les injecter dans des applications non Biopython mais qui utilisent ces matrices (ce qui évite de les recoder).

3. Switcher n’est pas forker

A partir de la version 1.84 l’utilisation du module pairwise2 est démodée (Biopython utilise maintenant PairwiseAligner) et la structure de données sur laquelle on travaille pour les alignements est différente. Il faut s’adapter  et cela commence par un import qui est capable de capter la version du paquetage via la chaine __version__ (l’intégrer dans un paquetage est une très bonne pratique) :

Si le code contient encore des appels à pairwise2, cela fonctionne mais Biopython envoie un message sur la sortie standard, mais que l’on peut désactiver grâce au paquetage warnings et à un filtre (fonction warnings.simplefilter) associé à un tag dans le message d’erreur (une autre très bonne pratique). Puis il nous faut une fonction pour définir un mode d’exécution :

La nous sommes calés, le code d’interface aura un switch pour se mettre en mode pairwise2 ou Pairwise Aligner et renvoyer les mêmes données dans les deux cas.

4. En mode pairwise2

Prenons le cas de l’appel en mode pairwise2 pour la fonction bio_interface_2str2aln_globalds qui effectue un alignement global avec quelques paramètres facilement reconnaissables (séquence ce référence, séquence query, matrice, gap penalties et mode de fonctionnement) :

bio_interface_2str2aln_globalds(r_seq, q_seq, matrix_name='blosum62', gap_open=-10, gap_extend=-0.5, opmode='')

La portion de code correspondante sera:

On exploite les résultats en sachant que alns est une liste car il peut y avoir plusieurs alignements possibles, chaque item est une liste de 4 chaines de caractères incluant des simples quotes pour encadrer les séquences et des = pour définir les égalités.

[Alignment(seqA='TGLLDGKRIL-SGIITDSSIAFHIARVAQEQGAQLVLTGFDRLQVLRRLTDRLPAKAPLLELDLQNEEHLASLAGRVTEAIGA', seqB='TGLLDGKRILVSGIITDSSIAFHIA-----QGAQLVLTGFDRLRLIQRITDRLPAKAPLLELDVQNEEHLASLAGRVTEAIGA', score=333.0, start=0, end=83)]

Mais le code de bio_interface_2str2aln_globalds ne renvoie que le premier élément de la liste. Si nous affichons les variables (r_aln, q_aln, score, begin, end) nous avons :

-> bio_interface_2str2aln_globalds Score [333.000000] begin [0] end [83] : for first alignment TGLLDGKRIL-SGIITDSSIAFHIARVAQEQGAQLVLTGFDRLQVLRRLTDRLPAKAPLLELDLQNEEHLASLAGRVTEAIGA TGLLDGKRILVSGIITDSSIAFHIA-----QGAQLVLTGFDRLRLIQRITDRLPAKAPLLELDVQNEEHLASLAGRVTEAIGA

5. En mode PairwiseAligner

Dans ce cas la portion de code est différente, car en entrée il faut paramétrer (matrice, gap penalties) un objet aligner qui va ensuite nous donner un résultat via la méthode aligner.align(r_seq, q_seq). Le résultat est soit itérable (construction d’une liste) soit directement transformé en liste pour être exploité, chaque item de cette liste est traité par une fonction d’interface bio_interface_get_PairwiseAligner :

Et l’ensemble reconstitue la même liste aln que dans le cas pairwise2, nous constatons également la modification dans l’entrée des paramètres (attributs de l’objet aligner) :

Cette liste aln sera traitée de la même manière par les modules en aval.

C’est le point clé, le code en aval doit rester stable, c’est les fonctions d’interface qui s’adaptent aux changement dans les codes/librairies métier.

Ce qui suppose les bons choix en amont, notamment pour les structures de données (listes, listes de dictionnaires, listes d’objets) qui vont transiter entre ces couches. La motivation du changement dans Biopython se comprennent tout de suite si on veut faire un dump de l’itérable alignments :

Ce qui va donner :

target      0 TGLLDGKRIL-SGIITDSSIAFHIARVAQEQGAQLVLTGFDRLQVLRRLTDRLPAKAPLL             0 ||||||||||-||||||||||||||-----|||||||||||||....|.||||||||||| query       0 TGLLDGKRILVSGIITDSSIAFHIA-----QGAQLVLTGFDRLRLIQRITDRLPAKAPLL target     59 ELDLQNEEHLASLAGRVTEAIGA 82            60 |||.||||||||||||||||||| 83 query      55 ELDVQNEEHLASLAGRVTEAIGA 78

Qui est plus proche de ce que l’on a l’habitude de voir en termes d’alignement. A noter que nous n’avons pas discuté sur les paramètres de l’alignement, en particulier les pénalités pour l’ouverture et pour l’extension d’un gap, qui peuvent aussi changer le résultat. Nous pourrions penser qu’utiliser un alignement de séquences, au sein de processus structure-based est anecdotique, en fait ce n’est pas le cas, il y a plusieurs applications, par exemple lorsqu’on veut calculer des fluctuations d’éléments structuraux d’une protéine.